Radial basis function network

A radial basis function (RBF) is a real-valued function whose value depends only on the distance from the origin. They are used in function approximation, time series prediction, and control. In artificial neural networks radial basis functions are utilized as activation functions.

Interpolation problem

Radial basis functions can be used to solve a function interpolation problem. A set of known input-output pairs values is needed:

$\left\{\left[\mathbf {x} (t),y(t)\right]:\left[\mathbb {R} ^{n},\mathbb {R} \right]\right\}_{t=1}^{K}$

where

$\mathbf {x} (t)$ is the input vector with index t (or the input at time t),

$y(t)$ is the output indexed with t,

$n$ is the dimension of the input space, and

$K$ is the number of points ($K$ can be infinite).

File:060804 architecture.png

Figure 1: Architecture of a radial basis function network. An input vector x is used as input to several radial basis functions, each with different parameters. The output of the network is a linear combination of the outputs from radial basis functions.

In the deterministic case the data is drawn from the set

$\left\{\left[\mathbf {x} (t),y(t)=f{\big (}\mathbf {x} (t){\big )}\right]\right\}_{t=1}^{K}$.

The data can be noisy, in which case is drawn from the set

$\left\{\left[\mathbf {x} (t),y(t)=f{\big (}\mathbf {x} (t){\big )}+\epsilon (t)\right]\right\}_{t=1}^{K}$

where $\epsilon (t)$ is a partially known random process.

In general stochastic case, data is drawn from the joint probability distribution

$P\left(\mathbf {x} \land y\right)$.

Architecture

RBF networks typically have 3 layers, the input layer, the hidden layer with the RBF non-linearity and a linear output layer. RBF networks have the advantage of not being locked into local minima as do the MLP networks. RBF architectures come in two forms, normalized and unnormalized. The forms can be expanded into a superposition of local linear models.

RBF types

The most popular choice for the non-linearity is the Gaussian.

• Gaussian: $\rho (r)=\exp(-\beta r^{2})$ for some $\beta >0$

Other forms, such as a multiquadratic are also used.

• Multiquadratic: $\rho (r)={\sqrt {r^{2}+\beta ^{2}}}$ for some $\beta >0$

File:060803 unnormalized radial basis functions.png

Figure 2: Two unnormalized radial basis functions in one input dimension. The basis function centers are located at $c_{1}=0.75$ and $c_{2}=3.25$.

Unnormalized

The unnormalized radial basis function architecture, $\varphi :\mathbb {R} ^{n}\to \mathbb {R}$ , is

$\varphi (\mathbf {x} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{i=1}^{N}a_{i}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

where $\varphi$ is the approximation to the data, $\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$, known as a "radial basis function," is a local function of the distance $\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert$ between the input vector $\mathbf {x}$ and a "basis function center"

$\mathbf {c} _{i}$ $(i=1,N)$,

and

$a_{i}$ $(i=1,N)$

are weights to be determined by data. Typically the distance is taken to be the Euclidean distance and the basis function is taken to be Gaussian

$\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}\propto \exp \left[-\beta \left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert ^{2}\right]$.

The weights $a_{i}$, $\mathbf {c} _{i}$, and $\beta$ are determined in a manner that optimizes the fit between $\varphi$ and the data.

Figure 3: Two normalized radial basis functions in one input dimension. The basis function centers are located at $c_{1}=0.75$ and $c_{2}=3.25$.

Normalized

Normalized architecture

The normalized RBF architecture is

$\varphi (\mathbf {x} )\ {\stackrel {\mathrm {def} }{=}}\ {\frac {\sum _{i=1}^{N}a_{i}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}}=\sum _{i=1}^{N}a_{i}u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

where

$u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}\ {\stackrel {\mathrm {def} }{=}}\ {\frac {\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}}$

is known as a "normalized radial basis function."

Figure 4: Three normalized radial basis functions in one input dimension. The additional basis function has center at $c_{3}=2.75$

Theoretical motivation for normalization

There is theoretical justification for this architecture in the case of stochastic data flow. Assume a Stochastic kernel approximation for the joint probability density

$P\left(\mathbf {x} \land y\right)=\sum _{i=1}^{N}\,\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}\,\sigma {\big (}\left\vert y-e_{i}\right\vert {\big )}$

where the weights $\mathbf {c} _{i}$ and $e_{i}$ are exemplars from the data and we require the kernels to be normalized

$\int \rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}\,d^{n}\mathbf {x} =1$

and

$\int \sigma {\big (}\left\vert y-e_{i}\right\vert {\big )}\,dy=1$.

Figure 5: Four normalized radial basis functions in one input dimension. The fourth basis function has center at $c_{4}=0$. Note that the first basis function (dark blue) has become localized.

The probability densities in the input and output spaces are

$P\left(\mathbf {x} \right)=\int P\left(\mathbf {x} \land y\right)\,dy=\sum _{i=1}^{N}\,\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

and

$P\left(y\right)=\int P\left(\mathbf {x} \land y\right)\,d^{n}\mathbf {x} =\sum _{i=1}^{N}\,\sigma {\big (}\left\vert y-e_{i}\right\vert {\big )}$

The expectation of y given an input $\mathbf {x}$ is

$\varphi \left(\mathbf {x} \right)\ {\stackrel {\mathrm {def} }{=}}\ E\left(y\mid \mathbf {x} \right)=\int y\,P\left(y\mid \mathbf {x} \right)dy$

where

$P\left(y\mid \mathbf {x} \right)$

is the conditional probability of y given $\mathbf {x}$. The conditional probability is related to the joint probability through Bayes theorem

$P\left(y\mid \mathbf {x} \right)={\frac {P\left(\mathbf {x} \land y\right)}{P\left(\mathbf {x} \right)}}$

which yields

$\varphi \left(\mathbf {x} \right)=\int y\,{\frac {P\left(\mathbf {x} \land y\right)}{P\left(\mathbf {x} \right)}}\,dy$.

This becomes

$\varphi \left(\mathbf {x} \right)={\frac {\sum _{i=1}^{N}a_{i}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}}=\sum _{i=1}^{N}a_{i}u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

when the integrations are performed.

Local linear models

It is sometimes convenient to expand the architecture to include local linear models. In that case the architectures become, to first order,

$\varphi \left(\mathbf {x} \right)=\sum _{i=1}^{N}\left(a_{i}+\mathbf {b} _{i}\cdot \left(\mathbf {x} -\mathbf {c} _{i}\right)\right)\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

and

$\varphi \left(\mathbf {x} \right)=\sum _{i=1}^{N}\left(a_{i}+\mathbf {b} _{i}\cdot \left(\mathbf {x} -\mathbf {c} _{i}\right)\right)u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

in the unnormalized and normalized cases, respectively. Here $\mathbf {b} _{i}$ are weights to be determined. Higher order linear terms are also possible.

This result can be written

$\varphi \left(\mathbf {x} \right)=\sum _{i=1}^{2N}\sum _{j=1}^{n}e_{ij}v_{ij}{\big (}\mathbf {x} -\mathbf {c} _{i}{\big )}$

where

$e_{ij}={\begin{cases}a_{i},&{\mbox{if }}i\in [1,N]\\b_{ij},&{\mbox{if }}i\in [N+1,2N]\end{cases}}$

and

$v_{ij}{\big (}\mathbf {x} -\mathbf {c} _{i}{\big )}\ {\stackrel {\mathrm {def} }{=}}\ {\begin{cases}\delta _{1j}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )},&{\mbox{if }}i\in [1,N]\\\left(x_{ij}-c_{ij}\right)\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )},&{\mbox{if }}i\in [N+1,2N]\end{cases}}$

in the unnormalized case and

$v_{ij}{\big (}\mathbf {x} -\mathbf {c} _{i}{\big )}\ {\stackrel {\mathrm {def} }{=}}\ {\begin{cases}\delta _{1j}u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )},&{\mbox{if }}i\in [1,N]\\\left(x_{ij}-c_{ij}\right)u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )},&{\mbox{if }}i\in [N+1,2N]\end{cases}}$

in the normalized case.

Here $\delta _{ij}$ is a Kronecker delta function defined as

$\delta _{ij}={\begin{cases}1,&{\mbox{if }}i=j\\0,&{\mbox{if }}i\neq j\end{cases}}$.

Objective functions

Main article: Optimization (mathematics)

The weights, which we signify by $\mathbf {w}$, in the RBF architecture are found through optimization of an objective function. The most common objective function is the least squares function

$K(\mathbf {w} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{t=1}^{\infty }K_{t}(\mathbf {w} )$

where

$K_{t}(\mathbf {w} )\ {\stackrel {\mathrm {def} }{=}}\ {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}^{2}$.

We have explicitly included the dependence on the weights. Minimization of the least squares objective function by optimal choice of weights optimizes accuracy of fit.

There are occasions in which multiple objectives, such as smoothness as well as accuracy, must be optimized. In that case it is useful to optimize a regularized objective function such as

$H(\mathbf {w} )\ {\stackrel {\mathrm {def} }{=}}\ K(\mathbf {w} )+\lambda S(\mathbf {w} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{t=1}^{\infty }H_{t}(\mathbf {w} )$

where

$S(\mathbf {w} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{t=1}^{\infty }S_{t}(\mathbf {w} )$

and

$H_{t}(\mathbf {w} )\ {\stackrel {\mathrm {def} }{=}}\ K_{t}(\mathbf {w} )+\lambda S_{t}(\mathbf {w} )$

where optimization of S maximizes smoothness and $\lambda$ is known as a regularization parameter.

Training

Training the centers and weights to optimize the objective function is typically done in hybrid fashion by first fixing the basis funtion centers and then optimizing the weights. In the sequential training the weights are updated at each time step as data streams in.

Training the basis function centers

Basis function centers can be either randomly sampled among the input instances or found by clustering the samples and choosing the the cluster means as the centers.

Gradient descent training of the linear weights

Main article: Gradient descent

The simplest training algorithm is Gradient descent. In gradient descent training the weights are adjusted at each time step by moving them in a direction opposite from the gradient of the objective function

$\mathbf {w} (t+1)=\mathbf {w} (t)-\nu {\frac {d}{d\mathbf {w} }}H_{t}(\mathbf {w} )$

where $\nu$ is a "learning parameter."

For the case of training the linear weights, $a_{i}$, the algorithm becomes

$a_{i}(t+1)=a_{i}(t)+\nu {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}\rho {\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}$

in the unnormalized case and

$a_{i}(t+1)=a_{i}(t)+\nu {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}u{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}$

in the normalized case.

For local-linear-architectures gradient-descent training is

$e_{ij}(t+1)=e_{ij}(t)+\nu {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}v_{ij}{\big (}\mathbf {x} (t)-\mathbf {c} _{i}{\big )}$

Projection operator training of the linear weights

For the case of training the linear weights, $a_{i}$ and $e_{ij}$, the algorithm becomes

$a_{i}(t+1)=a_{i}(t)+\nu {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}{\frac {\rho {\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho ^{2}{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}}$

in the unnormalized case and

$a_{i}(t+1)=a_{i}(t)+\nu {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}{\frac {u{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}u^{2}{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}}$

in the normalized case and

$e_{ij}(t+1)=e_{ij}(t)+\nu {\big [}y(t)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}{\frac {v_{ij}{\big (}\mathbf {x} (t)-\mathbf {c} _{i}{\big )}}{\sum _{i=1}^{N}\sum _{j=1}^{n}v_{ij}^{2}{\big (}\mathbf {x} (t)-\mathbf {c} _{i}{\big )}}}$

in the local-linear case.

For one basis function, projection operator training reduces to Newton's method.

Figure 6: Logistic map time series. Repeated iteration of the logistic map generates a chaotic time series. The values lie between zero and one. Displayed here are the 100 training points used to train the examples in this section. The weights c are the first five points from this time series.

Examples

Logistic map

The basic properties of radial basis functions can be illustrated with a simple mathematical map, the logistic map, which maps the unit interval onto itself. It can be used to generate a convenient prototype data stream. The logistic map can be used to explore function approximation, time series prediction, and control theory. The map originated from the field of population dynamics and became the prototype chaotic time series. The map, in the fully chaotic regime, is given by

$x(t+1)\ {\stackrel {\mathrm {def} }{=}}\ f\left[x(t)\right]=4x(t)\left[1-x(t)\right]$

where t is a time index. The value of x at time t+1 is a parabolic function of x at time t. This equation represents the underlying geometry of the chaotic time series generated by the logistic map.

Generation of the time series from this equation is the forward problem. The examples here illustrate the inverse problem; identification of the the underlying dynamics, or fundamental equation, of the logistic map from exemplars of the time series. The goal is to find an estimate

$x(t+1)=f\left[x(t)\right]\approx \varphi (t)=\varphi \left[x(t)\right]$

for f.

Function approximation

Unnormalized radial basis functions

The architecture is

Figure 7: Unnormalized basis functions. The Logistic map (blue) and the approximation to the logistic map (red) after one pass through the training set.

$\varphi (\mathbf {x} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{i=1}^{N}a_{i}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

where

$\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}=\exp \left[-\beta \left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert ^{2}\right]=\exp \left[-\beta \left(x(t)-c_{i}\right)^{2}\right]$.

Since the input is a scalar rather than a vector, the input dimension is one. We choose the number of basis functions as N=5 and the size of the training set to be 100 exemplars generated by the chaotic time series. The weight $\beta$ is taken to be a constant equal to 5. The weights $c_{i}$ are five exemplars from the time series. The weights $a_{i}$ are trained with projection operator training:

$a_{i}(t+1)=a_{i}(t)+\nu {\big [}x(t+1)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}{\frac {\rho {\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho ^{2}{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}}$

where the learning rate $\nu$ is taken to be 0.3. The training is performed with one pass through the 100 training points. The rms error is 0.15.

File:060731c Normalized basis functions.png

Figure 8: Normalized basis functions. The Logistic map (blue) and the approximation to the logistic map (red) after one pass through the training set. Note the improvement over the unnormalized case.

Normalized radial basis functions

The normalized RBF architecture is

$\varphi (\mathbf {x} )\ {\stackrel {\mathrm {def} }{=}}\ {\frac {\sum _{i=1}^{N}a_{i}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}}=\sum _{i=1}^{N}a_{i}u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}$

where

$u{\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}\ {\stackrel {\mathrm {def} }{=}}\ {\frac {\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}}}$.

Again:

$\rho {\big (}\left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert {\big )}=\exp \left[-\beta \left\Vert \mathbf {x} -\mathbf {c} _{i}\right\Vert ^{2}\right]=\exp \left[-\beta \left(x(t)-c_{i}\right)^{2}\right]$.

Again, we choose the number of basis functions as five and the size of the training set to be 100 exemplars generated by the chaotic time series. The weight $\beta$ is taken to be a constant equal to 6. The weights $c_{i}$ are five exemplars from the time series. The weights $a_{i}$ are trained with projection operator training:

$a_{i}(t+1)=a_{i}(t)+\nu {\big [}x(t+1)-\varphi {\big (}\mathbf {x} (t),\mathbf {w} {\big )}{\big ]}{\frac {u{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}u^{2}{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}}$

where the learning rate $\nu$ is again taken to be 0.3. The training is performed with one pass through the 100 training points. The rms error on a test set of 100 exemplars is 0.084, smaller than the unnormalized error. Normalization yields accuracy improvement. Typically accuracy with normalized basis functions increases even more over unnormalized functions as input dimensionality increases.

File:060803b chaotic time series prediction.png

Figure 9: Normalized basis functions. The Logistic map (blue) and the approximation to the logistic map (red) as a function of time. Note that the approximation is good for only a few time steps. This is a general characterisitc of chaotic time series.

Time series prediction

Once the underlyting geometry of the time series is estimated as in the previous examples, a prediction for the time series can be made by iteration:

$\varphi (0)=x(1)$

${x}(t)\approx \varphi (t-1)$

${x}(t+1)\approx \varphi (t)=\varphi [\varphi (t-1)]$.

A comparison of the actual and estimated time series is displayed in the figure. The estimated times series starts out at time zero with an exact knowledge of x(0). It then uses the estimate of the dynamics to update the the time series estimate for several time steps.

Note that the estimate is accurate for only a few time steps. This is a general characteristic of chaotic time series. This is a property of the sensitive dependence on initial conditions common to chaotic time series. A small initial error is amplified with time. A measure of the divergence of time series with nearly identical initial conditions is known as the Lyapunov exponent.

Control of a chaotic time series

Figure 10: Control of the logistic map. The system is allowed to evolve naturally for 49 time steps. At time 50 control is turned on. The desired trajectory for the time series is red. The system under control learns the underlying dynamics and drives the time series to the desired output. The architecture is the same as for the time series prediction example.

We assume the output of the logistic map can be manipulated through a control parameter $c[x(t),t]$ such that

${x}_{}^{}(t+1)=4x(t)[1-x(t)]+c[x(t),t]$.

The goal is to choose the control parameter in such a way as to drive the time series to a desired output $d(t)$. This can be done if we choose the control paramer to be

$c_{}^{}[x(t),t]\ {\stackrel {\mathrm {def} }{=}}\ -\varphi [x(t)]+d(t+1)$

where

$\varphi [x(t)]\approx f[x(t)]=x(t+1)-c[x(t),t]$

is an approximation to the underlying natural dynamics of the system.

The learning algorithm is given by

$a_{i}(t+1)=a_{i}(t)+\nu \varepsilon {\frac {u{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}{\sum _{i=1}^{N}u^{2}{\big (}\left\Vert \mathbf {x} (t)-\mathbf {c} _{i}\right\Vert {\big )}}}$

where

$\varepsilon \ {\stackrel {\mathrm {def} }{=}}\ f[x(t)]-\varphi [x(t)]=x(t+1)-c[x(t),t]-\varphi [x(t)]=x(t+1)-d(t+1)$.

See also

• Artificial neural networks
• Predictive analytics
• Chaos theory
• Autoregressive moving average model
• Autoregressive integrated moving average
• Autoregressive conditional heteroskedasticity
• Mitchell Feigenbaum

External links

• Model Predictive Control with radial basis functions
• Control of a negative ion source

References

• J. Moody and C. J. Darken, "Fast learning in networks of locally tuned processing units," Neural Computation, 1, 281-294 (1989). Also see Radial basis function networks according to Moody and Darken
• T. Poggio and F. Girosi, "Networks for approximation and learning," Proc. IEEE 78(9), 1484-1487 (1990).
• Roger D. Jones, Y. C. Lee, C. W. Barnes, G. W. Flake, K. Lee, P. S. Lewis, and S. Qian, “Function approximation and time series prediction with neural networks,” Proceedings of the International Joint Conference on Neural Networks, June 17-21, p. I-649 (1990).
• Martin D. Buhmann, M. J. Ablowitz (2003). Radial Basis Functions : Theory and Implementations. Cambridge University. ISBN 0-521-63338-9.
• Yee, Paul V. and Haykin, Simon (2001). Regularized Radial Basis Function Networks: Theory and Applications. John Wiley. ISBN 0-471-35349-3.{{cite book}}: CS1 maint: multiple names: authors list (link)
• John R. Davies, Stephen V. Coggeshall, Roger D. Jones, and Daniel Schutzer, "Intelligent Security Systems," in Freedman, Roy S., Flein, Robert A., and Lederman, Jess, Editors (1995). Artificial Intelligence in the Capital Markets. Chicago: Irwin. ISBN 1-55738-811-3. {{cite book}}: |author= has generic name (help)CS1 maint: multiple names: authors list (link)