ArviZ Quickstart

[1]:
%matplotlib inline
import arviz as az
import numpy as np

# ArviZ ships with style sheets!
az.style.use('arviz-darkgrid')

Get started with plotting

ArviZ is designed to be used with libraries like PyStan and PyMC3, but works fine with raw numpy arrays.

[2]:
az.plot_posterior(np.random.randn(100_000));
../_images/notebooks_Introduction_3_0.png

Plotting a dictionary of arrays, ArviZ will interpret each key as the name of a different random variable. Each row of an array is treated as an independent series of draws from the variable, called a chain. Below, we have 10 chains of 50 draws each for four different distributions.

[3]:
size = (10, 50)
az.plot_forest({
    'normal': np.random.randn(*size),
    'gumbel': np.random.gumbel(size=size),
    'student t': np.random.standard_t(df=6, size=size),
    'exponential': np.random.exponential(size=size)
});
../_images/notebooks_Introduction_5_0.png

Plotting with PyMC3 objects

ArviZ is designed to work well with high dimensional, labelled data. Consider the eight schools model, which roughly tries to measure the effectiveness of SAT classes at eight different schools. To show off ArviZ’s labelling, I give the schools the names of a different eight schools.

This model is small enough to write down, is hierarchical, uses labelling, and a centered parameterization causes divergences (which are interesting for illustration):

[4]:
import pymc3 as pm

J = 8
y = np.array([28.,  8., -3.,  7., -1.,  1., 18., 12.])
sigma = np.array([15., 10., 16., 11.,  9., 11., 10., 18.])
schools = np.array(['Choate', 'Deerfield', 'Phillips Andover', 'Phillips Exeter',
                    'Hotchkiss', 'Lawrenceville', "St. Paul's", 'Mt. Hermon'])


with pm.Model() as centered_eight:
    mu = pm.Normal('mu', mu=0, sd=5)
    tau = pm.HalfCauchy('tau', beta=5)
    theta = pm.Normal('theta', mu=mu, sd=tau, shape=J)
    obs = pm.Normal('obs', mu=theta, sd=sigma, observed=y)

    # This pattern is useful in PyMC3
    prior = pm.sample_prior_predictive()
    centered_eight_trace = pm.sample()
    posterior_predictive = pm.sample_posterior_predictive(centered_eight_trace)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [theta, tau, mu]
Sampling 4 chains: 100%|██████████| 4000/4000 [00:02<00:00, 1659.82draws/s]
There were 9 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.679274648861103, but should be close to 0.8. Try to increase the number of tuning steps.
There were 15 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.6222461372253658, but should be close to 0.8. Try to increase the number of tuning steps.
There were 32 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.5960725891337998, but should be close to 0.8. Try to increase the number of tuning steps.
There were 31 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.7091049759437956, but should be close to 0.8. Try to increase the number of tuning steps.
The estimated number of effective samples is smaller than 200 for some parameters.
100%|██████████| 2000/2000 [00:00<00:00, 4128.69it/s]

Most ArviZ functions work fine with trace objects from PyMC3:

[5]:
az.plot_autocorr(centered_eight_trace, var_names=['mu', 'tau']);
../_images/notebooks_Introduction_9_0.png

Convert to InferenceData

For much more powerful querying, analysis and plotting, we can use built-in ArviZ utilities to convert PyMC3 objects to xarray datasets. Note we are also giving some information about labelling.

ArviZ is built to work with InferenceData, and the more groups it has access to, the more powerful analyses it can perform. Here is a plot of the trace, which is common in PyMC3 workflows. Note the intelligent labels.

[6]:
data = az.from_pymc3(trace=centered_eight_trace,
                     prior=prior,
                     posterior_predictive=posterior_predictive,
                     coords={'school': schools},
                     dims={'theta': ['school'], 'obs': ['school']})
data
[6]:
Inference data with groups:
    > posterior
    > sample_stats
    > posterior_predictive
    > prior
    > observed_data
[7]:
az.plot_trace(data);




















../_images/notebooks_Introduction_12_1.png

Plotting with PyStan objects

ArviZ is built with first class support for PyStan objects, and can plot raw fit objects in a reasonable manner. Here is the same centered eight schools model:

[8]:
import pystan


schools_code = """
data {
  int<lower=0> J;
  real y[J];
  real<lower=0> sigma[J];
}

parameters {
  real mu;
  real<lower=0> tau;
  real theta[J];
}

model {
  mu ~ normal(0, 5);
  tau ~ cauchy(0, 5);
  theta ~ normal(mu, tau);
  y ~ normal(theta, sigma);
}
generated quantities {
    vector[J] log_lik;
    vector[J] y_hat;
    for (j in 1:J) {
        log_lik[j] = normal_lpdf(y[j] | theta[j], sigma[j]);
        y_hat[j] = normal_rng(theta[j], sigma[j]);
    }
}
"""

schools_dat = {'J': 8,
               'y': [28,  8, -3,  7, -1,  1, 18, 12],
               'sigma': [15, 10, 16, 11,  9, 11, 10, 18]}

sm = pystan.StanModel(model_code=schools_code, verbose=False)
fit = sm.sampling(data=schools_dat, iter=1000, chains=4)
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_b0eada194ad938e30ed47926216f98ca NOW.
/home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/Cython/Compiler/Main.py:367: FutureWarning: Cython directive 'language_level' not set, using 2 for now (Py2). This will change in a later release! File: /tmp/tmpwz_kifss/stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889.pyx
  tree = Parsing.p_module(s, pxd, full_module_name)
In file included from /home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1821:0,
                 from /home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,
                 from /home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpwz_kifss/stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889.cpp:688:
/home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
In file included from /tmp/tmpwz_kifss/stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889.cpp:692:0:
/home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/pystan/stan_fit.hpp: In function ‘int pystan::{anonymous}::command(pystan::StanArgs&, Model&, pystan::StanHolder&, const std::vector<long unsigned int>&, const std::vector<std::__cxx11::basic_string<char> >&, RNG_t&)’:
/home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/pystan/stan_fit.hpp:830:12: warning: ‘template<class> class std::auto_ptr’ is deprecated [-Wdeprecated-declarations]
       std::auto_ptr<stan::io::var_context> init_context_ptr;
            ^
In file included from /usr/include/c++/5/bits/locale_conv.h:41:0,
                 from /usr/include/c++/5/locale:43,
                 from /usr/include/c++/5/iomanip:43,
                 from /home/tbayes/miniconda3/envs/arviz3.6/lib/python3.6/site-packages/pystan/stan_fit.hpp:5,
                 from /tmp/tmpwz_kifss/stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889.cpp:692:
/usr/include/c++/5/bits/unique_ptr.h:49:28: note: declared here
   template<typename> class auto_ptr;
                            ^
/tmp/tmpwz_kifss/stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889.cpp: In function ‘PyObject* __pyx_pf_71stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889_2_call_sampler(PyObject*, PyObject*, PyObject*, PyObject*)’:
/tmp/tmpwz_kifss/stanfit4anon_model_b0eada194ad938e30ed47926216f98ca_1526190359138879889.cpp:9242:30: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
     __pyx_t_12 = ((__pyx_t_9 != __pyx_v_fitptr->param_names_oi().size()) != 0);
                              ^
WARNING:pystan:Rhat for parameter lp__ is 1.1101320455621984!
WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed
WARNING:pystan:103 of 2000 iterations ended with a divergence (5.15%).
WARNING:pystan:Try running with adapt_delta larger than 0.8 to remove the divergences.
WARNING:pystan:Chain 2: E-BFMI = 0.10311272757204448
WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model
[9]:
az.plot_density(fit, var_names=['mu', 'tau']);
../_images/notebooks_Introduction_15_0.png

Again, converting to InferenceData (a netcdf datastore that loads data into xarray datasets), we can get much richer labelling and mixing of data. Here is a plot showing where the Hamiltonian sampler had divergences

[11]:
data = az.from_pystan(posterior=fit,
                      posterior_predictive='y_hat',
                      observed_data=['y'],
                      log_likelihood='log_lik',
                      coords={'school': schools},
                      dims={'theta': ['school'], 'y': ['school'], 'log_lik': ['school'], 'y_hat': ['school'], 'theta_tilde': ['school']})
data
[11]:
Inference data with groups:
    > posterior
    > sample_stats
    > posterior_predictive
    > observed_data
[12]:
az.plot_pair(data, coords={'school': ['Choate', 'Deerfield', 'Phillips Andover']}, divergences=True);
../_images/notebooks_Introduction_18_0.png