Two coins discrete parameter model

[Ben Lambert]

Hi,

I am trying to implement a discrete parameter model to include in a document for my students. The model I have is the following,

I have two coins in my pocket - one has a probability of heads of 0.1, and the other is 0.9. In each experiment I pull one out at random and record the result for N flips. I want to implement a model that estimates the coin type in each experiment, as well as the probabilities of heads for each coin.

At the moment I have implemented this as below -- it seems to be working -- it estimates out the states of fake data very well. However I wanted to double check this implementation as I don't feel 100% comfortable with discrete stuff in Stan. 

Does this look ok? I know I could have used a hierarchical model where I pick coin one with probability phi, and coin two with probability 1-phi, but I wanted to keep it simple for now.

data { 

  int<lower=1> nStudy; // num studies

  int<lower=1> N; // samples per study

  int<lower=2> K; // num coin types

  int<lower=0, upper=1> X[N,nStudy];

} 

transformed data{

  real log_unif;

  log_unif = -log(K);

}

parameters { 

  ordered[K] alpha; 

}

transformed parameters{

  real<lower=0,upper=1> theta[K];

  matrix[nStudy,K] lp;

  

  for(i in 1:K){

    theta[i] = inv_logit(alpha[i]);

  }

  lp = rep_matrix(log_unif,nStudy,K);

  

  for(n in 1:nStudy){

    for(i in 1:N){

      for(s in 1:K){

        lp[n,s] = lp[n,s] + bernoulli_lpmf(X[i,n]|s==1? theta[1]:theta[2]);

      }

    }

  }

}

model { 

  for(n in 1:nStudy){

    target += log_sum_exp(lp[n]);

  }

}

generated quantities{

  int state[nStudy];

  for(n in 1:nStudy){

    state[n] = categorical_rng(softmax(to_vector(lp[n])));

  }

}

[Bob Carpenter]

> On Aug 19, 2016, at 5:33 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Hi,
>
> I am trying to implement a discrete parameter model to include in a document for my students. The model I have is the following,
>
> I have two coins in my pocket

Your model is partly designed for K coins but then hard codes K = 2.
I'd just go with a hard coded K = 2 for now and then generalize
when that's working. Our general advice is to start simple and then
extend.

> - one has a probability of heads of 0.1, and the other is 0.9.

Don't tell Andrew :-)

http://www.stat.columbia.edu/~gelman/research/published/diceRev2.pdf

> In each experiment I pull one out at random and record the result for N flips. I want to implement a model that estimates the coin type in each experiment, as well as the probabilities of heads for each coin.

The thing to do is start with the discrete indicator z[m] in { 1, 2 } indicating
which coin is drawn for study m in 1:M. I'm also going to let y[m] in 0:N
be the number of successes (er, heads) It's much more efficient than the
Bernoulli formulation. If theta[1] and theta[2] in (0,1) are the probabilities
of success (heads), then assuming the coin is selected randomly, the likelihood's just:

p(z, y | theta) = PROD_{m in 1:M} p(z[m], y[m] | theta)

Then you can marginalize inside the product

p(y[m] | theta) = 0.5 * binomial(y[m] | theta[1]) + 0.5 * binomial(y[m] | theta[2])

or

log p(y[m] | theta) = log_sum_exp(log(0.5) + log binomial(y[m] | theta[1]),
log(0.5) + log binomial(y[m] | theta[2]))

And that's just how you write the model in Stan, along with any prior
you might put on theta.

You can assume logit_theta[k] in (-inf, inf) and use binomial_logit
if you prefer that parameterization

> At the moment I have implemented this as below -- it seems to be working -- it estimates out the states of fake data very well. However I wanted to double check this implementation as I don't feel 100% comfortable with discrete stuff in Stan.
>
> Does this look ok?

Looks OK to me. The ordering you impose for identifiability can
be a problem because it can affect the posterior uncertainty (for
instance, the first coin will always have a lower probability of
heads than the second) but won't be an issue if the probabilities
are 0.1 and 0.9 and you have a reasonable amount of data.

Your sum over lp is just doing the binomial by hand.

> I know I could have used a hierarchical model where I pick coin one with probability phi, and coin two with probability 1-phi, but I wanted to keep it simple for now.

That's not what we'd call a hierarchical model --- it's just a mixture
with not necessarily symmetric mixing proportions.

The generated quantities should just be coded in expectation rather
than using a categorical_rng and you also don't need all that softmax
stuff for K = 2, which is where I'd start.

Pr[z[m] = 1 | data] = binomial(y[m] | theta[1] ) / (binomial(y[m] | theta[1]) + binomial(y[m] | theta[2]))

or on the log scale just takes the log.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Ben Lambert]

Bob thanks for your response. Agreed about the binomial - I will change to that! I also agree about your hard-coding idea.

A couple of questions. I have tried this,

for(n in 1:nStudy){

    for(s in 1:K){

      lp[n,s] = lp[n,s] + binomial_lpmf(Z[n],N|theta[s]);

    }

  }

But Stan seems to throw an error here, which doesn't really have any explanation. Any idea why this might be?

Finally, is this what you had in mind for the expectation? (Also why should I use this opposed to using categorical_rng?)

generated quantities{

  int pstate[nStudy];

  for(n in 1:nStudy){

    pstate[n] = exp(lp[n,1])/(exp(lp[n,1]) + exp(lp[n,2]));

  }

}

Thanks again for your help here.

Best,

Ben

[Bob Carpenter]

Yes, the vertical bar should be after the first argument.

But given the form, you don't need to do any incrementing.

You want to do this for each n, not sum over all the n in nStudy.
Your approach with the binomial would correspond to where you
draw a single coin and then use it for each study.

- Bob

[Ben Lambert]

Sorry am being slow here -  I don't follow what you mean by 'incrementing' and why what I have corresponds to using one coin across all studies. I have simplified the code now to hopefully make clear what's happening,

transformed parameters{

  real<lower=0,upper=1> theta[2];

  matrix[nStudy,2] lp;

  

  for(i in 1:2){

    theta[i] = inv_logit(alpha[i]);

  }

  

  for(n in 1:nStudy){

    for(s in 1:2){

      lp[n,s] = log(0.5) + binomial_lpmf(Z[n]|N,theta[s]);

    }

  }

}

Is this ok? So I am storing the log probability for each of the states in s for each of the individual experiment. There's no incrementing here. I am drawing one coin per study. 

In the model I do the following

In the generated quantities I then do as above,

generated quantities{

  real pstate[nStudy];

  for(n in 1:nStudy){

    pstate[n] = exp(lp[n,1])/(exp(lp[n,1]) + exp(lp[n,2]));

  }

}

where I use the log probs of each state (converted into probs) to infer the expectation a la Rao-Blackwell. Is this bit right?

Best,

Ben

[Bob Carpenter]

I think I misread the loop because of the non-standard
(for us) loop variables (still my fault---I try to
answer these questions too quickly so I don't spend all
day every day on our email list).

Yes, the generated quantities is right now.
And assuming the likelihood is

for (n in 1:nStudy)
target += log_sum_exp(lp[n, 1], lp[n, 2]);

then I think you're good to go.

Though no need to inv_logit the theta---you can just
use our binomial_logit or bernoulli_logit.

What people sometimes do is this:

p(y, theta) = 0.5 * PROD_n p(y[n] | theta[1])
+ 0.5 * PROD_n p(y[n] | theta[2]);

which is a mixture over the whole data set rather than the individual
observations. And it usually winds up looking sort of
like the loop you had.

- Bob

> On Aug 19, 2016, at 8:43 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Sorry am being slow here - I don't follow what you mean by 'incrementing' and why what I have corresponds to using one coin across all studies. I have simplified the code now to hopefully make clear what's happening,
>
> transformed parameters{
> real<lower=0,upper=1> theta[2];
> matrix[nStudy,2] lp;
>
> for(i in 1:2){
> theta[i] = inv_logit(alpha[i]);
> }
>
> for(n in 1:nStudy){
> for(s in 1:2){
> lp[n,s] = log(0.5) + binomial_lpmf(Z[n]|N,theta[s]);
> }
> }
> }
>
> Is this ok? So I am storing the log probability for each of the states in s for each of the individual experiment. There's no incrementing here. I am drawing one coin per study.
>
> In the model I do the following
>
> In the generated quantities I then do as above,
>
> generated quantities{

> int pstate[nStudy];

> for(n in 1:nStudy){
> pstate[n] = exp(lp[n,1])/(exp(lp[n,1]) + exp(lp[n,2]));
> }
> }
>
> where I use the log probs of each state (converted into probs) to infer the expectation a la Rao-Blackwell. Is this bit right?
>
> Best,
>
> Ben
>

[Bob Carpenter]

I should've also added that you probably want to keep
the Pr[z[n] = 1] calculation on the log scale, so
instead of

exp(l1) / (exp(l1) + exp(l2)),

where l1 and l2 are unnormalized log probabilities, you just
use

l1 - log_sump_exp(l1, l2).

- Bob