A simple example based on estimating the parameters of a linear regression model using
* Data annealing sequential Monte Carlo (LinReg).
* Likelihood annealing sequential Monte Carlo (LinRegLA).
* Likelihood annealing sequential Monte Carlo with the temperature schedule,
number of MCMC repeats and random walk covariance matrices adapted online (LinRegLA_adapt).
LinReg(model, particles = 1000, plot = FALSE)LinRegLA(model, particles = 1000, temperatures = seq(0, 1, 0.05)^5)
LinRegLA_adapt(model, particles = 1000, resampTol = 0.5, tempTol = 0.9)
The LinReg function returns a list containing the final particle
approximation to the target (\(\theta\) and the corresponding weights) as well as the logarithm
of the estimated model evidence.
The LinRegLA function returns a list containing the population
of particles and their associates log likelihoods, log priors and weights at each iteration. The
effective sample size at each of the iterations and several different
estimates of the logarithm of the model evidence are also returned.
The LinRegLA_adapt function returns a list containing all of the same
output as LinRegLA, in addition to the adaptively chosen temperature schedule
and number of MCMC repeats.
Choice of regression model (1 for density as the predictor and 2 for adjusted density as the predictor).
An integer specifying the number of particles.
A boolean variable to determine whether to plot the posterior estimates.
In likelihood annealing SMC the targets are defined as \(P(y|\theta)^{\gamma_t}P(\theta)\) where \(0=\gamma_0\le \ldots \le \gamma_T = 1\) can be referred to as the temperatures, \(P(y|\theta)\) is the likelihood and \(P(\theta)\) is the prior.
The adaptive implementation of likelihood annealing SMC allows for the resampling tolerance to be specified. This parameter can be set to a value in the range [0,1) corresponding to a fraction of the size of the particle set or it may be set to an integer corresponding to an actual effective sample size.
A tolerance for adaptive choice of the temperature schedule such that the conditional ESS is maintained at tempTol*particles.
Adam M. Johansen, Dirk Eddelbuettel and Leah F. South
Williams (1959) considers two competing linear regression models for the maximum compression strength parallel to the grain for radiata pine. Both models are of the form
\(y_i = \alpha + \beta (x_i - \bar{x}) + \epsilon_i\),
where \(\epsilon_i ~ N(0,\sigma^2)\) and \(i=1,\ldots,42\). Here \(y\) is the maximum compression strength in pounds per square inch. The density (in pounds per cubic foot) of the radiata pine is considered a useful predictor, so model 1 uses density for \(x\). Model 2 instead considers the density adjusted for resin content, which is associated with high density but not with strength.
This example is frequently used as a test problem in model choice (see for example Carlin and Chib (1995) and Friel and Pettitt (2008)). We use the standard uninformative normal and inverse gamma priors for this example along with the transformation \(\phi=log(\sigma^2)\) so that all parameters are on the real line and \(\theta = [\alpha,\beta,\phi]\). The evidence can be computed using numerical estimation for both of the competing models. The log evidence is -309.9 for model 1 and -301.4 for model 2.
The LinReg function implements a data annealing approach to this example.
The LinRegLA function implements a likelihood annealing approach to this example.
The LinRegLA_adapt function implements a likelihood annealing approach to this example
with adaptation of the temperature schedule, number of MCMC repeats and random walk covariance
matrices.
B. P. Carlin and S. Chib. Bayesian model choice via Markov chain Monte Carlo. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 57(3):473-484, 1995.
N. Friel and A. N. Pettitt. Marginal likelihood estimation via power posteriors. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 70(3):589-607, 2008.
Williams, E. (1959). Regression analysis. Wiley.
res <- LinReg(model=1, particles=1000, plot=TRUE)
res <- LinRegLA(model=1, particles=1000)
res <- LinRegLA_adapt(model=1, particles=1000)
Run the code above in your browser using DataLab