Cannot fit bayesian model with rethinking::quap for small (1e-4) parameter values

user_6326148 0
Viewed 232

I'm trying to calculate the posterior for p using rethinking package version 2.13:

rethinking::quap(
  alist(
    n ~ dbinom(N, p),
    p ~ dnorm(1e-4, 0.5e-4)    
  ),
  data = list(n = 10, N = 10e4),
  start = list(p = 10 / 10e4)
)

It throws an error:

Error in rethinking::quap(alist(n ~ dbinom(N, p), p ~ dnorm(1e-04, 5e-05)), : non-finite finite-difference value [1] Start values for parameters may be too far from MAP. Try better priors or use explicit start values. If you sampled random start values, just trying again may work. Start values used in this attempt: p = 1e-04

start is actually the MAP.
It works for bigger p values.
How to solve the problem?

bayesian numerical-methods r
1 Answers

tl;dr you're running into numerical approximation problems; adding control=list(ndeps=1e-6) to adjust the finite-difference tolerance is the fastest way to solve the problem.

Alternatively, you could change the model to fit p on the log scale:

alist(
    n ~ dbinom2(N, exp(logp)),
    logp ~ dnorm(-4, 2)
)

More generally, if you are working on a problem where p might be larger than, say, 0.5, you should probably fit on the log-odds or logit scale rather than the log scale, using plogis(logit_p) instead of exp(log_p). (I used log/exp instead of logit (qlogis)/plogis because I thought of it first, and because it's more familiar to most people.)

Two methods I tried that did not work (that I thought should have; I would need to dig in further to find out whynot): (1) specify the expected scale of the parameter via control=list(parscale=1e-4); (2) use method="L-BFGS-B", lower=0 to set a lower bound on the probability parameter.

The sad reality is that the tools for numerical evaluation of likelihoods, posteriors, etc. are fragile enough that you very quickly get to the point where you need to know something about what's going on under the hood.

We can diagnose what's going on by writing a dbinom2() function that's simply a "noisy" wrapper for dbinom:

dbinom2 <- function(x, size, prob, log) {
    r <- dbinom(x, size, prob, log)
    cat(x,size,prob,r,"\n")
    return(r)
}

Then re-run the estimation with dbinom2(...) in place of dbinom(...), and we get:

10 1e+05 1e-04 -2.078512 
10 1e+05 0.0011 -78.1496 
10 1e+05 -9e-04 NaN

The x and size values stay constant in this case; the search algorithm first tries prob=1e-4, then prob=0.0011, then prob=-9e-04 (a negative probability value), which is where it runs into trouble.

  • The error message gives more clues as to what is going on (if you already know quite a bit about how things work under the hood); quap() is using the "BFGS" method by default (it passes this straight to optim(), which means that it is trying to use the derivatives or gradients of the function to find the MAP estimate.
  • Because we don't have an explicit expression for the gradients, R does the next best thing (which is pretty crappy, but works OK for simple problems), which is to try to estimate the gradients by finite differences — basically replicating the definition of the derivative as lim(dx→0) (f(x+dx)-f(x))/dx, but using a small ("finite") dx value rather than trying to take the limit.
  • the problem arises in deciding on dx. Going way down into the help page for optim() (?optim()), you can see:

‘ndeps’ A vector of step sizes for the finite-difference approximation to the gradient, on ‘par/parscale’ scale. Defaults to ‘1e-3’.

This means that optim() is using dx=0.001, which is where your problem comes from. It tries f(x-dx) as well as f(x+dx), and since x-dx is <0, this breaks the calculation. Using ndeps=1e-6 specifies that dx should be smaller (small enough so that optim() doesn't try to evaluate the log-likelihood for a negative probability).

user_1902770
Related