nimble-dev / nimble

The base NIMBLE package for R
http://R-nimble.org
BSD 3-Clause "New" or "Revised" License
160 stars 24 forks source link

bug with `adaptInterval=1` for `RW_block` #1493

Closed paciorek closed 6 days ago

paciorek commented 1 month ago

It runs (which I was thinking it might not because of 1-d vs 2-d allocation for the samples saved for calculating the covariance), but in a simple example, there are no acceptances. @danielturek or I should look into for 1.3.0 release.

code <- nimbleCode({
    y[1:2] ~ dmnorm(z[1:2],pr[1:2,1:2])
})
m <- nimbleModel(code, inits = list(z=rep(0,2), pr=diag(2)))
m$y = rnorm(2)

conf  <- configureMCMC(m, nodes=NULL)
conf$addSampler(target='y', type='RW_block', control = list(adaptInterval=1))
mcmc <- buildMCMC(conf)
cm <- compileNimble(m)
cmcmc <- compileNimble(mcmc,project=m)
cmcmc$run(1000)
tmp=as.matrix(cmcmc$mvSamples)
danielturek commented 1 month ago

@paciorek The issue is essentially, when RW_block sampler estimates the empirical covariance from the recent samples (specifically, using the samples collected during the most recent adaptation period), that estimation of a 2nd order statistic (variance-covariance matrix) will fail when there's only 1 sample in each adaptation period (when adaptInterval = 1).

Looking the the code for the RW_block sampler, this problem manifests itself around line 745 or so, when first the "centering" of the empircal samples (here, a single sample), centers them to become all zero:

empirSamp[, i] <<- empirSamp[, i] - mean(empirSamp[, i])

and then worse, when it divides by $n-1$, or here timesRan-1, that's dividing by 0, which seems to give NaN, or something similar, and then the chol fails, and I guess it's downhill from there.

Maybe we could restrict adaptInterval to be at least 2? Or, set adaptScaleOnly = TRUE in the case when adaptInterval = 1, so that no adaptation of the proposal covariance happens? Or, modify the RW_block sampler it always uses more than 1 sample for this estimation? Just ideas..

paciorek commented 3 weeks ago

Simplest thing given this hasn't come up in practice is that we restrict adaptInterval to be at least 2. I'm going to add that error trap.