Dear all,

I am trying to fit a model with a Gaussian optimum and ZI Poisson distribution. I’m trying to see whether I can fit autocorrelation between years in said optimum.

Right now, my brm call looks like this:

```
prior <-
prior(gamma(0.01, 0.01), class = "b", nlpar = "Wmax", lb = 0) +
prior(gamma(0.01, 0.01), class = "b", nlpar = "omega", lb = 0) +
prior(normal(0, 10), class = "b", nlpar = "theta") +
prior(cauchy(0,5), class = "sd", group = "Year", nlpar = "theta") +
prior(cauchy(0,5), class = "sd", group = "Year", nlpar = "Wmax")
formula <-
bf(Fitness ~ Wmax * exp(-(((Pheno.scale - theta) / omega)^2)),
omega ~ 1,
theta + Wmax ~ (1|Year),
nl = TRUE)
mod_paper <-
brm(formula = formula,
data = data,
prior = prior,
family = zero_inflated_poisson(link = "identity"),
warmup = 500,
iter = 1500,
thin = 1)
```

Based on this response from Paul Bürkner, it seems I need to define the auto-correlation process from within a function and perform a call like this:

```
formula <-
bf(Fitness ~ nlfun(Wmax, Pheno.scale, theta, omega),
omega ~ 1,
theta + Wmax ~ (1|Year),
nl = TRUE,
loop = FALSE)
```

What I am not sure to understand is how I can put a stochastic process within such a non-linear “formula” call. Would for example the following thing work (setting up sigma as new parameter of course):

```
theta[n] <- rnorm(1, theta[n-1], sigma)
```

It fells somehow weird to include such a random call within the non-linear function itself. What would be best practice here?