#include "MCMCLib.h" // Define the log likelihood function double log_likelihood(const gsl_vector *x, void *params) { // Extract the parameters double mu = gsl_vector_get(x, 0); double sigma = gsl_vector_get(x, 1); // Compute the log likelihood double logL = -0.5 * log(2 * M_PI * sigma * sigma) - 0.5 * (y - mu) * (y - mu) / (sigma * sigma); return logL; } // Define the function to perform MCMC sampling void mcmc_sampling() { // Set the initial values for the parameters gsl_vector *x = gsl_vector_alloc(2); gsl_vector_set(x, 0, mu0); gsl_vector_set(x, 1, sigma0); // Set the proposal distribution gsl_matrix *cov = gsl_matrix_alloc(2, 2); gsl_matrix_set_identity(cov); gsl_matrix_scale(cov, sigma_prop); // Set the MCMC options MCMCOptions options; options.n_burnin = 1000; options.n_iter = 10000; options.thinning = 10; // Set the MCMC data MCMCData data; data.n_params = 2; data.logL = log_likelihood; data.params = NULL; data.x = x; data.cov = cov; // Perform the MCMC sampling MCMCSample *sample = MCMCSample_alloc(options.n_iter); MCMC(options, data, sample); // Extract the samples for (int i = 0; i < sample->n; i++) { mu_samples[i] = gsl_vector_get(sample->x[i], 0); sigma_samples[i] = gsl_vector_get(sample->x[i], 1); } // Free the allocated memory MCMCSample_free(sample); gsl_vector_free(x); gsl_matrix_free(cov); }
Preview:
downloadDownload PNG
downloadDownload JPEG
downloadDownload SVG
Tip: You can change the style, width & colours of the snippet with the inspect tool before clicking Download!
Click to optimize width for Twitter