calculation of loglikelihood value in RStan

[Sedat]

Dear list,

As suggested I just started to use Stan via R due to speed problems with WinBUGS. I read the manual and examples and saw that there are some differences between BUGS syntax and Stan syntax. I converted my BUGS syntax to Stan syntax and I could get IRT 1PL results and confirmed their correctness. When I try to calculate AIC and BIC values through loglikelihood value, I am getting problems. The thing is I cannot calculate loglik value. I may be making a silly mistake. Could you please tell me where or what I am doing wrong in my Stan syntax.

Here's my BUGS code:

#1-group 1PL model

 model 

{

 

#  1group model

for (j in 1:NE) {

for (k in 1:NI) {     

     tt1[j,k]<- exp(theta1[j] - beta1[k])

     p1[j,k]<-tt1[j,k]/(1 + tt1[j,k])

     r1[j,k]~ dbern(p1[j,k])

     l1[j,k]<-log(p1[j,k])*r1[j,k]+log(1-p1[j,k])*(1-r1[j,k])

}

}

     loglik<-sum(l1[1:NE,1:NI])

#  Priors for 1group

for (j in 1:NE) {

theta1[j] ~ dnorm(0, 1)

}

 for (k in 1:NI) { 

    beta1[k]~dnorm(0,1)

 }

par<-NI*1+2*1-1

AIC<--2*loglik+2*par

BIC<--2*loglik+par*log(NE)

}

Here's my Stan code:

{{{

library(Rcpp)

library(inline)

library(rstan) 

}}}

fcat<-'

data {

  int<lower=0> NE; // 2000, number of examinees

   int<lower=0> NI; // 28, number of items

  int<lower=0> r1[NE,NI];

}

transformed data {

real par;

par<-NI+2-1;

}

parameters {

real beta1[NI];

vector[NE] theta1;

}

model {

 matrix[NE,NI] tt1;

matrix[NE,NI] p1;

matrix[NE,NI] l1;

real loglik;

real AIC;

real BIC;

theta1 ~ normal(0, 1);

beta1 ~ normal(0,1);

for (j in 1:NE) {

for (k in 1:NI) {  

     tt1[j,k]<- exp(theta1[j] - beta1[k]);

     p1[j,k]<-tt1[j,k]/(1 + tt1[j,k]);

     r1[j,k]~ bernoulli(p1[j,k]);

     l1[j,k]<-log(p1[j,k])*r1[j,k]+log(1-p1[j,k])*(1-r1 [j,k]);

 }

}

for (j in 1:NE) {

for (k in 1:NI) {  

loglik<-sum(l1); // this is the problematic part ( corresponding to this part  loglik<-sum(l1[1:NE,1:NI]) from BUGS, I tried every combination in parenthesis //like l1[1:j,1:k],  l1[1:NE,1NI], and l1[]. None of these worked. 

}

}

AIC<--2*loglik+2*par;

BIC<--2*loglik+par*log(NE);

     

}

'

{{{

schools_dat <- list(NI = 28, NE = 2000, r1=r1)

fit <- stan(model_code = fcat, data = schools_dat, 

            iter = 1000, chains = 1)

}}} 

After I tried this code in R, it automatically closed program. I could not see any error message. This looks unusual. Could you please help me with that? I am also adding my data file in case you may wanna use it. I will be very appreciated if you can help me with that.

Thanks for your time.

Sedat

[Bob Carpenter]

On 3/9/13 6:28 PM, Sedat wrote:
> ...
> real loglik;
> matrix[NE,NI] l1;
> loglik<-sum(l1);
> this is the problematic part ( corresponding to this part * loglik<-sum(l1[1:NE,1:NI]) *from BUGS,

> I tried every combination in parenthesis //like l1[1:j,1:k], l1[1:NE,1NI], and l1[]. None of these worked.

> ...

Why do you think there's something wrong with sum(l1)? That's
the right syntax in Stan and your model compiles.

Conceptually, you can increment. So if you set

loglike <- 0;

eather than setting entries of l1, you can just increment
loglike.

> After I tried this code in R, it automatically closed program. I could not see any error message. This looks unusual.
> Could you please help me with that? I am also adding my data file in case you may wanna use it. I will be very
> appreciated if you can help me with that.

What do you mean by "automatically closed program"? Which R are
you running? Which version of Stan? Which operating system?

- Bob

[Ben Goodrich]

On Sunday, March 10, 2013 5:25:27 PM UTC-4, Bob Carpenter wrote:

 > After I tried this code in R, it automatically closed program. I could not see any error message. This looks unusual.
 > Could you please help me with that? I am also adding my data file in case you may wanna use it. I will be very
 > appreciated if you can help me with that.

What do you mean by "automatically closed program"?  Which R are
you running?  Which version of Stan?  Which operating system?

I think closed = segfault. If so, then follow

https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/FXHj6fn0SC0J

Ben

[Sedat Sen]

Bob and Ben,

Thanks for your replies.

Let me answer your questions Bob;

What do you mean by "automatically closed program"? 

When I run the code with the data file I sent you, after a while running(around 4-5 minutes) a small screen shows up and says: R for Windows GUI front-end has stopped working. A problem caused the profram to stop working correctly. Windows will close the program and notify you if solution is available.

Which R are you running?

I am using R.2.15.3

Which version of Stan?

rstan (Version 1.2.0, packaged: 2013-03-06 21:47:56 UTC, GitRev: 1cde932d15a3)  

Which operating system? 

I am using Windows 7.

Ben,

I also tried your suggestion and got the same problem again. I was using the level 3optimization in order to make it faster before your suggestion. Following the post, in rstan, I executed set_cppo(mode = "debug") before calling stan(). Nothing has changed. My R is still being closed automatically without showing any error. On the FAQ list you sent, it says model should be as simple as possible. I can get this model run without loglik part. But when I add loglik syntax it gives me this shut down error. I need to be able to calculate this loglik and AIC values for my research. I can already do this with WinBUGS but very slowly. I am attaching my syntax and data separately in case you may wanna try it.

Thanks for your time and help,

Sedat

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

--
SS

[Ben Goodrich]

On Sunday, March 10, 2013 10:42:36 PM UTC-4, Sedat wrote:

Ben,

I also tried your suggestion and got the same problem again. I was using the level 3optimization in order to make it faster before your suggestion. Following the post, in rstan, I executed set_cppo(mode = "debug") before calling stan(). Nothing has changed. My R is still being closed automatically without showing any error. On the FAQ list you sent, it says model should be as simple as possible. I can get this model run without loglik part. But when I add loglik syntax it gives me this shut down error. I need to be able to calculate this loglik and AIC values for my research. I can already do this with WinBUGS but very slowly. I am attaching my syntax and data separately in case you may wanna try it.

The problem is these lines

for (j in 1:NE) {
  for (k in 1:NI) {
    loglik <- sum(l1);
  }
}

It doesn't really make sense to me why you would do a double loop to sum over the l1 matrix. If you just do

loglik <- sum(l1);


it works for me. But with the double loop somehow we get a double free, which is a bug

TRANSLATING MODEL 'fcat' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'fcat' NOW.
SAMPLING FOR MODEL 'fcat' NOW (CHAIN 1).
*** glibc detected *** /usr/lib/R/bin/exec/R: double free or corruption (out): 0x00007fff6fae62d0 ***

...

Program received signal SIGABRT, Aborted.
0x00007ffff70bcdd5 in raise () from /lib/x86_64-linux-gnu/libc.so.6
(gdb) bt
#0  0x00007ffff70bcdd5 in raise () from /lib/x86_64-linux-gnu/libc.so.6
#1  0x00007ffff70bfefb in abort () from /lib/x86_64-linux-gnu/libc.so.6
#2  0x00007ffff70f825e in ?? () from /lib/x86_64-linux-gnu/libc.so.6
#3  0x00007ffff71021b6 in ?? () from /lib/x86_64-linux-gnu/libc.so.6
#4  0x00007ffff1738be2 in stan::agrad::sum<-1, -1> (m=...) at /opt/site-library/rstan/include//stansrc/stan/agrad/matrix.hpp:923
#5  0x00007ffff1737b26 in model61347b7c2652_fcat_namespace::model61347b7c2652_fcat::log_prob (this=0x2326418, params_r__=..., params_i__=..., pstream__=0x7ffff1a70198)
    at file61347e611dd1.cpp:189

Ben

[Sedat Sen]

Ben,

Thanks for your message. It fixed the problem. 

I have another question. I will be very appreciated if you or some one else can give an answer.

I can obtain results for only objects listed in parameters section of Stan syntax. How can I get the results for something, that is transformed parameters, I put in the model section of the syntax? When I put these into parameters section It gives an error which says you cannot declare the same thing twice. Let me be more specific. For example in my syntax I want to get results for AIC and BIC. I put these and their formula into model statement. I could not find a way to see these on my output. I tried to put them into transformed parameters statement just after model statement, RStan did not let me do that. When I put transfromed parameters before model statement with AIC and BIC declaration, it did not calculate the desired statistics because these 2 declarations require some parameters from model and to be declared before these statements.

I hope my question is clear to you. I just want to know how to see estimated or calculated values in the output. What do I have to do make them listed in the output?

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

--
SS

[Bob Carpenter]

On 3/11/13 9:05 PM, Sedat Sen wrote:
> Ben,
>
> Thanks for your message. It fixed the problem.
> I have another question. I will be very appreciated if you or some one else can give an answer.
>
> I can obtain results for only objects listed in parameters section of Stan syntax. How can I get the results for
> something, that is transformed parameters, I put in the model section of the syntax? When I put these into parameters
> section It gives an error which says you cannot declare the same thing twice. Let me be more specific. For example in my
> syntax I want to get results for AIC and BIC. I put these and their formula into model statement. I could not find a way
> to see these on my output. I tried to put them into transformed parameters statement just after model statement, RStan
> did not let me do that. When I put transfromed parameters before model statement with AIC and BIC declaration, it did
> not calculate the desired statistics because these 2 declarations require some parameters from model and to be declared
> before these statements.
>
>
> I hope my question is clear to you. I just want to know how to see estimated or calculated values in the output. What do
> I have to do make them listed in the output?

If you want to print a quantity, put it in
the transformed parameters block if it's used
in the model or in the generated quantities block if
it's not.

This may require you to move other definitions up,
but anything that is calculated as a local variable
in the model block can be calculated in the transformed
parameter block and saved.

If you're still having problems and need specifics, you'll
hae to send us an example.

You should define in transformed data if it only depends
on data, but those values don't get printed.

You should also read Andrew's notes to the list :-)

If you want to send the specific model and data that caused
the problem, just mail them to me as attachments (or to the
list as attachments with suffix .txt).

- Bob

>
>
> On Mon, Mar 11, 2013 at 12:12 AM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
> On Sunday, March 10, 2013 10:42:36 PM UTC-4, Sedat wrote:
>
> Ben,
>
> I also tried your suggestion and got the same problem again. I was using the level 3optimization in order to

> make it faster before your suggestion. Following the post, in rstan, I executed _set_cppo(mode =
> "debug")_ before calling _stan()._ Nothing has changed. My R is still being closed automatically without showing

> /
> /
>
> //

>
> --
> 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 <mailto:stan-users%2Bunsu...@googlegroups.com>.

> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
>
>
> --

> /SS
> /

[Ben Goodrich]

On Monday, March 11, 2013 9:05:45 PM UTC-4, Sedat wrote:

I can obtain results for only objects listed in parameters section of Stan syntax. How can I get the results for something, that is transformed parameters, I put in the model section of the syntax? When I put these into parameters section It gives an error which says you cannot declare the same thing twice. Let me be more specific. For example in my syntax I want to get results for AIC and BIC. I put these and their formula into model statement. I could not find a way to see these on my output. I tried to put them into transformed parameters statement just after model statement, RStan did not let me do that. When I put transfromed parameters before model statement with AIC and BIC declaration, it did not calculate the desired statistics because these 2 declarations require some parameters from model and to be declared before these statements.

It doesn't work to declare something in the transformed parameters block and define it in the model block? If not, you could possibly move the pieces that you need to calculate the AIC or BIC with from model to transformed parameters.

Ben

[Bob Carpenter]

> could possibly move the pieces that you need to calculate the AIC or BIC with from model to transformed parameters.s

Nope -- variables can only be set in the block in
which they're defined.

But you can always move any computation you want to
save into transformed parameters itself, because
even in the model, everything must be defined in
terms of the parameters and data.

- Bob