try RStan on Mac

[Jiqiang Guo]

For Mac users, please try install RStan as 

install.packages('rstan', repos = 'http://wiki.stan.googlecode.com/git/R')

If this is successful, we could try to use RStan as 

library(rstan)

schools.code <- '

 data {

   int[0,] J; // number of schools

    real y[J]; // estimated treatment effects

   real[0,] sigma_y[J]; // s.e. of effect estimates

  }

  parameters {

   real theta[J];

    real mu;

    real[0,] tau;

  }

  model {

   theta ~ normal(mu, tau);

    y ~ normal(theta, sigma_y);

 }

'

dat <- list(J = 8,

            y = c(28,  8, -3,  7, -1,  1, 18, 12),

           sigma_y = c(15, 10, 16, 11,  9, 11, 10, 18))

fit <- stan(model.code = schools.code, data = dat,

            n.iter = 2012, n.chains = 4)

The above should work, but as I have not got my Mac to try it.  If it goes well, we are very close to let Andrew try it. 

Thanks,

Jiqiang 

[Bob Carpenter]

Looks like they need more path info. Here's the error
message I get:

> install.packages('rstan', repos = 'http://wiki.stan.googlecode.com/git/R')

Warning: unable to access index for repository http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.14
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
package ‘rstan’ is not available (for R version 2.14.2)

I fired up a browser and see that the following path (with a 2.15
rather than 2.14) exists:

http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.15/

I'm downloading the 2.15 version of R now from CRAN and will let
you know how that goes.

Is there some way to make this work without requiring an update to
R itself?

- Bob

[Bob Carpenter]

It seemed to download and install fine this time (with R 2.15.1).
But then it crashes:

> library(rstan)
Error: package ‘Rcpp’ required by ‘rstan’ could not be found

Then I had the same problem with:

Rcpp
RcppEigen
inline

I've seen some packages auto-install other packages. Can/should
we do this?

Then, when they were all installed, the compiler crashed on
the model. Here's the message I got:

> fit <- stan(model.code = schools.code, data = dat,

+ n.iter = 2012, n.chains = 4)

COMPILING MODEL 'anon_model' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'anon_model' NOW.
filef8ab3cb9ddf2.cpp:6:30: error: rstan/rstaninc.hpp: No such file or directory
filef8ab3cb9ddf2.cpp:9:39: error: stan/model/model_header.hpp: No such file or directory
filef8ab3cb9ddf2.cpp:13: error: ‘std::vector’ has not been declared
filef8ab3cb9ddf2.cpp:14: error: ‘std::string’ has not been declared
filef8ab3cb9ddf2.cpp:15: error: ‘std::stringstream’ has not been declared
filef8ab3cb9ddf2.cpp:16: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:17: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:18: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:19: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:20: error: ‘std::istream’ has not been declared
filef8ab3cb9ddf2.cpp:21: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:21: error: ‘math’ is not a namespace-name
filef8ab3cb9ddf2.cpp:21: error: expected namespace-name before ‘;’ token
filef8ab3cb9ddf2.cpp:22: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:22: error: ‘prob’ is not a namespace-name
filef8ab3cb9ddf2.cpp:22: error: expected namespace-name before ‘;’ token
filef8ab3cb9ddf2.cpp:23: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:23: error: ‘agrad’ is not a namespace-name
filef8ab3cb9ddf2.cpp:23: error: expected namespace-name before ‘;’ token
filef8ab3cb9ddf2.cpp:25: error: ‘Eigen’ has not been declared
filef8ab3cb9ddf2.cpp:25: error: expected initializer before ‘<’ token
filef8ab3cb9ddf2.cpp:26: error: ‘Eigen’ has not been declared
filef8ab3cb9ddf2.cpp:26: error: expected initializer before ‘<’ token
filef8ab3cb9ddf2.cpp:27: error: ‘Eigen’ has not been declared
filef8ab3cb9ddf2.cpp:27: error: expected initializer before ‘<’ token
filef8ab3cb9ddf2.cpp:28: error: ‘Eigen’ has not been declared
filef8ab3cb9ddf2.cpp:28: error: expected initializer before ‘<’ token
filef8ab3cb9ddf2.cpp:29: error: ‘Eigen’ has not been declared
filef8ab3cb9ddf2.cpp:29: error: expected initializer before ‘<’ token
filef8ab3cb9ddf2.cpp:30: error: ‘Eigen’ has not been declared
filef8ab3cb9ddf2.cpp:30: error: expected initializer before ‘<’ token
filef8ab3cb9ddf2.cpp:32: error: expected class-name before ‘{’ token
filef8ab3cb9ddf2.cpp:35: error: ISO C++ forbids declaration of ‘vector’ with no type
filef8ab3cb9ddf2.cpp:35: error: expected ‘;’ before ‘<’ token
filef8ab3cb9ddf2.cpp:36: error: ISO C++ forbids declaration of ‘vector’ with no type
filef8ab3cb9ddf2.cpp:36: error: expected ‘;’ before ‘<’ token
filef8ab3cb9ddf2.cpp:38: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:38: error: expected `)' before ‘&’ token
filef8ab3cb9ddf2.cpp:311: error: expected `}' at end of input
filef8ab3cb9ddf2.cpp:311: error: expected unqualified-id at end of input
filef8ab3cb9ddf2.cpp:311: error: expected `}' at end of input
make: *** [filef8ab3cb9ddf2.o] Error 1

ERROR(s) during compilation: source code errors or compiler configuration errors!

Program source:
1:
2: // includes from the plugin
3:
4:
5: // user includes
6: #include <rstan/rstaninc.hpp>
7: // Code generated by Stan version 1.0
8:
9: #include <stan/model/model_header.hpp>
10:
11: namespace stan_model_namespace {
12:
13: using std::vector;
14: using std::string;
15: using std::stringstream;
16: using stan::agrad::var;
17: using stan::model::prob_grad_ad;
18: using stan::math::get_base1;
19: using stan::io::dump;
20: using std::istream;
21: using namespace stan::math;
22: using namespace stan::prob;
23: using namespace stan::agrad;
24:
25: typedef Eigen::Matrix<double,Eigen::Dynamic,1> vector_d;
26: typedef Eigen::Matrix<double,1,Eigen::Dynamic> row_vector_d;
27: typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> matrix_d;
28: typedef Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,1> vector_v;
29: typedef Eigen::Matrix<stan::agrad::var,1,Eigen::Dynamic> row_vector_v;
30: typedef Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,Eigen::Dynamic> matrix_v;
31:
32: class stan_model : public prob_grad_ad {
33: private:
34: int J;
35: vector<double> y;
36: vector<double> sigma_y;
37: public:
38: stan_model(stan::io::var_context& context__)
39: : prob_grad_ad::prob_grad_ad(0) {
40: static const char* function__ = "stan_model_namespace::stan_model(%1%)";
41: size_t pos__;
42: std::vector<int> vals_i__;
43: std::vector<double> vals_r__;
44: J = int(0);
45: if (!context__.contains_i("J"))
46: throw std::runtime_error("variable J not found.");
47: vals_i__ = context__.vals_i("J");
48: pos__ = 0;
49: J = vals_i__[pos__++];
50: y = std::vector<double>(J,double(0));
51: if (!context__.contains_r("y"))
52: throw std::runtime_error("variable y not found.");
53: vals_r__ = context__.vals_r("y");
54: pos__ = 0;
55: size_t y_limit_0__ = J;
56: for (size_t i_0__ = 0; i_0__ < y_limit_0__; ++i_0__) {
57: y[i_0__] = vals_r__[pos__++];
58: }
59: sigma_y = std::vector<double>(J,double(0));
60: if (!context__.contains_r("sigma_y"))
61: throw std::runtime_error("variable sigma_y not found.");
62: vals_r__ = context__.vals_r("sigma_y");
63: pos__ = 0;
64: size_t sigma_y_limit_0__ = J;
65: for (size_t i_0__ = 0; i_0__ < sigma_y_limit_0__; ++i_0__) {
66: sigma_y[i_0__] = vals_r__[pos__++];
67: }
68: // validate data
69: check_greater_or_equal(function__,J,0,"J");
70: for (int k0__ = 0; k0__ < J; ++k0__) {
71: check_greater_or_equal(function__,sigma_y[k0__],0,"sigma_y[k0__]");
72: }
73:
74: // validate transformed data
75:
76: set_param_ranges();
77: } // dump ctor
78:
79: void set_param_ranges() {
80: num_params_r__ = 0U;
81: param_ranges_i__.clear();
82: num_params_r__ += J;
83: ++num_params_r__;
84: ++num_params_r__;
85: }
86:
87: void transform_inits(const stan::io::var_context& var_context__,
88: std::vector<int>& params_i__,
89: std::vector<double>& params_r__) {
90: params_r__.clear();
91: params_i__.clear();
92: stan::io::writer<double> writer__(params_r__,params_i__);
93: size_t pos__;
94: std::vector<double> vals_r__;
95: std::vector<int> vals_i__;
96:
97:
98: if (!(var_context__.contains_r("theta")))
99: throw std::runtime_error("variable theta missing");
100: if (var_context__.dims_r("theta").size() != 1)
101: throw std::runtime_error("require 1 dimensions for variable theta");
102: vals_r__ = var_context__.vals_r("theta");
103: pos__ = 0U;
104: std::vector<double> theta(J,double(0));
105: for (int i0__ = 0U; i0__ < J; ++i0__)
106: theta[i0__] = vals_r__[pos__++];
107: for (int i0__ = 0U; i0__ < J; ++i0__)
108: writer__.scalar_unconstrain(theta[i0__]);
109:
110: if (!(var_context__.contains_r("mu")))
111: throw std::runtime_error("variable mu missing");
112: if (var_context__.dims_r("mu").size() != 0)
113: throw std::runtime_error("require 0 dimensions for variable mu");
114: vals_r__ = var_context__.vals_r("mu");
115: pos__ = 0U;
116: double mu(0);
117: mu = vals_r__[pos__++];
118: writer__.scalar_unconstrain(mu);
119:
120: if (!(var_context__.contains_r("tau")))
121: throw std::runtime_error("variable tau missing");
122: if (var_context__.dims_r("tau").size() != 0)
123: throw std::runtime_error("require 0 dimensions for variable tau");
124: vals_r__ = var_context__.vals_r("tau");
125: pos__ = 0U;
126: double tau(0);
127: tau = vals_r__[pos__++];
128: writer__.scalar_lb_unconstrain(0,tau);
129: params_r__ = writer__.data_r();
130: params_i__ = writer__.data_i();
131: }
132:
133: var log_prob(vector<var>& params_r__,
134: vector<int>& params_i__) {
135:
136: var lp__(0.0);
137: // model parameters
138: stan::io::reader<var> in__(params_r__,params_i__);
139:
140: vector<var> theta;
141: size_t dim_theta_0__ = J;
142: theta.reserve(dim_theta_0__);
143: for (size_t k_0__ = 0; k_0__ < dim_theta_0__; ++k_0__) {
144: theta.push_back(in__.scalar_constrain(lp__));
145: }
146: var mu = in__.scalar_constrain(lp__);
147: var tau = in__.scalar_lb_constrain(0,lp__);
148:
149: // transformed parameters
150:
151:
152: // validate transformed parameters
153:
154: // model body
155: lp__ += stan::prob::normal_log<true>(theta, mu, tau);
156: lp__ += stan::prob::normal_log<true>(y, theta, sigma_y);
157:
158: return lp__;
159:
160: } // log_prob()
161:
162:
163: void get_param_names(std::vector<std::string>& names__) {
164: names__.resize(0);
165: names__.push_back("theta");
166: names__.push_back("mu");
167: names__.push_back("tau");
168: }
169:
170:
171: void get_dims(std::vector<std::vector<size_t> >& dimss__) {
172: dimss__.resize(0);
173: std::vector<size_t> dims__;
174: dims__.resize(0);
175: dims__.push_back(J);
176: dimss__.push_back(dims__);
177: dims__.resize(0);
178: dimss__.push_back(dims__);
179: dims__.resize(0);
180: dimss__.push_back(dims__);
181: }
182:
183: void write_array(std::vector<double>& params_r__,
184: std::vector<int>& params_i__,
185: std::vector<double>& vars__) {
186: vars__.resize(0);
187: stan::io::reader<double> in__(params_r__,params_i__);
188: static const char* function__ = "stan_model_namespace::write_array(%1%)";
189: // read-transform, write parameters
190: vector<double> theta;
191: size_t dim_theta_0__ = J;
192: for (size_t k_0__ = 0; k_0__ < dim_theta_0__; ++k_0__) {
193: theta.push_back(in__.scalar_constrain());
194: }
195: double mu = in__.scalar_constrain();
196: double tau = in__.scalar_lb_constrain(0);
197: for (int k_0__ = 0; k_0__ < J; ++k_0__) {
198: vars__.push_back(theta[k_0__]);
199: }
200: vars__.push_back(mu);
201: vars__.push_back(tau);
202:
203: // declare and define transformed parameters
204: double lp__ = 0.0;
205:
206:
207: // validate transformed parameters
208:
209: // write transformed parameters
210:
211: // declare and define generated quantities
212:
213:
214: // validate generated quantities
215:
216: // write generated quantities
217: }
218:
219:
220: void write_csv_header(std::ostream& o__) {
221: stan::io::csv_writer writer__(o__);
222: for (int k_0__ = 1; k_0__ <= J; ++k_0__) {
223: writer__.comma();
224: o__ << "theta" << '.' << k_0__;
225: }
226: writer__.comma();
227: o__ << "mu";
228: writer__.comma();
229: o__ << "tau";
230: writer__.newline();
231: }
232:
233: void write_csv(std::vector<double>& params_r__,
234: std::vector<int>& params_i__,
235: std::ostream& o__) {
236: stan::io::reader<double> in__(params_r__,params_i__);
237: stan::io::csv_writer writer__(o__);
238: static const char* function__ = "stan_model_namespace::write_csv(%1%)";
239: // read-transform, write parameters
240: vector<double> theta;
241: size_t dim_theta_0__ = J;
242: for (size_t k_0__ = 0; k_0__ < dim_theta_0__; ++k_0__) {
243: theta.push_back(in__.scalar_constrain());
244: writer__.write(theta[k_0__]);
245: }
246: double mu = in__.scalar_constrain();
247: writer__.write(mu);
248: double tau = in__.scalar_lb_constrain(0);
249: writer__.write(tau);
250:
251: // declare, define and validate transformed parameters
252: double lp__ = 0.0;
253:
254:
255:
256: // write transformed parameters
257:
258: // declare and define generated quantities
259:
260:
261: // validate generated quantities
262:
263: // write generated quantities
264: writer__.newline();
265: }
266:
267: }; // model
268:
269: } // namespace
270:
271: /**
272: * Define Rcpp Module to expose stan_fit's functions to R.
273: */
274: RCPP_MODULE(stan_model){
275: Rcpp::class_<rstan::stan_fit<stan_model_namespace::stan_model,
276: boost::random::ecuyer1988> >("stan_model")
277: // .constructor<Rcpp::List>()
278: .constructor<SEXP>()
279: .constructor<SEXP, SEXP>()
280: .method("call_sampler",
281: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::call_sampler)
282: .method("param_names",
283: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::param_names)
284: .method("param_names_oi",
285: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::param_names_oi)
286: .method("param_fnames_oi",
287: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::param_fnames_oi)
288: .method("param_dims",
289: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::param_dims)
290: .method("param_dims_oi",
291: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::param_dims_oi)
292: .method("permutation",
293: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::permutation)
294: .method("update_param_oi",
295: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::update_param_oi)
296: .method("param_oi_tidx",
297: &rstan::stan_fit<stan_model_namespace::stan_model, boost::random::ecuyer1988>::param_oi_tidx)
298: ;
299: }
300:
301:
302: // declarations
303: extern "C" {
304: SEXP filef8ab3cb9ddf2( ) ;
305: }
306:
307: // definition
308:
309: SEXP filef8ab3cb9ddf2( ){
310: return R_NilValue;
311: }
312:
313:
Error in compileCode(f, code, language = language, verbose = verbose) :
Compilation ERROR, function(s)/method(s) not created! filef8ab3cb9ddf2.cpp:6:30: error: rstan/rstaninc.hpp: No such
file or directory
filef8ab3cb9ddf2.cpp:9:39: error: stan/model/model_header.hpp: No such file or directory
filef8ab3cb9ddf2.cpp:13: error: ‘std::vector’ has not been declared
filef8ab3cb9ddf2.cpp:14: error: ‘std::string’ has not been declared
filef8ab3cb9ddf2.cpp:15: error: ‘std::stringstream’ has not been declared
filef8ab3cb9ddf2.cpp:16: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:17: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:18: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:19: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:20: error: ‘std::istream’ has not been declared
filef8ab3cb9ddf2.cpp:21: error: ‘stan’ has not been declared
filef8ab3cb9ddf2.cpp:21: error: ‘math’ is not a namespace-name
filef8ab3cb9ddf2.cpp:21: error: expected namespace-name before ‘;’ token
fil
In addition: Warning message:
running command '/Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB filef8ab3cb9ddf2.cpp 2>
filef8ab3cb9ddf2.cpp.err.txt' had status 1

>

- Bob


On 8/8/12 6:09 PM, Jiqiang Guo wrote:
> |
>

[Jiqiang Guo]

Thanks very much.

From the error dump, I could not tell what went wrong.  I guess we need to step back: let's see if the following works. 

1. compile c++ code on the fly, please see 

http://code.google.com/p/stan/wiki/RStanGettingStarted#C++_compiler

2. install rstan from source (this will take more time though)

install.packages('rstan', repos = 'http://wiki.stan.googlecode.com/git/R', type = 'source')

The above is at 

http://code.google.com/p/stan/wiki/RStanGettingStarted

(Did not figure out how to create a link at http://code.google.com/p/stan/ to it, but I guess we need to wait until RStan works.)

Answers to previous questions: 

Then I had the same problem with:
   Rcpp
   RcppEigen
   inline
I've seen some packages auto-install other packages.  Can/should
we do this?

It seems that when the repos is CRAN, it will do that automatically.  Not sure how to make it install other packages needed in our case now.  

Q: about what if we have R version other than 2.15.x.  

I do not want to check many a binary packages to the wiki repository on google code.  But we can use the following though not sure if there is any other problem as I tested only on the latest version of R. 

install.packages('rstan', contrib = 'http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.15/')

Thanks,

Jiqiang 

[Michael Malecki]

I tried and failed at roughly the same points but have a couple possibly helpful things to add…

On Wed, Aug 8, 2012 at 7:08 PM, Jiqiang Guo <guo...@gmail.com> wrote:

It seems that when the repos is CRAN, it will do that automatically.  Not sure how to make it install other packages needed in our case now.  

You can set multiple repositories through options():

options(repos=c(getOption("repos"), Stan="'http://wiki.stan.googlecode.com/git/R"))

Then simply install.packages("rstan") should be able to find it.

Q: about what if we have R version other than 2.15.x.  

You can just make /2.14 a symlink to the 2.15 tree. You don't want to do that too far back, but I would say it can't hurt for 14 to 15. 

[Jiqiang Guo]

Thanks to Malecki and Bob for trying it. 

I am trying to recreate the problem though no clue for it.  For now, I hope the best that installing from source would work.  So on Mac and windows, please try 

install.packages('rstan', repos = 'http://wiki.stan.googlecode.com/git/R', type = 'source')

On Wed, Aug 8, 2012 at 10:47 PM, Michael Malecki <mal...@gmail.com> wrote:

I tried and failed at roughly the same points but have a couple possibly helpful things to add…

On Wed, Aug 8, 2012 at 7:08 PM, Jiqiang Guo <guo...@gmail.com> wrote:

It seems that when the repos is CRAN, it will do that automatically.  Not sure how to make it install other packages needed in our case now.  

You can set multiple repositories through options():

options(repos=c(getOption("repos"), Stan="'http://wiki.stan.googlecode.com/git/R"))

Good idea.  I have added this to 

http://code.google.com/p/stan/wiki/RStanGettingStarted 

Then simply install.packages("rstan") should be able to find it.

Q: about what if we have R version other than 2.15.x.  

You can just make /2.14 a symlink to the 2.15 tree. You don't want to do that too far back, but I would say it can't hurt for 14 to 15.

This is not working since symbol link in git is actually just a text file on the server.  

Jiqiang 

[Andrew Gelman]

J

This is what I got (see below) when I tried the line in R.  It seems to have given an ERROR.  Also, if it does work, could you post some test files so I can try to run stan() from within R?

A

-----

downloaded 2.6 Mb

ERROR: dependencies ‘Rcpp’, ‘inline’, ‘RcppEigen’ are not available for package ‘rstan’

* removing ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan’

Warning in install.packages :

  installation of package ‘rstan’ had non-zero exit status

The downloaded source packages are in

‘/private/var/folders/t1/2vc2f6t52w19hx_m166scfkc0000gn/T/RtmpooPUsH/downloaded_packages’

[Ben Goodrich]

I don't have a Mac. But I would suggest that if you want people to test the RStan package from Google Code that you just make an install script somewhere that works around various problems and platform-specific stuff. As an example, here is the one I give to people for mi:

source("http://www.columbia.edu/~bg2382/mi/install_mi.R", echo = TRUE)

You can open that in a web browser to see the code logic.

Ben

[Jiqiang Guo]

Thanks,

I have posted some code here:

http://code.google.com/p/stan/wiki/RStanGettingStarted, 

However, I need to figure out why the installation fails.  

Jiqiang 

[Jiqiang Guo]

Maybe you are right.  My thought is to make it as close as to installing it from CRAN.  And I thought people would be able to figure out if the problem is from the dependencies and solve that by themselves.  Then I suppose the call of install.packages would work.  But apparently not now. 

Thanks,

Jiqiang 

[Andrew Gelman]

J

I get the following error (see below).

A

> install.packages('Rcpp')

Installing package(s) into ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’

(as ‘lib’ is unspecified)

Warning in install.packages :

  package ‘Rcpp’ is not available (for R version 2.15.0)

[Jiqiang Guo]

Andrew:

Thanks. There is a problem with the binary version now though your error is because you do not have package Rcpp. 

Jiqiang 

[Andrew Gelman]

I'm confused.  Why could I not install Rcpp?  Does it not exist?  Do I need to upgrade to R 2.15.1? ??

[Jiqiang Guo]

I took a look from http://cran.r-project.org/web/packages/Rcpp/index.html,

the current version of Rcpp does need R 2.15.1.  That might be the reason preventing you installing Rcpp. 

I did not notice that.  This just complicates things. 

Jiqiang 

[Andrew Gelman]

That's a little scary!

[Andrew Gelman]

J
Thanks--it really works!!!!!
The only thing is, there are some problems because I was getting "404 Not Found" errors. I just ignored the errors and Stan ran just fine, but I could see that this sort of error message could bother people.
Below is the story.
A



I went to the webpage:
http://code.google.com/p/stan/wiki/RStanGettingStarted
and I followed the instructions. It all went well until this:

> install.packages('rstan')

Installing package(s) into ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’
(as ‘lib’ is unspecified)
Warning in install.packages :

cannot open: HTTP status was '404 Not Found'
trying URL 'http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.15/rstan_1.0.tgz'
Content type 'application/x-gzip; charset=binary' length unknown
opened URL
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........



I don't know if this was a problem, it said "cannot open, 404 Not Found" (this sounds bad) but then it did stuff (this sounds good). Once that was done, I tried the next command:

> install.packages('rstan', contrib = 'http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.15/')

Installing package(s) into ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’
(as ‘lib’ is unspecified)
Warning in install.packages :

cannot open: HTTP status was '404 Not Found'
trying URL 'http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.15//rstan_1.0.tgz'
Content type 'application/x-gzip; charset=binary' length unknown
opened URL
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........



Same thing. "404 Not Found". So I tried the next command:

> install.packages('rstan', repos = 'http://wiki.stan.googlecode.com/git/R', type = 'source')

Installing package(s) into ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’
(as ‘lib’ is unspecified)
Warning in install.packages :

cannot open: HTTP status was '404 Not Found'
trying URL 'http://wiki.stan.googlecode.com/git/R/src/contrib/rstan_1.0.tar.gz'
Content type 'application/x-gzip' length unknown
opened URL
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........
.



Also "404 Not Found". But it kept running and gave tons of output from the compilation. So I kept going. I loaded in stan and the schools example. It ran!

> fit <- stan(model.code = schools.code, data = dat,

+ n.iter = 2012, n.chains = 4)

COMPILING MODEL 'anon_model' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'anon_model' NOW.

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Iteration: 2000 / 2012 [ 99%] (Sampling)

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Iteration: 2000 / 2012 [ 99%] (Sampling)

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Iteration: 2000 / 2012 [ 99%] (Sampling)

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Iteration: 2000 / 2012 [ 99%] (Sampling)



It wasn't lighting-fast but it ran. Why is n.iter=2012, that's an odd number?

Anyway, I kept going:

> print (fit)
Inference for Stan model: anon_model.
4 chains: each with n.iter=2012; n.warmup=1006; n.thin=1; 2012 samples saved.

Mean SD 2.5% 25% 50% 75% 97.5% ESS Rhat
theta[1] 11.355 8.187 -1.747 6.220 9.999 15.666 30.802 254.532 1.036
theta[2] 7.790 5.926 -3.954 4.052 7.669 11.274 20.579 315.240 1.025
theta[3] 5.960 7.460 -10.332 1.932 6.591 10.417 20.719 816.329 1.009
theta[4] 7.386 6.507 -5.459 3.538 7.238 11.146 20.590 478.463 1.020
theta[5] 4.977 6.084 -8.308 1.542 5.516 8.524 16.913 382.239 1.017
theta[6] 5.761 6.402 -8.169 2.131 5.719 9.638 18.925 517.848 1.013
theta[7] 10.379 6.760 -1.375 5.537 9.827 14.581 25.025 123.649 1.047
theta[8] 8.163 7.820 -7.688 3.624 7.930 12.406 24.993 618.263 1.014
mu 7.844 4.981 -1.835 4.627 7.790 10.700 17.997 179.363 1.037
tau 6.716 5.065 1.364 3.085 5.189 8.872 19.837 61.132 1.069

Samples were drawn using NUTS2 at Thu Aug 9 18:26:11 2012.
For each parameters, ESS is a crude measure of effective samples size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
> plot (fit)

The plot is not the prettiest thing in the world but, hey, it worked!

[Andrew Gelman]

Hi all. Now that Stan works, I can see some new issues. Just for laffs I ran the 8-schaools with n.iter=100. Here's what we got:

Inference for Stan model: anon_model.

4 chains: each with n.iter=100; n.warmup=50; n.thin=1; 100 samples saved.

Mean SD 2.5% 25% 50% 75% 97.5% ESS Rhat

theta[1] 14.449 9.808 -1.448 8.403 13.700 21.419 33.724 96.486 1.018
theta[2] 8.365 7.128 -8.461 4.163 8.900 13.550 22.554 86.949 1.026
theta[3] 6.297 9.243 -9.866 0.308 7.860 12.529 22.728 135.647 1.025
theta[4] 8.039 7.697 -7.575 3.383 7.917 12.096 23.431 36.078 1.117
theta[5] 5.350 7.392 -10.747 0.073 6.124 10.524 17.986 31.981 1.140
theta[6] 5.925 7.823 -10.287 0.330 6.226 11.271 19.198 82.929 1.039
theta[7] 12.964 8.284 -3.970 8.884 13.282 17.443 31.436 57.826 1.045
theta[8] 10.138 10.123 -11.399 4.509 10.392 16.458 29.590 31.452 1.086
mu 9.155 5.600 -1.215 5.168 9.213 12.237 19.884 52.539 1.060
tau 9.462 6.327 2.569 5.198 8.281 12.017 21.786 47.363 1.046

Samples were drawn using NUTS2 at Thu Aug 9 18:31:37 2012.

For each parameters, ESS is a crude measure of effective samples size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

Then I tried it again:

4 chains: each with n.iter=100; n.warmup=50; n.thin=1; 100 samples saved.

Mean SD 2.5% 25% 50% 75% 97.5% ESS Rhat

theta[1] 10.543 6.960 -0.959 5.872 10.231 14.435 26.046 34.411 1.120
theta[2] 8.236 5.873 -2.728 4.397 8.527 11.917 18.902 101.412 1.008
theta[3] 6.910 6.554 -5.845 2.141 7.595 10.565 20.521 46.186 1.084
theta[4] 8.351 5.739 -3.476 4.498 8.529 11.890 18.385 58.148 1.079
theta[5] 5.643 5.601 -6.129 2.300 6.985 9.225 17.091 29.131 1.092
theta[6] 7.523 6.345 -6.929 3.674 8.230 11.110 20.970 88.508 1.002
theta[7] 10.823 5.186 2.402 7.682 9.968 13.406 21.518 40.083 1.072
theta[8] 8.030 6.633 -6.384 3.071 8.908 11.706 19.859 80.245 1.000
mu 8.412 4.494 0.210 5.082 8.649 10.785 17.620 29.995 1.086
tau 5.652 3.357 1.445 3.178 4.586 8.037 13.208 6.211 1.289




Here are some problems:
1. In the note at the bottom, "samples" should be "sample"
2. I'd like the default output to be 1 digit after the decimal place, not 3 (e.g., "10.5" rather than "10.543") and we should round ESS to 0 decimal places.
3. More worryingly, check out the first simulation above. R.hat < 1.1 (rounded) for all parameters and is less than 1.05 for tau. Yet the inference for tau is all wrong, the posterior mean is 9.5 rather than 7.2 (the result we get after running for 10000 iterations). Doing a quick calculation based on ESS: the mean should be estimated at an accuracy of sd/sqrt(n.eff), that would be 6.3/sqrt(47)=0.9. So the result is over 2 se's from the true value. Maybe this is ok, it's just bad luck? Perhaps the message is that, despite what David Mackay says, an ESS of 12 is not enough!
Are we using our current preferred ESS formula here?
4. Given point 3 above, I'm thinking that our output should include, right next to the mean, the se of the mean (calculated as above).

Now the good news: Stan is amazingly fast. I ran it for 100,000 iterations and, once compiling was done, it took aboout 15 seconds for all 4 chains together!

Compile time is still a bit slow (at least compared to runtime for a simple model). Here are 2 questions:
1. After compiling a model, can stan save the object code so that it doesn't need to recompile when we run it again? I ask because in practice I'm always running models over and over and over.
2. Even better, could stan automatically search the appropriate directory for the object code and then automatically just use the already-compiled model if it exists? That would be a huge help in practice.


Finally, one more thing that's interesting, regarding ESS. Here are some values of n.iter and ESS for tau (from second halves of 4 chains):

n.iter ESS
1000 55
10000 450
100000 990

Huh? Shouldn't that last number be about 5000? I tried the three again and got the folloing:

n.iter ESS
1000 86
10000 430
100000 850

What up?

[Andrew Gelman]

Hi. I just ran 8-schools for a million iterations (4 chains, it took about 2 or 3 minutes!). One thing I noticed by looking at the running counter is that Stan runs fast for the first half and the goes slow after the 50% point. That suggests to me that the warmup phase is taking us to a point where L is large in the Nuts trajectories.

4 chains: each with n.iter=1e+06; n.warmup=5e+05; n.thin=1; 1e+06 samples saved.

Mean SD 2.5% 25% 50% 75% 97.5% ESS Rhat

theta[1] 11.421 8.421 -2.231 5.931 10.260 15.640 31.759 21321.463 1.000
theta[2] 7.845 6.318 -4.785 3.892 7.796 11.762 20.719 66244.857 1.000
theta[3] 6.047 7.822 -11.646 1.869 6.554 10.804 20.565 103234.103 1.000
theta[4] 7.587 6.597 -5.808 3.537 7.582 11.636 20.981 74981.300 1.000
theta[5] 5.018 6.360 -8.983 1.225 5.478 9.277 16.401 48989.973 1.000
theta[6] 6.049 6.757 -8.682 2.099 6.399 10.380 18.673 73640.433 1.000
theta[7] 10.686 6.844 -1.223 6.089 10.066 14.647 26.066 18844.360 1.000
theta[8] 8.416 7.954 -7.193 3.768 8.143 12.738 25.671 87646.977 1.000
mu 7.883 5.207 -2.169 4.589 7.845 11.139 18.279 39607.885 1.000
tau 6.722 5.627 0.462 2.633 5.379 9.245 20.846 7571.285 1.001


And here's what we have now for ESS for tau:

n.iter ESS
1000 55
10000 450
100000 990

1000000 7500

The inferences all look basically reasonable

Does the sublinear growth of ESS as a function of n.iter suggest that the Nuts tuning is actually getting _worse_ when we have more time in the warmup phase? Giving Nuts enough rope to hang itself??

A

[Jiqiang Guo]

Andrew:

Thanks for your patience.  I guess for now install rstan from source seems to work and your patience paid off as it takes some time (but it might be not bad on your computer).  That's encouraging.  I am investigating why the binary version does not work.  If that works, it will be fast for the installation. 

And I will see if I have the 404 error. 

Why is n.iter=2012, that's an odd number? 

It is a number other than the default (2000), a way to demonstrate not to use the default.  It is just arbitrary though it is this year. 

Jiqiang  

[Jiqiang Guo]

I will fix some of the problems you raised.  Thanks for that. 

On Thu, Aug 9, 2012 at 12:51 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Here are some problems:
1.  In the note at the bottom, "samples" should be "sample"
2.  I'd like the default output to be 1 digit after the decimal place, not 3 (e.g., "10.5" rather than "10.543") and we should round ESS to 0 decimal places.
3.  More worryingly, check out the first simulation above.  R.hat < 1.1 (rounded) for all parameters and is less than 1.05 for tau.  Yet the inference for tau is all wrong, the posterior mean is 9.5 rather than 7.2 (the result we get after running for 10000 iterations).  Doing a quick calculation based on ESS:  the mean should be estimated at an accuracy of sd/sqrt(n.eff), that would be 6.3/sqrt(47)=0.9.  So the result is over 2 se's from the true value.  Maybe this is ok, it's just bad luck?  Perhaps the message is that, despite what David Mackay says, an ESS of 12 is not enough!
Are we using our current preferred ESS formula here?
4.  Given point 3 above, I'm thinking that our output should include, right next to the mean, the se of the mean (calculated as above).

Now the good news:  Stan is amazingly fast.  I ran it for 100,000 iterations and, once compiling was done, it took aboout 15 seconds for all 4 chains together!

Compile time is still a bit slow (at least compared to runtime for a simple model).  Here are 2 questions:
1.  After compiling a model, can stan save the object code so that it doesn't need to recompile when we run it again?  I ask because in practice I'm always running models over and over and over.
2.  Even better, could stan automatically search the appropriate directory for the object code and then automatically just use the already-compiled model if it exists?  That would be a huge help in practice.

For the above Q1, in one R session, it is easy to do that. I will see what interfaces (functions) I should provide. 

For Q2, that is a good idea.  It's tricky, but I believe we can do given more time to investigate.  If we can do this, Q1 is also possible for using R in another session with loading previous saved stuff. 

--

Jiqiang 

[Jiqiang Guo]

I forgot that actually for now we can compile model once and sampling differently (even with different data) as follows. 

For the eight schools in the wiki web page. 

http://code.google.com/p/stan/wiki/RStanGettingStarted

We can use function to stan.model to just compile the model.  Then we could use function sampling to draw samples. 

This is a demo:

m <- stan.model(model.code = schools.code)

fit1 <- sampling(m1, data = dat, n.iter = 100)

fit2 <- sampling(m2, data = dat, n.iter = 1000)

--

Jiqiang 

[Jiqiang Guo]

My bad, 

It should be 

m <- stan.model(model.code = schools.code)

fit1 <- sampling(m, data = dat, n.iter = 100)

fit2 <- sampling(m, data = dat, n.iter = 1000)

--

Jiqiang 

[Bob Carpenter]

I'm catching up on a bunch of messages inline.

On 8/9/12 11:30 AM, Jiqiang Guo wrote:

> I have posted some code here:
> http://code.google.com/p/stan/wiki/RStanGettingStarted,

This looks great. It's just what Andrew asked for.

We should probably do something similar for the command-line
Stan. I just hate to duplicate instructions in multiple
places because of the ongoing cost of keeping them consistent
and up to date.

I've seen other package installers automate the dependencies,
so we might be able to get rid of that step.

On 8/9/12 12:29 PM, Andrew Gelman wrote:> J

> Thanks--it really works!!!!!
> The only thing is, there are some problems because I was getting "404 Not Found" errors.

That's puzzling, because it then found the code. I wonder
if it's looking for some higher-level directory or something
we don't have?

>> install.packages('rstan')
> Installing package(s) into �/Library/Frameworks/R.framework/Versions/2.15/Resources/library�
> (as �lib� is unspecified)

> Warning in install.packages :
> cannot open: HTTP status was '404 Not Found'
> trying URL 'http://wiki.stan.googlecode.com/git/R/bin/macosx/leopard/contrib/2.15/rstan_1.0.tgz'
> Content type 'application/x-gzip; charset=binary' length unknown
> opened URL
> .......... .......... .......... .......... ..........
> .......... .......... .......... .......... ..........

On 8/9/12 12:59 PM, Andrew Gelman wrote:> Hi. I just ran 8-schools for a million iterations (4 chains, it took about 2

That's an interesting possibility. In the command-line
version, that information gets printed out to the CSV
file. Depth of tree (log2 # of steps) gets dumped as does
the value of epsilon at each iteration (to account for
when it varies). Also, the tuning parameters get dumped
as comments.

You can also track it down by setting the number of
warmup steps to 50

The way I get that info should be clear from the command
information -- I pushed most of it into an abstract base
class for HMC and then implemented versions in each of the
cases to write out alternative parameters. Rather than
writing to an output stream, they could just be returned as
numbers, too.

- Bob

[Andrew Gelman]

This is cool but I hope we can do it all using the stan() command.  That way, the user (e.g., me) can just do stan, stan, stan and it will do the right thing, rather than me having to do these specialized commands.

a

[Andrew Gelman]

We shoudl be able to get this from the R output using extractor functions.

[Bob Carpenter]

On 8/9/12 4:04 PM, Andrew Gelman wrote:
> We shoudl be able to get this [running multiple chains w/o

> from the R output using extractor functions.

Right -- we want to

1. make the HMC/NUTS debugging available via
extractor functions from the output, and

2. run multiple chains w/o recompiling from a single
call to the stan() function.

Does anyone know how to make the calls to run Stan parallel
under R?

mcmc::chains is thread safe for writing samples from different
chains -- I don't know about Jiqiang's extensions.

My post below was on where to find the info for (1) in
the C++ code.

- Bob

[Jiqiang Guo]

Please see my comments below inline. 

--

Jiqiang

On Thu, Aug 9, 2012 at 4:18 PM, Bob Carpenter <ca...@alias-i.com> wrote:

On 8/9/12 4:04 PM, Andrew Gelman wrote:
> We shoudl be able to get this [running multiple chains w/o

> from the R output using extractor functions.

Right -- we want to

1.  make the HMC/NUTS debugging available via
    extractor functions from the output, and

Currently, I do not collect those information such as step sizes back to R.  I will do that, which would be just mimicking Bob's code.  

2.  run multiple chains w/o recompiling from a single
    call to the stan() function.

Does anyone know how to make the calls to run Stan parallel
under R?

mcmc::chains is thread safe for writing samples from different
chains -- I don't know about Jiqiang's extensions.

For now, run multiple chains inside a call of stan only need to compile the c++ code once. 

But each of the following call of stan need to compile the c++ code 

stan(model.code = , data = , ...)

stan(model.code = , data = , ...)

 

I think we need some cache mechanism to make it possible to compile the code once if model.code (or file) is the same for another call.  This is doable, it is a bit like http://ccache.samba.org/ then. 

Another simpler appoach not to compile c++ code a second time for calling stan for the same model is to feed the result of the first time to later call such as the following (it is just a thought, that is, not implemented).  How does this look?

fit <- stan(model.code = , data = ,....)

fit2 <- stan(fit, .....) 

As I am not using mcmc::chains, actually I could use any parallel approaches based on multiple processes R has such as multicore (http://cran.r-project.org/web/packages/multicore/index.html), but even multicore may not work on all platforms.  I am not sure how easy it is to make it just an option. 

[Andrew Gelman]

The approach below (where I assume the object code (or a link to the location of the object code?) is saved in the fit object) would be fine, in some ways that is preferable in that it is a reminder that we're refitting the model, thus for example we fit the model once, then try again, etc.

I think it would be also good, if it were possible, for the regular stan() call to search for the object code and use it if possible. But maybe that is not so necessary if your code below works. I'd have to see how it goes when running.

I do think that some version of not-requiring-recompiling is important, expecially for little models where the compile time is longer than runtime.

A

[Bob Carpenter]

You can just add a flag to stan that
has something like "compile=true" as
a defualt.

Or, you could try to set up make-like
behavior that checks the time stamps on
the file and only recompiles if the
source file is newer than the compiled file.
I wrote this, but there's probably a better
way to do it:

time.of.file = function(file) {
return(file.info(file)$"ctime");
}

more.recently.modified = function(file1, file2) {
return(time.of.file(file1) > time.of.file(file2));
}

- Bob

[Jiqiang Guo]

Yeah, it is similar.  Just there will be more work in R as it is not that straightforward to get the compiled code and sometimes I am not dealing with files of stan code. 

--

Jiqiang 

[Jiqiang Guo]

I will do this and we will see how it goes. 

Jiqiang 

[Ben Goodrich]

On Thursday, August 9, 2012 4:18:02 PM UTC-4, Bob Carpenter wrote:

Does anyone know how to make the calls to run Stan parallel
under R?

In general, the R code would have to distinguish between !Windows and Windows, which is to say POSIX and !POSIX. If you call the compiled code via system(), there is a "wait" argument that defaults to TRUE but can be set to FALSE to async to system call. Conceivably, you could do multiple system calls that way and sleep() until they were done.

Ben

[Bob Carpenter]

I should've thought more about the overall setup.
I'm pretty sure Jiqiang integrated Stan using dynamic
linking so that it's running Stan in the same process
as R. That would require us sorting out multi-platform
threading on the C++ side and using thread-safe versions
of auto-dif.

If we are just executing Stan in a new process, then
we just want to be able to call four processes in parallel
and wait until they all finish (or two by two, etc. -- max
concurrency and total number of chains are independent).

- Bob