Getting Started

##Prerequisites MatlabStan has the following dependencies:

1. C++ compiler
2. CmdStan: 2.0.1 or greater
3. MatlabProcessManager: 0.3.0 or greater

The first two are covered in detail for different platforms in the Stan Manual (Section 2 and Appendix B).

MatlabProcessManager is a single Matlab file; installation simply add it to your Matlab path.

##Installation MatlabStan is written in pure Matlab. To install:

1. Obtain a copy here or clone the repo.
2. Add the resulting folder to your Matlab path.
3. Edit the file stan_home.m in the +mstan directory to point to the parent folder of your CmdStan installation.

###Optional Installing Steve Eddins's linewrap function is useful for dealing with unwrapped messages. His xUnit test framework is required if you want to run the unit tests.

##Using MatlabStan

Example: eight schools

This is a classic example from Section 5.5 of Gelman et al (2003). The following can be compared to the Rstan and Pystan versions.

schools_code = {
   'data {'
   '    int<lower=0> J; // number of schools '
   '    real y[J]; // estimated treatment effects'
   '    real<lower=0> sigma[J]; // s.e. of effect estimates '
   '}'
   'parameters {'
   '    real mu; '
   '    real<lower=0> tau;'
   '    real eta[J];'
   '}'
   'transformed parameters {'
   '    real theta[J];'
   '    for (j in 1:J)'
   '    theta[j] <- mu + tau * eta[j];'
   '}'
   'model {'
   '    eta ~ normal(0, 1);'
   '    y ~ normal(theta, sigma);'
   '}'
};
  
schools_dat = struct('J',8,...
                     'y',[28 8 -3 7 -1 1 18 12],...
                     'sigma',[15 10 16 11 9 11 10 18]);

fit = stan('model_code',schools_code,'data',schools_dat);

print(fit);

eta = fit.extract('permuted',true).eta;
mean(eta)

Stan models can also be defined using a file. For example, download the file eight_schools.stan into your working directory and use the following call:

fit1 = stan('file','eight_schools.stan','data',schools_dat,'iter',1000,'chains',4);

Once a model is fitted, we can reuse the result as an input to stan with other data or settings. This saves us from having to compile the C++ code again (see also ). For example, if we want to sample more iterations:

fit2 = stan('fit',fit1,'data',schools_dat,'iter',10000,'chains',4);

The stan function returns a StanFit object, which contains samples from the posterior distribution. StanFit objects possess a number of methods, including print, traceplot and extract. For example, a summary of the posterior samples as well as the log-posterior (which has the name lp__) is obtained using

print(fit2);

which should look something like this:

Inference for Stan model: eight_schools_model
4 chains: each with iter=(5000,5000,5000,5000); warmup=(0,0,0,0); thin=(1,1,1,1); 20000 iterations saved.
Warmup took (0.16, 0.18, 0.22, 0.20) seconds, 0.77 seconds total
Sampling took (0.29, 0.30, 0.30, 0.25) seconds, 1.1 seconds total
                    Mean     MCSE  StdDev        5%       50%    95%  N_Eff  N_Eff/s    R_hat
lp__            -4.8e+00  4.0e-02     2.6  -9.4e+00  -4.6e+00  -0.92   4364     3828  1.0e+00
accept_stat__    7.2e-01  7.0e-02    0.30   3.4e-02   8.5e-01    1.0     19       16  1.1e+00
stepsize__       4.7e-01  5.3e-02   0.075   3.7e-01   4.9e-01   0.58    2.0      1.8  1.5e+13
treedepth__      1.7e+00  1.8e-01    0.61   0.0e+00   2.0e+00    2.0     11     10.0  1.1e+00
n_divergent__    8.7e-03  1.7e-03   0.093   0.0e+00   0.0e+00   0.00   3095     2715  1.0e+00
mu               8.0e+00  9.4e-02     5.1  -9.5e-02   7.9e+00     16   2970     2605  1.0e+00
tau              6.7e+00  1.0e-01     5.5   5.3e-01   5.5e+00     17   2873     2520  1.0e+00
eta[1]           4.0e-01  9.1e-03    0.93  -1.2e+00   4.1e-01    1.9  10496     9206  1.0e+00
eta[2]          -2.7e-03  8.9e-03    0.87  -1.4e+00  -7.3e-03    1.4   9545     8372  1.0e+00
eta[3]          -2.0e-01  9.0e-03    0.93  -1.7e+00  -2.2e-01    1.4  10578     9279  1.0e+00
eta[4]          -3.2e-02  8.4e-03    0.89  -1.5e+00  -3.3e-02    1.4  11090     9727  1.0e+00
eta[5]          -3.5e-01  8.8e-03    0.87  -1.8e+00  -3.7e-01    1.1   9694     8503  1.0e+00
eta[6]          -2.2e-01  9.1e-03    0.90  -1.6e+00  -2.4e-01    1.3   9764     8564  1.0e+00
eta[7]           3.5e-01  8.7e-03    0.87  -1.1e+00   3.7e-01    1.7  10018     8787  1.0e+00
eta[8]           5.5e-02  8.9e-03    0.93  -1.5e+00   5.6e-02    1.6  10972     9624  1.0e+00
theta[1]         1.1e+01  1.0e-01     8.3   2.2e-01   1.0e+01     27   6320     5544  1.0e+00
theta[2]         7.9e+00  6.3e-02     6.3  -2.2e+00   7.8e+00     18  10076     8838  1.0e+00
theta[3]         6.1e+00  8.9e-02     7.8  -7.3e+00   6.6e+00     18   7582     6651  1.0e+00
theta[4]         7.7e+00  6.6e-02     6.6  -3.2e+00   7.7e+00     19   9987     8760  1.0e+00
theta[5]         5.1e+00  6.4e-02     6.4  -6.2e+00   5.6e+00     15  10122     8878  1.0e+00
theta[6]         6.1e+00  6.9e-02     6.8  -5.7e+00   6.4e+00     17   9765     8566  1.0e+00
theta[7]         1.1e+01  7.7e-02     6.8   8.4e-01   1.0e+01     23   7782     6826  1.0e+00
theta[8]         8.6e+00  1.0e-01     8.2  -3.9e+00   8.3e+00     22   6297     5523  1.0e+00
Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at 
convergence, R_hat=1).

The extract method returns a struct or struct arrayfor parameters of interest

% return a struct with all parameters
la = fit.extract('permuted',true);  
mu = la.mu

% return an array with selected parameter
mu2 = fit.extract('pars','mu').mu;

% returns individual chains (each array element is a chain)
a = fit.extract('permuted',false);

Plotting traces is pretty basic at the moment

fit.traceplot;