Fitting time-series issues

Modeling
1

I am trying to fit a consumer-resource model old and common in ecology in Stan, but struggling to get any of the parameter values returned. I have never done time-series in Stan so fear I am doing something quite obvious wrong, but have struggled to figure out what.

Below is the R and Stan code, any help appreciated.

Lizzie

library(rstan)
library(nimble)
# Set up the data
p <- list(r = .05, K = 2, Q = 5, H = .38, sigma = .01, a = 0.023, N = 1e4)
growth <- function(x, p) x * p$r * (1 - x / p$K)
consumption <- function(x,p) p$a * x ^ p$Q / (x^p$Q + p$H^p$Q)
# Define stochastic model in BUGS notation
may  <- nimble::nimbleCode({
  x[1] <- x0
  for(t in 1:(N-1)){
    mu[t] <- x[t] + x[t] * r * (1 - x[t] / K)  - a * x[t] ^ Q / (x[t] ^ Q + H ^ Q)
    y[t+1] ~ dnorm(mu[t], sd = sigma)
    x[t+1] <- max(y[t+1],0)
  }
  
})
model <- nimbleModel(may,constants = p, inits = list(x0 = 0.2))
cmodel <- compileNimble(model)
set.seed(123)
simulate(cmodel)

x=cmodel$x
dat <- list(x=x, N = length(x))

mod1 = stan('ghostfit.stan', data = dat,
               chains=4, iter = 2000, cores=4)

And the Stan code:

// Fitting mechanistic model (May 1976)
data {
  int<lower=1> N; 
	vector[N] x; // response
	}

parameters {
  real<lower=0.049, upper=0.051> r; // I realize this is bad, but not the problem I think ...
  real<lower=1.9, upper=2.1> K; // ibid
  real<lower=0> Q;  
  real<lower=0> H;   
  real<lower=0> a;   
  real<lower=0> x0;
  real<lower=0> sigma;   
}

model {	
   real mu[N];
   // x[1] = x0;
   for(t in 1:(N-1)){
            mu[t] =  x[t] + x[t] * r * (1 - x[t] / K)  - a * x[t] ^ Q / (x[t] ^ Q + H ^ Q);
            x[t+1] ~ lognormal(mu[t], sigma);
   }
            
    Q ~ normal(5, 0.1);
    H ~ normal(0.38, 0.1);
    a ~ normal(0.0233, 0.1);
    sigma ~ normal(0.01, 0.001);
    x0 ~ normal(0, 1);

}
2

I don’t see how the Stan script would fail to produce parameter results(though x0 is independent of the data).

3

I’m not familiar with nimble, but isn’t dnorm in BUGS the normal distribution? In the nimble code you have

y[t+1] ~ dnorm(mu[t], sd = sigma)

but in the stancode you have the lognormal distribution…

And, as @yizhang wrote, why is

// x[1] = x0;

commented out?