#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);
}