Extracting samples

MCMC samples obtained by calling the stan function or the sampling method of a StanModel object are stored in a StanFit object. Using the schools example from the Readme:

% This defaults to 1000 samples, following 1000 warmup samples
% Note that if 'inc_warmup' is not set true, the default is not to keep warmup samples
fit = stan('model_code',schools_code,'data',schools_dat,'inc_warmup',true)

fit = 

  StanFit with properties:

          model: [1×1 StanModel]
      processes: [1×4 processManager]
    output_file: {1×4 cell}
        verbose: 0
     exit_value: [0 0 0 0]
         loaded: [1 1 1 1]
            sim: [1×1 mcmc]
        version: '0.8.0'

The samples are stored in the sim property, where warmup samples are stored separately if inc_warmup was true when the model was sampled:

fit.sim = 

  mcmc with properties:

            names: {11×1 cell}
         n_warmup: [1×4 struct]
        n_samples: [1×4 struct]
    permute_index: [1×4000 double]
        user_data: {1×4 cell}
           warmup: [1×4 struct]
          samples: [1×4 struct]
          version: '0.3.1'

The samples are stored as a struct array (one element for each chain), with fields corresponding to parameters defined in the Stan model, as well as details about each chain:

fit.sim.samples = 

  1×4 struct array with fields:

    lp__
    accept_stat__
    stepsize__
    treedepth__
    n_leapfrog__
    divergent__
    energy__
    mu
    tau
    eta
    theta

Samples can be conveniently obtained using the extract method:

% Calling without arguments returns samples from all chains concatenated and randomly permuted
samples = fit.extract()

samples = 

  struct with fields:

             lp__: [4000×1 double]
    accept_stat__: [4000×1 double]
       stepsize__: [4000×1 double]
      treedepth__: [4000×1 double]
     n_leapfrog__: [4000×1 double]
      divergent__: [4000×1 double]
         energy__: [4000×1 double]
               mu: [4000×1 double]
              tau: [4000×1 double]
              eta: [4000×8 double]
            theta: [4000×8 double]

% Return separate chains without permutation
samples = fit.extract('permuted',false)
samples = 

  1×4 struct array with fields:

    lp__
    accept_stat__
    stepsize__
    treedepth__
    n_leapfrog__
    divergent__
    energy__
    mu
    tau
    eta
    theta

% Prepend the warmup samples
samples = fit.extract('permuted',false,'inc_warmup',true);

samples(1)

ans = 

  struct with fields:

             lp__: [2000×1 double]
    accept_stat__: [2000×1 double]
       stepsize__: [2000×1 double]
      treedepth__: [2000×1 double]
     n_leapfrog__: [2000×1 double]
      divergent__: [2000×1 double]
         energy__: [2000×1 double]
               mu: [2000×1 double]
              tau: [2000×1 double]
              eta: [2000×8 double]
            theta: [2000×8 double]

To retrieve samples for a subset of variables, use the par input with a cell array of names:

samples = fit.extract('par',{'eta' 'theta'})

samples = 

  struct with fields:

      eta: [4000×8 double]
    theta: [4000×8 double]