Implement a beta-possion in rstan

[Jin-Rui Xu]

Hello,

I am a first time user of rstan. I implemented the following model using rstan, and got sampling results. 

The implementation is so easy that I suspect its correctness.....

Could anyone read my code and give me some pointers?

Thanks in advance!

Jinrui

The model: 

x|s,p ~ Poisson(s*p)

p|kon,koff ~ Beta(kon,koff)

hyper-priors

koff ~ Gamma(akoff, bkoff)

kon ~ Gamma(akon, bkon)

s ~ Gamma(as,bs)

My code:

dat <- read.table(fin)

x <- dat; x <- as.vector(x); x <- as.numeric(x);

x <- round(x,digits=0)

dataList <- list(

         N = length(x),

         x = as.vector(x),

         as = 1, bs = 0.001,

         akon = 1, bkon = 0.001,

         akoff = 1, bkoff = 0.001

         )

excode <- '

data {

  int<lower=0> N;  # number of cells

  int<lower=0> x[N]; # expression

  real<lower=0> as; # expression

  real<lower=0> bs; # expression

  real<lower=0> akon; # expression

  real<lower=0> bkon; # expression

  real<lower=0> akoff; # expression

  real<lower=0> bkoff; # expression

#  real spj[N];

}

parameters {

  real<lower=1e-100> s;

  real<lower=1e-100> kon;

  real<lower=1e-100> koff;

  real<lower=1e-100,upper=1> pj[N];

#  vector[N] mu;

#  vector[N] delta;

#  real delta_new;

}

transformed parameters {

  real spj[N];

  for (j in 1:N){

    spj[j] <- s*pj[j];

  }

}

model {

  kon ~ gamma(akon, bkon);

  koff ~ gamma(akoff, bkoff);

  s ~ gamma(as,bs);

  for( i in 1:N){

    pj[i] ~ beta(kon, koff);

    x[i] ~ poisson(spj[i]);

  }

}

'

[Jonah Gabry]

Hi Jinrui,

In terms of translating that model into a Stan program your code looks pretty good. You can (and should) replace all the <lower=1e-100> constraints on parameters with <lower=0> unless 1e-100 is a logically necessary constraint (but I doubt it). You can also vectorize some of the code in the model block, which will be more efficient. In the model block you can replace 

for (i in 1:N) {
 pj[i] ~ beta(kon,koff)

 x[i] ~ poisson(spj[i])
}
 

with the vectorized version (no loop)

pj ~ beta(kon, koff)
x ~ poisson(spj)

which is equivalent to what you have but more efficient.

Also, I would double check that Stan's parameterization of the Gamma distribution is what you're expecting it to be (I can never remember off the top of my head but I think the second parameter of the Gamma in Stan's parameterization is a rate parameter whereas in other places you might see a scale parameter, i.e. 1/rate) . 

Jonah

[Jonah Gabry]

Also, a good way to check for problems is to simulate fake data from the assumed data generating process, fit the model to the fake data, and see if you can approximately recover the parameters used to generate the data.

[Bob Carpenter]

Teh thing to do is simulate data and make sure the model
is calibrated in the sense of recovering roughly p% of
the parameters in the posterior p% intervals.

You can speed this up a lot by using vectors, removing
the transformed parameters block altogether, and replacing
this:

for( i in 1:N){

pj[i] ~ beta(kon, koff);

x[i] ~ poisson(spj[i]);

}

with just these two lines:

pj ~ beta(kon, koff);
x ~ poisson(s * pj);

And while you can put independent priors on kon and
koff, we have a suggestion to parameterize in terms of

phi = kon / (kon + koff) [prior mean]

and

kappa = (kon + koff) [prior count + 2]

with phi in (0,1) and kappa in (0, inf),
and then use

pj ~ beta(phi * kappa, (1 - phi) * kkappa);

- 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.

[Jin-Rui Xu]

Thank you Jonah and Bob for your tips!

Jinrui