I’m implementing an IRT model and have followed the relevant chapter in the user guide to the letter (except for changing the `kk`

index to `ii`

, to match edstan’s format). Hence the model block is simply

```
model {
alpha ~ std_normal();
beta ~ normal(0, sigma_beta);
gamma ~ lognormal(0, sigma_gamma);
mu_beta ~ cauchy(0, 5);
sigma_beta ~ cauchy(0, 5);
sigma_gamma ~ cauchy(0, 5);
y ~ bernoulli_logit(gamma[ii] .* (alpha[jj] - (beta[ii] + mu_beta)));
}
```

I want to use `vb()`

in R because my data is rather large, and the model above works fine. But when I rewrite it as

```
model {
// ...
y ~ bernoulli(inv_logit(gamma[ii] .* (alpha[jj] - (beta[ii] + mu_beta))));
}
```

Stan balks and outputs `Error: stan::variational::advi::adapt_eta: Cannot compute ELBO using the initial variational distribution. Your model may be either severely ill-conditioned or misspecified.`

Now before you say why not just stick to `bernoulli_logit`

? It’s because eventually I want to generalize the two-parameter logistic to a three-parameter logistic, i.e.

\operatorname{Pr}[y_n = 1] = \delta_{ii} + (1 - \delta_{ii})\operatorname{logit}^{-1}(\gamma_{ii} [ \alpha_{jj} - (\beta_{ii} + \mu_{\beta}) ] )

which means having to split up `bernoulli_logit`

inevitably, into something like

```
model {
// ...
y ~ bernoulli(delta[ii] + (1 - delta[ii]) .* inv_logit(gamma[ii] .* (alpha[jj] - (beta[ii] + mu_beta))));
}
```

So I’m curious why `bernoulli(inv_logit(...))`

doesn’t work with `vb()`

. Because, surprisingly, `sampling()`

(with a small random subset of the data) works with *either* model (although the later runs slower). But `vb()`

does not. Any idea why? Is this a known issue?

P.S.: I suggest using edstan’s “aggression” data if you’re trying to replicate the error.

```
library(edstan)
data("aggression")
agg <- irt_data(y=aggression$dich,ii=aggression$item, jj=aggression$person)
vbsample <- vb(irt_model, data=agg)
```