Stan ( http://mc-stan.org ) is one of several tools available to perform Bayesian sampling. Statistical inference based on Stan's result (the posterior samples in 1 or more chains) are typically done in R, Python or Matlab. Stan.jl uses Stan as an external program and simplifies running Stan models and subsequent inference, all from within Julia, thus providing aother 'interface' to Stan.
Stan is written in C++ but still attractive to Julia users because it is a very good tool for users interested in Bayesian models that can be efficiently sampled with Stan. It is also supported by a very capable team and community willing to help out with modeling and installation issues. Reasons enough for the existance of Stan.jl.
For Julia users who prefer a Notebook-like interface (over a Julia REPL + Editor interface), Jupyter is available and this Notebook demonstrates one way of using Stan.jl in a Notebook, including several options for the graphical display of results.
JuliaBox contributes another major benefit in that it removes the burden of installing Stan from the end users. The user still has a 'private' permanent file system to store Julia scripts that include the Stan model and available data.
Finally, pushing Stan's envelope can require knowledge of C++, C++ templates, the C++ Boost libraries, R, the Stan modeling language, and its extensions (e.g. functions). This can be a daunting task. Stan.jl can't help here but Julia's ecosystem provides options entirely written in Julia, like Mamba.jl ( http://mambajl.readthedocs.org/en/latest/index.html ) and Lora.jl which, if access to lower level MCMC functionality is required, might be a better route. Using Stan.jl will pave the way to such a migration.
Running below cell is occasionally useful to make sure all packages are up to date.
In [1]:
Pkg.update()
Stan.jl relies on Mamba.jl's very broad collection of functions for posterior inference and display after the MCMC sampling has been performed by Stan.
In [7]:
using Mamba, Stan;
In this example we'll use one of several examples that are included in Stan.jl, i.e. estimating the distribution of the success probability (theta) based on 4 sets of 10 observed outcomes.
I like to control where I store intermediate and result files. In this case I 'borrow' the Bernoulli directory inside Stan.jl.
In [8]:
old = pwd()
ProjDir = Pkg.dir("Stan", "Examples", "Bernoulli")
cd(ProjDir)
Define a simple Stan model. On JuliaBox it is best to define the model inline, but it could also be stored in an external file and uploaded to JuliaBox.
In [9]:
bernoullimodel = "
data {
int<lower=0> N;
int<lower=0,upper=1> y[N];
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1,1);
y ~ bernoulli(theta);
}
";
Define the observed data, in this case the number of observations in the 4 experiments and the outcomes of these observations. This is an array of dictionaries. By default Stan will run 4 chains and in this setup each chain is assigned a dictionary with the results of 1 experiment. If that is an acceptable way to analyze these results depends on the purpose of the 4 experiments. It's used here purely as an illustration, e.g. one could increase the probability of success in one of the experiment by updating the y vector. This will show up in the overall value of theta below and in the plots for that chain.
In [10]:
bernoullidata = [
Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
Dict("N" => 10, "y" => [0, 0, 0, 1, 1, 0, 0, 1, 0, 1])
];
Stanmodel() and stan() are the main functions provided by Stan.jl. For details on the arguments of these functions see the CmdStan documentation and the README in Stan.jl ( https://github.com/goedman/Stan.jl )
In [11]:
stanmodel = Stanmodel(update=1200, thin=2, name="bernoulli", model=bernoullimodel);
In [13]:
sim1 = stan(stanmodel, bernoullidata, ProjDir, diagnostics=false);
From above summary you can see that by default Stan.jl runs 4 chains. Stan will not perform any thinning, this is done later on in fact as we defined a thinning factor of 2 in the Stanmodel() call.
Notice that if we were to run the simulation again the Stan Model would not be compiled unless the model has been changed.
Stan records 6 variables by default augmented by model specific variables (only theta in this example). Not all of these are suitable for posterior inference. In this example we'll only perform posterior analysis on "theta", "lp" and "accept_stat".
In some cases Stan might record 10s, 100s even 1000s of model specific variables. In those cases a variable like 'posterior' can be used to constrain the number of variables used in posterior analysis.
In [14]:
posterior = ["theta", "lp__", "accept_stat__"];
Apply the subset 'posterior' before calling describe().
In [15]:
sim = sim1[:, posterior, :]
describe(sim)
Mamba.jl provides a range of function to diagnose and analyze the simulation result, e.g. gelmandiag() is shown below. See the Mamba documentation ( http://mambajl.readthedocs.org/en/latest/index.html ) or the Bernoulli example in Stan.jl for other options.
In [16]:
gelmandiag(sim, mpsrf=true, transform=true)
Out[16]:
In [17]:
p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
In [18]:
draw(p, ncol=4)
Intermediate results are all in the 'tmp' directory. Results that you might want to keep around are left in the current working directory (ProjDir as defined above).
In [19]:
;ls
What's in the 'tmp' directory?
In [20]:
;ls tmp
In [21]:
isdir("tmp") && rm("tmp", recursive=true)
In [16]: