How are the values included in gs.Summary exactly computed?

[jeli2103]

Dear Liesels,
The code for the tutorial "Parameter transformation" has some functionalities from older Liesel versions. Here is a small update:

import jax
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt

# We use distributions and bijectors from tensorflow probability
import tensorflow_probability.substrates.jax.distributions as tfd
import tensorflow_probability.substrates.jax.bijectors as tfb

import liesel.goose as gs
import liesel.model as lsl

rng = np.random.default_rng(42)

# data-generating process
n = 500
true_beta = np.array([1.0, 2.0])
true_sigma = 1.0
x0 = rng.uniform(size=n)
X_mat = np.c_[np.ones(n), x0]
y_vec = X_mat @ true_beta + rng.normal(scale=true_sigma, size=n)

# Model
# Part 1: Model for the mean
beta_prior = lsl.Dist(tfd.Normal, loc=0.0, scale=100.0)
beta = lsl.Param(value=np.array([0.0, 0.0]), distribution=beta_prior,name="beta")

X = lsl.Obs(X_mat, name="X")
mu = lsl.Var(lsl.Calc(jnp.dot, X, beta), name="mu")

# Part 2: Model for the standard deviation
a = lsl.Var(0.01, name="a")
b = lsl.Var(0.01, name="b")
sigma_sq_prior = lsl.Dist(tfd.InverseGamma, concentration=a, scale=b)
sigma_sq = lsl.Param(value=10.0, distribution=sigma_sq_prior, name="sigma_sq")

sigma = lsl.Var(lsl.Calc(jnp.sqrt, sigma_sq), name="sigma")

# Observation model
y_dist = lsl.Dist(tfd.Normal, loc=mu, scale=sigma)
y = lsl.Var(y_vec, distribution=y_dist, name="y")

gb = lsl.GraphBuilder().add(y)

gb.transform(sigma_sq, tfb.Exp)

model = gb.build_model()

builder = gs.EngineBuilder(seed=1339, num_chains=4)

builder.set_model(lsl.GooseModel(model))
builder.set_initial_values(model.state)

builder.add_kernel(gs.NUTSKernel(["beta", "sigma_sq_transformed"]))

builder.set_duration(warmup_duration=1000, posterior_duration=1000)

# by default, goose only stores the parameters specified in the kernels.
# let's also store the standard deviation on the original scale.
builder.positions_included = ["sigma_sq"]

engine = builder.build()

engine.sample_all_epochs()

results = engine.get_results()

# what is the output here? 
summary = gs.Summary(results).to_dataframe().reset_index()

However, I cannot make use of all summary statistics, as I am not sure what they mean. Could you maybe give some explanations in the documentary or in this tutorial? I would appreciate it.

[jobrachem]

Hi @jeli2103, thanks for your question! You are right that more documentation of the summary object would be nice and important. I opened a corresponding issue: #158

In the meantime, there's an interim solution. Under the hood, we use arviz to calculate the summary statistics. It will probably already help you if you consult their API documentation, for example:

• arviz.hdi
• arviz.ess
• arviz.rhat