![]() |
|
Macros | |
#define | markit(test, value) |
Functions | |
int | sigma_adapt (apop_mcmc_proposal_s *ps, apop_mcmc_settings *ms) |
Apop_settings_init (Apop_settings_copy(apop_mcmc, Apop_varad_set(periods, 6e3);Apop_varad_set(burnin, 0.05);Apop_varad_set(target_accept_rate, 0.35);Apop_varad_set(gibbs_chunks, 'b');Apop_varad_set(start_at, '1');Apop_varad_set(base_step_fn, step_to_vector);Apop_varad_set(base_adapt_fn, sigma_adapt);) | |
Apop_settings_free (apop_mcmc, if(in->proposal_is_cp){for(int i=0;i< in->block_count;i++) apop_model_free(in->proposals[i].proposal);free(in->proposals);}) | |
apop_model * | maybe_prep (apop_data *d, apop_model *m, bool *is_a_copy) |
int | apop_model_metropolis_draw (double *out, gsl_rng *rng, apop_model *model) |
void | main_mcmc_loop (apop_data *d, apop_model *m, apop_data *out, gsl_vector *draw, apop_mcmc_settings *s, gsl_rng *rng, int *constraint_fails) |
apop_model * | apop_model_metropolis (apop_data *d, gsl_rng *rng, apop_model *m) |
Markov Chain Monte Carlo.
#define markit | ( | test, | |
value | |||
) |
apop_model* apop_model_metropolis | ( | apop_data * | d, |
gsl_rng * | rng, | ||
apop_model * | m | ||
) |
Use Metropolis-Hastings Markov chain Monte Carlo to make draws from the given model.
The basic storyline is that draws are made from a proposal distribution, and the likelihood of your model given your data and the drawn parameters evaluated. At each step, a new set of proposal parameters are drawn, and if either they are more likely than the previous set the new proposal is accepted as the next step, else with probability (prob of new params)/(prob of old params), they are accepted as the next step anyway. Otherwise the last accepted proposal is repeated.
The output is an apop_pmf model with a data set listing the draws that were accepted, including those repetitions. The output model is modified so that subsequent draws are one more step from the Markov chain, via apop_model_metropolis_draw.
constraint
element of the model you input, then the proposal is thrown out and a new one selected. By the default proposal distribution, this is not mathematically correct (it breaks detailed balance), and values near the constraint will be oversampled. The output model will have outmodel->error=='c'
. It is up to you to decide whether the resulting distribution is good enough for your purposes or whether to take the time to write a custom proposal and step function to accommodate the constraint.Attach an apop_mcmc_settings group to your model to specify the proposal distribution, burnin, and other details of the search. See the apop_mcmc_settings documentation for details.
base_adapt_fn
in the apop_mcmc_settings group to a do-nothing function, or one that damps its adaptation as gibbs_chunks
element of the apop_mcmc_settings group. If you set gibbs_chunks='a'
, all parameters are drawn as a set, and accepted/rejected as a set. The variances are adapted at an identical rate. If you set gibbs_chunks='i'
, then each scalar parameter is assigned its own proposal distribution, which is adapted at its own pace. With gibbs_chunks='b'
(the default), then each of the vector, matrix, and weights of your model's parameters are drawn/accepted/adapted as a group (and so on to additional chunks if your model has ->more
pages). This works well for complex models which naturally break down into subsets of parameters.Each chunk counts as a step in the Markov chain. Therefore, if there are several chunks, you can expect chunks to repeat from step to step. If you want a draw after cycling through all chunks, try using apop_model_metropolis_draw, which has that behavior.
d | The apop_data set used for evaluating the likelihood of a proposed parameter set. |
rng | A gsl_rng , probably allocated via apop_rng_alloc. (Default: an RNG from apop_rng_get_thread) |
m | The apop_model from which parameters are being drawn. (No default; must not be NULL ) |
parameters
element of your likelihood model has the last accepted parameter proposal.apop_opts.verbose=2
or greater, I will report the accept rate of the M-H sampler. It is a common rule of thumb to select a proposal so that this is between 20% and 50%. Set apop_opts.verbose=3
to see the stream of proposal points, their likelihoods, and the acceptance odds. You may want to set apop_opts.log_file=fopen("yourlog", "w")
first.draw
method that returns another step from the Markov chain with each draw.out->error='c' | Proposal was outside of a constraint; see above. |
int apop_model_metropolis_draw | ( | double * | out, |
gsl_rng * | rng, | ||
apop_model * | model | ||
) |
The draw method for models estimated via apop_model_metropolis.
That method produces an apop_pmf, typically with a few thousand draws from the model in a batch. If you want to get a single next step from the Markov chain, use this.
A Markov chain works by making a new draw and then accepting or rejecting the draw. If the draw is rejected, the last value is reported as the next step in the chain. Users sometimes mitigate this repetition by making a batch of draws (say, ten at a time) and using only the last.
If you run this without first running apop_model_metropolis, I will run it for you, meaning that there will be an initial burn-in period before the first draw that can be reported to you. That run is done using model->data
as input.
out | An array of doubles , which will hold the draw, in the style of apop_draw. |
rng | A gsl_rng , already initialized, probably via apop_rng_alloc. |
model | A model which was probably already run through apop_model_metropolis. |
out
is filled with the next step in the Markov chain. The ->data
element of the PMF model is extended to include the additional steps in the chain. If a proposal failed the model constraints, then return 1; else return 0. See the notes in the documentation for apop_model_metropolis.base_model
in the mcmc settings group == the parent model.