binomial-binomial model with stan

[Soga Shota]

Hi all,

   I just wrote a binomial-binomial model using stan via brute-force marginalization. I think this would be a good example for those who want to marginalize discrete variables, so I post it here.

In BUGS notation, the following stan code is equivalent to 

model

{

   theta1 ~ dunif(0,1)

   theta2 ~ dunif(0,1)

   x ~ dbin(theta1,50)

   y ~ dbin(theta2,x)

}

DATA

list(

y = 5

)

In stan:

transformed data

{

   int y;

   int N;

   int K;

   y <- 5;

   N <- 50;

   K <- N-y+1;

}

parameters

{

   real<lower=0,upper=1> theta1;

   real<lower=0,upper=1> theta2;

}

model

{

   

   {

      real tmp_vec[K];

      for(k in y:N)

      {

         real tmp_factorial; #some terms are redundant, but I leave it for those who how I derive this 

         tmp_factorial <- (tgamma(N+1.0) /tgamma(k+1.0) / tgamma(N-k+1.0)) * (tgamma(k+1.0)/ tgamma(y+1.0) / tgamma(k-y+1.0));

         tmp_vec[k-y+1] <- tmp_factorial * pow(theta1,k) *pow(1.0-theta1,(N-k)) * pow(theta2,y)  * pow(1.0-theta2,(k-y));

      }

      increment_log_prob(log(sum(tmp_vec)));

   }

}

Enjoy!

[Soga Shota]

I just made it via user defined function block.

functions

{

   void binbin_lp(int y_,int N_,real theta1_,real theta2_)

   {

      int K_;

      K_ <- N_ - y_ + 1;

      {

      real tmp_vec[K_];

      for(k_ in y_:N_)

      {

         real tmp_factorial;

         tmp_factorial <- (tgamma(N_+1.0) /tgamma(k_+1.0) / tgamma(N_-k_+1.0)) * (tgamma(k_+1.0)/ tgamma(y_+1.0) / tgamma(k_-y_+1.0));

         tmp_vec[k_-y_+1] <- tmp_factorial * pow(theta1_,k_) *pow(1.0-theta1_,(N_-k_)) * pow(theta2_,y_)  * pow(1.0-theta2_,(k_-y_));

      }

      increment_log_prob(log(sum(tmp_vec)));

      }

      return;

   }

}

Now, you can use binbin_jp in model block (in the manual, _lp signature is used for transformed parameters block, but we can use it in model block, too).

2014年7月15日火曜日 14時25分51秒 UTC+9 Soga Shota:

[Bob Carpenter]

Thanks for sharing. It's the same kind of thing
we'd have to do for binomial missing data. Was there
an application driving this?

I would think this would all be more arithmetically
stable on the log scale (less chance of underflow).
So my knee-jerk reaction is to convert the linear scale
formulation

> {

> real tmp_vec[K];
> for(k in y:N)
> {
> real tmp_factorial; #some terms are redundant, but I leave it for those who how I derive this
> tmp_factorial <- (tgamma(N+1.0) /tgamma(k+1.0) / tgamma(N-k+1.0)) * (tgamma(k+1.0)/ tgamma(y+1.0) / tgamma(k-y+1.0));
> tmp_vec[k-y+1] <- tmp_factorial * pow(theta1,k) *pow(1.0-theta1,(N-k)) * pow(theta2,y) * pow(1.0-theta2,(k-y));
> }
> increment_log_prob(log(sum(tmp_vec)));
> }

into the following log-scale version (which won't underflow or overflow
like the tgamma andprod are prone to do with large k or N, especially
wen multiplied together):

{
real log_tmp_vec[K];
for (k in y:N)
log_tmp_vec[k - y + 1]
<- lgamma(N + 1)
- lgamma(k + 1)
- lgamma(N - k + 1)
+ lgamma(k + 1)
- lgamma(y + 1)
- lgamma(k - y + 1)
+ k * log(theta1)
+ (N - k) * log1m(theta1)
+ y * log(theta2)
+ (k - y) * logm(theta2);
increment_log_prob(log_sum_exp(log_tmp_vec));
}

But now I see the structure is just two binomials,
so I'm pretty sure everything reduces to this:

{
real tmp[K];
for (k in y:N)
tmp[k - y + 1]
<- binomial_log(k, N, theta1)
+ binomial_log(y, k, theta2);
increment_log_prob(log_sum_exp(tmp));
}

At least I think that'll give you the same model. It makes
sense given the original BUGS formulation, too.

And speaking of arithmetic, it's usually better to just work
on the int scale (both faster and more accurate) rather than forcing
early double conversion. That is, (k + 1) will evaluate the addition as
integers, then convert to double as the argument to lgamma. But doing
(k + 1.0) will first cast k to a double, then do the addition.

- Bob

On Jul 15, 2014, at 1:25 AM, Soga Shota <soga....@gmail.com> wrote:

> Hi all,
>
> I just wrote a binomial-binomial model using stan via brute-force marginalization. I think this would be a good example for those who want to marginalize discrete variables, so I post it here.
>
> In BUGS notation, the following stan code is equivalent to
>
> model
> {
> theta1 ~ dunif(0,1)
> theta2 ~ dunif(0,1)
> x ~ dbin(theta1,50)
> y ~ dbin(theta2,x)
> }
> DATA
> list(
> y = 5
> )
>

> transformed data
> {
> int y;
> int N;
> int K;
> y <- 5;
> N <- 50;
> K <- N-y+1;
> }
> parameters
> {
> real<lower=0,upper=1> theta1;
> real<lower=0,upper=1> theta2;
> }
> model
> {
>
> {
> real tmp_vec[K];
> for(k in y:N)
> {
> real tmp_factorial; #some terms are redundant, but I leave it for those who how I derive this
> tmp_factorial <- (tgamma(N+1.0) /tgamma(k+1.0) / tgamma(N-k+1.0)) * (tgamma(k+1.0)/ tgamma(y+1.0) / tgamma(k-y+1.0));
> tmp_vec[k-y+1] <- tmp_factorial * pow(theta1,k) *pow(1.0-theta1,(N-k)) * pow(theta2,y) * pow(1.0-theta2,(k-y));
> }
> increment_log_prob(log(sum(tmp_vec)));
> }
> }
>
> Enjoy!
>

> --
> 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.
> For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

On Jul 15, 2014, at 1:56 AM, Soga Shota <soga....@gmail.com> wrote:

> I just made it via user defined function block.
>
> functions
> {
> void binbin_lp(int y_,int N_,real theta1_,real theta2_)
> {
> int K_;
>
> K_ <- N_ - y_ + 1;
> {
> real tmp_vec[K_];
>
> for(k_ in y_:N_)
> {
> real tmp_factorial;
> tmp_factorial <- (tgamma(N_+1.0) /tgamma(k_+1.0) / tgamma(N_-k_+1.0)) * (tgamma(k_+1.0)/ tgamma(y_+1.0) / tgamma(k_-y_+1.0));
> tmp_vec[k_-y_+1] <- tmp_factorial * pow(theta1_,k_) *pow(1.0-theta1_,(N_-k_)) * pow(theta2_,y_) * pow(1.0-theta2_,(k_-y_));
> }
> increment_log_prob(log(sum(tmp_vec)));
> }
> return;
> }
> }

If you define instead (and I'd suggest not using trailing underscores
for arguments --- those are conventionally used for member variables
in classes, at least in C++):

real binbin_log(int y, int N, real theta1, real theta2) {
... as in your code ...
return log(sum(tmp_vec));
}

The "_log" suffix is important --- it's what lets you call it as:

y ~ binbin(N, theta1, theta2);

in the model block. You can also use binbin_log() as a function,
and just wrap it in increment_log_prob() externally.

It's a very clunky name-mangling and identifier-based syntax
we're using --- very old-school Fortran.

- Bob

[Soga Shota]

Hi Bob,

   The binomial-binomial model will be used in industrial application. Actually, my model is more complicated: binomial-binomial-logitnormal model.

   log_sum_exp is a very interesting function. so log_sum_exp(vec) basically does log(sum(exp(vec))).

It is quite useful in simplifying marginalization.

Shota

2014年7月15日火曜日 15時18分08秒 UTC+9 Bob Carpenter:

[Luc Coffeng]

I think in this case the structure can be further simplified with the statement k ~ binomial(N, theta1 * theta2), which also touches on the theta's currently being unidentified (and defeats the purpose of the exercise).

[Luc Coffeng]

Whoops, I meant y ~ binomial(N, theta1 * theta2), of course!

[Bob Carpenter]

On Jul 15, 2014, at 2:44 AM, Soga Shota <soga....@gmail.com> wrote:

> Hi Bob,
>
> The binomial-binomial model will be used in industrial application. Actually, my model is more complicated: binomial-binomial-logitnormal model.
>
> log_sum_exp is a very interesting function. so log_sum_exp(vec) basically does log(sum(exp(vec))).
> It is quite useful in simplifying marginalization.

The important feature of log_sum_exp is that it won't overflow with
the exp() operation in either values or derivatives. It also has
analytic derivatives so it's faster than calling log(sum(exp(vec)))
and takes less memory.

- Bob

[Bob Carpenter]

I was thinking there might be such a reduction.

I can see that binomial(N, theta1 * theta2) will have
the same mean, but will it have the same variance? I
would think the compound version would be overdispersed,
but I don't know how to work out its variance other than
by simulation.

Presumably in the richer model, these won't be simple
parameters but the results of regressions on a log-odds
scale.

- Bob

[Luc Coffeng]

Hi Bob,

Hm, I see your point.

I've experimented with binomial-binomial models in a setting where observations can be misclassified, so there are unobserved true values and observed values that are a function of the true values and the specificity (% true negative detected) and sensitivity (% true positive detected) of some test that determines positiveness. Because there are many sums of false and true positives that can lead to the observed values (and sums of true and false negatives that can lead to observed negatives), it is significantly more complicated to marginalize over all possible combinations, but the concept of marginalizing over all possible sums is similar to the current example (I think). Anyway, I found that a stan model based on marginalization behaved exactly the same as a simple model with: y ~ binomial(N, p_true * sensitivity + (1-p_true) * (1-specificity), including the uncertainty estimates. I tried to find the code this morning so that I could post it here, but couldn't find it (and had to head out to the office). Will post it later if I can find it.

Best,
Luc

[Luc Coffeng]

In terms of my misclassification model, the current example would be a specific case where specificity of the test is 100%, and p_true and sensitivity are theta1 and theta2, respectively:

y ~ binomial(N, p_true * sensitivity + (1-p_true) * (1-specificity))  -->  y ~ binomial(N, p_true * sensitivity)

-Luc

[Soga Shota]

Hi Luc,

   Thank you for your information. I am very curious in your suggestion because neither my version nor Bob's version of binbin is applicable due to its speed... :(

   I would greatly appreciate it if you would provide any reference that justify y~binomial(N,theta1*theta2) since I have a trouble to show it is analytically correct.

Shota

2014年7月16日水曜日 5時07分56秒 UTC+9 Luc Coffeng:

[Soga Shota]

Hi Luc,

   I just found it in wikipedia(it does not have a citation, tho)... http://en.wikipedia.org/wiki/Binomial_distribution#Conditional_binomials

EDIT:

http://math.stackexchange.com/questions/626457/conditional-binomials

Shota

2014年7月16日水曜日 10時21分07秒 UTC+9 Soga Shota:

[Luc Coffeng]

Hi Shota,

Thanks for providing that link! I tried the analytical way before, and quickly decided that I couldn't deal with all the gamma functions flying around :). So I just tested it empirically, as I mentioned before. Will post my script if I manage to dig it up!

Luc

[Soga Shota]

Hi Luc,

   when I edit the mailing list on web, I am not sure it will be sent it to mailing list again, so I post new link again.

http://math.stackexchange.com/questions/626457/conditional-binomials

I haven't check the math, but this guy showed it analytically.

Shota

2014年7月16日水曜日 10時55分11秒 UTC+9 Luc Coffeng: