Bayesian MARS.
A Bayesian approach to multivariate adaptive regression spline (MARS) fitting (Friedman, 1991) is proposed. This takes the form of a probability distribution over the space of possible MARS models which is explored using reversible jump Markov chain Monte Carlo methods (Green, 1995). The generated sample of MARS models produced is shown to have good predictive power when averaged and allows easy interpretation of the relative importance of predictors to the overall fit.
Introduction
A common problem in statistics, and other disciplines, is to approximate adequately
a function of several variables. In a statistical setting this is known
as multiple regression and the task can be performed either parametrically
by global modelling (e.g. linear regression), or nonparametrically (see below
for examples of these methods).
The aim is to model the dependence of a response variable Y on one
or more predictor variables the data is given by
The data is assumed to come
from a relationship described by
where f is an unknown regression function that we wish to estimate, ffl is a
zero-mean error distribution, most commonly assumed to be Gaus-
sian, and D is the domain of interest, usually taken to be the convex hull
defined by the predictor variables. The regression function f gives the predictive
relationship of Y on X, i.e. the conditional expectation of Y given X.
Thus we may use f to predict future values of Y at previously unseen points
in the domain D. The aim of the regression analysis is to use the data to
construct an estimate b
f (X) to the true regression function which can serve
as a reasonable approximation over the domain of interest, D.
Many methods exist to model the function of interest f [e.g. Additive
models, Hastie and Tibshirani (1990); CART, Breiman et al. (1984); Projection
pursuit regression, Friedman and Stuetzle (1981); Alternating conditional
expectation, Breiman and Friedman (1985)]. We, however, concentrate
on the multivariate adaptive regression spline (MARS) methodology
proposed by Friedman (1991). This method seems to be both highly flexible
and easily interpretable. It was motivated by the recursive partitioning
approach to regression (Morgan and Sonquist, 1963; Breiman et al., 1984)
but produces a continuous model which can be made to have continuous
derivatives and has greater flexibility to model relationships that are nearly
additive or involve at most a few variables. The model can be represented in
such a way that the additive contributions of each predictor variable and the
interactions between variables can be easily identified which helps to identify
variables which are important in the model.
To highlight the progression from recursive partition regression to MARS
we start by giving the partition regression model,
and the a i are the suitably chosen coefficients of
the basis functions B i and k is the number of basis functions in the model.
These basis functions are such that B i I is the indicator
function which is one where the argument is true, zero elsewhere and
the R i a partition of D.
The usual MARS model is the same as that given in (1) except that the
basis functions are different. Instead the B i are given by
(2)
is the degree of the interaction of basis B i , the s ji ,
which we shall call the sign indicators, equal \Sigma1, the v(j; i) give the index
of the predictor variable which is being split on and the t ji (known as knot
points) give the position of the splits. The v(j; \Delta) (j = are constrained
to be distinct so each predictor only appears once in each interaction term.
See Section 3, Friedman (1991) for a comprehensive illustration of the model.
To illustrate the initially confusing notation we present an example. Suppose
a MARS model contains the basis function B i given by
We can immediately see that there are two factors in the interaction term
so 2. The sign indictors are s with the knot points
given by t \Gamma3:7 and the labels for the predictors split on are
The MARS model is continuous in D and can be made to have continuous
first derivatives by replacing the truncated linear basis functions B i by
truncated cubic basis functions. This has the effect of "rounding" the basis
function at the split points.
The MARS algorithm proceeds as follows. A forward stepwise search
for basis functions takes place with the constant basis function the only one
present initially. At each step the split which minimises some "lack-of-fit"
criterion from all the possible splits on each basis function is chosen. Splits
are only permissable at the marginal predictor values. If the split was on
predictor x   at t   this corresponds to the two new basis func-
tions, henceforth referred to as a . Note that unlike the
recursive partitioning algorithm the basis function with which the split was
made is not removed. This continues until the model reaches some predetermined
maximum number of basis functions, which should be about twice
the number expected in the model to aid the subsequent backwards stepwise
deletion of basis functions. This just involves removing basis functions one
at a time until the lack-of-fit criterion is at a minimum. The basis which
improves the fit the most or degrades it the least is removed at each step.
Finally the resulting model can be made to have a continuous first derivative
by "rounding" at the split points as mentioned above. The lack-of-fit
measure used by Friedman (1991) is the generalised cross-validation criterion
which was originally proposed by Craven and Wahba (1979).
The aim of this paper is to provide a Bayesian algorithm which mimics the
MARS procedure. This is done by considering the number of basis functions,
along with their type (see Section 2.1), their coefficients and their form (the
positions of the split points and the sign indicators) random. We treat these
as additional parameters in the problem and make inference about them
using the data.
The problem of routine calculation of the posterior distribution of the
models is addressed by designing a suitable Markov chain Monte Carlo (MCMC)
reversible jump simulation algorithm as set out by Green (1995). The simulated
sample contains many different MARS models with corresponding
posterior weights but if a estimate for f with high predictive power is all
that is required then pointwise averaging over all the models in the sample
is suggested.
This work is an extension to the Bayesian approach to curve fitting in
one dimension given by Denison et al. (1998b) and is related to the Bayesian
CART algorithms proposed by Denison et al. (1998a) and Chipman et al.
(1998).
In Section 2 we outline the Bayesian MARS method and show examples
using the method on simulated data in Section 3 and real data in Section 4.
Section 5 contains a discussion of the proposed methodology.
2 The Bayesian MARS Method
2.1 Basic Ideas
We must first define what we mean by the type of a basis function. Using
the notation in (1) and (2), we consider basis functions to be
of the same type if is identical to some permutation of
Hence with m predictor variables there are N different
types of basis function where N is given by
Note that the sum does not include the constant basis function term [for
which i would equal 0 in (3)] as this basis B 1 is always the sole constant
basis function in the model so it cannot be chosen as a candidate basis.
Frequently some maximum order of interaction I is assigned to the model
such that J i - I in which case the sum in (3) only runs from
1 up to I. We let the type of basis function B i
thus T i , in effect, just tells us which predictor variables we are splitting on,
ie what the values of v(1; are.
As an example suppose we have a problem with just two predictors
2). Then the types of basis functions that could be in the model
(not including the constant one) are [\Sigma(x
(say type 2) and [\Sigma(x (say type 3). So for all those basis
functions which split only on predictor x 2 their types T i are all equal to 2.
We propose a model which can be used to set up a probability distribution
over the space of possible MARS structures. Any MARS model can be
uniquely defined by the number of basis functions present, the coefficients
and the types of the basis functions, together with the knot points and the
sign indicators associated with each interaction term. This means that we
make the k; a random with the J i uniquely defined via the T i .
We change the dimension of the model when we change k and so, for
Bayesian computation, we use a reversible jump MCMC approach (Green,
1995; Richardson and Green, 1997) when we are considering changes in the
number of basis functions in the model.
Inference is carried out assuming that the "true" model is unknown but
comes from the class of models denotes the model with
exactly k basis functions. The overall parameter space \Theta can then be written
as a countable union of subspaces
is a subspace of the
Euclidean space R n(k) , where R n(k) denotes the
dimensional
parameter space corresponding to model M k . Here '
where each B i is the 2(1+J i )-dimensional vector (a
which corresponds to basis function B i .
There is a natural hierarchical structure to this setup, which, denoting
a generic element of \Theta k by ' (k) and the data vector by y, we formalise by
modelling the joint distribution of (k; '
that is, as the product of model probability, parameter prior and likelihood.
Bayesian inference about k and ' (k) will be based on the joint posterior
which we shall explore and summarise by regarding it as the
target distribution for tailored MCMC computations. It will often be useful
to consider this in the factorised form
We will generate samples from the joint posterior of (k; ' (k) ) by using a
class of reversible jump Metropolis-Hastings algorithms (Green, 1995). Full
details of the method can be found in the reference cited. Here, we focus on
the essence of the methodology and the particular forms of the algorithms in
our current context.
2.2 The Bayesian Model
We assume ffl in (1) follows a N(0; oe 2 ) distribution where oe 2 is unknown. As
a result we extend the parameter vector ' to include this new unknown. This
leads to the log-likelihood of the model, l k ('jy), being given by
l k \Gamman log oe \Gamma2oe 2
f is of the form given in (1) and (2).
We use a vague, but proper, prior for the variance of the error distribution
ie -(oe \Gamma2 and the T i are assumed to be uniformly
distributed on Ng. The sign indicators s ji and knot points t ji are
also assumed uniform on the sets f\Gamma1; 1g and ng respectively. We
use another vague, but proper prior, for the coefficients of basis functions so
that we assume the a i - N(0; - 2 ) where we take the variance - These
priors may be chosen differently but this formulation is used to let the data
dictate the form of the model and leads to a proper posterior distribution as
all the priors are themselves proper. In particular, the prior on the type of the
basis functions could be chosen more carefully so that, a priori, main effects
are favoured over interaction terms. This could be useful when interpretation
of the model, rather than prediction, is more important. This refinement is
used in Mallick et al. (1997, 1998) which concentrate on situations where
interpretability of the model is of great interest.
A Poisson distribution (with parameter -) is used to specify the prior
probabilities for the number of basis functions, giving
In practice, a Poisson distribution truncated to k ! k max , for a suitable
choice of k max , is adopted. However, to help the sampler to mix better and
to limit prior influence we put a gamma hyperprior (with both parameters
equal to 10) over -. This reflects knowledge that we expect just a few basis
functions will fit the data well which controls overfitting.
2.3 Computational Strategy
Our aim is to simulate samples from the joint posterior distribution of p('
since analytic or numerical analyses are totally intractable in this situation.
For this purpose we design a reversible jump algorithm of the general type
discussed by Green (1995), to which the reader is referred for details.
In the context of our problem, with multiple parameter subspaces of different
dimensionality, it will be necessary to devise different types of moves
between the subspaces. These will be combined to form what Tierney (1994)
calls a hybrid sampler, making random choice between available moves at
each transition, in order to traverse freely around the combined parameter
space.
We use the following move types: (a) a change in a knot location; (b) the
addition of a basis function; (c) the deletion of a basis function. Note that in
steps (b) and (c) we are changing the dimension of the model and that we do
not add basis functions in pairs as in the standard MARS forward-stepwise
procedure: in fact, we depart completely from the any sort of recursive partitioning
approach. We have found that adding basis functions singly makes
our procedure more flexible and the reversibility condition easier to satisfy.
When we change the MARS structures, as described below, the coefficients
of the basis functions a in the new model have little
relationship to those in the current model so inference about them is difficult
and can lead to labelling problems, as in Richardson and Green (1997). Instead
we choose to integrate out the coefficients a i and the error variance oe 2
from the parameter vector ' (k) , which now only contains the model param-
eters. This is straightforward because we chose to use conjugate priors for
both the coefficients and the error variance. However, to make predictions
using the generated sample of models we draw coefficients for each model
in the sample from their full conditional distributions given the other model
parameters.
Given the current model step (a) is straight-forward. First we pick a basis
uniformly at random and then we pick one of the
factors we alter the current knot location t ji and with
probability 1reverse the sign indicator. We choose a new knot location from
the marginal predictor values of variable x v(j;i) and set this to the new t ji .
This move type is then undertaken using a Metropolis step (Metropolis et
al., 1953; Hastings, 1970) to accept or reject the proposed new state.
The addition of a basis function (BIRTH), step (b), is carried out by
choosing uniformly a type of basis function, say T i , to add to the model.
Then we uniformly choose a knot location and sign indicator for each of the
factors in this new basis.
Step (c) (DEATH), the deletion of a basis function is constructed in such
a way as to make the jump step reversible. This is easily done by choosing
a basis function uniformly from those present (except the constant basis
removing it.
At the end of each iteration after the move step has been performed we
use Gibbs steps (Gelfand and Smith, 1990) to generate a new -. This is
straightforward as its full conditional is simple to calculate and is given by
at each full cycle of the algorithm we obtain
a sample of (k; ' (k) ).
2.4 Algorithm
In the reversible jump algorithm we use the three move types described above
so that we can write the set of moves as refers to
changing a knot location and refers to increasing the number of
terminal nodes from m to m+1 or decreasing it from m+1 to m. Independent
move types are randomly chosen with probabilities ae
k and d k for k. In this problem
we took b
with the constant c, a parameter of the sampler, taken to be 0.4.
For
We find that by marginalising over the coefficients and the error variance
the acceptance probability given in Green (1995) simplifies to
where ' denotes the current model parameters, ' 0 the proposed model pa-
rameters, the probability of proposing a move to ' 0 from ' and D
the data. Thus the acceptance probability is just a Bayes factor (Kass and
multiplied by prior and proposal odds terms. This approach
is common in many fixed dimensional parameter inference problems and has
recently been used in other contexts where there are an unknown number of
parameters (Chipman et al., 1998; Holmes and Mallick, 1997).
The use of conjugate prior distributions allows simple evaluation of the
integrals in the Bayes factor term as demonstrated in O'Hagan (1994). This
leads to the Bayes factor being given by
d
where 0 refers to the proposed model, I is the identity matrix,
a a. Note that X and Y refer to the usual design
and data matrices of the current regression with b a
being the Bayesian least squares regression estimates of the coefficients. The
constants are the parameters of the gamma prior distribution over
oe \Gamma2 and so, from Section 2.2, were both taken to be 0.01. Note that an
ordinary least squares approach to estimating the regression coefficients was
undertaken in Denison et al. (1998a,b) and this is equivalent to the above
method when using reference priors for the coefficients and error variance.
We now show how the prior and proposal odds terms are calculated using
the birth step as an example. A birth step (b) adds basis B k+1 when there
are currently k basis functions in the model. The prior odds are given by
prior for the functions and ' (k+1)
prior for the k basis functions and ' (k)
where the terms in the numerator of (6) are given by the prior for the number
of basis functions, the prior for the type of the basis functions and the priors
for the knot positions and sign indicators: similarly for the denominator.
The prior probability for the set of bases, when there are k in the model,
can be thought of as the probability of choosing items from a set of
N where the ordering does not matter. The constant basis function is fixed
so it does not need to be chosen. The prior for the model parameters is the
product of the probability of each factor having a certain sign indicator (i.e.2
) and a certain knot point
) to the power of the number of factor terms
in the model (
The corresponding proposal odds is given by
p(propose death '
p(propose birth '
d k+1 =k
where we propose a death by randomly choosing one of the basis functions
not including the constant one, and a birth by randomly choosing a type of
basis function, with probability 1=N , together with a sign indicator and knot
point for each factor in the new basis fN(2n)g \GammaJ k+1 . Hence, using equations
(4-5) and (7-8), we can find the acceptance probability for a birth step. The
acceptance probability for a death step is worked out similarly.
The algorithm we use is straightforward and works quickly. The BIRTH
step is described in detail in the appendix and the other steps are very similar.
Algorithm
1. Start with just the constant basis function present.
2. Set k equal to the number of basis function in the current structure.
3. Generate u uniformly on [0,1].
4. Goto move type determined by u.
ffl else if (b
ffl else goto CHANGE step.
5. Draw - using Gibbs steps.
6. Repeat 2 for a suitable number of iterations once there is evidence of
convergence.
3 Simulated Examples
3.1 Bivariate Predictors
We first of all test our methodology on the examples given by Hwang et al.
(1994) and studied by Roosen and Hastie (1994). Following their approaches
we generate 225 pairs of predictors uniformly on the unit square and the
response is f(x is the true value of the test
function and the ffl i are drawn from a N(0; 0:25 2 ) distribution. The test
functions are
ffl Simple interaction function
ffl Radial function
f (2)
ffl Harmonic function
ffl Additive function
ffl Complicated interaction function
These functions are scaled and translated to have a standard deviation of
one and a non-negative range. We use the fraction of variance unexplained
to test the fits of the models to the data, given by
where f(x i ) is the true value of the function, b
is the fitted value and -
f
is the mean of the true values. We use FVU as it is helpful in comparing fits
of the model with differently generated datasets. To evaluate the FVU for a
fit we replace the expectations by averages over a test set of 10,000 points.
For these bivariate examples this is simply a 100 by 100 grid on the unit
square ie (1=200; 199=200). For the
higher dimensional examples which follows the test set is composed of 10,000
random uniform values over the domain of interest D. For the training set
we calculate the FVU over the training sample treating the observed y values
as f(x).
We took the results from the last 30,000 iterations of the sampler after
an initial burn-in period which was long enough for convergence to have
occurred by the end of it. Convergence was assumed when the mean squared
error of the fit had been settled for some time, as in the similar curve fitting
algorithm of Denison et al. (1998b).
In

Table

1 we display the FVU for the training set of data and the test
set for the standard MARS algorithm together with the number of basis
functions it found. For the BMARS model we give the average FVU for the
posterior mean model (obtained by pointwise averaging) from 10 runs of the
algorithm and we display the average number of basis functions in the samples
produced by the 10 repetitions of the algorithm. The standard MARS models
are referred to as LMARS (piecewise-linear MARS) and CMARS (piecewise-
cubic MARS) and we give results for both these MARS models. Table 1
shows that the BMARS model gives comparable, and often better, results
than both LMARS and CMARS for this wide variety of examples. Also, the
average number of basis functions found by BMARS is commonly less than
the number found using standard MARS.
The true surfaces for these examples are shown in Fig. 1 and the corresponding
BMARS estimates are given in Fig. 2.
Figs. 1 and 2 about here

Table

about here
3.2 High Dimensional Predictors
We take this example from Friedman et al. (1983). The basic function is
Following Friedman et al. we generate 200 random uniform predictors from
the six-dimensional unit hypercube and take the response to be f(x
where the ffl i are independent and identically distributed N(0; 1) random
variables. The extra predictor is spurious and does not affect the response.
Friedman (1991) also uses this function but in this paper a ten-dimensional
predictor is used (five being spurious) and the sample size is reduced to 100.
As is commonly the case in Friedman (1991) we do not allow more than
two-way interactions in the MARS models we use. Higher-order interaction
terms do not generally improve the fit and make the model unnecessarily
complex even though they could be incorporated if required.
In

Table

2 we display the FVU for the training and test set for these examples
using the standard MARS and BMARS methods. As before we took
the results from the last 30,000 iterations of 10 runs of the algorithm after
an initial burn-in period. The results shown demonstrate how the BMARS
model parsimoniously fits the data when compared to the standard MARS
fits.

Table

about here
4 Real Data Example
We now illustrate our methodolgy using a real dataset. We use data from a
study by Bruntz et al. (1974) of the dependence of ozone on some meteorological
variables on 111 days from May to September 1973 at sites in the New
York metropolitan area. As in Cleveland et al. (1988) we work with the cube
root of ozone. This dataset is known as air and is available in Splus (Becker
et al., 1988). There are three predictor variables, radiation, temperature and
wind speed but because of the vastly differing ranges of the response and
predictors we standardised the data beforehand, i.e. we linearly transformed
each variable so that it had zero mean and unit variance. Again we allow
only main-effect and two-way interaction terms in the MARS models.
The MARS fit had 6 basis functions and the residual sum of squares
(RSS) given by
was 18.41 with the linear approximation and 20.19 with the cubic one. The
final model was of the form f(wind)+f(temperature)+f(wind; temperature)+
f(radiation; temperature) with one basis function for each of these terms except
the final one for which two basis functions existed.
Over 5 runs of the algorithm, using the same priors as in the previous
section, the average RSS given by the BMARS model was 18.32 with 4.06
basis functions. This is lower than both the RSS's given by the MARS model
with fewer basis functions. As shown by Fig. 3 the BMARS estimate in the
(temperature, wind) plane is smooth whereas the piecewise-linear estimate
using MARS is not (Fig. 4). The BMARS estimate also has a smaller RSS
with less degrees of freedom.
In

Table

3 we display the estimated posterior probabilities of the possible
terms in the models, the RSS and the average number of basis functions in
each of the 5 runs. Immediately it can be seen that they each produced
similar results which suggests that convergence had occurred by the end of
the burn-in period. Also, it appears that the most important basis functions
in the fit were the main effect terms for radiation and temperature and the
(wind,temperature) interaction term. This is borne out by the fact that the
model which included only the terms just referred to had easily the largest
posterior probability and made up well over 50% of the generated samples.
The use of BMARS as a tool for performing a stochastic search for variable
selection comes "for free" when performing the analysis for prediction.
Variable selection, itself, is a difficult problem and one that has received
much attention in the literature. The BMARS method could be used in
a similar way to the Gibbs sampling-based method outlined in George and
McCulloch (1993) which is shown to identify good models using a stochastic
search procedure.
Figs. 3 and 4 about here

Table

3 about here
We have presented a Bayesian approach to finding regression surfaces which
uses truncated linear basis functions to build up the surface. We use the
data to find the knot points, the main effect and/or interaction terms and
the number of basis functions that are required to adequately approximate
the required surface. We simulate a random sample of models using the
reversible jump MCMC approach of Green (1995).
This Bayesian approach to multiple regression produces a model with high
predictive power due to the posterior averaging over all the models. This not
only leads to a good overall fit for the data but can also help to combat
overfitting problems. We can however choose a single model as our estimate
by picking out single models from the sample by various criteria. These
models have less predictive power but have more interpretability as we can
easily produce their ANOVA decompositions. We can only produce expected
numbers and posterior probabilities of basis functions for the posterior mean
estimate but this too can be used to identify variables that are important to
the fit.
In some applications we may not wish to draw the coefficients but use the
mean of their sampling distributions instead. If identifying good models, and
not prediction, is the main aim then not drawing the coefficients is slightly
quicker and produces models with smaller (squared) error. This may be
useful in the related BMARS algorithms for analysing financial time series
(Denison, 1997) and failure time data (Mallick et al., 1997). This is the
approach Denison et al. (1998a) undertake to find "good" CART models via
a stochastic search as in this example prediction using the posterior mean is
not helpful.
The structure of the algorithm which allows the sampler to move fluidly
about the probability space means that Bayesian MARS does not inherit the
problems of Bayesian CART algorithms (Denison et al., 1998a; Chipman et
al., 1998) by only searching restricted portions of the entire space. This is
greatly helped by not "splitting" on current basis functions thus inducing no
hierarchical structure to the form of the bases. This point is demonstrated
in full in Chapter 5 of Denison (1997). If the space of predictor variables is
substantial (over 50, say) then convergence of the sampler would seem to be
unlikely but still the BMARS approach should yield "good" models.
The BMARS algorithm we used in this paper was written in ANSI C and
is available, together with the datasets, from the World Wide Wed address
http://ma.ic.ac.uk/-dgtd/. The algorithm takes around 5 minutes to
run on a DEC Alpha workstation when there are 200 datapoints.

Acknowledgements

The work of the first author was supported by an EPSRC research stu-
dentship. We acknowledge the helpful comments made by the anonymous
referees, the associate editor and Mr C.C. Holmes.


Appendix


The move types CHANGE, BIRTH and DEATH given in the algorithm in
Section 2.3 are undertaken similarly so we just describe the BIRTH step in
pseudo-code.The notation follows that in Section 2.
1. Uniformly choose a type of basis function T i to add from the N possible
ones.
2. Uniformly choose the knot positions, predictors to split on and the sign
indicators in this new basis remembering that each predictor may only
occur once in each basis function.
3. Generate u uniformly on [0,1].
4. Work out the acceptance probability, ff.
5. IF accept the proposed model.
ELSE keep the current model.
6. Return to main algorithm.



--R

The New S Language
Estimating optimal transformations for multiple regression and correlation (with discussion).
Classification and Regression Trees
The dependence of ambient ozone on solar radiation

Bayesian CART model search (with discussion).

Smoothing noisy data with spline func- tions
Simulation based Bayesian nonparametric regression methods.


Multivariate adaptive regression splines (with dis- cussion)
Multidimensional additive spline approximation.
Projection pursuit regression.

Variable selection via Gibbs sam- pling
Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.
Generalized additive models
Monte Carlo sampling methods using Markov chains and their applications.
Bayesian wavelet networks for non-parametric regression
Regression modeling in back-propagation and projection pursuit learn- ing
Bayes factors.


In Generalized Linear Models: A Bayesian Perpective
Equations of state calculations by fast computing machines.
Problems in the analysis of survey data and a proposal.

Bayesian analysis of mixtures with an unknown number of components (with discussion).
Automatic smoothing spline projection pursuit.
Markov chains for exploring posterior distributions (with discussion).

Figure 2 Hwang
Figure 4 Linear-piecewise (left) and cubic-piecewise (right) MARS estimates in the (temperature

--TR

--CTR
C. C. Holmes , B. K. Mallick, Bayesian radial basis functions of variable dimension, Neural Computation, v.10 n.5, p.1217-1233, July 1, 1998
David J. Nott , Anthony Y. Kuk , Hiep Duc, Efficient sampling schemes for Bayesian MARS models with many predictors, Statistics and Computing, v.15 n.2, p.93-101, April     2005
Christophe Andrieu , Nando De Freitas , Arnaud Doucet, Robust Full Bayesian Learning for Radial Basis Networks, Neural Computation, v.13 n.10, p.2359-2407, October 2001
multivariate adaptive regression splines;reversible jump Markov Chain Monte Carlo;bayesian methods;multiple regression