bat / BAT.jl

A Bayesian Analysis Toolkit in Julia
Other
198 stars 30 forks source link

unshape result_trafo and add initval_alg to MCMCChainPoolInit #369

Closed Cornelius-G closed 2 years ago

Cornelius-G commented 2 years ago

As discussed in #364 all samplers now return unshaped samples in the field result_trafo, also for DoNotTransform() This one should replace #364 (sorry for opening another PR, had to make a new branch)

@oschulz, there is something strange going on with the trafos in "mcmc_sample.jl" line 91 following. If trafo is something like identity ∘ identity broadcasting with inverse(trafo).(samples_trafo) fails and is not needed since inverse(trafo)(samples_trafo) works. But for the other trafos the broadcasting is needed (see stacktrace.txt).

Strangely, in "importance_sampler.jl" and "ellipsoidal_nested_sampling.jl" the same code is used but works with broadcasting for all trafos. It seems to me there is an issue if DensitySampleVector contains something in info. At least, that's the only difference in the samples for mcmc and the importance sampler. Do you have any idea what might be going on there?

codecov-commenter commented 2 years ago

Codecov Report

Merging #369 (4e7da10) into main (61ae26f) will increase coverage by 0.00%. The diff coverage is 94.73%.

@@           Coverage Diff           @@
##             main     #369   +/-   ##
=======================================
  Coverage   55.59%   55.60%           
=======================================
  Files         114      114           
  Lines        5439     5442    +3     
=======================================
+ Hits         3024     3026    +2     
- Misses       2415     2416    +1     
Impacted Files Coverage Δ
src/samplers/mcmc/mcmc_algorithm.jl 65.21% <0.00%> (-1.45%) :arrow_down:
src/samplers/mcmc/mh/mh_tuner.jl 94.11% <ø> (-1.48%) :arrow_down:
src/samplers/importance/importance_sampler.jl 100.00% <100.00%> (ø)
src/samplers/mcmc/chain_pool_init.jl 89.02% <100.00%> (+0.27%) :arrow_up:
src/samplers/mcmc/mcmc_sample.jl 100.00% <100.00%> (ø)
...dal_nested_sampling/ellipsoidal_nested_sampling.jl 100.00% <100.00%> (ø)
src/samplers/nested_sampling/ultranest.jl 97.14% <100.00%> (+0.08%) :arrow_up:
src/variates/density_sample.jl 85.47% <0.00%> (-0.86%) :arrow_down:
src/transforms/distribution_transform.jl 86.86% <0.00%> (+0.29%) :arrow_up:
... and 1 more

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 61ae26f...4e7da10. Read the comment docs.

oschulz commented 2 years ago

Do you have any idea what might be going on there?

No, that's weird. Can you try to make a minimal code snippet to reproduce?

Cornelius-G commented 2 years ago

No, that's weird. Can you try to make a minimal code snippet to reproduce?

Sure. This is the same example posterior as in "test/densities/test_truncate_density.jl"

using BAT, Distributions, InverseFunctions, DensityInterface, ValueShapes

mvec = [-0.3, 0.3]
cmat = [1.0 1.5; 1.5 4.0]
mv_dist = MvNormal(mvec, cmat)

likelihood = logfuncdensity(logdensityof(BAT.DistMeasure(mv_dist)))
prior = product_distribution(Uniform.([-5, -8], [5, 8]))

posterior = PosteriorMeasure(likelihood, prior)

trafo_alg = DoNotTransform()
density, trafo = BAT.transform_and_unshape(trafo_alg, posterior)

# for importance & nested samples it works:
samples_sobol = bat_sample(posterior, SobolSampler(nsamples=10^2, trafo=trafo_alg))
samples_notrafo_sobol = inverse(trafo).(samples_sobol.result_trafo) #this works

# for mcmc samples it fails:
samples_mh = bat_sample(posterior, MCMCSampling(mcalg = MetropolisHastings(), trafo = trafo_alg, nchains = 2, nsteps = 10^2))
samples_notrafo_mh = inverse(trafo).(samples_mh.result_trafo) #this fails
oschulz commented 2 years ago

This should fix it, I hope: #370

Cornelius-G commented 2 years ago

@oschulz Thanks for the fix of the broadcast. Now it works as expected. The one CI test failing is a stochastic thing and ok, I guess:

RandResampling: Test Failed at /Users/runner/work/BAT.jl/BAT.jl/test/samplers/test_bat_sample.jl:44 Expression: isapprox(mean(samples_rdm), mean(result), rtol = 10 ^ -1) Evaluated: isapprox(-0.06715567259585313, -0.06032598807659495; rtol = 0.1)

So from my side, the PR would be ready to merge.

However, for using this to start sampling in/very close to the mode, we found that we need to decrease tuner.scale in tuning_init! in line 96 of "mh_tuner.jl" since otherwise, the chain won't start running as the initial proposal covariance is too wide and no proposed point will be accepted and the chain will be marked "non-viable" before the tuning cycles are even started. Currently, there is no option for changing this initial tuner scale parameter (it is fixed to scale = 2.38^2 / m). But I would make another PR for that if we want to provide an interface to change this scale.

oschulz commented 2 years ago

However, for using this to start sampling in/very close to the mode ... tuner scale ...

We need to make tuner scale adaption more dynamic anyway.

oschulz commented 2 years ago

Thanks!