<jupyter_start><jupyter_text>Inference Data Cookbook`InferenceData` is the central data format for ArviZ. `InferenceData` itself is just a container that maintains references to one or more `xarray.Dataset`. Below are various ways to generate an `InferenceData` object. See [here](XarrayforArviZ.ipynb) for more on xarray.<jupyter_code>import arviz as az
import numpy as np<jupyter_output><empty_output><jupyter_text>From 1d numpy array<jupyter_code>size = 100
dataset = az.convert_to_inference_data(np.random.randn(size))
print(dataset)
dataset.posterior<jupyter_output>Inference data with groups:
	> posterior<jupyter_text>From nd numpy array<jupyter_code>shape = (1, 2, 3, 4, 5)
dataset = az.convert_to_inference_data(np.random.randn(*shape))
print(dataset)
dataset.posterior<jupyter_output>Inference data with groups:
	> posterior<jupyter_text>From a dictionary<jupyter_code>datadict = {
    'a': np.random.randn(100),
    'b': np.random.randn(1, 100, 10),
    'c': np.random.randn(1, 100, 3, 4),
}
dataset = az.convert_to_inference_data(datadict)
print(dataset)
dataset.posterior<jupyter_output>Inference data with groups:
	> posterior<jupyter_text>From dictionary with coords and dims<jupyter_code>datadict = {
    'a': np.random.randn(100),
    'b': np.random.randn(1, 100, 10),
    'c': np.random.randn(1, 100, 3, 4),
}
coords = {'c1' : np.arange(3), 'c2' : np.arange(4), 'b1' : np.arange(10)}
dims = {'b' : ['b1'], 'c' : ['c1', 'c2']}

dataset = az.convert_to_inference_data(datadict, coords=coords, dims=dims)
print(dataset)
dataset.posterior<jupyter_output>Inference data with groups:
	> posterior<jupyter_text>From pymc3<jupyter_code>import pymc3 as pm
draws = 500
chains = 2

eight_school_data = {'J': 8,
                     'y': np.array([28., 8., -3., 7., -1., 1., 18., 12.]),
                     'sigma': np.array([15., 10., 16., 11., 9., 11., 10., 18.])
                    }

with pm.Model() as model:
    mu = pm.Normal('mu', mu=0, sd=5)
    tau = pm.HalfCauchy('tau', beta=5)
    theta_tilde = pm.Normal('theta_tilde', mu=0, sd=1, shape=eight_school_data['J'])
    theta = pm.Deterministic('theta', mu + tau * theta_tilde)
    pm.Normal('obs', mu=theta, sd=eight_school_data['sigma'], observed=eight_school_data['y'])
    
    trace = pm.sample(draws, chains=chains)
    prior = pm.sample_prior_predictive()
    posterior_predictive = pm.sample_posterior_predictive(trace, 500, model)
    
    data = az.from_pymc3(
            trace=trace,
            prior=prior,
            posterior_predictive=posterior_predictive,
            coords={'school': np.arange(eight_school_data['J'])},
            dims={'theta': ['school'], 'theta_tilde': ['school']},
        )
data<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 4 jobs)
NUTS: [theta_tilde, tau, mu]
Sampling 2 chains: 100%|██████████| 2000/2000 [00:01<00:00, 1407.23draws/s]
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
100%|██████████| 500/500 [00:00<00:00, 1951.70it/s]<jupyter_text>From pystan<jupyter_code>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_tilde[J];
        }

        transformed parameters {
            real theta[J];
            for (j in 1:J)
                theta[j] = mu + tau * theta_tilde[j];
        }

        model {
            mu ~ normal(0, 5);
            tau ~ cauchy(0, 5);
            theta_tilde ~ normal(0, 1);
            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]);
            }
        }
    '''
stan_model = pystan.StanModel(model_code=schools_code)
fit = stan_model.sampling(data=eight_school_data,
                          iter=draws,
                          warmup=0,
                          chains=chains)

data = az.from_pystan(posterior=fit,
                      posterior_predictive='y_hat',
                      observed_data=['y'],
                      log_likelihood='log_lik',
                      coords={'school': np.arange(eight_school_data['J'])},
                      dims={'theta': ['school'],
                             'y': ['school'],
                             'log_lik': ['school'],
                             'y_hat': ['school'],
                             'theta_tilde': ['school']
                            }
                     )

data<jupyter_output>INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_e7aeb8b836685a923269e6171e7377cd 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/tmpehtg_iwy/stanfit4anon_model_e7aeb8b836685a923269e6171e7377cd_1721201981547910013.pyx
  tree = Parsing.p_module(s, pxd, full_module_name)<jupyter_text>From pyro<jupyter_code>import torch 

import pyro
import pyro.distributions as dist
import pyro.poutine as poutine
from pyro.infer.mcmc import MCMC, NUTS

pyro.enable_validation(True)
pyro.set_rng_seed(0)

draws = 1000
warmup_steps = 0
eight_school_data = {'J' : 8,
                     'y' : torch.tensor([28,  8, -3,  7, -1,  1, 18, 12]).type(torch.Tensor),
                     'sigma' : torch.tensor([15, 10, 16, 11, 9, 11, 10, 18]).type(torch.Tensor)
                    }


def model(sigma):
    eta = pyro.sample('eta', dist.Normal(torch.zeros(eight_school_data['J']), torch.ones(eight_school_data['J'])))
    mu = pyro.sample('mu', dist.Normal(torch.zeros(1), 10 * torch.ones(1)))
    tau = pyro.sample('tau', dist.HalfCauchy(scale=25 * torch.ones(1)))

    theta = mu + tau * eta

    return pyro.sample("obs", dist.Normal(theta, sigma))


def conditioned_model(model, sigma, y):
    return poutine.condition(model, data={"obs": y})(sigma)



nuts_kernel = NUTS(conditioned_model, adapt_step_size=True)
posterior = MCMC(nuts_kernel,
                 num_samples=draws,
                 warmup_steps=warmup_steps).run(model, eight_school_data['sigma'], eight_school_data['y'])

pyro_data = az.from_pyro(posterior)
pyro_data<jupyter_output>INFO:pyro.infer.mcmc.mcmc:Starting MCMC using kernel - NUTS ...
INFO:pyro.infer.mcmc.mcmc:Iteration: 50 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.125390 	 Acceptance rate: 0.860000
INFO:pyro.infer.mcmc.mcmc:Iteration: 100 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.789698 	 Acceptance rate: 0.910000
INFO:pyro.infer.mcmc.mcmc:Iteration: 150 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.282786 	 Acceptance rate: 0.920000
INFO:pyro.infer.mcmc.mcmc:Iteration: 200 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.105557 	 Acceptance rate: 0.930000
INFO:pyro.infer.mcmc.mcmc:Iteration: 250 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.527265 	 Acceptance rate: 0.944000
INFO:pyro.infer.mcmc.mcmc:Iteration: 300 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.200404 	 Acceptance rate: 0.940000
INFO:pyro.infer.mcmc.mcmc:Iteration: 350 [SAMPLE]
INFO:pyro.infer.mcmc.mcmc:Step size: 0.286618 	 Acceptance rate: 0.945714
INFO:pyro.infer.mcmc.mcmc:Iteration: 400 [SAMPLE]
INFO:pyro.infer.mcm[...]<jupyter_text>From emcee<jupyter_code>import emcee

eight_school_data = {'J': 8,
                     'y': np.array([28., 8., -3., 7., -1., 1., 18., 12.]),
                     'sigma': np.array([15., 10., 16., 11., 9., 11., 10., 18.])
                    }

def log_prior_8school(theta,J):
    mu = theta[0]
    tau = theta[1]
    eta = theta[2:]
    # Half-cauchy prior
    if tau<0:
        return -np.inf
    hwhm = 25
    prior_tau = -np.log(tau**2+hwhm**2)
    prior_mu = -(mu/10)**2  # normal prior, loc=0, scale=10
    prior_eta = -np.sum(eta**2)  # normal prior, loc=0, scale=1
    return prior_mu + prior_tau + prior_eta
    
def log_likelihood_8school(theta,y,sigma):
    mu = theta[0]
    tau = theta[1]
    eta = theta[2:]
    return -np.sum(((mu + tau * eta - y) / sigma)**2)
    
def lnprob_8school(theta,J,y,sigma):
    prior = log_prior_8school(theta,J)
    if prior <= -np.inf:
        return -np.inf
    like = log_likelihood_8school(theta,y,sigma)
    return like+prior

nwalkers = 40
ndim = eight_school_data['J']+2
draws = 1500
pos = np.random.normal(size=(nwalkers,ndim))
pos[:,1] = np.absolute(pos[:,1])
sampler = emcee.EnsembleSampler(nwalkers, 
                                ndim, 
                                lnprob_8school, 
                                args=(eight_school_data['J'], 
                                      eight_school_data['y'], 
                                      eight_school_data['sigma']
                                     )
                               )
sampler.run_mcmc(pos, draws)

# define variable names, it cannot be inferred from emcee
var_names = ['mu','tau']+['eta{}'.format(i) for i in range(eight_school_data['J'])]
emcee_data = az.from_emcee(sampler, var_names = var_names)
emcee_data<jupyter_output><empty_output>