Hi All:

I recently mimicked a Richard McElreath blog entry regarding marginalizing an unobserved discrete random variable. His work is in Stan, so I reproduced the workflow in Greta and found it much easier/simpler/elegant.

My biggest fear is that it looks so easy that maybe I am missing something. Going in to the exercise, I thought I would use a mixture() distribution with weights determined by the data. If you take a look, you will see that is not what I had to do. It was really quite astounding how it just worked. Anyway, check out my write-up at:

and feel free to compare to McElreath’s orginal: http://elevanth.org/blog/2018/01/29/algebra-and-missingness/

I appreciate any feedback.

The minimum code to reproduce the work extracted from the blog entry is pasted below for convenience.

This code generates some sample data:

```
library(dplyr)
library(greta)
# Actual Values That Models Will Try To Recover
actual_p_u = 0.5; actual_p_s = 1; actual_sigma = 0.75
set.seed(1)
N_children <- 51
s <- rbinom( N_children , size=1 , prob=actual_sigma )
s_obs <- s # copy to get censored obs
s_obs[ sample( 1:N_children , size=21 ) ] <- NA ##censored obs NA
tea <- rbinom( N_children , size=1 , prob= s*actual_p_s + (1-s)*actual_p_u )
dataDF = tibble(s_obs,tea)
nonCensoredDataDF = dataDF %>% filter(!is.na(s_obs))
censoredDataDF = dataDF %>% filter(is.na(s_obs))
```

And this code, runs the greta model:

```
library(greta)
s_i = as_data(nonCensoredDataDF$s_obs) #data
tea_i = as_data(nonCensoredDataDF$tea) #data
tea_j = as_data(censoredDataDF$tea) #data
p_s = beta(2,2) #prior
p_u = beta(2,2) #prior
sigma = beta(2,2) # hyperprior
pi_i = s_i * p_s + (1 - s_i) * p_u #operation for i
pi_j = sigma * p_s + (1 - sigma) * p_u # operation for j
distribution(s_i) = bernoulli(sigma) #likelihood i
distribution(tea_i) = bernoulli(pi_i) #likelihood i
distribution(tea_j) = bernoulli(pi_j) #likelihood j
m = model(p_s,p_u,sigma) #model
```

And here is the model block of the analogous Stan code:

```
model{
// priors
p_cheat ~ beta(2,2);
p_drink ~ beta(2,2);
sigma ~ beta(2,2);
// probability of tea
for ( i in 1:N_children ) {
if ( s[i] == -1 ) {
// ox unobserved
target += log_mix(
sigma ,
bernoulli_lpmf( tea[i] | p_drink ) ,
bernoulli_lpmf( tea[i] | p_cheat ) );
} else {
// ox observed
tea[i] ~ bernoulli( s[i]*p_drink + (1-s[i])*p_cheat );
s[i] ~ bernoulli( sigma );
}
}//i
}
```