Description Usage Arguments Details Note Author(s) Examples

Sequential Monte Carlo Sampler

1 2 3 | ```
smcSampler(bayesianSetup, initialParticles = 1000, iterations = 10,
resampling = T, resamplingSteps = 2, proposal = NULL,
adaptive = T, proposalScale = 0.5)
``` |

`bayesianSetup` |
either an object of class bayesianSetup created by |

`initialParticles` |
initial particles - either a draw from the prior, provided as a matrix with the single parameters as columns and each row being one particle (parameter vector), or a numeric value with the number of desired particles. In this case, the sampling option must be provided in the prior of the BayesianSetup. |

`iterations` |
number of iterations |

`resampling` |
if new particles should be created at each iteration |

`resamplingSteps` |
how many resampling (MCMC) steps between the iterations |

`proposal` |
optional proposal class |

`adaptive` |
should the covariance of the proposal be adapted during sampling |

`proposalScale` |
scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection |

The sampler can be used for rejection sampling as well as for sequential Monte Carlo. For the former case set the iterations to one.

The SMC currently assumes that the initial particle is sampled from the prior. If a better initial estimate of the posterior distribution is available, this the sampler should be modified to include this. Currently, however, this is not included in the code, so the appropriate adjustments have to be done by hand.

Florian Hartig

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 | ```
## Example for the use of SMC
# First we need a bayesianSetup - SMC makes most sense if we can for demonstration,
# we'll write a function that puts out the number of model calls
MultiNomialNoCor <- generateTestDensityMultiNormal(sigma = "no correlation")
parallelLL <- function(parMatrix){
print(paste("Calling likelihood with", nrow(parMatrix), "parameter combinations"))
out = apply(parMatrix, 1, MultiNomialNoCor)
return(out)
}
bayesianSetup <- createBayesianSetup(likelihood = parallelLL, lower = rep(-10, 3),
upper = rep(10, 3), parallel = "external")
# Defining settings for the sampler
# First we use the sampler for rejection sampling
settings <- list(initialParticles = 1000, iterations = 1, resampling = FALSE)
# Running the sampler
out1 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
#plot(out1)
# Now for sequential Monte Carlo
settings <- list(initialParticles = 100, iterations = 5, resamplingSteps = 1)
out2 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
#plot(out2)
## Not run:
## Example for starting a new SMC run with results from a previous SMC run
# Generate example data (time series)
# x1 and x2 are predictory, yObs is the response
t <- seq(1, 365)
x1 <- (sin( 1 / 160 * 2 * pi * t) + pi) * 5
x2 <- cos( 1 / 182.5 * 1.25 * pi * t) * 12
# the model
mod <- function(par, t1 = 1, tn = 365) {
par[1] * x1[t1:tn] + par[2] * x2[t1:tn]
}
# the true parameters
par1 <- 1.65
par2 <- 0.75
yObs <- mod(c(par1, par2)) + rnorm(length(x1), 0, 2)
# split the time series in half
plot(yObs ~ t)
abline(v = 182, col = "red", lty = 2)
# First half of the data
ll_1 <- function(x, sum = TRUE) {
out <- dnorm(mod(x, 1, 182) - yObs[1:182], 0, 2, log = TRUE)
if (sum == TRUE) sum(out) else out
}
# Fit the first half of the time series
# (e.g. fit the model to the data soon as you collect the data)
setup_1 <- createBayesianSetup(ll_1, lower = c(-10, -10), upper = c(10, 10))
settings_1 <- list(initialParticles = 1000)
out_1 <- runMCMC(setup_1, "SMC", settings_1)
summary(out_1)
# Second half of the data
ll_2 <- function(x, sum = TRUE) {
out <- dnorm(mod(x, 183, 365) - yObs[183:365], 0, 2, log = TRUE)
if (sum == TRUE) sum(out) else out
}
# Fit the second half of the time series
# (e.g. fit the model to the data soon as you collect the data)
setup_2 <- createBayesianSetup(ll_2, lower = c(-10, -10), upper = c(10, 10))
# This is the important step, we use the final particles from the
# previous SMC run to initialize the new SMC run
settings_2 <- list(initialParticles = out_1$particles)
out_2 <- runMCMC(setup_2, "SMC", settings_2)
summary(out_2)
par_pred <- apply(out_2$particles, 2, median)
pred <- mod(par_pred)
plotTimeSeries(yObs, pred)
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.