How to use conditionals?

#1

Imagine a model with two binary data vectors a and 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)

So whenever 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?

#2

I was working on a similar problem. It looks like the following solution will work,by separating the “b” values. The only “b” values that are important to the model are when a != 0. Let me know if you agree:

library(greta)

ar <- rep(c(0,1), each = 1000)
br <- rep(c(0,1), 1000)
br[ar == 0] = 0
table(ar,br)

a <- as_data(ar)
b <- as_data(br[ar != 0])


b0 <- normal(0,1)
p <- iprobit( b0 )
distribution(b) <- bernoulli( p )

m <- model(b0)
plot(m)

inits <- initials( b0 = 0)

draws <- mcmc(m, n_samples = 1000,warmup = 1000, initial_values = inits)
summary(draws)
#3

I’ll have to check to see whether this will work for my model, but I suspect not. The index is lost, and so if these parameters are then used again, they are no longer matched to other potential data points.

#4

Yes, in your prediction phase, you would have to manually track the index and build in the logic in R equivalent to
p= ifelse(a == 0, 0 , pnorm( b ))

where b is the simulated chain of parameters

#5

That’s right, there isn’t a nice syntax for this type of conditional statement in greta. I’m not sure how commonly used that would be or what the syntax would look like, but if you have an idea for that as a feature please do post as an issue!

The reason you got an error saying you can’t evaluate at the initial values is that you’re passing a probability of 0 to the Bernoulli distribution (whenever a is not 0), and that’s outside the support of the parameter for that distribution. Ie. if the parameter of the Bernoulli distribution is 0, the log-probability of observing any value of b is -Inf.

As @howard says, the way forward here is to omit those elements where a is not equal to 0 from the data.