change point model

[Song]

I went through a few BUGS examples in the Stan home page to understand the language better. The change point problem in the Vol2 "stagnant" example is of special interest to me.  This example was labeled as an example of "how not to do MCMC." In a previous data analysis work, I used the change point model to analyze similar data for 4 subjects (see attached PDF). I can run the change point model on these subjects individually using Stan (modelfile1.txt and dumpdata1.Rdata).  But when I tried to fit a multilevel model (assuming exchangeability on model coefficients) (modelfile2.txt and dumpdata2.Rdata), I got the following error. Thank you in advance for your help.

SAMPLING FOR MODEL 'input.to.bugs$model' NOW (CHAIN 1).

 *** caught segfault ***

address 0xb7e92f6c, cause 'memory not mapped'

Traceback:

 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x2806ea0>,     dll = list(name = "Rcpp", path = "/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/libs/i386/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x28106c0>,         info = <pointer: 0x2ccad0>), numParameters = -1L), <pointer: 0x53a8a40>,     <pointer: 0x5385490>, .pointer, ...)

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

 3: .local(object, ...)

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

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

 6: stan(model_code = input.to.bugs$model, data = input.to.bugs$data,     pars = input.to.bugs$para, iter = n.iter, thin = thin, chains = n.chains)

Possible actions:

1: abort (with core dump, if enabled)

2: normal R exit

3: exit R without saving workspace

4: exit R saving workspace

Selection: 

Song Qian

Department of Environmental Sciences

University of Toledo

[Ben Goodrich]

As of rstan 1.0.2, there is a set_cppo() function whose argument can be set to mode = "debug" to compile with less optimization, debug symbols, and more range-checking from Eigen, which would point you toward the problem

In this case, when I run it under gdb without optimization, there is an out-of-range message referring to stid

goodrich@CYBERPOWERPC:/tmp/Song$ gdb --args ./modelfile2 --data=dumpdata2.Rdata
GNU gdb (GDB) 7.4.1-debian
Copyright (C) 2012 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>...
Reading symbols from /tmp/Song/modelfile2...done.
(gdb) run
Starting program: /tmp/Song/modelfile2 --data=dumpdata2.Rdata
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".

Exception: INDEX OPERATOR [] OUT OF BOUNDS; index=5; lower bound=1; upper bound=4; index position=1; stid

Diagnostic information:
Dynamic exception type: std::out_of_range
std::exception::what: INDEX OPERATOR [] OUT OF BOUNDS; index=5; lower bound=1; upper bound=4; index position=1; stid


[Inferior 1 (process 14886) exited with code 0377]
(gdb) q

Your modelfile2.txt file has

data {
  int<lower=0> n;
  int<lower=0> nst;
  real y[n];
  int<lower=0> stid[nst]; // station ID
  ...
}

but I think you meant int<lower=0> stid[n] or at least that is what is consistent with dumpdata2.Rdata .

Other notes as code comments

parameters{
  real beta[nst];
  real<lower=0> delta;
  real<lower=minYr, upper=maxYr> phi[nst];
  real mu0;
  real mu1;
  real<lower=0> sigmay;  // real<lower=0,upper=100> sigmay;
  real<lower=0> sigma[2]; // real<lower=0,upper=100> sigma[2];
}
model{
        mu0 ~ normal(0, 100);
        mu1 ~ normal(0, 100);
        for (j in 1:nst){ // better to vectorize these instead of a for loop
            beta[j] ~ normal(mu0, sigma[1]); // beta ~ normal(mu0, sigma[1]);
            phi[j] ~ normal(mu1,sigma[2]); // phi ~ normal(mu1, sigma[2]);
        }
        delta ~ normal(0,100);
        for (k in 1:2) { // not necessary when support of sigma is bounded on both sides
            sigma[k] ~ uniform(0, 100);
        }
        sigmay ~ uniform(0, 100); // not necessary when support of sigmay is bounded on both sides
        for (i in 1:n) { // sdid needs to be of length n
            y[i] ~ normal(beta[stid[i]] - delta * int_step(x[i] - phi[stid[i]])*(x[i]-phi[stid[i]]), sigmay);
        }
}

Ben

[Song S. Qian]

Thank you very much. stid[nst] was the problem.  

Song