Patterns in static

Apophenia

Macros | Functions
apop_mcmc.c File Reference

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_modelmaybe_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_modelapop_model_metropolis (apop_data *d, gsl_rng *rng, apop_model *m)
 

Detailed Description

Markov Chain Monte Carlo.

Macro Definition Documentation

#define markit (   test,
  value 
)
Value:
if (test) \
s->block_starts[this++] = ctr += value;

Function Documentation

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.

  • If a proposal fails to meet the 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.

  • The default proposal includes an adaptive step: you specify a target accept rate (default: .35), and if the accept rate is currently higher the variance of the proposals is widened to explore more of the space; if the accept rate is currently lower the variance is narrowed to stay closer to the last accepted proposal. Technically, this breaks ergodicity of the Markov chain, but the consensus seems to be that this is not a serious problem. If it does concern you, you can set the base_adapt_fn in the apop_mcmc_settings group to a do-nothing function, or one that damps its adaptation as $n\to\infty$.
  • Note the 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.

Parameters
dThe apop_data set used for evaluating the likelihood of a proposed parameter set.
rngA gsl_rng, probably allocated via apop_rng_alloc. (Default: an RNG from apop_rng_get_thread)
mThe apop_model from which parameters are being drawn. (No default; must not be NULL)
  • If the likelihood model no parameters, I will allocate them. That means you can use one of the stock models that ship with Apophenia. If I need to run the model's prep routine to get the size of the parameters, then I'll make a copy of the likelihood model, run prep, and then allocate parameters for that copy of a model.
  • On exit, the parameters element of your likelihood model has the last accepted parameter proposal.
  • If you set 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.
Returns
A modified apop_pmf model representing the results of the search. It has a specialized draw method that returns another step from the Markov chain with each draw.
Exceptions
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.

Parameters
outAn array of doubles, which will hold the draw, in the style of apop_draw.
rngA gsl_rng, already initialized, probably via apop_rng_alloc.
modelA model which was probably already run through apop_model_metropolis.
Returns
On return, 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.
  • After pulling the attached settings group, the parent model is ignored. One expects that base_model in the mcmc settings group == the parent model.
  • If your settings break the model parameters into several chunks, this function returns after stepping through all chunks.

Autogenerated by doxygen on Wed Oct 15 2014 (Debian 0.999b+ds3-2).