── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Making sure initially that we can plot the mystery function:
xs = seq(-30,30, 0.01)
data = tibble(x = xs, y = mystery(xs))
ggplot(data=data, aes(x=x, y=y)) + geom_point()
Setting up our sampler:
mcmc_sampler <- function(func, start, step_size, num_steps) {
samples = numeric(num_steps)
samples[1] = start
accepted = 0
for (i in 2:num_steps) {
current = samples[i-1]
proposed = rnorm(1, current, step_size)
ratio = func(proposed) / func(current)
if (ratio > runif(1)) {
samples[i] = proposed
accepted = accepted + 1
} else {
samples[i] = current
}
}
return(list(samples=samples, accept_rate=accepted / length(samples)))
}
Running the sampler
out = mcmc_sampler(mystery, -10, 4, 50000)
samples = out$samples
accept_rate = out$accept_rate
accept_rate
Looking at the trace plot:
data = tibble(iter=1:length(samples), x=samples)
ggplot(data=data, aes(x=iter, y=x)) + geom_line()
Looking at the lag plot:
burned = data %>% filter(iter > 1000)
ggplot(data=burned, aes(x=x, y=lag(x))) + geom_point(alpha=0.1)
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).
And looking at the overall histogram:
ggplot(data=burned, aes(x=x)) + geom_histogram(bins=100)
