Hello,

I am fitting a non-linear model of the form

y(x) = p_2 + (p_4 - p_2) * e^{-e^{p_1} * (x - p_3) * g((x - p_3) * 10)}

where g is the inverse logit function. This is being fitted to a beta-distributed data (proportion):

```
data_str <- "y,x
0.838313611,0.09162907
0.879521184,0.09162907
0.850763463,0.09162907
0.865867666,0.09162907
0.881867882,0.09162907
0.849611448,0.09162907
0.999000000,0.91629073
0.920652405,0.91629073
0.919116384,0.91629073
0.947916772,1.60943791
0.948684782,1.60943791
0.893691199,1.60943791
0.890315996,2.30258509
0.926028477,2.30258509
0.916428348,2.30258509
0.662600928,2.70805020
0.628424467,2.70805020
0.590023950,2.70805020
0.446406015,3.21887582
0.366148933,3.21887582
0.335064725,3.21887582
0.201794719,3.55534806
0.287427873,3.55534806
0.298685288,3.55534806
0.076329756,3.91202301
0.114032082,3.91202301
0.104257405,3.91202301
0.007632976,4.60517019
0.007632976,4.60517019
0.007632976,4.60517019"
data <- read.table(text = data_str, sep = ",", header = TRUE)
```

I am trying to fit it in brms using my own initial values, but I am getting lots of the classic rejection â€śRejecting initial valueâ€ť sampling statement:

```
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: beta_lpdf: First shape parameter[1] is -0.071096, but must be > 0! (in 'model4936186f04f5_844d16a56a3a8ac20fc7d8c9ff39d375' at line 47)
```

Line 47 in the generated stan code is:

```
target += beta_lpdf(Y | mu * phi, (1 - mu) * phi);
```

this means that `mu < 1`

. I have checked for that though, and I am 100% positive that my initial values yield `mu`

values that conform with the beta distribution expectation of shape parameters > 0 and < 1.

```
my_inits <- list(
list(b_p1 = -0.07, b_p2 = 0.11, b_p3 = 4.52, b_p4 = 0.96, phi = 0.67),
list(b_p1 = 0.19, b_p2 = 0.09, b_p3 = 3.36, b_p4 = 0.97, phi = 0.07),
list(b_p1 = 1.50, b_p2 = 0.07, b_p3 = 3.58, b_p4 = 0.93, phi = 0.30),
list(b_p1 = 0.49, b_p2 = 0.30, b_p3 = 3.83, b_p4 = 0.97, phi = 0.64)
)
mu_formula <- function(p1, p2, p3, p4, x) {
p2 + (p4 - p2) * exp(-exp(p1) * (x - p3) *
plogis((x - p3) * 10))
}
for (i in seq_along(my_inits)) {
print(range(mu_formula(my_inits[[i]]$b_p1[1], my_inits[[i]]$b_p2[1],
my_inits[[i]]$b_p3[1], my_inits[[i]]$b_p4[1],
data$x)))
}
[1] 0.9139801 0.9611010
[1] 0.2852354 0.9999382
[1] 0.07869372 0.97270772
[1] 0.4891421 0.9883545
```

Yet, brms cannot fit the model. I do not quite understand what could be going on. Does anyone have an idea? Shouldnâ€™t those values ensure that the very first calculation of mu is valid? Or are the initial values intended to do different thing?

Full reprex:

```
library(brms)
priors <- prior("normal(0, 1", nlpar = "p1") +
prior("beta(1, 5)", nlpar = "p2") +
prior("gamma(5, 1)", nlpar = "p3") +
prior("beta(5, 1)", nlpar = "p4") +
prior("gamma(0.1, 0.1)", class = "phi")
bform <- bf(y ~ p2 + (p4 - p2) * exp(-exp(p1) *
(x - p3) * inv_logit((x - p3) * 10)),
p2 + p4 + p1 + p3 ~ 1,
nl = TRUE)
test <- brm(bform, data = data, family = Beta(link = "identity"),
prior = priors, iter = 2e4, warmup = 1.6e4,
chains = 4, cores = 4, inits = my_inits)
```