Trouble passing initial values

Interfaces RStan
1

I think I’m missing something obvious here, but it appears that my initial values are being ignored. I’m trying to initialize a nonlinear state space model where the state (biomass) must be positive. On the scale of depletion P_t (fraction of carrying capacity K), with catch C_t the system evolves as

\tilde{P}_t = P_{t-1} + r P_{t-1} (1 - P_{t-1}) - C_t / K,\quad t = 2,\dots,23
P_t \sim \operatorname{log\ Normal}(\log(\tilde{P}_t), \sigma^2)

I’ve declared the state variables P_t to be positive, but the median prediction \tilde{P}_t can go negative depending on the previous state, growth rate, carrying capacity, and catch (passed as data) values. Rather than setting \tilde{P}_t to some small positive number, I’d like to fit it by rejecting parameter values that give these predictions. We know a priori that this population does not go extinct, so these values are not possible. I understand that this violates the Reference Manual warning that “Rejection is not for constraints”, but given the nonlinear multiparametric nature of these constraints, I’m not sure there’s a better way.

This obviously requires carefully setting initial values. Using fits where I did set \tilde{F}_t to a small number if it went negative, I have plenty of reasonable starting values. For simplicity, here I just try to set them to the rough means of each parameter. Following the examples in the stan help file, I defined

init_fn <- function() {
  list(r = 0.3,
       K = 250,
       q = 0.28,
       sigma2 = 0.0015,
       tau2 = 0.02,
       P = seq(1, 0.35, length.out = 22))
}

When I pass either of these to the stan function I get

> fit <- stan("src/models_reject/01_centered.stan", data = tuna_data, chains = 1, init = init_fn)

SAMPLING FOR MODEL '01_centered' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: Negative depletion, P_med[2] = -4.58502 P[t - 1] = 0.756282 r = 3.39455 K = 2.66466  (in 'model66024d7dfcb6_01_centered' at line 28)

.
.
.

Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
error occurred during calling the sampler; sampling not done

Which appears that it’s using init = "random" rather than the values from my passed function. The practical problem here is that K is on the order of 200, so init = "random" gives values that are far too small, and result in infeasible values for \tilde{P}_t. A workaround is to declare K_scaled as a parameter, then K = K * 200 as a transformed parameter (I can get this to sample), but I’d prefer to get the initial values to pass correctly.

I’ve attached a (slightly) simplified copy of my model and an R script that reproduces this behavior.

example.R (1.0 KB)
example.stan (1.7 KB)

Is there some reason that the stan function would ignore my inits?

2

I think it does recognize your inits. If you specify init = "0", you will see that it stops after one function evaluation. In your case, it gets past the first function evaluation but the leapfrog process goes into inadmissible territory before completing the first trajectory. Then it tries 100 more trajectories (with random values) before giving up.

Is P_t also supposed to have an upper bound of 1 that is not being accounted for? Usually the trick to these multivariate inequality constraints is to define a function in the functions block that calculates the minimum or maximum value for r that makes it admissible for a given value of K (and the other stuff).

3

That is what I’m seeing, and that makes more sense. Thanks for the explanation.

P_t doesn’t have an upper bound of 1; P_1 \sim \operatorname{log\ Normal}(\log(1), \sigma^2), so it often ends up (a bit) larger than one at the beginning of the trajectory.

Edit edit:

The constraints on each P_t can be expressed as

\frac{(1 + r) - \sqrt{(1 + r)^2 - 4 r C_t / K}}{2 r} < P_t < \frac{(1 + r) - \sqrt{(1 + r)^2 - 4 r C_t / K}}{2 r}

Looks like tomorrow’s project is to follow 16.4 in the User’s Guide to implement varying bounds for each P_t.