Rejection of Metropolis proposal

General
1

Dear all, Could someone please have a look at my code which always complains of the rejection of Metropolis proposal? The priors seem reasonable. Here is the error:

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: lognormal_lpdf: Location parameter is -nan, but must be finite! (in ‘/tmp/RtmpcpVYpz/model-26327a33e7fb.stan’, line 40, column 4 to column 54)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

and here is the code:

library(rethinking)
data("Lynx_Hare")
d <- Lynx_Hare
Lynx.dat = list(N =21, L = d$Lynx, H = d$Hare, Y = d$Year)
Lynx.c = '
functions {
     vector MyODE( real t,                 // time
                  vector pop_init,          // initial state {lynx, hares}
                  vector theta){         // parameters
    real L = pop_init[1];
    real H = pop_init[2];
    real bH = theta[1];
    real mH = theta[2];
    real bL = theta[3];
    real mL = theta[4];
    // differential equations
    vector[2] dif;
    dif[1] = (bH - mH * L) * H;
    dif[2] = (bL * H - mL) * L;
    return dif; 
  }
}
    data{
    int N; 
    vector[N] L; 
    vector[N] H; 
    }
    parameters{
    vector<lower=0>[4] theta; //{bH, mH, bL, mL}
    vector<lower=0>[2] sigma;
    vector<lower=0,upper=1>[2] p;
    }
    model{
    array[N-1] real ts;
    for(i in 1:(N-1)) ts[i] = i;
    array[N] vector[2] pop;
    pop[1,1] = H[1];
    pop[1,2] = L[1];
    pop[2:N, 1:2] = ode_rk45(MyODE, to_vector(pop[1]), 0, ts, theta);
    for ( t in 1:N ){
    target += lognormal_lpdf(H[t] | log(pop[t , 1]*p[1]), sigma[1]); 
    target += lognormal_lpdf(L[t] | log(pop[t , 2]*p[2]), sigma[2]); 

    }
    p ~ beta(40, 200);
    sigma ~ exponential(1);
    theta[{1,4}] ~ normal(1, 0.5); // {bH, mL}
    theta[{2, 3}] ~ normal(0.1, 0.05); // {mH, bL}
    }
'
Lynx.f <- stan( model_code=Lynx.c,data=Lynx.dat,chains=4, iter = 3000, cores=1 , control=list( adapt_delta=0.95 ) )

Thanks :)

2

The error message tells you where the problem is. This is why we recommend putting Stan programs in their own files—it helps in figuring out the line number.

Given that it says it’s a lognormal and you only have two, it’s one of these:

    target += lognormal_lpdf(H[t] | log(pop[t , 1]*p[1]), sigma[1]); 
    target += lognormal_lpdf(L[t] | log(pop[t , 2]*p[2]), sigma[2]);

Given that it says the location parameter is -nan (i.e., not-a-number), it means either log(pop[t , 1]*p[1]) or log(pop[t , 2]*p[2]) produced a not-a-number. If you look at the doc for lognormal_lpdf, you’ll see the first argument (after the vertical bar) is the location parameter. So there’s a problem with pop[t, 1] * p[1] or pop[t, 2] * p[2] not being finite positive numbers.

The value p is constrained to be in (0, 1), but there are no constraints on pop as it’s just set to L and H, which are unconstrained data variables. Is it possible that L and H have negative entries? I’d declare them with <lower=0> just to make sure. Or that ode_rk45 is returning negative values?

I also think there’s a problem at H[1] and L[1] in that the variables show up on the right and left hand sides. What you normally want to do here is not set the population to be H[1] and L[1] but to have it be an unknown centered around H[1] and L[1]. Otherwise the model won’t be properly generative.

3

Dear Bob,
Thank you for your time and comprehensive reply.
I tried to set a lower bound on p but I got this error(I am using rstan 2.26.15)

31:      array[N-1] real ts;
32:      for(i in 1:(N-1)) ts[i] = i;
33:      array[N] vector<lower=0>[2] pop;
                        ^
34:      pop[1,1] = pop_init[1];
35:      pop[1,2] = pop_init[2];

Expected “[” expression “]” for vector size.

You suggested to set a constraint on L and H (which are data) and also on pop(which is calculated deterministically from data and parameters). Isn’t it true that you can constrain only when you sample a random number? What does it mean to set a constraint on data?

Given a time series data that you want to fit a dynamics to, the first datapoint (t=0) is not informative as you are interested in change with time. Right? Why not forcing the fit to start at the first datapoint.
Thanks

4

I should have been more specific. You can’t put a constraint on a local variable. You have to put the constraint on parameters that will make sure your derived quantities obey the constraint. If you’re getting these as solutions to a differential equation, the equation probably needs to be set up to ensure you can’t get a negative answer. One way to do that is to reparameterize on the log scale. In some cases, you might get negative solutions from bad initializations out in the tail—if you can initialize closer to an expected answer you might be able to avoid the negative values.