From JAGS code to STAN code (Categorical distribution Problems)

[Newbie (:]

I have initially use JAGS to run MCMC multinomial ordered probit model. There is a variable named pi1 which is the probability of a person getting infected with a disease. In my JAGS code, I have a vector k: k[1]<- 1-pi1 and k[2]<- pi1 and inf[i] ~ dcat(k[]). I want to let inf[i] to be 1 for non-infected and 2 for infected so that I can let mu[i]<- alpha when inf[i]==2. Btw i is the number of subjects.

 

So I wanted to code these in STAN but I will always get an error. I have used another method to replace those abovementioned by generating random number uniformly and determine the value of mu when this generated number is less than pi1. But I cannot include the code u[i] ~ uniform(0,1) in the model section of the STAN code so I create a vector of uniformly generated numbers first in R and let the vector to be in the data section.

 

So I hope someone can enlighten me on how to include the dcat code successfully in the model part in STAN. Attached named model_1 is the JAGS code. Thank you! :)

 

 

[Bob Carpenter]

On 9/24/13 4:05 AM, Newbie (: wrote:
> I have initially use JAGS to run MCMC multinomial ordered probit model. There is a variable named pi1 which is the
> probability of a person getting infected with a disease. In my JAGS code, I have a vector k: k[1]<- 1-pi1 and k[2]<- pi1
> and inf[i] ~ dcat(k[]). I want to let inf[i] to be 1 for non-infected and 2 for infected so that I can let mu[i]<- alpha
> when inf[i]==2. Btw i is the number of subjects.
> So I wanted to code these in STAN but I will always get an error.

It helps if you report the error you get. Otherwise, we have
to guess (see below).

> I have used another method to replace those
> abovementioned by generating random number uniformly and determine the value of mu when this generated number is less
> than pi1. But I cannot include the code u[i] ~ uniform(0,1) in the model section of the STAN code

Why not? All you have to do is declare u as
a parameter instead of data; the same declaration will
work:

real<lower=0,upper=1> u[I];

In fact, with that, you don't even need the u[i] ~ uniform(0,1)
statement because it's implicit in the constraints (all variables
are uniform on constraints unless they have additional sampling
statements or manipulations of lp__).

> so I create a vector
> of uniformly generated numbers first in R and let the vector to be in the data section.

Maybe if you explained what the model's trying to do, it would
be easier to help.

> So I hope someone can enlighten me on how to include the dcat code successfully in the model part in STAN. Attached
> named model_1 is the JAGS code. Thank you! :)

There are other problems you're going to need to
address first, mainly that the following declarations
and sampling statements are inconsistent.

parameters {
...
ordered[J-1] tau;
}
model {
...
tau[1]~normal(0,0.01);
tau[2]~uniform(tau[1],tau[3]);
tau[3]~uniform(tau[2],tau[4]);
tau[4]~uniform(tau[3],tau[5]);
tau[5]~uniform(tau[4],tau[6]);
tau[6]~uniform(tau[5],tau[7]);
tau[7]~uniform(tau[6],tau[8]);
tau[8]~uniform(tau[7],100);

Issue 1: You should either hard-code both J and 8 or
neither. So the sampling could be:

tau[1] ~ normal(0,0.01);
tau[J] ~ uniform(tau[J-1],100);
for (j in 2:(J-1))
tau[j] ~ uniform(tau[J-1],tau[J+1]);

Issue 2: It still won't work, becuase the declared type for
tau places no constraints on the values other than that they're
ordered. So [-8,-7,...,0] is valid according to the declared
constraint. On the other hand, your sampling is more constrained.

The problem that will arise is that you'll go outside of support
for the sampling distributions and everything will get rejected,
or if you're very careful with initialization, you'll just see greatly
reduced sampling times.

The right way to do this is as follows:

parameters {
simplex[J-1] tau_raw;
}
transformed parameters {
ordered[J-1] tau;
tau <- 100 * tau_raw;
}

What this produces is a uniform distribution over tau
with values constrained to [0,100]. You can add
back the

tau[1] ~ normal(0, 0.01);

but it looks like it was probably a hack anyway.

Issue 3: I have no idea what you're trying to accomplish with this
block in the model:

for (i in 1:I) {
if (u[i]<pi1)
mu[i]<- alpha;
else
mu[i]<- -alpha;
p[1] <- 1-Phi(mu[i]-tau[1]);
for(j in 2:(J-1))
p[j] <- Phi(mu[i]-tau[j-1])-Phi(mu[i]-tau[j]);
p[J] <- Phi(mu[i]-tau[J-1]);
y[i] ~ categorical(p);
}

But the solution to Issue 2 should get around any errors that may
have arisen due to automatic rejection.

Did you see that we have an ordered logistic distribution? It's
not probit, but it's close, and handy to use and much faster to
compute (those Phi() calls are expensive compared to inv_logit).

If you're having problem with categorical, then you need to
make sure that the p in categorical(p) is in fact a simplex
(all values non-negative, values sum to 1). You can try printing
it.

- Bob

[Andrew Gelman]

Regarding the code below, I think it is a peculiarity of Bugs that you have to code constrained parameters in that way. So I can see how this confusion arose when translating the model to Stan.  One refreshing thing about Stan is that the lines in the model code can be interpreted directly as augmentations to the log posterior density.

A

--
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/groups/opt_out.

[Newbie (:]

Hi Bob,

 

Thank you so much for helping! :) I have attached the stan R code with inf~categorical(k) in the model block. What I want to do is to let k[1]=1-pi1 (probability of a person being non-infected) and k[2]=pi1 (probability of someone being infected. int inf[I] and inf~categorical(k) so that inf[i] can be either equal 1 or 2. When inf[i]=2, I let mu[i]=alpha and when inf[i]=1, I let mu[i]=-alpha.

 

The abovementioned code works in JAGS: inf[i]~dcat(k[]).

 

But I will get an error after running the R stan code:

 

inf~categorical(k);
 ^-- here

DIAGNOSTIC(S) FROM PARSER:
no matches for function name="categorical_log"
    arg 0 type=real
    arg 1 type=vector
available function signatures for categorical_log:
0.  categorical_log(int, vector) : real
unknown distribution=categorical
Parser expecting: "}"

 

I am not sure whether is it because inf is unknown that cause this error. But I want to let inf to be equal 1 or 2 depending on pi1 (probability of infected).

Btw the issue 3 you have mentioned, the tau in my model actually can be negative and they do not have to sum up to 1. So I let tau[1] to be normally distributed and the rest of the tau to follow a uniform distribution. So tau[1] can actually be negative too.

 

Thank you. =D

 

 

[Bob Carpenter]

On 9/25/13 8:48 AM, Newbie (: wrote:
> Hi Bob,
> Thank you so much for helping! :) I have attached the stan R code with inf~categorical(k) in the model block. What I
> want to do is to let k[1]=1-pi1 (probability of a person being non-infected) and k[2]=pi1 (probability of someone being
> infected. int inf[I] and inf~categorical(k) so that inf[i] can be either equal 1 or 2. When inf[i]=2, I let mu[i]=alpha
> and when inf[i]=1, I let mu[i]=-alpha.
> The abovementioned code works in JAGS: inf[i]~dcat(k[]).
> But I will get an error after running the R stan code:
> inf~categorical(k);
> ^-- here
>
> DIAGNOSTIC(S) FROM PARSER:
> no matches for function name="categorical_log"
> arg 0 type=real
> arg 1 type=vector
> available function signatures for categorical_log:
> 0. categorical_log(int, vector) : real
> unknown distribution=categorical
> Parser expecting: "}"
> I am not sure whether is it because inf is unknown that cause this error. But I want to let inf to be equal 1 or 2
> depending on pi1 (probability of infected).

You just need to define the variable "inf" as an int rather
than a real.

This is a subtle issue in our parser. The name "inf" shouldn't
be used for a variable because the parser wants to parse it as
a numeric literal denoting +infinity. That's why it thinks
the type of "inf" in:

inf~categorical(k);

is real.

We can solve this problem by declaring a different variable:

int foo[I];

And now

for (i in 1:I) {

foo ~ categorical(k);

gives you a more interpretable error message, namely:

EXPECTATION FAILURE LOCATION: file=../temp2/newbie.stan; line=29, column=5

foo ~ categorical(k);

^-- here


DIAGNOSTIC(S) FROM PARSER:
no matches for function name="categorical_log"

arg 0 type=int[1]

arg 1 type=vector
available function signatures for categorical_log:
0. categorical_log(int, vector) : real

What you really wanted to do was say:

for (i in 1:I) {

foo[i] ~ categorical(k);

We like to stick to Gelman's notation where Greek letters
are used for parameters (e.g., theta), lower-case roman letters (e.g., k)
are used for loop variables and indexes and upper-case roman
letters (e.g., K) for un-modeled constants.

> Btw the issue 3 you have mentioned, the tau in my model actually can be negative and they do not have to sum up to 1. So
> I let tau[1] to be normally distributed and the rest of the tau to follow a uniform distribution. So tau[1] can actually
> be negative too.

The way you've declared these variables won't work because
the declared support of the variables still does not match
the sampling constraints. Everything's OK but the upper bound
of 100.

parameters {
...
ordered[8] tau;

}

model {
...
tau[1]~normal(0,0.01);
tau[2]~uniform(tau[1],tau[3]);
tau[3]~uniform(tau[2],tau[4]);
tau[4]~uniform(tau[3],tau[5]);
tau[5]~uniform(tau[4],tau[6]);
tau[6]~uniform(tau[5],tau[7]);
tau[7]~uniform(tau[6],tau[8]);
tau[8]~uniform(tau[7],100);

The particular issue is that tau[8]'s sampling distribution
puts an upper limit of support at 100 whereas there's no such
constraint on the declaration. There's not a good way to code
this constraint into Stan without declaring each variable tau[1]
through tau[8] separately with the appropriate constraints, which
has to be done in reverse order because tau[8] is the key value.

parameters {
real<upper=100> tau8;
real<upper=tau8> tau7;
...
real<upper=tau2> tau1;


This keeps tau8 under 100, and each successive value after that
under the previous one.

I think this will work, but it's ugly.

An alternative is to do the transforms (and Jacobians) yourself,
which is also clunky (not to mention error prone).

- Bob

P.S. I'm opening an issue to put "inf" on our list of reserved
words so you get a more sensible error message up front when
trying to define the variable:

https://github.com/stan-dev/stan/issues/239

You can use the above link to track progress. We'll probably get
this fix in before the official 2.0 release. Beta will be real soon.
We're feature frozen and just finalizing some doc at this point.

P.P.S. Could everyone please start indenting their code so it's more
readable? My standard operating procedure is to cut out of R strings,
move to a file, then let emacs reset it all. But it's for your own
good to do this the "right" way to begin with --- it'll save you huge
headaches in missing where blocks start and end if the visual code
layout matches the logical layout. The result for this model is:

data {
int<lower=0> I; //number of subjects
int<lower=2> J;
int<lower=1,upper=J> y[I];
}

parameters {
real<lower=0,upper=10> alpha;
real<lower=0,upper=1> pi1;
ordered[8] tau;
}

model {
vector[J] p;
real mu[I];
vector[2] k;
int inf[I];

tau[1]~normal(0,0.01);
tau[2]~uniform(tau[1],tau[3]);
tau[3]~uniform(tau[2],tau[4]);
tau[4]~uniform(tau[3],tau[5]);
tau[5]~uniform(tau[4],tau[6]);
tau[6]~uniform(tau[5],tau[7]);
tau[7]~uniform(tau[6],tau[8]);
tau[8]~uniform(tau[7],100);

k[1]<- 1-pi1;
k[2]<- pi1;

for (i in 1:I) {

inf~categorical(k);
if (inf==2) {

[Newbie (:]

Hi Bob,

 

I will still get an error for foo[i]~categorical(k) and I think this is because foo is unknown. It works for dcat in JAGS but I think in Stan, foo has to be known right? The error is as follows:

 

SAMPLING FOR MODEL 'c_code' NOW (CHAIN 1).
Error : Error in function stan::prob::categorical_log(i): Number of categories is 0, but must be between (1, 2)
error occurred during calling the sampler; sampling not done

 

 

For the tau issue, if I remove tau[8]~uniform(tau[7],100) and since tau[8] will always be bigger than the rest so it actually doesn't matter if I don't put a constraint to it.

 

Thank you. =)

[Bob Carpenter]

Well, you need to declare "foo".

On 9/26/13 10:41 AM, Newbie (: wrote:
> Hi Bob,
> I will still get an error for foo[i]~categorical(k) and I think this is because foo is unknown.

You don't need to literally use the name "foo" --- you
just need to change the variable name to something other
than "inf" in both the declaration and model.

- Bob

[Newbie (:]

Hi Bob,

 

I have attached my JAGS code and the Stan code.

 

My problem here is that, running the JAGS code will give my pi1 (probability of infected) to be around 0.2. But running the Stan code, my pi1 will become very very small, about 0.002. Why is this so? The reason that I change to Stan is because JAGS gives me a very wide quantile value ranging from values somewhat 0 to 1 which is impossible. But Stan gives quite a reasonable quantile but now the probability is too small. If hopefully, I can get the same probability that I gotten from JAGS around 0.2, and have a narrow quantile, it will be great! :)

 

 

[Daniel Lee]

Newbie, take a look at the manual and make sure things are parameterized properly. For example, the equivalent prior on the tau's is:

tau[j] ~ normal(0, 10);

not:

tau[j] ~ normal(0, 0.1);

Stan parameterizes the normal distribution with mean and standard deviation.

I'm not sure of the other differences in parameterization between Stan and JAGS. After double-checking the model, can you report whether you get similar results?


Daniel

--

[Newbie (:]

Hi,

 

I realised the normal distribution in JAGS and STAN are different. I have decided to have standard deviation 0.1 means precision 100. In other words, in JAGS, tau~dnorm(0,100)and in STAN, tau~normal(0,0.1). This is actually good because all the parameters are converging very well with good mixing in the traceplots. BUT pi1, the probability of infected become very small abt 2%. I realised that by increasing the standard deviation (reducing precision), the probability (pi1) will increase to 20% which is good but the traceplots are not mixing means not converging very well. Is it possible to have the standard deviation remain as 0.1 to maintain convergence but increase the probability at the same time? Like maybe change to other distribution instead of Uniform(0,1) for the probability?? Or is there any other idea?

 

Btw don't know whether is this infor useful, for pi1, the plot is more to right-skewed instead of bell shape. So maybe there is some ways to do sth to pi1?

[Bob Carpenter]

On 10/6/13 4:21 AM, Newbie (: wrote:
> Hi,
> I realised the normal distribution in JAGS and STAN are different. I have decided to have standard deviation 0.1 means
> precision 100. In other words, in JAGS, tau~dnorm(0,100)and in STAN, tau~normal(0,0.1). This is actually good because
> all the parameters are converging very well with good mixing in the traceplots. BUT pi1, the probability of infected
> become very small abt 2%. I realised that by increasing the standard deviation (reducing precision), the probability
> (pi1) will increase to 20% which is good but the traceplots are not mixing means not converging very well. Is it
> possible to have the standard deviation remain as 0.1 to maintain convergence but increase the probability at the same
> time? Like maybe change to other distribution instead of Uniform(0,1) for the probability?? Or is there any other idea?

I'm not exactly sure what you're asking. The model
is just the model, so you can't have it use one probability
for sampling and another for the probability evaluations.

Can you try running simulated data to see if your
model's doing the right thing?

- Bob

> On Saturday, October 5, 2013 12:16:40 AM UTC+8, Daniel Lee wrote:
>

> stan-users+...@googlegroups.com <javascript:>.
> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.