Normal model

Bob Verity and Pete Winskill

2024-06-26

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

Purpose: to compare drjacoby infered posterior distribution against the exact analytical solution for a very simple model.

Model

A vector of data x drawn from a normal distribution with unknown mean mu and known variance of 1.0.

n <- 1e2
data_list <- list(x = rnorm(n))

Likelihood and prior:

Parameters dataframe:

df_params <- define_params(name = "mu", min = -1, max = 1)

MCMC

mcmc <- run_mcmc(data = data_list,
                 df_params = df_params,
                 loglike = "loglike",
                 logprior = "logprior",
                 burnin = 1e3,
                 samples = 1e5,
                 chains = 10,
                 silent = TRUE)

Posterior plots

# extract sampling draws
mu_draws <- subset(mcmc$output, phase = "sampling")$mu

## Warning: In subset.data.frame(mcmc$output, phase = "sampling") :
##  extra argument 'phase' will be disregarded

# calculate analytical solution
x <- seq(-1, 1, l = 1001)
fx <- dnorm(x, mean = mean(data_list$x), sd = sqrt(1/n))

# histogram and overlay analytical
hist(mu_draws, breaks = seq(-1, 1, 0.01), probability = TRUE, col = "black")
lines(x, fx, col = 2, lwd = 4)