Imagine a model with two binary data vectors
b. Both represent events that can either occur or not occur. However,
b is impossible when
a does not occur. I’d like to model the probability that
b occurs given that
a occurs with a probit model.
In stan, I can do this with a conditional statement, e.g.
p[i] = (a[i] == 0) ? 0 : Phi( b0 );
which just makes the probability that
b occurs equal to 0 whenever
a does not occur; when
a does occur, I get my probit model. This seems to work well.
In greta, I don’t see that there are conditions. I’m trying to emulate this with a multiplication:
library(greta) ar <- rep(c(0,1), each = 100) br <- rep(c(0,1), 100) br[a0 == 0] = 0 table(ar,br) a <- as_data(ar) b <- as_data(br) b0 <- normal(0,1) p <- (a == 0) * iprobit( b0 ) distribution(b) <- bernoulli( p ) m <- model(b0) plot(m) inits <- initials( b0 = 0) draws <- mcmc(m, n_samples = 1000, initial_values = inits)
a is 0, the probability is 0; else, it is Phi(b0). However, I cannot get this model to run. I get
Error: The log density could not be evaluated at these initial values. Try using calculate() see whether they lead to values of other greta arrays in the model.
When I remove the
(a == 0) * bit, the model runs fine, but gives me the marginal estimate of the probability (which is not what I want, of course).
What’s the best way of accomplishing this in greta?