Take in a prior and likelihood distribution, and output a posterior distribution.
This function first checks a table of conjugate distributions for the pair you sent in. If the names match the table, then the function returns a closed-form model with updated parameters. If the parameters aren't in the table of conjugate priors/likelihoods, then it uses Markov Chain Monte Carlo to sample from the posterior distribution, and then outputs a histogram model for further analysis. Notably, the histogram can be used as the input to this function, so you can chain Bayesian updating procedures.
- If the prior distribution has a
p
or log_likelihood
element, then I use apop_model_metropolis to generate the posterior.
- If the prior does not have a
p
or log_likelihood
but does have a draw
element, then I make draws from the prior and weight them by the p
given by the likelihood distribution. This is not a rejection sampling method, so the burnin is ignored.
Here are the conjugate distributions currently defined:
- The conjugate table is stored using a vtable; see Registering new methods in vtables for details. The typedef new functions must conform to and the hash used for lookups are:
1 typedef apop_model *(*apop_update_type)(apop_data *, apop_model , apop_model);
2 #define apop_update_hash(m1, m2) ((size_t)(m1).draw + (size_t)((m2).log_likelihood ? (m2).log_likelihood : (m2).p)*33)
- Parameters
-
data | The input data, that will be used by the likelihood function (default = NULL .) |
prior | The prior apop_model. If the system needs to estimate the posterior via MCMC, this needs to have a log_likelihood or p method. (No default, must not be NULL .) |
likelihood | The likelihood apop_model. If the system needs to estimate the posterior via MCMC, this needs to have a log_likelihood or p method (ll preferred). (No default, must not be NULL .) |
rng | A gsl_rng , already initialized (e.g., via apop_rng_alloc). (default: an RNG from apop_rng_get_thread) |
- Returns
- an apop_model struct representing the posterior, with updated parameters.
Here is a test function that compares the output via conjugate table and via Metropolis-Hastings sampling:
void distances(gsl_vector *v1, gsl_vector *v2, double tol){
Apop_stopif(error/updated_size > tol, exit(1), 0,
"The error is %g, which is too big.", error/updated_size);
}
int main(){
double binom_start = 0.6;
double beta_start_a = 0.3;
double beta_start_b = 0.5;
double n = 4000;
int i, draws = 1.3e5;
for(i=0; i < draws; i ++)
}