Interpolation models with multiple hyperparameters.
A traditional interpolation model is characterized by the choice of regularizer applied to the interpolant, and the choice of noise model. Typically, the regularizer has a single regularization constant , and the noise model has a single parameter . The ratio / alone is responsible for determining globally all these attributes of the interpolant: its complexity, flexibility, smoothness, characteristic scale length, and characteristic amplitude. We suggest that interpolation models should be able to capture more than just one flavour of simplicity and complexity. We describe Bayesian models in which the interpolant has a smoothness that varies spatially. We emphasize the importance, in practical implementation, of the concept of conditional convexity when designing models with many hyperparameters. We apply the new models to the interpolation of neuronal spike data and demonstrate a substantial improvement in generalization error.
Introduction
In this paper our philosophy of generalization is as follows: the best generalization
will be achieved by a Bayesian model that is well-matched to the problem
and that is accurately implemented. The aim of obtaining the best generalization
is thus subsumed under the aim of searching for good models. In this
paper we expand the space of interpolation models by introducing additional
hyperparameters, and demonstrate that the generalization performance on a
real problem is substantially improved.
A traditional linear interpolation model 'H 1 ' is characterized by the choice
of the regularizer R, or prior probability distribution, that is applied to the
interpolant; and the choice of noise model N . The choice of basis functions A
used to represent the interpolant may also be important if only a small number
of basis functions are used. Typically the regularizer is a quadratic functional
of the interpolant and has a single associated regularization constant ff, and the
noise model is also quadratic and has a single parameter fi. For example, the
splines prior for the function y(x) (Kimeldorf and Wahba 1970) is: 1
log
ff
Z
dx [y (p)
where y (p) denotes the pth derivative of y. The probability of the data measurements
assuming independent Gaussian noise is:
log
(The constants in equations (1) and (2) are functions of ff and fi respectively.)
When we use these distributions with find the most probable y(x) we
obtain the cubic splines interpolant. For any quadratic regularizer and quadratic
1 Strictly this prior is improper since addition of an arbitrary polynomial of degree
y(x) is not constrained. It can be made proper by adding terms corresponding to boundary
conditions to (1). In the present implementations of the models, we enforce the boundary
conditions
Figure

1: An inferred spike signal from a zebra finch neuron.
Courtesy of M. Lewicki and A. Doupe, California Institute of Technology.
log likelihood, the most probable interpolant depends linearly on the data values.
This is the property by which we define a 'linear' interpolation model.
Such models may be optimized and compared using Bayesian methods as
reviewed in MacKay (1992). In such models, for fixed fi, the ratio ff=fi alone
determines globally all the following attributes of the interpolant: its complexity,
flexibility, smoothness, characteristic scale length, and characteristic amplitude.
whilst some of these terms may be synonyms, surely others describe distinct
properties. Should not our models be able to capture more than just one
flavour of simplicity and complexity? And should not the interpolant's smooth-
ness, for example, be able to vary spatially?
1.1 Example: Neural spike modelling
An example of a function from a real system is shown in figure 1; this is the action
potential of a neuron deduced from recordings of 40 distinct events (Lewicki
1994). The graph was created by fitting a simple spline model (with
the data. This function has one 'spiky' region with large characteristic amplitude
and short spatial scale. Elsewhere the true function is smooth. However the
fitted function shown in figure 1, controlled by only one regularization constant
ff, overfits the noise on the right, having a rough appearance where it should
plausibly be smooth. The value of ff appropriate for fitting the spiky region is
too small for the rest of the curve. It would be useful here to have a model
capable of capturing the concepts of local smoothness, because such a model,
having a prior better matched to the real world, would require less data to
yield information of the same quality. Furthermore, when different hypotheses
are compared, broad priors introduce a bias toward simpler hypotheses. For
example, if we ask whether one or two distinct spike functions are present in
a data set, the traditional model's prior with small ff will bias the conclusion
in favour of the single spike function. Only with well-matched priors can the
results of Bayesian hypothesis comparison be trusted.
In this paper we discuss methods for introducing multiple flavours of simplicity
and complexity into a hierarchical probabilistic model in a computationally
tractable way, and demonstrate new interpolation models with multiple hyper-parameters
that capture a spatially varying smoothness.
Prior work making use of variable hyperparameters includes the modelling of
data with non-Gaussian innovations or observation noise (see, e.g., (West 1984;
Carter and Kohn 1994; Shephard 1994)). The interpolation models we propose
might be viewed as Bayesian versions of the 'variable bandwidth' kernel regression
technique (Muller and Stadtmuller 1987). The aim of our new model is
also similar to the goal of inferring the locations of discontinuities in a function,
studied by Blake and Zisserman (1987). Traditional interpolation models have
difficulty with discontinuities: if the value of ff=fi is set high, then edges are
blurred out in the model; if ff=fi is lowered the edge is captured, but ringing
appears near the edge, and noise is overfitted everywhere. Blake and Zisserman
introduce additional hyperparameters defining the locations of edges. The
models they use are computationally non-convex, so that finding good representatives
of the posterior distribution is challenging. They use 'graduated
non-convexity' techniques to find good solutions. By contrast we attempt to
create new hierarchical models that are, for practical purposes, convex.
Tractable hierarchical modelling: Convexity
Bayesian statistical inference is often implemented either by Gaussian approximations
about modes of distributions, or by Markov Chain Monte Carlo methods
(Smith 1991). Both methods clearly have a better chance of success if the
posterior probability distribution over the model parameters and hyperparameters
is not dominated by multiple distinct optima. If we know that most of the
probability mass is in just one 'hump', then we know that we need not engage in
a time-consuming search for the more probable optima, and we might hope that
some approximating distribution (e.g., involving the mode of the distribution)
might be able to capture the key properties of that hump. Furthermore, convex
conditional distributions may be easier to sample from with, say, Gibbs sampling
methods (Gilks and Wild 1992). It would be useful if all the conditional
and marginal probability distributions of our models were log convex:
probability distribution is log convex if there is a representation
x of the variables such that the matrix M defined by
log P (x) (3)
is everywhere positive definite.
It is hard, however, to make interesting hierarchical models such that all
conditional and marginal distributions are log convex. We introduce a weaker
criterion:
model is conditionally convex if its variables can be divided
into groups such that, for every group, their distribution conditioned on any
values for the other variables is log convex.
An example of a conditionally convex model is the traditional interpolation
model with three groups of variables: D (data), w (parameters), and ff (one
hyperparameter). The probability distribution P (Djw;
convex over D (it is Gaussian). The distribution P (wjD; ff) is log convex over
w (it is Gaussian). And the distribution P (ffjw;
over ff (it is a Gamma distribution).
That a model is conditionally convex does not guarantee that marginal distributions
of all variables are unimodal. For example the traditional model's
posterior marginals P (wjD) and P (ffjD) are not necessarily unimodal; but good
unimodal approximations to them can often be made (MacKay 1996). So we
conjecture that conditional convexity is a desirable property for a tractable
model.
We now generalize the spline model of equation (1) to a model with multiple
hyperparameters that is conditionally convex, and demonstrate it on the neural
spike data. We then discuss general principles for hierarchical modelling with
multiple hyperparameters.
3 A new interpolation model
We replace the regularizer of equation (1) by:
log
Z
dx ff(x)[y (p)
where ff(x) is written in terms of hyperparameters thus:
and the constant of equation (4) is a function of ff(x; u) which becomes important
when ff(x; u) is inferred. The exponentiated quantity has the form
of a linear interpolant using basis functions / h (x). In the special case
we obtain the traditional single alpha model. This representation
is chosen because (1) it embodies our prior belief that ff(x) should be a
smooth function of x, and (2) the model is conditionally convex (a partial proof
is given in section 4).
When implementing this model we optimize the hyperparameters u and fi
by maximizing the marginal likelihood or 'evidence',
Z
where k is the dimensionality of our representation y of y(x). Some authors view
this 'empirical Bayes' approach as controversial and inaccurate (Wolpert 1993),
but it is widely used under various names such as 'ML-II', and is closely related
to 'generalized maximum likelihood' (Gu and Wahba 1991). The ideal Bayesian
method would put a proper prior on the hyperparameters and marginalize over
them, but optimization of the hyperparameters is computationally more convenient
and often gives predictive distributions that are indistinguishable (MacKay
1996).
We use a discrete representation of y(x) and ff(x) on a finely spaced grid,
y, ff(x; u) ! fff c j ff(x c ; u)g and / In
this representation the Hessian of the log posterior is a sum of band-diagonal
terms from the log prior and a diagonal matrix from the log likelihood, A j
\Gammarr logP (yjD; fffg;
fiI. The gradient of the log evi-
dence, which we use for the optimization, is then:
@
log
@
log P (Djfff c g) (7)
where
@
log P (Djfff c
Trace
-2000
-2000

Figure

2: Traditional models: 2. The diamond-shaped points in
the upper plots are the artifical data. The solid line shows the most probable
interpolant found using the traditional single alpha model. The predictive error
bars (dotted lines) are one-standard-deviation error bars. The lower row shows
the errors between the interpolant and the original function to which the noise
was added to make the artificial data. The predictive error bars are also shown.
Contrast with figure 3.
-2000
-2000

Figure

3: New models with multiple hyperparameters: 2. Top
row: The diamond-shaped points are the artifical data. The solid line shows
the most probable interpolant and the predictive error bars (dotted lines) are
one-standard-deviation error bars. Second row: the inferred ff(x) on a log
scale (contrast with the values of 5.9 \Theta10 \Gamma7 and 2:0 \Theta 10 \Gamma6 inferred for the
traditional models). The third row shows the nine basis functions / used to
represent ff(x). The bottom row shows the errors between the interpolant and
the original function to which the noise was added to make the artificial data.
The predictive error bars are also shown. The top and bottom graphs should
be compared with those of figure 2. 9

Table

1: Comparison of models on artificial data.
The first three columns give the evidence, the effective number of parameters, and the
RMS error for each model when applied to the data shown in figures 2-3. The fourth
column gives the RMS error averaged over four similar data sets.
Model log fl RMS avg. RMS
Evidence error error
3.1 Demonstration
We made an artificial data set by adding Gaussian noise of standard deviation
1000 to the function depicted in figure 1. [This function plays the role, in
these experiments, of a true underlying function; the presence of some actual
roughness in this function is believed to be unimportant since our chosen noise
level is substantially greater than the apparent size of the roughness.] Figure
2 shows the data, interpolated using the traditional single alpha models with
2. The hyperparameter ff was optimized by maximizing the
evidence, as in Lewicki (1994). The noise level oe - was set to the known noise
level. In order for the spiky part of the data to be fitted, ff has to be set to
a small value, and the most probable interpolant is able in both models to go
very close to all the data points. There is considerable overfitting everywhere,
and the predictive error bars are large everywhere.
We then interpolated the data with two new models defined by equations (4)
and (5), with 2. We set the basis functions / to the hump-shaped
functions shown in figure 3. These functions define a scale length on which the
smoothness is permitted to vary. This scale length was optimized roughly by
maximizing the evidence. The new models had nine hyperparameters u. These
hyperparameters were set by maximizing the evidence using conjugate gradi-
ents. Because the new models are conditionally convex, we had hoped that the
maximization of the evidence would lead to a unique optimum uMP . However,
there were multiple optima in the evidence as a function of the hyperparam-
eters; but these did not cause insurmountable problems. We found different
optima by using different initial conditions u for the optimization. The best
evidence optima were found by initializing u in a way that corresponded to our
prior knowledge that neuronal spike functions start and end with a smooth re-
gion; we set u initially to fu h 0g. This
prior knowledge was not formulated into an informative prior over u during
the optimization, though doing so would probably be a good idea for practical
purposes.

Figure

3 shows the solutions found using the new interpolation models with
2. The inferred value of ff is small in the region of the spike, but
elsewhere a larger value of ff is inferred, and the interpolant is correspondingly
smoother.
The log evidence for the four models is shown in table 1. The reported evidence
values are log e P (Djff MP we were to make a
proper model comparison we would integrate over the hyperparameters; this integration
would introduce additional small subjective Occam factors penalizing
the extra hyperparameters in H 2 , c.f. MacKay (1992). The root mean square
errors between the interpolant and the original function to which the noise was
added to make the artificial data are given in table 1, and the errors themselves
are displayed at bottoms of figures 2-3.
By both the evidence value and the RMS error values, the new models are
significantly superior to the traditional model. Table 1 also displays the value of
the 'effective number of well-determined parameters' (Gull 1989; MacKay 1992),
fl, which, when the hyperparameters are optimized, is given by:
Z
dx ff(x)y (p)
The smaller the effective number of parameters, the less overfitting of noise
there is, and the smaller the error bars on the interpolant become. The total
number of parameters used to represent the interpolant was in all cases 100.
3.2 Model criticism
It is interesting to assess whether the observed errors with respect to the original
function are compatible with the one-standard-deviation error bars that were
obtained. These are shown together at the bottom of figure 3. The errors are
only significantly larger than the error bars at the leftmost five data points,
where the small amount of noise in the original function is incompatible with
the assumed boundary conditions Omitting those five
data points, we find for the new model that the other 95 errors have
expectation 95 \Sigma 14), and for the
of the 95 errors in either case exceed 2.5 standard deviations. We therefore
see no significant evidence for the observed errors to be incompatible with the
predictive error bars.
3.3 Discussion
These new models offer two practical benefits. First, while the new models
still fit the spiky region well (indeed the errors are slightly reduced there),
they give a smoother interpolant elsewhere. This reduction in overfitting allows
more information to be extracted from any given quantity of experimental data;
neuronal spikes will be distinguishable given fewer samples. To quantify the
potential savings in data we fitted the four models to fake data equivalent to
independent observations of the function shown in figure 1, that is,
data points with noise level oe (we we did this
Traditional p=1
Traditional p=2

Figure

4: Average RMS error of the traditional and new models as a function
of amount of data
by decreasing the actual noise level in the artificial data). The figures and
tables shown thus far correspond to the case of one observation, In
figure 4 we show the RMS error of each model as a function of the number of
data points, averaged over four runs with different artificial noise. To achieve
the same performance (RMS error) as the new models, the traditional models
require about three times as much data.
Second, the new models have greater values of the evidence. This does not
only mean that they are more probable models (assuming that the omitted
Occam factors for the hyperparameters are smaller than these evidence differ-
ences). It also means that model comparison questions can be answered in a
more reliable way. For example, if we wish to ask 'are two distinct spike types
present in several data sets or just one?' then we must compare two hypotheses:
HB , which explains the data in terms of two spike functions, and HA , which
just uses one function. In such model comparisons, the 'Occam factors' that
penalize the extra parameters of HB are important. If we used the traditional
interpolation model, we would obtain Occam factors about e 20 bigger than those
obtained using the new interpolation model. Broad priors bias model comparisons
toward simpler models. The new interpolation model, when optimized,
produces a prior in which the effective number of degrees of freedom of the
interpolant is reduced so that the prior is less broad.
Of course, inference is open-ended, and we expect that these models will in
turn be superceded by even better ones. Close inspection of figure 3 reveals
that the smoothness assumption on the regularizer may be imperfect - we
know from prior experience that the true function's spikiness is confined to a
very small time interval, but the new model gives a jagged interpolant in the
time interval before the spike too because the function ff(x) is assumed to vary
smoothly. Future models might include a continuum of alternative values of p
(non-integer values of p can be implemented in a Fourier representation). It
might also make sense for the characteristic length scale of the basis functions
/ with which ff(x) is represented to be shorter where ff is small.
The advantages conferred by the new models are not accompanied by a significant
increase in computational cost. The optimization of the hyperparameters
requires that the Hessian matrix be inverted a small number of times.
Other approaches to the implementation of models with multiple hyperparameters
could be considered. The confidence intervals in the present approach,
in which the hyperparameters are optimized, are likely to be too small. One
could use Markov chain Monte Carlo methods such as Gibbs sampling or hybrid
Monte Carlo, both of which would involve a similar computational load (see
Neal (1993) for an excellent review). We have used the Gibbs sampling software
'BUGS' (Thomas et al. 1992) to implement a similar interpolation model
in which the Gaussian noise level is a spatially varying function fi(x) (MacKay
1995).
4 Some Generalizations
4.1 Strategies for making models with multiple hyperparameter

We now discuss more generally the construction of hierarchical models with
multiple hyperparameters.
Consider a Gaussian prior on some parameters w, equivalent to the function
y(x) in the earlier example. There are various ways of defining a model with
multiple hyperparameters such that each hyperparameter controls a different
flavour of simplicity or complexity in w.
Sum Model
Firstly, one might define the inverse covariance matrix as a sum:
exp
\Gamma2
where fC c g are arbitrary positive semi-definite matrices and ff c - 0; 8c.
Covariance Sum Model
Secondly, one might define the covariance matrix as a sum:
Z
exp@ \Gamma2 w T
with hyperparameters ' c - 0; 8c.
Exponential Sum Model
Thirdly, we can take a sum model of the form (10) (though not necessarily using
the same matrices fC c g) and rewrite the coefficients as an exponential sum:
with hyperparameters u h 2 (\Gamma1; 1), so that
exp
\Gamma2
exp
These models have very different capabilities.
The sum model implements the paradigm of starting from a flexible distri-
bution, then adding in extra terms ff c C c so as to kill degrees of freedom. This
model has no way of introducing selective flexibility. If one hyperparameter ff c
is large, there is no way that other hyperparameters can be set to undo the
stiffness introduced.
The covariance sum model uses an alternative paradigm, starting from a stiff
distribution, and introducing lacunae of flexibility into it.
The important difference between these two paradigms is that whereas the
sum model is conditionally convex, the covariance sum model is not; it is possible
for there to be multiple optima over the hyperparameters even in the limit of
perfect data. This will be demonstrated and explained subsequently.
The exponential sum model, of which the interpolation model of section 3 is
an example, is intended to combine the best of both worlds. Consider the case
where the matrix elements / ch are non-negative. As one hyperparameter u h is
increased, it introduces selective stiffness, and as it is decreased, it introduces
selective flexibility. The model, being a reparameterization of the sum model, is
still conditionally convex (as long as / does not have pathological properties).
4.2 Convexity of the sum model
We give a partial proof of conditional convexity for the sum model. It is
straightforward to confirm that the conditional distributions P (Djw; fffg) and
are log convex. The non-trivial property is that P (fffgjw; D) /
(fffg)P (wjfffg) is convex. We assume that the prior over fffg is defined to be
a) Sum model b) Covariance sum c) Exponential sum

Figure

5: Toy problem probability contours.
Each figure shows the likelihood of two hyperparameters given
Hyperparameters ff a ; ca and ua are on the horizontal axes, and ff b ; cb and ub on the
vertical axes. In all figures the top e 8 of the function is shown and the contours are
equally spaced in log probability.
convex and examine the second factor. Defining
@ff c @ff d
log
This second derivative is negative definite.
Proof
For arbitrary x,
c;d
d
4.3 A toy illustration
As an illustration, we examine the conditional convexity of a model that assigns
a zero-mean Gaussian distribution to a three component vector w. This
distribution is to be parameterized by two hyperparameters. For simplicity, we
assume w is directly observed: This choice of w favours
priors that give flexibility to component 2. Components 1 and 3 do not call for
such flexibility.
Sum model: We build M as a sum of two matrices, diag(1,1,0) and
diag(0,1,1).
Figure 5a shows the log probability log P (wjfffg) as a function of log ff a and
log ff b . The function is convex.
Covariance sum model: We now build M \Gamma1 as a sum of diag(1,1,0) and
diag(0,1,1), letting:
diag
c a
Figure 5b shows the log probability log P (wjfcg) as a function of log c a and
log c b . The function is not convex. The two alternative flavours of flexibility
compete with each other to give the required variance for component 2 of w.
Either we may switch on c a to a large value, or we may switch on c b - but we
may not switch on both to an intermediate degree.
Exponential sum model: We build M as a sum of three matrices, diag(1,0,0),
diag(0,1,0). and diag(0,0,1), with the aid of basis functions /
diag
This model has the same number of hyperparameters as the previous two models
but uses them differently. Figure 5c shows the log probability log P (wjfug) as
a function of u a and u b . The function is convex. Two alternative flavours of
flexibility are embodied, but (just) do not compete with each other destructively.
The sum model starts from flexibility and adds in constraints of stiffness that
kill degrees of freedom in w. The covariance sum representation starts from stiffness
and adds in selective flexibility to create required degrees of freedom. The
covariance sum model is not convex because different forms of flexibility compete
to account for the data. There is a struggle for existence, because any potential
piece of flexibility is penalized by Occam factors in the det M term, encouraging
it to stay switched off. In contrast, alternative ways of introducing stiffness (as
in the sum model and the exponential sum model) do not compete. If two sorts
of stiffness are compatible with the data, they can both be switched on without
incurring any penalty. This is why the sum model is convex. The exponential
sum model, we conjecture, pushes flexibility to the limits of convexity. We
believe these ideas may be relevant to the design of computationally tractable
Gaussian process models for non-linear regression (Williams and Rasmussen
1996).
4.4 How to represent a covariance matrix
In this paper we have used interpolation of neural spike data as a test bed for
the new models. We now discuss another task to which the general principles
we have discussed may apply.
Imagine that we wish to model correlations between k variables
that are assumed to be Gaussian with a covariance matrix V that varies with
other variables x. How should this varying covariance matrix V(x) be param-
eterized? We assume that a representation V(U(x)) is to be used. We would
like the parameterization V(U) to satisfy the following desiderata.
1. Any setting of the parameters U should produce a valid positive definite
matrix V.
2. Any positive definite matrix V should be realizable by a unique value of
the parameters U.
3. The parameterization and its inverse should be continuous and differentiable

4. The representation should treat all indices of the covariance matrix sym-
for example, the first row of V should not be treated differently
from the second row.
5. U should have k(k degrees of freedom, that being the number of
independent elements in the symmetric matrix V.
6. Finally we would like the representation to be conditionally convex; that
is, given one or more vectors y, the conditional probability of U should
be log convex.
These desiderata rule out most obvious representations of V. The raw matrix
V is not permitted because it violates desideratum 1. A triangular decomposition
violates 4. An eigenvector / eigenvalue representation violates 2,3,5. The
'variance component model' representation used for example in Gu and Wahba
(1991) is a covariance sum representation and violates desiderata 5 and 6.
The ideas of this paper motivate the following representation, which is conditionally
convex. Let y be k dimensional, and let R k\Gamma1 be the unit spherical
surface, with v being a unit vector in that space. As parameters we introduce
a symmetric matrix U that is not constrained to be positive definite. Then we
as the inverse of a sum of outer products thus:
d
This representation satisfies all the desiderata. Since this may not be self-
evident, we include a sketch of a proof of half of property 2, namely, that the
mapping from U to V is one to one. We first transform into the eigenvector
basis of U (by an orthogonal transformation that leaves R k\Gamma1 invariant) and
prove that the eigenvectors feg of U are also eigenvectors of V. Let fw i g be
the components of v in the eigenvector basis so that
where the
eigenvectors and eigenvalues of U satisfy Ue
. Then from equation
(20) we can write
Z
d
The integrand, for i 6= j, is antisymmetric in w i and w j , so the integral is zero
in these cases. Thus
e (i) e T
Z
d
that is, V has the same eigenvectors as U, and its eigenvalues are given by:
Z
d
Then the mapping from U to V is one to one if the above mapping from the
eigenvalues of U, f- U g, to the eigenvalues of V, f- V g, is one to one. We
differentiate equation (23) to obtain the Jacobian; if the Jacobian is full-rank
then the mapping is one to one.
@- U
Z
d
This Jacobian is a sum of outer products of positive vectors z given by z
so it either defines a positive semi-definite or a positive definite matrix. The
matrix can only be positive semi-definite if there is some direction h such that
having non-zero measure under the integral over R k\Gamma1 .
Because the integral is over all of R k\Gamma1 , there is no such vector h. Thus the
matrix is full rank, and the mapping from U to V is one to one.
The only problem with this representation is that it involves a high-dimensional
integral. We propose for practical purposes the following approximation:
"C
exp(v T
c
where fv c g C
are fixed unit vectors lying in R k\Gamma1 , selected either at random
or systematically. This representation is conditionally convex and is able to
in the limit C !1.
Conclusions
This work builds on a data modelling philosophy previously illustrated by
work on the 'Automatic Relevance Determination' model for neural networks
(MacKay 1994; Neal 1996): use a huge, flexible model with an essentially infinite
number of parameters; and control the complexity of the model with
sophisticated regularizers. Models with large numbers of hyperparameters can,
if carefully designed, be practically implemented. The hyperparameters reduce
the effective number of degrees of freedom of the model in a manner appropriate
to the the properties of the data, leading to substantial improvements in
generalization error.

Acknowledgements

D.J.C.M. thanks the Isaac Newton Institute and T. Matsumoto, Waseda Uni-
versity, for hospitality, and Radford Neal, Mike Lewicki, David Mumford and
Brian Ripley for helpful discussions. R.T. thanks T. Matsumoto for his support.
We also thank the referees for helpful feedback.



--R

Visual Reconstruction.
On Gibbs sampling for state-space models

Applied Statistics 41:
Minimizing GCV/GML scores with multiple smoothing parameters via the Newton method.
Developments in maximum entropy data analysis.
A correspondence between Bayesian estimation of stochastic processes and smoothing by splines.
Bayesian modeling and classification of neural signals.
Neural Computation 6 (5):
Bayesian interpolation.
Bayesian non-linear modelling for the prediction competition
Probabilistic networks: New models and new meth- ods
Hyperparameters: Optimize
Variable bandwidth kernel estimators of regression-curves

Bayesian Learning for Neural Networks.
Partial non-Gaussian state-space
Bayesian computational methods.
BUGS: A program to perform Bayesian inference using Gibbs sampling.
Outlier models and prior distributions in Bayesian linear- regression
Gaussian processes for regression.
On the use of evidence in neural networks.
--TR
convexity;regularizer;bayesian inference;neuronal spike;spline;regression