Choice model

[M.F. Jonker]

Hi Daniel/Bob,

It took a while but I completed a simple version of my model. It works well (see the attached code and data) but I do have several questions about the code. Can I ask them on the forum?

a) it appears that print commands in the "generated quantities" section still do not work; is this correct?
   I'd like to use them for debugging and for printing a few key parameters during estimation.
b) the compiled command-line print.exe sometimes crashes when I try to print the samples.csv. I don't use R, so this function is GREAT! I can attach one of the problematic files?
c) my pilot model currently uses a subset of ~600 respondents; even without respondent-specific beta parameters, there are ~1,200 parameters that need not be monitored (i.e. the s[] and s_star[] parameters). With respondent-specific parameters, there are ~ 14,000 parameters that need not be monitored. And with a real-life dataset with many more respondents, even many many more. So, is a new STAN feature that allows several parameters to be excluded from the samples.csv output a possibility? I've briefly mentioned this before, apologies for the duplication.
d) using an array in the likelihood declaration does not work (yet)?
e) I've tried to replace the poisson with a categorical distribution. In BUGS this is coded as follows:

  for (n in 1:N){
   for (t in 1:T){
    Y[n,t] ~ dcat(prob[n, t, 1:A])   
  }}
 
   I've tried several options, but always got an error in the calculation of the probabilities when I initialized them as: vector[A] prob[N,T], which is -I think- the required form(?). A suggestion would be much appreciated!

f) the multiplication of s[n] and beta[] is not fully identified, hence the mean-of-1 constraint on the s[n]. I don't understand, however, why the posterior mean of the s_star[n] are so much smaller than 1. Did I specify the constraint on mu_s incorrectly, or perhaps any other idea why?

Thank you very much,
Marcel

[Jiqiang Guo]

It would be to nice to give the detailed error messages when you say it is not working, and/or crashes. 

Now Stan would output all the parameters, derived parameters, and generated quantities in the csv files. 

I think you should look at your "s" (not s_star) to see if they look like what you might expect. 

--

Jiqiang 

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

[Bob Carpenter]

On 3/23/13 8:44 AM, M.F. Jonker wrote:

> It took a while but I completed a simple version of my model. It works well (see the attached code and data) but I do
> have several questions about the code. Can I ask them on the forum?

Ask away...

> a) it appears that print commands in the "generated quantities" section still do not work; is this correct?
> I'd like to use them for debugging and for printing a few key parameters during estimation.

They should work in version 1.2.0.

> b) the compiled command-line print.exe sometimes crashes when I try to print the samples.csv. I don't use R, so this
> function is GREAT! I can attach one of the problematic files?

Yes, please send us something that crashes.

> c) my pilot model currently uses a subset of ~600 respondents; even without respondent-specific beta parameters, there
> are ~1,200 parameters that need not be monitored (i.e. the s[] and s_star[] parameters). With respondent-specific
> parameters, there are ~ 14,000 parameters that need not be monitored. And with a real-life dataset with many more
> respondents, even many many more. So, is a new STAN feature that allows several parameters to be excluded from the
> samples.csv output a possibility? I've briefly mentioned this before, apologies for the duplication.

As I believe we've replied before, we're reluctant to do this
because of the implications for convergence monitoring.

> d) using an array in the likelihood declaration does not work (yet)?
> e) I've tried to replace the poisson with a categorical distribution. In BUGS this is coded as follows:
>
> for (n in 1:N){
> for (t in 1:T){
> Y[n,t] ~ dcat(prob[n, t, 1:A])
> }}
>
> I've tried several options,

You can use

Y[n,t] ~ categorical(prob[n,t]);

in Stan.

but always got an error in the calculation of the probabilities when I initialized them as:
> vector[A] prob[N,T], which is -I think- the required form(?). A suggestion would be much appreciated!

You want a simplex if it's categorical probability.
Stan requires the constraints to be specified. So that's

simplex[A] prob[N,T];

I'd stay away from using variables named "prob" in a model!

> f) the multiplication of s[n] and beta[] is not fully identified, hence the mean-of-1 constraint on the s[n]. I don't
> understand, however, why the posterior mean of the s_star[n] are so much smaller than 1. Did I specify the constraint on
> mu_s incorrectly, or perhaps any other idea why?

I don't think the following does what you intend:

parameters {
...
vector[N] s_star;
}
transformed parameters{
vector[N] s;
s <- s_star/mean(s_star);
}

model {
...
s_star ~ lognormal(-0.5 * square(sigma_s) , sigma_s);

First, you probably want a prior on sigma_s -- it's unbounded here.
I'd then expect sigma_s to go to zero anyway given the way you
set up the sampling distribution on s_star. But what you really
want is something like a prior on the transformed parameter. But
to do that you'll need to adjust for the Jacobian.

Or find some other way to identify the model.

- Bob

[M.F. Jonker]

Dear Bob,

Thank you very much :) I'll give it a try and also attach one of the problematic sample.csv files.

All the best,
Marcel

[M.F. Jonker]

Dear list,

I've been trying to use the softmax function (section 28.5) to create the required input for the categorical distribution.  Somehow it doesn't work (the parser error says: "no matches for function name="softmax" - arg 0 type = real[1] - expression is ill formed - expecting <eps>). Perhaps someone can explain what I did wrong /how it is done properly? :)

data {
 ...
  int<lower=1,upper=A>   Y[N,T];
  vector<lower=0>[V]   X[N,T,A];
}

model {
  real lambda[N,T,A];

 
  for (n in 1:N){  
   for (t in 1:T){  

     for (a in 1:A)
       lambda[n,t,a] <- alpha[a] + s[n] * dot_product(X[n,t,a], beta);
     Y[n,t] ~ categorical( softmax(lambda[n,t]) );
  }}
  ...
}

[Bob Carpenter]

softmax requires a vector argument. Here's the signature
in the index:

softmax
(vector x):vector, 209

You gave it an array. If you declare

vector[A] lambda[N,T];

then softmax(lambda[n,t]) will be well-formed and
return a simplex.

Almost all of the multivariate probability functionality
is defined in terms of vectors. We may extend this to
arrays in the future, because the type signatures are
confusing people.

- Bob

[Marcel Jonker]

Hi Bob,

With the correct specification of lambda, the model works! Now I've been trying to vectorize the calculations.
In Matlab, I can specify a dot-product for 2 vectors:

A = [1 2 3];
B = [4 5 6];
dot(A,B)
-> 32

but also for 2 matrices (the result is a vector):

A = [1 2 3; 1 2 3; 1 2 3];
B = [4 5 6; 4 5 6; 4 5 6];
dot(A,B,2)
-> 32
    32
    32

The latter doesn't work in STAN. Is there perhaps another way to vectorize the below loops?

data {
  ...
  int<lower=1,upper=A> Y[N,T];
  vector<lower=0>[V] X[N,T,A];
}

parameters {
  vector[V] beta[N];
  vector[N] s;
}

model {
  vector[A] lambda[N,T];
 
  # likelihood

  for (n in 1:N){
   for (t in 1:T){
    for (a in 1:A)

      lambda[n,t,a] <- s[n] * dot_product(X[n,t,a], beta[n]);

    Y[n,t] ~ categorical( softmax(lambda[n,t]) );
  }}

Thank you very much (again)!

Op dinsdag 26 maart 2013 00:27:00 UTC+1 schreef Bob Carpenter het volgende:

[Bob Carpenter]

I put a dot-product of matrices on our to-do list.
It's easy to implement.

We haven't yet vectorized categorical. We really haven't
vectorized any of the distributions that have multivariate
parameters or outcomes. Once we get everything fully
vectorized, you'll be able to write

Y[n,t] ~ categorical( softmax(lambda[n,t]) );

as

Y ~ categorical(softmax(lambda));


Right now, you could vectorize this further:

> lambda[n,t,a] <- s[n] * dot_product(X[n,t,a], beta[n]);

but I don't really think it'll be worth doing.

In general, you can convert categorical to multinomial,
but it's not going to help you here because you have
varying distributions going into each call to categorical().

- Bob

[M.F. Jonker]

Thanks! I'll wait for the factorized truncated log-normal and categorical distributions then ;)
A little off-topic but almost unbelievable how fast you guys keep developing STAN.. it's brilliant and keeps getting faster and more versatile!

Best,
Marcel

> Date: Tue, 26 Mar 2013 13:57:17 -0400
> From: ca...@alias-i.com
> To: stan-...@googlegroups.com
> Subject: Re: [stan-users] Choice model