Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Proof-of-concept parallel tempering #106

Draft
wants to merge 7 commits into
base: main
Choose a base branch
from
Draft

Proof-of-concept parallel tempering #106

wants to merge 7 commits into from

Conversation

richfitz
Copy link
Member

@richfitz richfitz commented Nov 8, 2024

Merge after #104, contains those commits

Closes #24

A basic parallel tempering sampler. There's much more to add (particularly embedding other samplers and sorting out sampling from the hot chain), but this is a start at least

Marc, a simple example to get you going. This actually looks pretty good now I think:

set.seed(1) # so we see the same thing
pkgload::load_all()
likelihood <- ex_mixture(5)
prior <- monty_dsl({
  x ~ Normal(0, 10)
})
posterior <- likelihood + prior
s <- monty_sampler_parallel_tempering(n_rungs = 10, vcv = matrix(0.1))
res <- monty_sample(posterior, s, 10000, n_chains = 4)

hist(c(res$pars), freq = FALSE)
x <- seq(min(res$pars), max(res$pars), length.out = 1001)
y <- exp(posterior$density(rbind(x)))
lines(x, y / sum(y) / diff(x)[[1]], col = "red")

(takes 17s on my machine for 40,000 steps, so 2k steps/s which is ok, especially as this is really 11 times that many samples).

I've run out a much longer chain and this is what I see:

set.seed(42) # so we see the same thing
pkgload::load_all()
likelihood <- ex_mixture(5)
prior <- monty_dsl({
  x ~ Normal(0, 10)
})
posterior <- likelihood + prior
s <- monty_sampler_parallel_tempering(n_rungs = 10, vcv = matrix(0.1))
r <- monty_runner_callr(dust2::dust_openmp_threads(10, "fix"))
res <- monty_sample(posterior, s, 50000, n_chains = 10, runner = r)

hist(c(res$pars), freq = FALSE)
x <- seq(min(res$pars), max(res$pars), length.out = 1001)
y <- exp(posterior$density(rbind(x)))
lines(x, y / sum(y) / diff(x)[[1]], col = "red")

(27s for 500,000 points, I can live with that)

image

I suspect we're out by a factor related to the beta swaps at the moment?

Copy link

codecov bot commented Nov 8, 2024

Codecov Report

Attention: Patch coverage is 99.00000% with 1 line in your changes missing coverage. Please review.

Project coverage is 99.64%. Comparing base (84df023) to head (027b66a).

Files with missing lines Patch % Lines
R/sampler-parallel-tempering.R 99.00% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main     #106      +/-   ##
==========================================
- Coverage   99.65%   99.64%   -0.02%     
==========================================
  Files          70       71       +1     
  Lines        5477     5577     +100     
==========================================
+ Hits         5458     5557      +99     
- Misses         19       20       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@richfitz richfitz force-pushed the mrc-5110 branch 2 times, most recently from 841914d to d9ecefb Compare November 8, 2024 10:30
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant