Re: STAN question?

[Andrew Gelman]

Hi--I'm passing this along to the stan-users mailing list!

A

On Sep 12, 2012, at 7:10 AM, Gianluca Baio wrote:

Dear Andrew,

I'm guessing you get this a lot right now (in addition to all the rest!), so you probably won't have much time to look at this. 

However, I'd really appreciate if you had a moment to possibly clear my doubts about this (I suppose it may even help you debug the function --- assuming I'm not creating the bug!).

I'm trying to run a hierarchical proportional odds model (which we've successfully coded and run in JAGS/BUGS). The dataset is relatively large (with 12000+ individual observations, clustered in 1900+ repeated measurements), so I'd like to test STAN on this one --- it looks like a perfect situation.

The model is like this:

data {

  int<lower=0> N; // number of observed data points (including repeated measures)

  int<lower=0> S; // number of possible scores

  int<lower=0> K; // number of observed combinations for (voter,performer,time)

  int<lower=0> V; // number of observed voters

  int<lower=0> P; // number of observed performers

  int<lower=0> R; // number of "cultural" regions

  int<lower=0,upper=S> y[N]; // observed score for (voter,performer,time)

  int comb[N]; // index for the observed combination (voter,performer,time)

  int x1[N]; // year in which the contest happens

  int x2[N]; // language in which a song is sung (3 categories)

  int x3[N]; // type of artist (3 categories)

  int voter[K];   // voter indicator

  int performer[K]; // performer indicator

  int w[K]; // indicator for geographical proximity between voter and performer

  int index_migr[K]; // indicator for existing data on migration stocks

  real z[K]; // migration stocks

  real phi_fix; // need this to deal with partially observed quantities

  real beta2_fix; // need this to deal with partially observed quantities

  real beta3_fix; // need this to deal with partially observed quantities

  real a[R]; // prior values for the Dirichlet parameters

  int region[K]; // fixed (and arbitrary) "cultural" region indicator

} 

transformed parameters {

  vector[2] phi;

  vector[3] beta2;

  vector[3] beta3;

  vector[R] prob;

  phi[1] <- phi_fix;           // need this trick because phi[1]=0 (if no migration stock data present) and phi[2]~Normal(0,100) (if data on migration stocks present)

  phi[2] <- phi_var;

  beta2[1] <- beta2_fix;       // need this trick because beta2[1]=0 (reference category) and beta2[2,3] ~ Normal(0,100)

  beta2[2] <- beta2_var[1];

  beta2[3] <- beta2_var[2];

  beta3[1] <- beta3_fix;       // need this trick because beta3[1]=0 (reference category) and beta3[2,3] ~ Normal(0,100)

  beta3[2] <- beta3_var[1];

  beta3[3] <- beta3_var[2];

  for (r in 1:R) {

      prob[r] <- d[r]/sum(d);  // trick to code the Dirichlet distribution starting from Gamma iid variables (d[1],...,d[R])

  }

}

model { 

  // Data model --- ordered logistic regression

  for (i in 1:N) {

      y[i] ~ ordered_logistic(alpha[comb[i]]+beta1*x1[i]+beta2[x2[i]]+beta3[x3[i]], lambda);

  }

  

  // Structured effects

  for (k in 1:K) {

      alpha[k] ~ normal(gamma+delta[T[voter[k]],performer[k]]+psi*w[k]+phi[index_migr[k]]*z[k], sd_alpha);

  }

 

  // Tries to infer the "region" to which each voter belongs using a categorical/Dirichlet model

  for (v in 1:V) {

      T[v] ~ categorical(prob);

  }

  for (r in 1:R) {

      d[r] ~ gamma(a[r],1);

  }

  // Structured "regional" effects

  for (r in 1:R) {

      for (p in 1:P)

          delta[r,p] ~ normal(0,sd_delta);

  }

  // Unstructured coefficients

  beta1 ~ normal(0,100);

  for (b in 1:2) {

       beta2_var[b] ~ normal(0,100);

       beta3_var[b] ~ normal(0,100);

  }

  gamma ~ normal(0,100);

  psi ~ normal(0,100);

  phi_var ~ normal(0,100);

  // Standard deviations for the structured effects in the model

  sd_alpha ~ uniform(0,10);

  sd_delta ~ uniform(0,10);

}

In other words, we have an ordered logistic regression, where the linear predictor is defined as the sum of some unstructured effects (scalar beta1 and vectors beta2 and beta3) and a structured effect (alpha) which accounts for clustering at the combination level (repeated measures of voter/performer at different time points). This structured effect alpha is then modelled using normal distributions with mean given by a residual gamma, two unstructured effects psi and phi and a structured effect delta. This is nested within the region to which voter k belongs and the performer. 

In the first version of the model, I assume that the region is arbitrarily fixed (and thus in the code above I modify T[voter[k]] with region[k], which is then an observed variable). This version of the model seems to do OK.

In the second version (more interesting and effectively the one reported above), I assume that the regions are not codified a priori, but rather use the (little and indirect) evidence from the voting patterns to cluster the voters in R groups. In JAGS/BUGS we've done that by including a variable T[v] ~ categorical(probs) with probs ~ Dirichlet(a[1:R]) and a[1:R] given --- which is what I was trying to do in the code above.

I'm not sure: a) whether you could in fact do something like that in STAN; and b) whether the strategy I've used is actually correct. However, when I run the second version of the model, it seems to compile OK, but when it comes to sampling I get an error message like 

 *** caught segfault ***

address 0x170000000e, cause 'memory not mapped'

Traceback:

 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x1c27030>,     dll = list(name = "Rcpp", path = "/usr/local/lib/R/site-library/Rcpp/libs/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x273e240>,         info = <pointer: 0x7faff193d090>), numParameters = -1L),     <pointer: 0x63fd4f0>, <pointer: 0x4371a40>, .pointer, ...)

 2: sampler$call_sampler(args_list[[i]])

 3: .local(object, ...)

 4: sampling(sm, data, pars, n_chains, iter, warmup, thin, seed,     init, sample_file = sample_file, verbose = verbose, check_data = FALSE,     ...)

 5: sampling(sm, data, pars, n_chains, iter, warmup, thin, seed,     init, sample_file = sample_file, verbose = verbose, check_data = FALSE,     ...)

 6: stan(file = "model.stan", data = EF, pars = pars, iter = 2000,     chains = 2)

The actual address of the segfault changes every time I run the code, but the Traceback seems to consistently produce the same message.

Am I missing something here, or is my code fundamentally wrong? Apologies if utterly trivial!

Thanks for your help,

Best

Gianluca

--
Gianluca Baio

Lecturer in Medical Statistics
Department of Statistical Science
University College London
1-19 Torrington Place, London WC1E 6BT, UK
Telephone +44(0)20 7679 1248
Fax +44 (0)20 7383 4703
Website: http://www.ucl.ac.uk/statistics/people/gianlucabaio

Blog: http://gianlubaio.blogspot.co.uk

[Jiqiang Guo]

I will take a look. This kind of segfault from Rcpp has been haunted me for a long time, but I though I have taken care of it. 

Gianluca, if you do not mind, could you please send us (or me) the data somehow?  Otherwise, I will simulate some data, which takes long. 

Jiqiang 

[Bob Carpenter]

The seg fault is probably coming from the generated
model. But that's not the root cause of the model's
problem (see below for at least one problem in the model
specification).

We should have eliminated all the seg faults from
the models by the next release. The two places they
can occur now (version 1.0.1) when there are indexes out of
bounds of the arrays/vectors/matrices on the left-hand side
of assignment statements (<-) or when matrix operations are
applied to the wrong size arguments (e.g., trying to multiply
a 3x2 matrix by a 5x7 matrix).

Gianluca -- the code in your transformed parameters
block is wrong. You have:

>> transformed parameters {
>> vector[2] phi;
>> vector[3] beta2;
>> vector[3] beta3;
>> vector[R] prob;
>>
>> phi[1] <- phi_fix; // need this trick because phi[1]=0 (if no migration stock data present) and
>> phi[2]~Normal(0,100) (if data on migration stocks present)
>> phi[2] <- phi_var;

The problem is that Stan does not work like BUGS (see
the new appendix for an overview of some differences).
When Stan execute the second statement,

phi[2]~Normal(0,100) (if data on migration stocks present)

phi[2] has not yet been defined (we initialize to NaN).
So the result is going to be

lp__ <- lp__ + normal_log(NaN,0,100);

so that the log probability will become NaN and never recover.

That's not the cause of the seg fault, but it's definitely
a bug in the model specification.

We try to describe how the transformed parameters work and
the parameters work in the manual.

Several of these user examples have been very
complex. I'm thinking we need some way of debugging
models step-wise. It's a big project, but it might
pay off with tracking down this kind of bug, which should
be pretty obvious under stepwise evaluation.

Jiqiang -- I don't think we need to go to the trouble
of writing simulations for user models to help them
with debugging. Great if you have the time, but we
also have other things to do!

- Bob

On 9/12/12 4:56 PM, Jiqiang Guo wrote:
> I will take a look. This kind of segfault from Rcpp has been haunted me for a long time, but I though I have taken care
> of it.
>
> Gianluca, if you do not mind, could you please send us (or me) the data somehow? Otherwise, I will simulate some data,
> which takes long.
>
> Jiqiang
>

[Jiqiang Guo]

Bob is right. The model is not correctly specified.  Actually, I do not even see a parameters block.  Please follow Bob's suggestions and see if you still have the segfault, which though seemingly is because of somehow the Rcpp cannot find the loaded dynamic shared object compiled from the C++ code for the model. 

Jiqiang

[Jiqiang Guo]

See comments below for one questions raised earlier.

Jiqiang

On Sep 12, 2012, at 7:10 AM, Gianluca Baio wrote:

In the second version (more interesting and effectively the one reported above), I assume that the regions are not codified a priori, but rather use the (little and indirect) evidence from the voting patterns to cluster the voters in R groups. In JAGS/BUGS we've done that by including a variable T[v] ~ categorical(probs) with probs ~ Dirichlet(a[1:R]) and a[1:R] given --- which is what I was trying to do in the code above.

Though the parameters block is not given in the mail, it seems that you want to specify T here as parameters, which then have a *discrete* distribution.  Unfortunately, for now Stan does not have a sampler for discrete parameters, so this second version of your model cannot be run in Stan.  

[Bob Carpenter]

Stan 1.0.1 has a test in the model parser that will
stop and report an error if the model attempts to
define discrete parameters.

We're planning to add a discrete sampler soon, but
not one that'll be particularly efficient. So we'll
still be recommending marginalizing out discrete
parameters wherever possible. There are examples in
the manual for the typical case of mixture models.

- Bob