GAMs and GPs: Flexibility and Calibration

import pandas as pd
import matplotlib.pyplot as plt
import bambi as bmb
import seaborn as sns
from pygam.datasets import mcycle
from pygam import LinearGAM, s, f, GAM, l, utils
import numpy as np
import arviz as az
import pymc as pm
import scipy as sp
from patsy import bs as bs_patsy, dmatrix
import pytensor.tensor as pt


random_seed = 100

import warnings

warnings.simplefilter('ignore')

GAMs and GPs: Flexibility and Calibration

Flexible Spline models risk overfit and need to be carefully calibrated against real data. Hierarchical models impose structure and aid calibration. We combine them.

Summary

Statistical models both encode structure in their definition and learn structure as they are calibrated against real data. There is a tension between how much structure we encode and how much we allow the model to discover. Too much discovery turns the model into a black-box and we lose interpretability, too much encoded structure and we bias the discovery process and risk losing sight of the data generating mechanism.

Generalised Additive spline models and Gaussian process models are a tractable set of mathematical abstractions fit for arbitrary function approximation. We demonstrate their flexibility and then continue to explore the interplay between such free wheeling function-approximation and theoretically imposed constraints in the case of hierarchical spline modelling of insurance losses.

Theoretical Posits and Scientific Commitments

The scientific enterprise must abstract away details of reality to provide insight into the structure of reality. Statistical models of real mechanisms are cobbled together from a variety of abstractions; probability distributions and functional relations cojoined, are then calibrated against observed data. We are building on disparate foundations. Good abstractions aim at modelling the structures or real patterns that matter. It doesn’t matter if those abstractions are themselves reflective of reality.

[T]here is no way to establish [a] fully secured, … starting point of the sciences. There is no tabula rasa. We are like sailors who have to rebuild their ship on the open sea, without ever being able to dismantle it in dry-dock and reconstruct it from its best components. - Otto Neurath

This is an apt metaphor for spline models in particular - like these ships built at sea, spline models are constructed from ad-hoc materials fished from the vast ocean of possible models. They are linear approxiations invoked to ride the next wave form.

In this blog post we’re going to dive into modelling of non-linear functions and explore some of the tooling available in the python eco-system. We’ll start by looking into Generalised Additive Models with splines in pyGAM before proceeding to look at how Bayesian versions of spline modelling compares to Gaussian processes in Bambi and PyMC. Finally we will show how hierarchical bayesian models avoid some of the issues of overfit in simpler spline models.

Our interest in these models stems from their flexibility to approximate functions of arbitrary complexity. We’ll see how the methods work in the case of relatively straightforward toy example and then we’ll apply each of the methods to deriving insights into the functional form of insurance loss curves. In this application we adapt a data set discussed in Mick Cooney’s Stan case study to demonstrate the power of hierarchical spline models. Throughout we’ll draw on the discussion of these methods in Osvaldo Martin’s “Bayesian Analysis with Python” for practical details implementing these models.

Generalised Additive models

Nonlinear functions can be approximated with linearly additive combinations of component features. Before delving into the practicalities we’ll quickly note some of the theoretical background. The canonical reference for these additive models is Simon Wood’s “Generalised Additive Models: An Introduction with R” which outlines in some detail the theoretical background of splines and univariate smoothers. The book stresses the trade-offs between the flexibility of splines and the need for cross-validation and penalised estimation methods for spline based modelling. Spline models ape the complexities of esoteric functions by regression-like formulas combining subordinate functions of observed variables to to predict another:

\[y \sim f(x_1) + f(x_2) ...\]

The structure of the individual \(f\) functions is approximated. So we are approximating \(y\) by adding a series of input approximations \(\sum f_{i}\). In R these penalised models fits can be achieved in mgcv which incorporates a Wilkinson like formula syntax for model specification: y ~ s(x) + s(x1). The closest implementation in python is available in PyGam and we will adopt this package to illustrate the spline based smoothing.

PyGAM and Penalised Fits

Let’s first look at an example data set on which to demonstrate univariate smoothing patterns using penalised splines. We’ll initially look simply at the function calls before going “under the hood”.

Spline models can be optimised by fitting differing strength penalities over a varying number of splines. Splines are used to construct non-linear functions of the input variable which are combined additively in our model equation. The penalities used in the optimisation routines constrain how “wiggly” these constructions will be. The number of splines determine how precise our will be.

X, y = mcycle(return_X_y=True)

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(X, y)
ax.set_ylabel("Acceleration")
ax.set_xlabel("Time Step")
ax.set_title("Crash Test Dummy Acceleration \n Simulated Motorcycle Crash", fontsize=20);

Now we fit a number of different models to account for the herky-jerky nature of the data generating processs. We vary the parameterisations to see how the numbers of splines and strength of the penalty help account for the variation in \(y\) over the support of \(X\).

gam1 = LinearGAM(s(0, n_splines=5)).fit(X, y)
gam2 = LinearGAM(s(0, n_splines=7)).fit(X, y)
gam3 = LinearGAM(s(0, n_splines=10)).fit(X, y)
gam4 = LinearGAM(s(0, n_splines=15)).fit(X, y)
gam5 = LinearGAM(s(0, lam=.1)).fit(X, y)
gam6 = LinearGAM(s(0, lam=.5)).fit(X, y)
gam7 = LinearGAM(s(0, lam=5)).fit(X, y)
gam8 = LinearGAM(s(0, lam=15)).fit(X, y)


def plot_fit(gam, X, y, ax, t, c1='b', c2='r'):
    XX = gam.generate_X_grid(term=0, n=500)

    ax.plot(XX, gam.predict(XX), color=c2, linestyle='--')
    ax.plot(XX, gam.prediction_intervals(XX, width=.95), color=c1, ls='--')

    ax.scatter(X, y, facecolor='gray', edgecolors='none')
    ax.set_title(f"""95% prediction interval with {t} \n LL: {gam.statistics_['loglikelihood']}""");

fig, axs = plt.subplots(4,2, figsize=(10, 20))
axs = axs.flatten()
titles = ['5_splines', '7_splines', '10_splines', '15_splines',
'lam=.1', 'lam=.5', 'lam=5', 'lam=15']
gs = [gam1, gam2, gam3, gam4, gam5, gam6, gam7, gam8]
for ax, g, t in zip(axs, gs, titles):
    plot_fit(g, X, y, ax, t)

Here we’ve seen the PyGAM package applied to fitting our model to the data. In the formula specification we see y ~ s(i) where i denotes the index of the column variable in the X data.

Over the range of of the x-axis we can see how the conditional expectation is more or less well fit to the data depending on how the penalities and complexity of the model is specified. The art and science of developing a GAM model is finding the constraints on the approximation that help capture the observed patterns but do not excessively overfit to the observed data.

Optimising The Parameter Setting

The manner of the the linear combination achieved by GAMs is constrained by the optimisation goal in their model fit. We can see from the model summary what is going on under the hood. For a given model specification the summary will report a number of model-fit statistics such as the log-likelihood and the AIC.

## Naive Model manually specified splines
gam_raw = LinearGAM(s(0, n_splines=10), fit_intercept=False).fit(X, y)
print("log_likelihood:", gam_raw.statistics_['loglikelihood'])
print("AIC:", gam_raw.statistics_['AIC'])

log_likelihood: -1042.2807399926558
AIC: 2097.7407642114244

The question then becomes what changes are induced in the model as we seek to optimise these model fit statistics.

## model optimised
gam = LinearGAM(s(0),  fit_intercept=False)
gam.gridsearch(X, y)
gam.summary()

  0% (0 of 11) |                         | Elapsed Time: 0:00:00 ETA:  --:--:--

  9% (1 of 11) |##                       | Elapsed Time: 0:00:00 ETA:   0:00:00

100% (11 of 11) |########################| Elapsed Time: 0:00:00 Time:  0:00:00

LinearGAM                                                                                                 
=============================================== ==========================================================
Distribution:                        NormalDist Effective DoF:                                     11.8135
Link Function:                     IdentityLink Log Likelihood:                                  -952.2409
Number of Samples:                          133 AIC:                                             1930.1088
                                                AICc:                                            1933.0789
                                                GCV:                                              609.3811
                                                Scale:                                            512.7965
                                                Pseudo R-Squared:                                   0.7984
==========================================================================================================
Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
================================= ==================== ============ ============ ============ ============
s(0)                              [0.2512]             20           11.8         1.11e-16     ***         
==========================================================================================================
Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

WARNING: Fitting splines and a linear function to a feature introduces a model identifiability problem
         which can cause p-values to appear significant when they are not.

WARNING: p-values calculated in this manner behave correctly for un-penalized models or models with
         known smoothing parameters, but when smoothing parameters have been estimated, the p-values
         are typically lower than they should be, meaning that the tests reject the null too readily.

Fortunately, this routine can be performed directly and results in the following differences between the naive and optimised model.

Plot the GAM fits

fig, ax = plt.subplots(figsize=(10, 6))
plot_fit(gam_raw, X, y, ax, "Unoptimised Fit", c1='orange', c2='green')
plot_fit(gam, X, y, ax, "Optimised Fit")

This is all well and good! We’ve seen an approach to modelling that can capture eccentric patterns in raw data. But how does it work and why should we care? If you’re familiar with the nomenclature of machine learning, you should think of spline modelling as a variety of feature creation. It generates “synthetic” features over the range of the observed variable. These synthetic features are the splines in question. Model building here too then relies on feature selection and principles for model-comparison.

Digression on the usage of “Splines”

The history of the term “spline” is related to the history of draftmanship. Historically splines were thin strips of flexible wood or plastic that could be bent or shaped around a weight or “knot” points to express a traceable curve over the space of a “numberline”. The elastic nature of the spline material allowed it to be bent around the knot points of curvature expressing a smooth or continuous bend.

The mathematical technique apes these properties by defining a curve over in an analogous way. We specify “knots” to carve up the support of the random variable \(X\) into portions that require different weighting schemes to represent the outcome \(y\) in each partition of the support variable.

Extracting the Splines

Returning to our model we can extract the spline features used in the PyGAM by invoking the following commands. We first identify the knot points and create the b-spline basis appropriate for the variable \(X\).

knot_edges=utils.gen_edge_knots(X,dtype='numerical')
knots=np.linspace(knot_edges[0],knot_edges[-1],len(gam.coef_))

splines = utils.b_spline_basis(X, edge_knots=knot_edges, sparse=False)

splines_df = pd.DataFrame(splines, columns=[f'basis_{i}' for i in range(len(gam.coef_))])

splines_df.head(10)

	basis_0	basis_1	basis_2	basis_3	basis_4	basis_5	basis_6	basis_7	basis_8	basis_9	basis_10	basis_11	basis_12	basis_13	basis_14	basis_15	basis_16	basis_17	basis_18	basis_19
0	0.166667	0.666667	0.166667	0.000000	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1	0.132997	0.661606	0.205334	0.000063	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2	0.059688	0.594827	0.341426	0.004059	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
3	0.030095	0.518726	0.437481	0.013698	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
4	0.012374	0.428013	0.527144	0.032470	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
5	0.000000	0.040337	0.551431	0.399315	0.008917	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
6	0.000000	0.018232	0.465467	0.492630	0.023671	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
7	0.000000	0.011137	0.418489	0.535407	0.034967	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
8	0.000000	0.000014	0.189310	0.664817	0.145859	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
9	0.000000	0.000000	0.120914	0.656897	0.222015	0.000174	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

These spline features range the extent of the of covariate space \(X\) defining “partitions” of the space. The model “learns” to capture the shape of the outcome variable \(y\) by figuring out how to weight the different portion of this spline basis matrix i.e. the linear combination of this basis matrix with the derived coefficients is a model of our outcome variable.

Note how each row is 0 everywhere except within the columns that represent a partition of \(X\). Additionally each row sums to unity. These properties are important because they ensure that any weighted combination of this basis represents the outcome variable in a controlled and quite granular manner. The more splines we use the more control we have of the representation, but we also risk overfit.

Plotting the Weighted Spline

ax = splines_df.dot(gam.coef_).plot(title='Weighted Splines', label='Weighted Combination of Spline Basis', figsize=(10, 6))
ax.set_ylabel("Acceleration")
ax.set_xlabel("Time Steps")
ax.legend();

In this manner we can see how the specification of a spline basis can help us model eccentric curves and waves in an outcome space. Spline features are like ad-hoc joists structuring our ship’s rigging. Features filling the gaps in our theory, they act like duct-tape binding the pieces together. If we hold this part down there, and jiggle another part here, our theory stays afloat. Spline features exemplify the contorted nature of empirically informed theory construction.

Next we’ll see how to more directly work with the specification of basis splines, passing these feature matrices into Bambi models.

Bayesian Splines with bambi

Spline models can be built and assessed using Bayesian approaches to quantify the uncertainty in the constructed function. Under the hood Bambi makes use of the ‘formulae’ package to allow for formula-like syntax to specify spline basis terms and mixed effect model terms. We leverage this to streamline the specification of spline models and demonstrate different kinds of spline-features.

knots_6 = np.linspace(0, np.max(X), 6+2)[1:-1]
knots_10 = np.linspace(0, np.max(X), 10+2)[1:-1]
knots_20 = np.linspace(0, np.max(X), 20+2)[1:-1]

df = pd.DataFrame({'X': X.flatten(), 'y': y})
formula1 = 'bs(X, degree=0, knots=knots_6)'
formula2 = 'bs(X, degree=1, knots=knots_6, intercept=False)'
formula3 = 'bs(X, degree=3, knots=knots_6, intercept=False)'
formula4 = 'bs(X, degree=3, knots=knots_10, intercept=False)'
formula5 = 'bs(X, degree=3, knots=knots_20, intercept=False)'
model_spline1 = bmb.Model(f"y ~ {formula1}", df)
model_spline2 = bmb.Model(f"y ~ {formula2}", df)
model_spline3 = bmb.Model(f"y ~ {formula3}", df)
model_spline4 = bmb.Model(f"y ~ {formula4}", df)
model_spline5 = bmb.Model(f"y ~ {formula5}", df)
model_spline5

       Formula: y ~ bs(X, degree=3, knots=knots_20, intercept=False)
        Family: gaussian
          Link: mu = identity
  Observations: 133
        Priors: 
    target = mu
        Common-level effects
            Intercept ~ Normal(mu: -25.5459, sigma: 222.6623)
            bs(X, degree=3, knots=knots_20, intercept=False) ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
                0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.], sigma: [1407.0172 1448.7263 1394.069   877.5692
                912.3295  679.8904  507.1974
              696.2455  860.7155  666.7815  665.863   932.178   883.0419  740.3925
             1010.3524  914.7406  923.9545 1218.535  1355.2211 1420.6241 1363.5314
             2403.9764 1393.2356])
        
        Auxiliary parameters
            sigma ~ HalfStudentT(nu: 4.0, sigma: 48.14)

As you can see here we are specifying a range of different degrees of spline basis. The different degrees corresspond to the smoothness of the overlapping splines. The degree=0 splines mean we specify a piecewise constant basis i.e. 0 or 1 within each region of the partition. But we can add more degrees to see more flexible representations of the space.

model_spline5.build()
model_spline5.graph()

The differences here will become clearer as we plot the various spline basis matrices below, but the main thought is that there are different degrees of “smoothness” to the linear basis of the spline features.

Plot the Spline Basis

The below functions extract the basis specification from each model and plots the basis design for an increasingly complex series of spline basis matrices.

def plot_spline_basis(basis, X, ax, title="Spline Basis"):
    df = (
        pd.DataFrame(basis)
        .assign(X=X)
        .melt("X", var_name="basis_idx", value_name="y")
    )


    for idx in df.basis_idx.unique():
        d = df[df.basis_idx == idx]
        ax.plot(d["X"], d["y"])
    
    ax.set_title(title)
    return ax

def plot_knots(knots, ax):
    for knot in knots:
        ax.axvline(knot, color="0.1", alpha=0.4)
    return ax



fig, axs = plt.subplots(5, 1, figsize=(9, 20))
axs = axs.flatten()
axs.flatten()
B1 = model_spline1.response_component.design.common[formula1]
plot_spline_basis(B1, df["X"].values, ax=axs[0], title="Piecewise Constant Basis")
plot_knots(knots_6, axs[0]);

B2 = model_spline2.response_component.design.common[formula2]
ax = plot_spline_basis(B2, df["X"].values, axs[1], 
title="Piecewise Linear Basis")
plot_knots(knots_6, axs[1]);

B3 = model_spline3.response_component.design.common[formula3]
ax = plot_spline_basis(B3, df["X"].values, axs[2], 
title="Cubic Spline Basis (6 Knots)")
plot_knots(knots_6, axs[2]);

B4 = model_spline4.response_component.design.common[formula4]
ax = plot_spline_basis(B4, df["X"].values, axs[3], 
title="Cubic Spline Basis (10 Knots)")
plot_knots(knots_10, axs[3]);


B5 = model_spline5.response_component.design.common[formula5]
ax = plot_spline_basis(B5, df["X"].values, axs[4], 
title="Cubic Spline Basis (20 Knots)")
plot_knots(knots_20, axs[4]);

Here we’ve seen the nature of the modelling splines we’ll fit to the data. By “fit” we mean estimate a set of linear coefficients that we can use to additively combine these spline features creating a smooth linear representation of the outcome.

Fit the Individual Spline Models

We now combine the spline components within a linear model fit, pulling and tying them together in such a way as to ape the shape of the observed sequence.

idata_spline1 = model_spline1.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline2 = model_spline2.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline3 = model_spline3.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline4 = model_spline4.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline5 = model_spline5.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

Finally we can plot fits achieved by each of the models and compare how the different spline smooths contribute to approxiating the outcome variable.

Plot the Weighted Mean

def plot_weighted_splines(B, idata, formula, ax, knots):
    posterior_stacked = az.extract(idata)
    wp = posterior_stacked[formula].mean("sample").values

    plot_spline_basis(B * wp.T, df["X"].values, ax)
    ax.plot(df.X.values, np.dot(B, wp.T), color="black", lw=3, label='Weighted Splines')
    plot_knots(knots, ax);
    ax.legend()



fig, axs = plt.subplots(5, 1, figsize=(10, 20))
axs = axs.flatten()
axs.flatten()

plot_weighted_splines(B1, idata_spline1, formula1, axs[0], knots_6)
plot_weighted_splines(B2, idata_spline2, formula2, axs[1], knots_6)
plot_weighted_splines(B3, idata_spline3, formula3, axs[2], knots_6)
plot_weighted_splines(B4, idata_spline4, formula4, axs[3], knots_10)
plot_weighted_splines(B5, idata_spline5, formula5, axs[4], knots_20)

Here we can see how the models with increasingly complex splines are more exactly able to fit the herky jerky trajectory of the outcome variable in each interval. The fewer the intervals, the less flexibility available to the model.

Compare Model Fits

As before we can evaluate these model fits and compare them based on leave-one-out cross validation scores and information theoretic complexity measures.

fig, axs = plt.subplots(1, 2, figsize=(10, 7))
axs = axs.flatten()
models_dict = {"piecewise_constant": idata_spline1, "piecewise_linear": idata_spline2, "cubic_bspline": idata_spline3, "cubic_bspline_10": idata_spline4, 
"cubic_bspline_20": idata_spline5}
df_compare = az.compare(models_dict)
az.plot_compare(df_compare, ax=axs[0])
axs[0].get_legend().remove()
az.plot_compare(az.compare(models_dict, 'waic'), ax=axs[1])
axs[1].set_yticklabels([])
df_compare

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
cubic_bspline_10	0	-612.572745	11.646898	0.000000	9.011087e-01	9.646678	0.000000	True	log
cubic_bspline_20	1	-620.709337	19.669201	8.136592	1.360171e-13	9.600116	2.074742	True	log
cubic_bspline	2	-634.647180	8.703519	22.074435	6.196578e-14	8.915728	6.642850	True	log
piecewise_constant	3	-643.781042	6.981098	31.208296	5.299668e-02	9.770740	8.728072	False	log
piecewise_linear	4	-647.016885	5.987587	34.444140	4.589467e-02	7.914787	8.116163	False	log

Here we see that the extra complexity of using 20 splines leads to slightly worse performance measures than the less complex but seemingly adequate 10 splines. In other words, it’s starting to overfit to the data at 20. This is the price of flexibility and a sign of a model unlikely to generalise well on out-of-sample data.

new_data = pd.DataFrame({"X": np.linspace(df.X.min() - 5, df.X.max() + 5, num=500)})
    
model_spline4.predict(idata_spline4, data=new_data, 
kind='pps', inplace=True)

idata_spline4

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:                                               (chain: 4,
                                                             draw: 1000,
                                                             bs(X, degree=3, knots=knots_10, intercept=False)_dim: 13,
                                                             y_obs: 500)
  Coordinates:
    * chain                                                 (chain) int64 0 1 2 3
    * draw                                                  (draw) int64 0 ... 999
    * bs(X, degree=3, knots=knots_10, intercept=False)_dim  (bs(X, degree=3, knots=knots_10, intercept=False)_dim) int64 ...
    * y_obs                                                 (y_obs) int64 0 ......
  Data variables:
      Intercept                                             (chain, draw) float64 ...
      bs(X, degree=3, knots=knots_10, intercept=False)      (chain, draw, bs(X, degree=3, knots=knots_10, intercept=False)_dim) float64 ...
      y_sigma                                               (chain, draw) float64 ...
      y_mean                                                (chain, draw, y_obs) float64 ...
  Attributes:
      created_at:                  2024-05-27T06:00:45.172558
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      sampling_time:               2.553478956222534
      tuning_steps:                1000
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • bs(X, degree=3, knots=knots_10, intercept=False)_dim: 13
    • y_obs: 500
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • bs(X, degree=3, knots=knots_10, intercept=False)_dim
      (bs(X, degree=3, knots=knots_10, intercept=False)_dim)
      int64
      0 1 2 3 4 5 6 7 8 9 10 11 12

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (4)
    • Intercept
      (chain, draw)
      float64
      -30.48 -9.768 17.5 ... -4.816 -3.18

      array([[-30.47836564,  -9.76829627,  17.49712677, ...,  17.81763441,
               13.95495119,  -6.10602699],
             [-16.23316086,   9.49899062,   3.70481638, ..., -25.03974245,
               -2.91235316, -18.76784286],
             [ -0.38968653, -23.55404923,  21.85558793, ..., -15.88911032,
              -12.17693056, -22.799846  ],
             [-23.78536816, -10.33769596, -18.72531118, ...,  10.30614657,
               -4.81605641,  -3.18037091]])

    • bs(X, degree=3, knots=knots_10, intercept=False)
      (chain, draw, bs(X, degree=3, knots=knots_10, intercept=False)_dim)
      float64
      60.08 -26.33 72.03 ... 26.62 -19.5

      array([[[ 60.08361849, -26.32632766,  72.03440963, ..., -24.69505575,
                41.29966173, 126.0957609 ],
              [ 42.01923744,  13.97976273, -10.66570712, ...,  37.69872884,
                32.66180174,  56.71948798],
              [-20.79003642, -60.86218268,   7.89670524, ..., -79.792871  ,
                26.13817326,   5.2715798 ],
              ...,
              [-53.05009482, -44.90737225,  -0.53936709, ...,  11.86630968,
               -30.56679357, -47.68641912],
              [-31.19449725, -22.73967732,  20.06162131, ..., -28.10157326,
                 4.07136337, -21.05300466],
              [ 35.06057706,  -9.59093606,  11.09479783, ...,  -3.88987727,
                -8.987879  ,  27.57356096]],
      
             [[ 33.08556413,  -1.11856386,  24.73543073, ...,  73.28686802,
               -24.37794241,  33.15302761],
              [-28.57730471,  -3.73441648,  -2.57506729, ..., -18.45867357,
                -4.60293539,  -5.79920437],
              [ -4.42415991, -12.80936972,  20.31341056, ...,   1.64870928,
                19.07410069,  26.02601531],
      ...
              [  9.81825793,  16.07126062,  40.10272142, ..., -21.65327653,
                34.82747494,  28.59603009],
              [ 50.49534197, -16.63350377,  36.81576572, ...,  36.77027802,
                 8.30573837,  15.51060296],
              [ 33.01164624,  -1.33302557,  38.41471255, ...,   3.94749231,
                46.54871666,  43.92361804]],
      
             [[ 18.45099645, -28.24611558,  49.26005839, ...,   4.79065054,
                50.07280024, -17.17531704],
              [ 24.27302084,  12.70125656,  13.48983162, ...,  23.33510231,
               -38.49057179, 101.69703639],
              [ 22.36020074, -13.16116681,  32.08258294, ...,  56.15103319,
                 4.60000597,  72.46523374],
              ...,
              [-71.35771382,  25.42919025, -19.3976876 , ..., -30.2130114 ,
                 5.95767059, -25.33043675],
              [ 10.9586426 , -46.45830404,  45.31732601, ...,  47.8378147 ,
               -28.12282789,  41.36265809],
              [-31.53039023,  20.83808587, -18.46035112, ...,  -4.39658487,
                26.61753145, -19.50485034]]])

    • y_sigma
      (chain, draw)
      float64
      22.77 23.59 21.69 ... 20.98 25.07

      array([[22.77247029, 23.5905779 , 21.68882473, ..., 21.91856404,
              23.1234507 , 20.69603421],
             [23.30953138, 21.78668016, 23.3377476 , ..., 23.554608  ,
              22.30139546, 23.25599855],
             [25.0826106 , 24.14700624, 23.17634169, ..., 22.55796555,
              22.00396287, 23.83528365],
             [22.97336709, 22.34507905, 24.09241758, ..., 23.84889601,
              20.97955574, 25.07045989]])

    • y_mean
      (chain, draw, y_obs)
      float64
      -1.786e+03 -1.688e+03 ... -384.5

      array([[[-1785.63598357, -1688.46091242, -1594.80861303, ...,
                 471.11786457,   486.11321041,   501.38658015],
              [-1008.56095247,  -955.43500036,  -904.13867636, ...,
                 209.04228219,   217.17041391,   225.5454136 ],
              [  171.46343997,   166.77923669,   162.10649158, ...,
                -298.9913081 ,  -317.81225915,  -337.30145789],
              ...,
              [ 1113.19763849,  1056.55963212,  1001.80503673, ...,
                 -46.36650365,   -45.39216019,   -44.30916241],
              [  672.75237366,   638.6237798 ,   605.62977443, ...,
                -210.41648402,  -220.83665251,  -231.5695011 ],
              [ -981.31984966,  -927.88664463,  -876.36327872, ...,
                 245.7530471 ,   256.32233432,   267.17337954]],
      
             [[ -892.49690627,  -844.79097031,  -798.77984615, ...,
                 543.89918243,   572.12293189,   601.25719405],
              [  729.10432472,   690.36248691,   652.97626248, ...,
                 -29.69218582,   -31.76353623,   -33.91297753],
              [   29.72027438,    29.25039646,    28.75480424, ...,
                  45.18641388,    45.54217006,    45.89560459],
      ...
              [ -186.24078895,  -177.60095429,  -169.24721219, ...,
                -129.71594267,  -138.38733622,  -147.37882584],
              [-1442.9878318 , -1364.25391604, -1288.35121975, ...,
                 100.15273049,   105.7273578 ,   111.49734282],
              [ -904.10355066,  -855.98742694,  -809.58838617, ...,
                 -72.7318446 ,   -78.59078957,   -84.67012043]],
      
             [[ -700.6554728 ,  -661.92143273,  -624.64606141, ...,
                -521.29646421,  -545.21735853,  -569.83464216],
              [ -566.28543846,  -536.82610575,  -508.37957578, ...,
                1032.55566447,  1078.47471162,  1125.69033827],
              [ -695.22027979,  -657.54938003,  -621.25358394, ...,
                 556.45783286,   581.91710165,   608.14063371],
              ...,
              [ 2030.9713329 ,  1919.9629329 ,  1812.93379418, ...,
                -265.67005783,  -278.5688838 ,  -291.86028771],
              [ -600.01947051,  -564.67068253,  -530.70395972, ...,
                 590.9906349 ,   619.60614144,   649.1001913 ],
              [  954.09836044,   900.90776771,   849.64997795, ...,
                -351.31548764,  -367.65439889,  -384.46730166]]])

  • Indexes: (4)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • bs(X, degree=3, knots=knots_10, intercept=False)_dim
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], dtype='int64', name='bs(X, degree=3, knots=knots_10, intercept=False)_dim'))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (8)

    created_at :

    2024-05-27T06:00:45.172558

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    2.553478956222534

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 500)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 492 493 494 495 496 497 498 499
  Data variables:
      y        (chain, draw, y_obs) float64 -1.769e+03 -1.707e+03 ... -438.7
  Attributes:
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 500
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      -1.769e+03 -1.707e+03 ... -438.7

      array([[[-1.76850707e+03, -1.70683103e+03, -1.58503478e+03, ...,
                4.75269803e+02,  4.71717129e+02,  4.93095980e+02],
              [-9.78963706e+02, -9.57350338e+02, -8.94242702e+02, ...,
                2.23703058e+02,  1.99989491e+02,  2.36217010e+02],
              [ 1.51236792e+02,  1.44573452e+02,  1.59874803e+02, ...,
               -2.92521983e+02, -3.39277905e+02, -3.39765298e+02],
              ...,
              [ 1.12245700e+03,  1.02901445e+03,  1.00163629e+03, ...,
               -5.76818185e+01, -5.79655659e+01, -6.95232453e+01],
              [ 7.03262674e+02,  6.89295681e+02,  5.72987411e+02, ...,
               -2.19231936e+02, -1.91846651e+02, -1.77780985e+02],
              [-1.00564975e+03, -9.26117757e+02, -8.32114624e+02, ...,
                2.43447636e+02,  2.46073703e+02,  2.61526238e+02]],
      
             [[-8.78027516e+02, -8.39940972e+02, -7.70858431e+02, ...,
                5.62919028e+02,  5.81119362e+02,  5.98368471e+02],
              [ 7.59822216e+02,  7.35672407e+02,  6.66036776e+02, ...,
               -3.86148939e+01, -2.62397986e+01, -7.83702403e+01],
              [ 6.01153527e+01, -8.09278342e-01, -4.14103639e+01, ...,
                1.77331558e+01,  3.35647548e+01,  8.67001318e+01],
      ...
              [-2.03929257e+02, -1.59688749e+02, -1.67829907e+02, ...,
               -1.26352656e+02, -1.58838619e+02, -1.33905356e+02],
              [-1.42981773e+03, -1.35031319e+03, -1.30844300e+03, ...,
                1.11580849e+02,  8.39822162e+01,  1.55903514e+02],
              [-9.65920008e+02, -8.62270119e+02, -7.92953633e+02, ...,
               -5.96365921e+01, -5.50070341e+01, -9.95204835e+01]],
      
             [[-7.28155125e+02, -6.77887003e+02, -5.89951035e+02, ...,
               -5.12375162e+02, -5.30806452e+02, -6.00515410e+02],
              [-5.86770932e+02, -5.53033193e+02, -4.89166246e+02, ...,
                9.92634196e+02,  1.04493361e+03,  1.12704971e+03],
              [-7.01773152e+02, -6.43685059e+02, -5.89255386e+02, ...,
                5.27179378e+02,  5.58208456e+02,  5.88341537e+02],
              ...,
              [ 2.02804805e+03,  1.90391514e+03,  1.80455354e+03, ...,
               -2.68343325e+02, -2.78799372e+02, -2.96707748e+02],
              [-5.94745406e+02, -5.62236707e+02, -5.62459190e+02, ...,
                6.09916668e+02,  6.51977577e+02,  6.51319039e+02],
              [ 9.61991994e+02,  9.05536001e+02,  8.37728249e+02, ...,
               -3.91949049e+02, -3.49475144e+02, -4.38689017e+02]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (2)

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 133)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (chain, draw, y_obs) float64 -4.94 -4.349 -4.059 ... -4.232 -5.027
  Attributes:
      created_at:                  2024-05-27T06:00:45.279798
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 133
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      -4.94 -4.349 ... -4.232 -5.027

      array([[[ -4.94013004,  -4.34926962,  -4.05923397, ...,  -4.07989994,
                -4.32526297, -10.99702367],
              [ -4.16551546,  -4.0798865 ,  -4.34890397, ...,  -4.27847476,
                -4.90199179,  -5.2604801 ],
              [ -4.32114608,  -4.21922923,  -4.01523436, ...,  -4.00171403,
                -4.46666592,  -4.15055295],
              ...,
              [ -4.33667694,  -4.08524041,  -4.16193486, ...,  -4.25075843,
                -4.0367804 ,  -5.71916378],
              [ -4.2418906 ,  -4.13815308,  -4.06475279, ...,  -4.08566333,
                -4.13276703,  -4.35600265],
              [ -3.99240307,  -3.95297982,  -4.22738498, ...,  -4.39489511,
                -3.96850328,  -4.08422128]],
      
             [[ -4.31029941,  -4.1356895 ,  -4.08943162, ...,  -4.39621813,
                -4.11988518,  -4.10340212],
              [ -4.09528528,  -4.02978374,  -4.0169432 , ...,  -4.09620393,
                -4.02276603,  -4.05185642],
              [ -4.08161109,  -4.08438013,  -4.07760521, ...,  -4.11690648,
                -4.53977755,  -4.40149229],
      ...
              [ -4.28309395,  -4.19158751,  -4.07661987, ...,  -4.08978805,
                -4.11622924,  -4.03898436],
              [ -4.16328531,  -4.01182659,  -4.32670151, ...,  -4.05121559,
                -4.0454269 ,  -4.0661975 ],
              [ -4.54760772,  -4.2926589 ,  -4.09280683, ...,  -4.09149036,
                -4.35227017,  -4.18573222]],
      
             [[ -4.58924409,  -4.39837658,  -4.20729252, ...,  -4.07955443,
                -4.09450053,  -6.58164643],
              [ -4.13256189,  -4.0441922 ,  -4.07194408, ...,  -4.64029899,
                -4.06235348, -10.54056155],
              [ -4.40287769,  -4.25092797,  -4.11800136, ...,  -4.11077154,
                -4.21757925,  -5.69653922],
              ...,
              [ -4.18405048,  -4.09405162,  -4.6372543 , ...,  -4.15126206,
                -4.11639733,  -4.67240483],
              [ -3.98883574,  -3.9655664 ,  -3.96310644, ...,  -4.21756776,
                -3.96543135,  -4.72138657],
              [ -4.14867518,  -4.1893545 ,  -4.35918185, ...,  -4.14286099,
                -4.23197849,  -5.02728226]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (6)

    created_at :

    2024-05-27T06:00:45.279798

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:                (chain: 4, draw: 1000)
  Coordinates:
    * chain                  (chain) int64 0 1 2 3
    * draw                   (draw) int64 0 1 2 3 4 5 ... 994 995 996 997 998 999
  Data variables: (12/17)
      perf_counter_diff      (chain, draw) float64 0.0007404 ... 0.0007515
      lp                     (chain, draw) float64 -721.3 -715.6 ... -710.3 -711.0
      step_size_bar          (chain, draw) float64 0.3142 0.3142 ... 0.3369 0.3369
      smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan
      acceptance_rate        (chain, draw) float64 0.9519 0.8589 ... 0.7224 0.9861
      reached_max_treedepth  (chain, draw) bool False False False ... False False
      ...                     ...
      largest_eigval         (chain, draw) float64 nan nan nan nan ... nan nan nan
      energy                 (chain, draw) float64 728.6 729.2 ... 725.7 714.8
      tree_depth             (chain, draw) int64 4 4 4 4 3 4 4 5 ... 4 5 4 3 4 4 4
      diverging              (chain, draw) bool False False False ... False False
      index_in_trajectory    (chain, draw) int64 -7 -13 -14 -3 6 ... -3 4 -4 7 8
      process_time_diff      (chain, draw) float64 0.000741 0.000737 ... 0.000751
  Attributes:
      created_at:                  2024-05-27T06:00:45.180249
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      sampling_time:               2.553478956222534
      tuning_steps:                1000
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (17)
    • perf_counter_diff
      (chain, draw)
      float64
      0.0007404 0.0007361 ... 0.0007515

      array([[0.00074042, 0.00073613, 0.00073075, ..., 0.00072508, 0.00071846,
              0.00071829],
             [0.00072588, 0.00143679, 0.00037763, ..., 0.00073321, 0.00073304,
              0.00073225],
             [0.00073129, 0.00073483, 0.00108604, ..., 0.00142812, 0.00037529,
              0.00073333],
             [0.00144608, 0.00037754, 0.0003765 , ..., 0.00075521, 0.00076121,
              0.0007515 ]])

    • lp
      (chain, draw)
      float64
      -721.3 -715.6 ... -710.3 -711.0

      array([[-721.32325798, -715.60693507, -712.81600649, ..., -711.63567075,
              -711.86251963, -711.25215936],
             [-711.83173082, -714.18292872, -708.88038704, ..., -710.81095225,
              -707.61049434, -710.43240125],
             [-712.7395976 , -713.40853693, -713.62383765, ..., -712.97340089,
              -709.41442383, -710.14460595],
             [-714.37055597, -716.70641237, -710.51320632, ..., -713.81412225,
              -710.25741961, -710.97040762]])

    • step_size_bar
      (chain, draw)
      float64
      0.3142 0.3142 ... 0.3369 0.3369

      array([[0.31415288, 0.31415288, 0.31415288, ..., 0.31415288, 0.31415288,
              0.31415288],
             [0.33735905, 0.33735905, 0.33735905, ..., 0.33735905, 0.33735905,
              0.33735905],
             [0.35905043, 0.35905043, 0.35905043, ..., 0.35905043, 0.35905043,
              0.35905043],
             [0.33686986, 0.33686986, 0.33686986, ..., 0.33686986, 0.33686986,
              0.33686986]])

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • acceptance_rate
      (chain, draw)
      float64
      0.9519 0.8589 ... 0.7224 0.9861

      array([[0.95188821, 0.85891772, 0.71451549, ..., 0.8702077 , 0.50215875,
              0.97026438],
             [0.71509145, 0.51444105, 0.90341721, ..., 0.79993904, 0.57204968,
              0.66470717],
             [0.53066916, 0.66611802, 0.72800025, ..., 0.61955158, 0.92591962,
              0.77854416],
             [0.99982809, 0.97668472, 0.87164137, ..., 0.74461051, 0.72243678,
              0.98614764]])

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • step_size
      (chain, draw)
      float64
      0.2833 0.2833 ... 0.3204 0.3204

      array([[0.28329617, 0.28329617, 0.28329617, ..., 0.28329617, 0.28329617,
              0.28329617],
             [0.23729018, 0.23729018, 0.23729018, ..., 0.23729018, 0.23729018,
              0.23729018],
             [0.2420878 , 0.2420878 , 0.2420878 , ..., 0.2420878 , 0.2420878 ,
              0.2420878 ],
             [0.32038673, 0.32038673, 0.32038673, ..., 0.32038673, 0.32038673,
              0.32038673]])

    • perf_counter_start
      (chain, draw)
      float64
      2.361e+06 2.361e+06 ... 2.361e+06

      array([[2360978.51092729, 2360978.51171625, 2360978.51249971, ...,
              2360979.30091483, 2360979.30168592, 2360979.30245238],
             [2360978.47167667, 2360978.47245279, 2360978.47393825, ...,
              2360979.25252446, 2360979.2533045 , 2360979.2540865 ],
             [2360978.51288592, 2360978.51366746, 2360978.51445129, ...,
              2360979.26986125, 2360979.27133717, 2360979.27175967],
             [2360978.46373525, 2360978.46523   , 2360978.46565383, ...,
              2360979.27246896, 2360979.27327358, 2360979.27464221]])

    • n_steps
      (chain, draw)
      float64
      15.0 15.0 15.0 ... 15.0 15.0 15.0

      array([[15., 15., 15., ..., 15., 15., 15.],
             [15., 31.,  7., ..., 15., 15., 15.],
             [15., 15., 23., ..., 31.,  7., 15.],
             [31.,  7.,  7., ..., 15., 15., 15.]])

    • max_energy_error
      (chain, draw)
      float64
      -1.01 0.7776 ... 0.8487 -0.3605

      array([[-1.00971736,  0.77757398,  1.28280502, ...,  0.37495702,
               1.17970654, -0.34836573],
             [ 1.3345751 ,  2.3915387 , -1.06730429, ...,  1.15497563,
               1.75653514,  0.89915678],
             [ 1.40682187,  1.52508573,  1.19349047, ...,  1.84259665,
              -1.22585309,  0.93739617],
             [-1.54429096, -0.75846954,  0.55295383, ...,  0.76575828,
               0.84873093, -0.36049323]])

    • energy_error
      (chain, draw)
      float64
      -0.2612 -0.1061 ... -0.1596 -0.1877

      array([[-0.26121615, -0.10611572, -0.69771087, ...,  0.06849994,
               0.4074744 , -0.15667011],
             [-0.29476717,  0.28462281, -1.06730429, ..., -0.313449  ,
              -0.58684002,  0.26092216],
             [ 0.37035487,  0.68575505,  0.05343487, ...,  0.74127564,
              -1.20632225,  0.29989649],
             [-0.85058098, -0.46592943, -0.24132289, ..., -0.12983655,
              -0.15964365, -0.18767355]])

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • energy
      (chain, draw)
      float64
      728.6 729.2 721.0 ... 725.7 714.8

      array([[728.6255778 , 729.24609198, 720.95261511, ..., 718.23864589,
              723.51162522, 718.37013269],
             [725.96636541, 725.37893664, 719.42266114, ..., 721.44682066,
              717.94143164, 714.93926686],
             [724.0278209 , 722.19477874, 721.36794365, ..., 720.83274011,
              716.33821598, 715.0968274 ],
             [723.0827726 , 724.98007703, 723.47496067, ..., 724.58398944,
              725.68406617, 714.81290627]])

    • tree_depth
      (chain, draw)
      int64
      4 4 4 4 3 4 4 5 ... 5 4 5 4 3 4 4 4

      array([[4, 4, 4, ..., 4, 4, 4],
             [4, 5, 3, ..., 4, 4, 4],
             [4, 4, 5, ..., 5, 3, 4],
             [5, 3, 3, ..., 4, 4, 4]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • index_in_trajectory
      (chain, draw)
      int64
      -7 -13 -14 -3 6 -8 ... -3 4 -4 7 8

      array([[ -7, -13, -14, ...,   5,   3,   5],
             [  6,  -4,  -4, ...,  -8,  14, -13],
             [ -7, -13, -14, ...,  10,  -5,   8],
             [ -8,  -5,   4, ...,  -4,   7,   8]])

    • process_time_diff
      (chain, draw)
      float64
      0.000741 0.000737 ... 0.000751

      array([[0.000741, 0.000737, 0.000731, ..., 0.000724, 0.000719, 0.000718],
             [0.000725, 0.001437, 0.000377, ..., 0.000734, 0.000732, 0.000732],
             [0.000731, 0.000735, 0.001086, ..., 0.001429, 0.000375, 0.000734],
             [0.001446, 0.000378, 0.000377, ..., 0.000756, 0.000762, 0.000751]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (8)

    created_at :

    2024-05-27T06:00:45.180249

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    2.553478956222534

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (y_obs: 133)
  Coordinates:
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (y_obs) float64 0.0 -1.3 -2.7 0.0 -2.7 ... -2.7 10.7 -2.7 10.7
  Attributes:
      created_at:                  2024-05-27T06:00:45.182737
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • y_obs: 133
  • Coordinates: (1)
    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (y_obs)
      float64
      0.0 -1.3 -2.7 ... 10.7 -2.7 10.7

      array([   0. ,   -1.3,   -2.7,    0. ,   -2.7,   -2.7,   -2.7,   -1.3,
               -2.7,   -2.7,   -1.3,   -2.7,   -2.7,   -2.7,   -5.4,   -2.7,
               -5.4,    0. ,   -2.7,   -2.7,    0. ,  -13.3,   -5.4,   -5.4,
               -9.3,  -16. ,  -22.8,   -2.7,  -22.8,  -32.1,  -53.5,  -54.9,
              -40.2,  -21.5,  -21.5,  -50.8,  -42.9,  -26.8,  -21.5,  -50.8,
              -61.7,   -5.4,  -80.4,  -59. ,  -71. ,  -91.1,  -77.7,  -37.5,
              -85.6, -123.1, -101.9,  -99.1, -104.4, -112.5,  -50.8, -123.1,
              -85.6,  -72.3, -127.2, -123.1, -117.9, -134. , -101.9, -108.4,
             -123.1, -123.1, -128.5, -112.5,  -95.1,  -81.8,  -53.5,  -64.4,
              -57.6,  -72.3,  -44.3,  -26.8,   -5.4, -107.1,  -21.5,  -65.6,
              -16. ,  -45.6,  -24.2,    9.5,    4. ,   12. ,  -21.5,   37.5,
               46.9,  -17.4,   36.2,   75. ,    8.1,   54.9,   48.2,   46.9,
               16. ,   45.6,    1.3,   75. ,  -16. ,  -54.9,   69.6,   34.8,
               32.1,  -37.5,   22.8,   46.9,   10.7,    5.4,   -1.3,  -21.5,
              -13.3,   30.8,  -10.7,   29.4,    0. ,  -10.7,   14.7,   -1.3,
                0. ,   10.7,   10.7,  -26.8,  -14.7,  -13.3,    0. ,   10.7,
              -14.7,   -2.7,   10.7,   -2.7,   10.7])

  • Indexes: (1)
    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (6)

    created_at :

    2024-05-27T06:00:45.182737

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  

Next we plot the posterior predictive distribution of our observed variable and compare against the observed data. Additionally we plot the 89th and 50% HDI.

ax = az.plot_hdi(new_data['X'], idata_spline4['posterior_predictive']['y'], fill_kwargs={'alpha': 0.2, 'color':'firebrick'}, hdi_prob=0.89, figsize=(10, 8))

az.plot_hdi(new_data['X'], idata_spline4['posterior_predictive']['y'], fill_kwargs={'alpha': 0.8, 'color':'firebrick'}, hdi_prob=0.5)

y_mean = idata_spline4['posterior_predictive']['y'].mean(dim=('chain', 'draw'))

ax.plot(new_data['X'], y_mean, label='Expected posterior predictive', color='k')
ax.set_xlabel("Time Point")
ax.set_ylabel("Acceleration")

ax.scatter(df['X'], df['y'], label='Observed Datapoints')

ax.legend()

ax.set_title("Posterior Predictive Distribution \n Based on 10 Knots");

This represents a good a clean model fit to the observed data using univariate spline smoothers. However, it’s clear that the uncertainty spikes massively with data points even slightly out of the training data. Next we’ll see another alternative approach to model this outcome variable using approximate gaussian processes.

Approximate Gaussian processes

Gaussian processes models allow us to express our beliefs about the structure of a function and calibrate these beliefs against the observed data. The topic of gaussian processes is rich and complex. Too rich to be fairly covered in this blog post, so we’ll just say that we’re using a method designed for function approximation that makes use of drawing samples from a multivariate normal distribution under a range of different covariance relationships. The properties of these covariance relations determine the shape and fluidity of the function realisations. We will use the HSGP approximation of gaussian processes models.

The gaussian process models can be somewhat interrogated by examining the different combinations of covariance relationships with priors over the parameters governing the covariance of a sequence of points. For example consider the following parameterisations.

fig, axs = plt.subplots(2, 2, figsize=(9, 10))
axs = axs.flatten()

def plot_cov_draws(ax1, ax2, lengthscale=3, sigma=13):
    cov = sigma**2 * pm.gp.cov.ExpQuad(1, lengthscale)
    X = np.linspace(0, 60, 200)[:, None]
    K = cov(X).eval()

    sns.heatmap(pd.DataFrame(K), center=0, xticklabels=[], yticklabels=[], ax=ax1, cmap='crest');
    ax1.set_title(f"Covariance Length Scale {lengthscale}")
    ax2.plot(
        X,
        pm.draw(
            pm.MvNormal.dist(mu=np.zeros(len(K)), cov=K, shape=K.shape[0]), draws=10, random_seed=random_seed
        ).T, color='blue', alpha=0.5
    )
    ax2.set_title(f"Samples from the GP prior \n lengthscale: {lengthscale}, sigma: {sigma}")
    plt.ylabel("y")
    plt.xlabel("X");

plot_cov_draws(axs[0], axs[1])
plot_cov_draws(axs[2], axs[3], lengthscale=10, sigma=13)

We’ve specified the range of X to reflect the support of the acceleration example and allowed the draws to be informed by a covariance function we have parameterised using the Exponentiated Quadratic kernel. We’ve additionally tweaked the lengthscale to demonstrate how proximity in the sequence informs the kernel and determines the shape of the covariance structure.

\[k(x, x') = \mathrm{exp}\left[ -\frac{(x - x')^2}{2 \ell^2} \right]\]

The plot on the left highlights the importance of the lengthscale parameter over the sequence. The wider the central shading is the more pronounced is the correlation among more points in the sequence. The patterns to the right show a good range of “wiggliness” that they should be flexible enough to capture the shape of the acceleration. These are the constraints of structure we can impose of the theoretical realisations of our model. If we can calibrate the model and derive posterior parameters against the observed data, we can learn the probable shape of our functional relationship.

Priors on Gaussian Processes

Consider the following specification for the priors

fig, axs = plt.subplots(1, 2, figsize=(9, 6))
axs = axs.flatten()
axs[0].hist(pm.draw(pm.InverseGamma.dist(mu=1, sigma=1), 1000), ec='black', bins=30);
axs[0].set_title("Priors for Lengthscale \n in ExpQuad Kernel")
axs[1].hist(pm.draw(pm.Exponential.dist(lam=1), 1000), ec='black', bins=30);
axs[1].set_title("Priors for Amplitude \n in ExpQuad Kernel")

Text(0.5, 1.0, 'Priors for Amplitude \n in ExpQuad Kernel')

We use these to specify priors on the Hilbert space approximation of gaussian priors available in the Bambi package.

prior_hsgp = {
    "sigma": bmb.Prior("Exponential", lam=1), # amplitude
    "ell": bmb.Prior("InverseGamma", mu=1, sigma=1) # lengthscale
}

# This is the dictionary we pass to Bambi
priors = {
    "hsgp(X, m=10, c=1)": prior_hsgp,
    "sigma": bmb.Prior("HalfNormal", sigma=4)
}
model_hsgp = bmb.Model("y ~ 0 + hsgp(X, m=10, c=1)", df, priors=priors)
model_hsgp

       Formula: y ~ 0 + hsgp(X, m=10, c=1)
        Family: gaussian
          Link: mu = identity
  Observations: 133
        Priors: 
    target = mu
        HSGP contributions
            hsgp(X, m=10, c=1)
                cov: ExpQuad
                sigma ~ Exponential(lam: 1.0)
                ell ~ InverseGamma(mu: 1.0, sigma: 1.0)
        
        Auxiliary parameters
            sigma ~ HalfNormal(sigma: 4.0)

Here we’ve set the m=10 to determine the number of basis vectors used in the Hilbert space approximation. The idea differs in detail from the spline based approximations we’ve seen, but it’s perhaps useful to think of the process in the same vein. Here again we have a theory of the world expressed as a function of opaque components jerry-rigged for modelling some phenomena.

idata_hsgp = model_hsgp.fit(inference_method="nuts_numpyro",target_accept=0.95, random_seed=121195, 
idata_kwargs={"log_likelihood": True})
print(idata_hsgp.sample_stats["diverging"].sum().to_numpy())

Compiling...

Compilation time = 0:00:01.138923

Sampling...

Sampling time = 0:00:01.870798

Transforming variables...

Transformation time = 0:00:00.176689

Computing Log Likelihood...

Log Likelihood time = 0:00:00.144039

0

This model fits and the sampling seems to have worked well.

az.plot_trace(idata_hsgp, backend_kwargs={"layout": "constrained"}, figsize=(9, 15));

The lengthscale and sigma parameters we have learned by calibrating our priors against the data. The degree to which these parameters are meaningful depend a little on how familar you are with covariance matrix kernels and their properties, so we won’t dwell on the point here.

az.summary(idata_hsgp, var_names=['hsgp(X, m=10, c=1)_ell', 'hsgp(X, m=10, c=1)_sigma', 'y_sigma', 'hsgp(X, m=10, c=1)_weights'])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
hsgp(X, m=10, c=1)_ell	3.301	0.666	2.099	4.567	0.016	0.011	1763.0	1815.0	1.0
hsgp(X, m=10, c=1)_sigma	23.675	2.679	18.925	28.847	0.069	0.049	1497.0	2167.0	1.0
y_sigma	20.589	1.121	18.590	22.751	0.019	0.014	3440.0	2485.0	1.0
hsgp(X, m=10, c=1)_weights[0]	-131.588	13.784	-157.019	-105.189	0.216	0.154	4071.0	3020.0	1.0
hsgp(X, m=10, c=1)_weights[1]	-94.391	18.010	-125.921	-58.020	0.363	0.256	2500.0	2876.0	1.0
hsgp(X, m=10, c=1)_weights[2]	106.121	20.717	69.627	147.247	0.458	0.324	2046.0	2563.0	1.0
hsgp(X, m=10, c=1)_weights[3]	149.289	22.154	106.508	188.107	0.521	0.368	1807.0	2467.0	1.0
hsgp(X, m=10, c=1)_weights[4]	-83.399	22.728	-124.812	-40.406	0.539	0.383	1776.0	2517.0	1.0
hsgp(X, m=10, c=1)_weights[5]	-125.322	23.273	-169.346	-82.187	0.531	0.377	1917.0	2940.0	1.0
hsgp(X, m=10, c=1)_weights[6]	37.835	21.153	-1.712	77.574	0.512	0.362	1713.0	2577.0	1.0
hsgp(X, m=10, c=1)_weights[7]	94.574	20.745	53.167	132.182	0.423	0.299	2393.0	3105.0	1.0
hsgp(X, m=10, c=1)_weights[8]	-1.494	17.184	-31.835	33.063	0.361	0.255	2264.0	2853.0	1.0
hsgp(X, m=10, c=1)_weights[9]	-45.442	15.955	-76.649	-16.636	0.289	0.208	3041.0	3183.0	1.0

But again we can sample from the posterior predictive distribution of the outcome variable

model_hsgp.predict(idata_hsgp, data=new_data, 
kind='pps', inplace=True)

idata_hsgp

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:                         (chain: 4, draw: 1000,
                                       hsgp(X, m=10, c=1)_weights_dim: 10,
                                       y_obs: 500)
  Coordinates:
    * chain                           (chain) int64 0 1 2 3
    * draw                            (draw) int64 0 1 2 3 4 ... 996 997 998 999
    * hsgp(X, m=10, c=1)_weights_dim  (hsgp(X, m=10, c=1)_weights_dim) int64 0 ...
    * y_obs                           (y_obs) int64 0 1 2 3 4 ... 496 497 498 499
  Data variables:
      hsgp(X, m=10, c=1)_weights_raw  (chain, draw, hsgp(X, m=10, c=1)_weights_dim) float64 ...
      y_sigma                         (chain, draw) float64 20.86 19.86 ... 21.17
      hsgp(X, m=10, c=1)_sigma        (chain, draw) float64 26.85 24.0 ... 18.06
      hsgp(X, m=10, c=1)_ell          (chain, draw) float64 3.861 3.509 ... 3.76
      hsgp(X, m=10, c=1)_weights      (chain, draw, hsgp(X, m=10, c=1)_weights_dim) float64 ...
      y_mean                          (chain, draw, y_obs) float64 15.74 ... -4...
      hsgp(X, m=10, c=1)              (chain, draw, y_obs) float64 15.74 ... -4...
  Attributes:
      created_at:                  2024-05-27T06:01:00.623094
      arviz_version:               0.17.0
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • hsgp(X, m=10, c=1)_weights_dim: 10
    • y_obs: 500
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • hsgp(X, m=10, c=1)_weights_dim
      (hsgp(X, m=10, c=1)_weights_dim)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (7)
    • hsgp(X, m=10, c=1)_weights_raw
      (chain, draw, hsgp(X, m=10, c=1)_weights_dim)
      float64
      -1.452 -0.974 ... 0.9659 -0.2453

      array([[[-1.45186697, -0.97401433,  1.51626692, ...,  1.84709969,
                0.92584892, -1.19616826],
              [-1.73799921, -1.55920285,  1.42827017, ...,  1.78787921,
               -1.14159607, -1.66042935],
              [-1.3059213 , -1.22687096,  2.0214052 , ...,  2.80079646,
                0.17462599, -0.06073935],
              ...,
              [-2.26781832, -1.21871881,  1.80249133, ...,  2.65551316,
                0.30040289, -0.7085128 ],
              [-2.23372877, -1.1832454 ,  1.74365597, ...,  2.71376469,
                0.41123845, -0.60671544],
              [-1.96117304, -1.13093178,  2.24600086, ...,  2.3750183 ,
               -0.00901474, -0.75147244]],
      
             [[-1.95999601, -1.27507311,  1.52467611, ...,  3.28813978,
               -0.1598232 , -1.10988783],
              [-1.59314189, -1.44864255,  1.95956808, ...,  2.01431167,
                1.0199955 , -0.8962025 ],
              [-1.8432561 , -1.32415496,  1.80874966, ...,  3.05464118,
                0.81105708, -1.02872434],
      ...
              [-1.94277245, -1.5433581 ,  1.78771407, ...,  2.11464836,
               -0.62470516, -1.01229181],
              [-1.96747626, -1.39055532,  1.24775496, ...,  2.93940091,
                0.46989555, -1.57735663],
              [-2.20654554, -1.042289  ,  1.46402154, ...,  3.37106106,
               -0.11061588, -2.15605928]],
      
             [[-1.68378706, -0.84035677,  1.41475357, ...,  2.46320953,
                0.88793346, -0.56144583],
              [-1.66959934, -1.01014124,  1.6874097 , ...,  2.71030445,
                0.68597924, -0.76198643],
              [-1.54563186, -0.89279988,  1.77115866, ...,  2.3842957 ,
                0.08392828, -0.78058783],
              ...,
              [-2.17148384, -1.59167361,  1.34960308, ...,  2.34802173,
                0.44912308, -1.0698665 ],
              [-1.31376138, -0.63355604,  1.97759132, ...,  1.8841587 ,
               -0.13453184, -1.36605679],
              [-2.0018976 , -1.36470426,  2.15352253, ...,  3.17100177,
                0.96590611, -0.24527631]]])

    • y_sigma
      (chain, draw)
      float64
      20.86 19.86 21.17 ... 20.12 21.17

      array([[20.85613959, 19.86298062, 21.16886643, ..., 19.48940449,
              19.31511941, 19.44304302],
             [18.5908299 , 21.71605063, 21.22255223, ..., 19.44547467,
              19.21042952, 22.0099345 ],
             [20.48425625, 21.17420065, 22.11318822, ..., 18.40825711,
              21.65013089, 22.05150211],
             [21.08143229, 20.7215663 , 22.60954219, ..., 21.15149371,
              20.11622356, 21.16568101]])

    • hsgp(X, m=10, c=1)_sigma
      (chain, draw)
      float64
      26.85 24.0 27.81 ... 26.09 18.06

      array([[26.84572148, 24.00472166, 27.80754218, ..., 23.34108117,
              23.72941745, 22.60241785],
             [22.01917097, 24.17709343, 26.16315265, ..., 22.03189338,
              22.17588023, 24.00832607],
             [24.12239909, 24.36393141, 25.79692715, ..., 24.33866729,
              20.63877848, 22.26540398],
             [25.41092492, 25.05093181, 28.3790263 , ..., 22.30157727,
              26.0901197 , 18.06420147]])

    • hsgp(X, m=10, c=1)_ell
      (chain, draw)
      float64
      3.861 3.509 3.695 ... 3.679 3.76

      array([[3.86057038, 3.50871714, 3.69475981, ..., 3.20961859, 3.07237787,
              2.99194079],
             [4.21093385, 4.11858175, 3.55744155, ..., 3.54737759, 2.43665447,
              3.67009957],
             [3.38790729, 3.56906294, 3.07026931, ..., 3.17401967, 4.77493695,
              4.1757221 ],
             [3.73124202, 3.62463311, 4.06803959, ..., 3.06030257, 3.67928011,
              3.75995625]])

    • hsgp(X, m=10, c=1)_weights
      (chain, draw, hsgp(X, m=10, c=1)_weights_dim)
      float64
      -120.2 -78.54 ... 27.36 -5.933

      array([[[-120.19159663,  -78.54474016,  117.04015633, ...,
                 88.13161831,   38.07210878,  -41.65695052],
              [-122.83658383, -107.83699596,   95.27680851, ...,
                 80.15736607,  -45.26655282,  -57.3944273 ],
              [-109.63228325, -100.55013819,  159.16260856, ...,
                141.94233647,    7.72311493,   -2.30700494],
              ...,
              [-149.23661164,  -78.75788623,  113.01499747, ...,
                119.40081633,   12.18790241,  -25.62639515],
              [-146.28304701,  -76.21158457,  109.23906639, ...,
                125.36332205,   17.28993874,  -22.96014529],
              [-120.75657952,  -68.54675039,  132.60289356, ...,
                105.03375274,   -0.36461133,  -27.50691292]],
      
             [[-138.76216692,  -87.49701007,   99.32062774, ...,
                120.85182044,   -4.92170357,  -28.04716014],
              [-122.53330814, -108.14123925,  139.17924148, ...,
                 82.74939574,   35.37871319,  -25.72818479],
              [-142.94342858, -100.42482592,  132.17632445, ...,
                148.36712045,   34.72137014,  -38.24384896],
      ...
              [-132.58698442, -103.47676242,  116.36868708, ...,
                 99.43921353,  -26.56716678,  -38.47612061],
              [-138.61535955,  -94.11467693,   78.98484564, ...,
                 89.13997713,   11.35087487,  -29.54951483],
              [-157.32953055,  -72.06983662,   96.1818003 , ...,
                126.1515611 ,   -3.47852942,  -55.82170316]],
      
             [[-129.78695588,  -63.20664967,  102.15003445, ...,
                113.4760498 ,   35.60103971,  -19.27405014],
              [-125.10213961,  -73.9587159 ,  118.87369923, ...,
                124.94562328,   27.73901676,  -26.61384578],
              [-138.71497899,  -77.82483202,  147.07295446, ...,
                116.05385128,    3.46342745,  -26.78443203],
              ...,
              [-133.39280566,  -96.17651826,   79.33895456, ...,
                102.02385007,   17.77408526,  -38.1416893 ],
              [-103.26872735,  -48.62806342,  145.87733997, ...,
                 89.78661037,   -5.60101795,  -48.90554932],
              [-110.1016679 ,  -73.21183754,  110.83494684, ...,
                103.40789151,   27.35517734,   -5.93337047]]])

    • y_mean
      (chain, draw, y_obs)
      float64
      15.74 15.74 15.71 ... -4.651 -4.606

      array([[[ 15.73717823,  15.74095015,  15.71442577, ...,  -4.67876513,
                -4.20969038,  -3.71327448],
              [-22.08564087, -21.76667383, -21.39772149, ...,   7.03009177,
                 7.12845194,   7.21535709],
              [ 56.02093901,  56.10232981,  56.08787761, ...,   1.40193661,
                 1.10510341,   0.79052751],
              ...,
              [ 27.64916801,  27.65115973,  27.60142293, ...,   1.4750735 ,
                 1.59074503,   1.70839567],
              [ 29.49505409,  29.44086994,  29.3291777 , ...,   3.89804528,
                 4.04221731,   4.18596665],
              [ 24.70348988,  24.9191706 ,  25.10055256, ...,  -1.75167732,
                -1.70197532,  -1.64814065]],
      
             [[ 13.09479226,  13.08797832,  13.05446772, ...,  11.0508337 ,
                11.0972292 ,  11.12522822],
              [ 19.45360075,  19.35303968,  19.21116795, ..., -17.11306852,
               -16.95029643, -16.75310125],
              [ 39.40925935,  39.40667898,  39.33002811, ...,  -2.12983751,
                -2.01625076,  -1.89843614],
      ...
              [-14.4322126 , -14.38059677, -14.30401941, ...,  11.55147703,
                11.55940569,  11.5447161 ],
              [  3.99064388,   3.88036701,   3.75439771, ...,   6.42065763,
                 6.71826106,   7.01812801],
              [ 11.85726202,  12.01889355,  12.16180986, ...,  12.90196389,
                13.3121751 ,  13.717131  ]],
      
             [[ 27.47986715,  27.31475262,  27.09142621, ...,   7.10292311,
                 7.47525977,   7.85365034],
              [ 34.24033269,  34.24902789,  34.19561616, ...,   5.51087296,
                 5.77912687,   6.04952041],
              [ 31.73947173,  31.94613931,  32.10509971, ...,   9.39499668,
                 9.48056436,   9.55191877],
              ...,
              [  9.46004312,   9.35384251,   9.22108586, ...,   3.61236621,
                 3.94961437,   4.29454081],
              [ 17.75725777,  18.27924573,  18.79439922, ...,   5.95699839,
                 6.24094799,   6.5246295 ],
              [ 33.41206521,  33.28138834,  33.08465457, ...,  -4.68577545,
                -4.65087068,  -4.60586526]]])

    • hsgp(X, m=10, c=1)
      (chain, draw, y_obs)
      float64
      15.74 15.74 15.71 ... -4.651 -4.606

      array([[[ 15.73717823,  15.74095015,  15.71442577, ...,  -4.67876513,
                -4.20969038,  -3.71327448],
              [-22.08564087, -21.76667383, -21.39772149, ...,   7.03009177,
                 7.12845194,   7.21535709],
              [ 56.02093901,  56.10232981,  56.08787761, ...,   1.40193661,
                 1.10510341,   0.79052751],
              ...,
              [ 27.64916801,  27.65115973,  27.60142293, ...,   1.4750735 ,
                 1.59074503,   1.70839567],
              [ 29.49505409,  29.44086994,  29.3291777 , ...,   3.89804528,
                 4.04221731,   4.18596665],
              [ 24.70348988,  24.9191706 ,  25.10055256, ...,  -1.75167732,
                -1.70197532,  -1.64814065]],
      
             [[ 13.09479226,  13.08797832,  13.05446772, ...,  11.0508337 ,
                11.0972292 ,  11.12522822],
              [ 19.45360075,  19.35303968,  19.21116795, ..., -17.11306852,
               -16.95029643, -16.75310125],
              [ 39.40925935,  39.40667898,  39.33002811, ...,  -2.12983751,
                -2.01625076,  -1.89843614],
      ...
              [-14.4322126 , -14.38059677, -14.30401941, ...,  11.55147703,
                11.55940569,  11.5447161 ],
              [  3.99064388,   3.88036701,   3.75439771, ...,   6.42065763,
                 6.71826106,   7.01812801],
              [ 11.85726202,  12.01889355,  12.16180986, ...,  12.90196389,
                13.3121751 ,  13.717131  ]],
      
             [[ 27.47986715,  27.31475262,  27.09142621, ...,   7.10292311,
                 7.47525977,   7.85365034],
              [ 34.24033269,  34.24902789,  34.19561616, ...,   5.51087296,
                 5.77912687,   6.04952041],
              [ 31.73947173,  31.94613931,  32.10509971, ...,   9.39499668,
                 9.48056436,   9.55191877],
              ...,
              [  9.46004312,   9.35384251,   9.22108586, ...,   3.61236621,
                 3.94961437,   4.29454081],
              [ 17.75725777,  18.27924573,  18.79439922, ...,   5.95699839,
                 6.24094799,   6.5246295 ],
              [ 33.41206521,  33.28138834,  33.08465457, ...,  -4.68577545,
                -4.65087068,  -4.60586526]]])

  • Indexes: (4)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • hsgp(X, m=10, c=1)_weights_dim
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='hsgp(X, m=10, c=1)_weights_dim'))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (4)

    created_at :

    2024-05-27T06:01:00.623094

    arviz_version :

    0.17.0

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 500)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 492 493 494 495 496 497 498 499
  Data variables:
      y        (chain, draw, y_obs) float64 76.04 21.77 64.12 ... -41.12 -26.02
  Attributes:
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 500
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      76.04 21.77 64.12 ... -41.12 -26.02

      array([[[ 7.60446220e+01,  2.17746429e+01,  6.41213079e+01, ...,
                2.88184766e+01,  5.35607407e+00,  8.28161473e+00],
              [-3.71916430e+01, -1.02105981e+01, -3.75556355e+01, ...,
                6.13483186e+00,  5.17152558e+00, -7.14534907e+00],
              [ 4.18037444e+01,  8.94677926e+01,  6.01499750e+01, ...,
                5.15362586e+00,  1.50220973e+01,  5.10623210e+00],
              ...,
              [ 1.17455216e+01,  7.18088614e+00,  2.27295340e+01, ...,
               -1.40562367e+01, -1.55645500e+01, -1.59305317e+01],
              [ 4.57996940e+01,  4.09450669e+01,  6.17203837e+01, ...,
                5.43347241e+00,  2.16723522e+01, -1.34623303e+01],
              [-1.99298911e+01, -7.31053097e+00,  2.33680338e+01, ...,
               -1.86693562e+01, -3.45100068e+01, -4.37085484e-02]],
      
             [[ 1.68744925e+01,  4.86269717e+01, -3.54535158e-01, ...,
                6.97092074e+00,  1.38401427e+01,  4.41897478e+01],
              [ 3.80673611e+01,  5.08088452e+01,  3.24135895e+01, ...,
               -1.64036397e+00,  9.20252673e-01, -2.28191503e+01],
              [ 2.45419863e+01,  3.67156591e+01,  8.06594166e+01, ...,
                2.35165419e+01, -2.27769611e+00,  3.19786177e+01],
      ...
              [ 7.65492507e+00, -2.52651612e+01, -1.14319746e+01, ...,
                2.10783876e+00,  6.51461474e+01,  1.39449220e+01],
              [-1.52884615e+01, -3.81858568e+01,  3.83976350e+01, ...,
                6.70649060e+00, -4.19015833e+00,  3.10102382e+01],
              [ 3.22122816e+01,  1.40541501e-01,  1.39868866e+01, ...,
                1.39995757e+00,  1.40413480e+01,  2.15327900e+01]],
      
             [[ 6.45006934e+01,  1.89982303e+01,  4.92943955e+01, ...,
                4.48989355e+01,  6.00626370e+00, -4.11335459e+00],
              [ 7.69567142e+01,  3.85217722e+01,  6.00470265e+01, ...,
               -2.85540970e+01,  1.17829549e+00,  2.88735960e+01],
              [ 3.82736929e+01,  1.35754696e+01,  3.18636391e+01, ...,
                1.30573340e+01,  2.58695675e+01,  3.63775721e+00],
              ...,
              [ 2.64333429e+01,  2.96956055e+01,  9.65696351e+00, ...,
               -8.66072179e+00, -1.22879325e+00, -1.74010423e+01],
              [-3.67585411e+01, -3.45711874e+00,  1.22087765e+01, ...,
                2.78805963e+01, -2.67630198e+01, -1.17549578e+01],
              [ 6.74946544e+01,  2.73016469e+01,  1.89342532e+01, ...,
               -1.80490133e+01, -4.11156817e+01, -2.60158793e+01]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (2)

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 133)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (chain, draw, y_obs) float64 -3.957 -3.957 -3.957 ... -4.013 -4.099
  Attributes:
      created_at:                  2024-05-27T06:01:00.626582
      arviz_version:               0.17.0
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 133
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      -3.957 -3.957 ... -4.013 -4.099

      array([[[-3.95662204, -3.95668036, -3.95697947, ..., -3.96786393,
               -4.06567925, -4.08819109],
              [-3.92804385, -3.95248402, -4.00438808, ..., -4.19279911,
               -3.90894979, -4.05289002],
              [-4.11236339, -4.08418887, -3.99366348, ..., -4.1948162 ,
               -3.97180293, -4.09921461],
              ...,
              [-3.88915792, -3.88893229, -3.89724674, ..., -4.04248536,
               -3.89847518, -4.03951892],
              [-3.88161419, -3.88254965, -3.9052254 , ..., -4.08668059,
               -3.88211318, -4.03326819],
              [-4.08257196, -4.09265666, -4.04648571, ..., -3.99108035,
               -3.91131912, -4.03785686]],
      
             [[-3.85025382, -3.84763163, -3.85608195, ..., -4.34165517,
               -3.86684608, -4.00723736],
              [-4.0090095 , -4.00707108, -4.017144  , ..., -4.01175855,
               -4.25268066, -4.11837821],
              [-3.97404168, -3.97456237, -4.00764182, ..., -4.0433394 ,
               -4.00511976, -4.10110199],
      ...
              [-3.85025004, -3.83654156, -3.832685  , ..., -4.37637177,
               -3.86484545, -4.00066997],
              [-4.06709934, -4.05434047, -4.05493257, ..., -4.17570336,
               -3.99476606, -4.11607839],
              [-4.01346366, -4.01406212, -4.01257167, ..., -4.30726434,
               -4.01839009, -4.13004222]],
      
             [[-3.98552821, -3.98561292, -4.0079376 , ..., -4.15815157,
               -3.96833515, -4.09613757],
              [-3.96540632, -3.96291088, -3.95021137, ..., -4.13476712,
               -3.95195981, -4.08343265],
              [-4.17463712, -4.17201421, -4.11695573, ..., -4.3197009 ,
               -4.04463544, -4.14929407],
              ...,
              [-4.0692212 , -4.06241625, -4.09034963, ..., -4.09469991,
               -3.98159516, -4.09860353],
              [-4.68482086, -4.73793167, -4.72387693, ..., -4.11913623,
               -3.92229916, -4.06192875],
              [-3.97170872, -3.97200671, -3.98370832, ..., -4.02417539,
               -4.01254991, -4.09910258]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (4)

    created_at :

    2024-05-27T06:01:00.626582

    arviz_version :

    0.17.0

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:          (chain: 4, draw: 1000)
  Coordinates:
    * chain            (chain) int64 0 1 2 3
    * draw             (draw) int64 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
  Data variables:
      acceptance_rate  (chain, draw) float64 0.9484 0.989 0.9772 ... 0.9905 0.9283
      step_size        (chain, draw) float64 0.1631 0.1631 ... 0.1827 0.1827
      diverging        (chain, draw) bool False False False ... False False False
      energy           (chain, draw) float64 671.4 665.4 671.0 ... 670.8 674.0
      n_steps          (chain, draw) int64 31 31 31 15 15 31 ... 15 31 31 15 15 31
      tree_depth       (chain, draw) int64 5 5 5 4 4 5 5 4 4 ... 4 5 5 4 5 5 4 4 5
      lp               (chain, draw) float64 661.7 661.3 665.7 ... 665.6 668.8
  Attributes:
      created_at:                  2024-05-27T06:01:00.625568
      arviz_version:               0.17.0
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (7)
    • acceptance_rate
      (chain, draw)
      float64
      0.9484 0.989 ... 0.9905 0.9283

      array([[0.94837217, 0.98900712, 0.97718268, ..., 0.99599719, 0.97974622,
              0.93727528],
             [0.9824305 , 0.87426419, 0.99008681, ..., 0.99942845, 0.97008375,
              0.99555447],
             [0.99910745, 0.98130743, 0.95550076, ..., 0.97069333, 0.9342073 ,
              0.95033605],
             [0.97363799, 0.97895034, 0.93808177, ..., 0.9514047 , 0.99054107,
              0.92830969]])

    • step_size
      (chain, draw)
      float64
      0.1631 0.1631 ... 0.1827 0.1827

      array([[0.16306409, 0.16306409, 0.16306409, ..., 0.16306409, 0.16306409,
              0.16306409],
             [0.15671706, 0.15671706, 0.15671706, ..., 0.15671706, 0.15671706,
              0.15671706],
             [0.15280013, 0.15280013, 0.15280013, ..., 0.15280013, 0.15280013,
              0.15280013],
             [0.1826574 , 0.1826574 , 0.1826574 , ..., 0.1826574 , 0.1826574 ,
              0.1826574 ]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • energy
      (chain, draw)
      float64
      671.4 665.4 671.0 ... 670.8 674.0

      array([[671.39512608, 665.35858794, 670.96295141, ..., 670.72923684,
              668.42111758, 669.24601286],
             [665.04518693, 671.40450254, 672.47055359, ..., 665.8504822 ,
              674.75907256, 672.99992931],
             [672.00375668, 668.16520377, 666.3166769 , ..., 671.79404347,
              671.22701287, 671.74308362],
             [669.46327702, 671.11966704, 667.28825365, ..., 665.11423361,
              670.79366043, 673.9890605 ]])

    • n_steps
      (chain, draw)
      int64
      31 31 31 15 15 ... 31 31 15 15 31

      array([[31, 31, 31, ..., 31, 15, 15],
             [15, 31, 31, ..., 31, 31, 31],
             [31, 31, 15, ..., 15, 31, 15],
             [31, 15, 15, ..., 15, 15, 31]])

    • tree_depth
      (chain, draw)
      int64
      5 5 5 4 4 5 5 4 ... 5 5 4 5 5 4 4 5

      array([[5, 5, 5, ..., 5, 4, 4],
             [4, 5, 5, ..., 5, 5, 5],
             [5, 5, 4, ..., 4, 5, 4],
             [5, 4, 4, ..., 4, 4, 5]])

    • lp
      (chain, draw)
      float64
      661.7 661.3 665.7 ... 665.6 668.8

      array([[661.70967519, 661.2875458 , 665.70729556, ..., 661.40453042,
              661.53660449, 661.51768248],
             [662.47186408, 665.25230302, 666.28391286, ..., 662.10849958,
              668.1822538 , 664.18500466],
             [662.60487008, 660.69369885, 662.99303071, ..., 663.00720061,
              664.60205279, 667.39059102],
             [664.1010232 , 660.27378818, 664.08879552, ..., 661.26338997,
              665.58839085, 668.77321037]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (4)

    created_at :

    2024-05-27T06:01:00.625568

    arviz_version :

    0.17.0

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (y_obs: 133)
  Coordinates:
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (y_obs) float64 0.0 -1.3 -2.7 0.0 -2.7 ... -2.7 10.7 -2.7 10.7
  Attributes:
      created_at:                  2024-05-27T06:01:00.626851
      arviz_version:               0.17.0
      inference_library:           numpyro
      inference_library_version:   0.13.2
      sampling_time:               1.870798
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • y_obs: 133
  • Coordinates: (1)
    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (y_obs)
      float64
      0.0 -1.3 -2.7 ... 10.7 -2.7 10.7

      array([   0. ,   -1.3,   -2.7,    0. ,   -2.7,   -2.7,   -2.7,   -1.3,
               -2.7,   -2.7,   -1.3,   -2.7,   -2.7,   -2.7,   -5.4,   -2.7,
               -5.4,    0. ,   -2.7,   -2.7,    0. ,  -13.3,   -5.4,   -5.4,
               -9.3,  -16. ,  -22.8,   -2.7,  -22.8,  -32.1,  -53.5,  -54.9,
              -40.2,  -21.5,  -21.5,  -50.8,  -42.9,  -26.8,  -21.5,  -50.8,
              -61.7,   -5.4,  -80.4,  -59. ,  -71. ,  -91.1,  -77.7,  -37.5,
              -85.6, -123.1, -101.9,  -99.1, -104.4, -112.5,  -50.8, -123.1,
              -85.6,  -72.3, -127.2, -123.1, -117.9, -134. , -101.9, -108.4,
             -123.1, -123.1, -128.5, -112.5,  -95.1,  -81.8,  -53.5,  -64.4,
              -57.6,  -72.3,  -44.3,  -26.8,   -5.4, -107.1,  -21.5,  -65.6,
              -16. ,  -45.6,  -24.2,    9.5,    4. ,   12. ,  -21.5,   37.5,
               46.9,  -17.4,   36.2,   75. ,    8.1,   54.9,   48.2,   46.9,
               16. ,   45.6,    1.3,   75. ,  -16. ,  -54.9,   69.6,   34.8,
               32.1,  -37.5,   22.8,   46.9,   10.7,    5.4,   -1.3,  -21.5,
              -13.3,   30.8,  -10.7,   29.4,    0. ,  -10.7,   14.7,   -1.3,
                0. ,   10.7,   10.7,  -26.8,  -14.7,  -13.3,    0. ,   10.7,
              -14.7,   -2.7,   10.7,   -2.7,   10.7])

  • Indexes: (1)
    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (7)

    created_at :

    2024-05-27T06:01:00.626851

    arviz_version :

    0.17.0

    inference_library :

    numpyro

    inference_library_version :

    0.13.2

    sampling_time :

    1.870798

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  

and plot the model fit to see if it can recover the observed data.

ax = az.plot_hdi(new_data['X'], idata_hsgp['posterior_predictive']['y'], fill_kwargs={'alpha': 0.2, 'color':'firebrick'}, figsize=(9, 8))

az.plot_hdi(new_data['X'], idata_hsgp['posterior_predictive']['y'], fill_kwargs={'alpha': 0.8, 'color':'firebrick'}, hdi_prob=0.5)

y_mean = idata_hsgp['posterior_predictive']['y'].mean(dim=('chain', 'draw'))

ax.plot(new_data['X'], y_mean, label='Expected posterior predictive', color='k')

ax.scatter(df['X'], df['y'], label='Observed Datapoints')

ax.legend()

ax.set_title("Posterior Predictive Distribution \n Based on HSGP approximation");

Here again we’ve allowed the model to extrapolate beyond the observed data but we achieve far more reasonable extrapolation than with the Spline model. And we can compare other performance metrics versus the spline models to see that by the aggregate performance measures our HSGP model seems to come out on top too.

models_dict = {"piecewise_constant": idata_spline1, "piecewise_linear": idata_spline2, "cubic_bspline": idata_spline3, "cubic_bspline_10": idata_spline4, 
"cubic_bspline_15": idata_spline5, 'hsgp': idata_hsgp}
df_compare = az.compare(models_dict)
df_compare

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
hsgp	0	-608.819807	9.958864	0.000000	0.872617	12.009745	0.000000	False	log
cubic_bspline_10	1	-612.578311	11.652463	3.758504	0.000000	9.646544	3.125929	True	log
cubic_bspline_15	2	-620.709337	19.669201	11.889531	0.000000	9.600116	3.836813	True	log
cubic_bspline	3	-634.647180	8.703519	25.827373	0.000000	8.915728	7.666315	True	log
piecewise_constant	4	-643.781042	6.981098	34.961235	0.069309	9.770740	9.490854	False	log
piecewise_linear	5	-647.016885	5.987587	38.197078	0.058073	7.914787	9.204140	False	log

For a deeper dive on HSGP you might consult Juan Orduz’s work here

Bayesian Penalised Splines (Smoothing Priors)

Priors can be used in the same way to that regularisation techniques are deployed in machine learning. They serve as a kind of extra model parameter used to effect a variable selection routine. In the case of splines we are creating these synthetic features and this can be seen as a generating a feature selection problem. So far we have discussed priors as a means of imposing structure based on our knowledge of the problem but we can also use priors to automate feature selection over spline features.

Here we fit a two spline models to the speed-test data set but we use 60 knot points to generate our spline features. This is wildy excessive for the data set in question. We show how to use non-centered parameterisation with a “smoothing” trick to induce a good fit despite excessive features. In this way we’re removing the burden of manually specifying the spline structure. This is akin to penalisation methods used in PyGam but with a Bayesian flavour.

max_dev= len(y)
num_knots=60
knot_list = np.linspace(0, np.max(X), num_knots)[2:-2]
dev_periods = np.linspace(0, np.max(X), len(y))

Bx = dmatrix(
    "bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=False)",
    {"dev_periods": dev_periods, "knots": knot_list},
)
Bx

def penalised_model(Bx, penalised=True):
    with pm.Model() as model:
        basis = pm.MutableData('basis', np.asfortranarray(Bx))

        if penalised:
            sigma_a = pm.HalfCauchy('sigma_a', 5.)
            a0 = pm.Normal('a0', 0., 10.)
            delta_a = pm.Normal('delta_a', 0., 1.,shape=num_knots)
            a = pm.Deterministic('a', a0 + (sigma_a * delta_a).cumsum())
        else: 
            b = pm.Normal('delta_a', 0., 25,
            shape=num_knots)
            a0 = pm.Normal('a0', 0, 10)
            a = pm.Deterministic('a', a0 + b)
        
        sigma = pm.HalfCauchy('sigma', 5.)
        obs = pm.Normal('obs', pm.math.dot(basis, a), sigma, observed=y)
        idata = pm.sample(target_accept=.95, idata_kwargs={"log_likelihood": True})
    
    return model, idata

model_raw, idata_raw = penalised_model(Bx, penalised=False)
model_penalised, idata_penalised = penalised_model(Bx, penalised=True)

The model structure is as follows:

pm.model_to_graphviz(model_penalised)

This parameterisation trick is owed to work by Austin Rochford and Adrian Seyboldt discussed here. We can see how the two parameterisations induce quite distinct values on the coefficients for splines.

fig, ax = plt.subplots(figsize=(9, 20))
az.plot_forest([idata_raw, idata_penalised], var_names=['a'], combined=True, ax=ax, model_names=['raw', 'penalised'])
ax.fill_betweenx(range(160), -25, 25, alpha=0.2, color='red');

This results in much different posterior predictive fits. The un-smoothed spline fit overfits to the small partitions of the spline features, whereas the smoothing spline model achieves a reasonable fit recovering the smooth structure of the data generating process.

B_df = pd.DataFrame(Bx)
weights = az.summary(idata_penalised, var_names='a')['mean']
fig, ax = plt.subplots(figsize=(8, 7))
processed = np.array([np.dot(B_df, az.extract(idata_penalised['posterior']['a'])['a'].to_numpy()[:, i]) for i in range(4000)])
az.plot_hdi(range(len(y)), processed, ax=ax,  fill_kwargs={'alpha': 0.4, 'color':'firebrick'})
az.plot_hdi(range(len(y)), processed, ax=ax,  fill_kwargs={'alpha': 0.8, 'color':'firebrick'}, hdi_prob=.50)
weights_raw = az.summary(idata_raw, var_names='a')['mean']
ax.plot(range(len(y)), np.dot(B_df, weights), color='black', label='Cubic_Bsplines_60_Penalised')
ax.plot(range(len(y)), np.dot(B_df, weights_raw), color='blue', linestyle='--', label='Cubic_Bsplines_60_unpenalised')
ax.set_title("Bayesian Automatic Regularisation \n Smoothing Priors over 60 Knots")
ax.legend();

This is reflected in the model performance measures as well.

models_dict = {"piecewise_constant": idata_spline1, "piecewise_linear": idata_spline2, "cubic_bspline": idata_spline3, "cubic_bspline_10": idata_spline4, 
"cubic_bspline_20": idata_spline5, 'hsgp': idata_hsgp, 'spline_penalised_bayes': idata_penalised, 'cubic_bspline_60': idata_raw}
df_compare = az.compare(models_dict)
df_compare

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
hsgp	0	-608.819807	9.958864	0.000000	0.764385	12.009745	0.000000	False	log
spline_penalised_bayes	1	-610.258217	19.573938	1.438411	0.115567	11.444307	2.822797	True	log
cubic_bspline_10	2	-612.578311	11.652463	3.758504	0.000000	9.646544	3.125929	True	log
cubic_bspline_20	3	-620.709337	19.669201	11.889531	0.000000	9.600116	3.836813	True	log
cubic_bspline	4	-634.647180	8.703519	25.827373	0.000000	8.915728	7.666315	True	log
cubic_bspline_60	5	-642.070674	24.951194	33.250867	0.000667	9.135971	8.313174	False	log
piecewise_constant	6	-643.781042	6.981098	34.961235	0.061876	9.770740	9.490854	False	log
piecewise_linear	7	-647.016885	5.987587	38.197078	0.057506	7.914787	9.204140	False	log

This functionality of Bayesian priors is not unique to spline models. We can mirror Lasso regression models with double-exponential prior specifications. Ridge regrssion is provably a species of regression models with Normal priors pulled toward 0. Other alternatives such horse-shoe and slab and spike specifications offer a rich range of bayesian approaches to variable selection methods with sparse data.

Recap

So far we’ve seen how we can use splines and gauassian processes to model highly eccentric functional relationships where the function could be approximated with univariate smoothing routine. These are two distinct abstractions which seem adequately fit to the world, but demand very different ways of thinking about the underlying reality. This is no fundamental contradiction in so far as the world admits many descriptions.

[K]nowledge develops in a multiplicity of theories, each with it’s limited utility … These theories overlap very considerable in their logical laws and in much else, but that they add up to an integrated and consistent whole is only a worthy ideal and not a pre-requistite of scientific progress … Let reconciliations proceed. - W.V.O Quine in Word and Object

Another observation in a similar vein is that penalised spline models are provably equivalent to hierarchical regression (random effects) models. This is striking because the character of these types of model seems diametrically opposed. With spline models you jerry-rig your features to an optimisation goal, with hierarchical model you tend to impose theoretical structure to induce shrinkage. It’s hard to see how this would work? With a penalised spline model are you inducing a hierarchy of latent features you can’t name? Should you even try to translate between the two!? How can we leverage our prior knowledge of the data generating process to add structure to spline model predictions?

The abstract components of our model are less graspable than their predictive performance yet the qualitative character of theory buiding differs markedly between these abstractions and aesthetics and explainability feature heavily in model preference. Next we’ll add another layer to our abstractions and show how to use hierarchical modelling over spline fits can be used to extract insight into the data generating process over a family of curves. In particular we’ll focus on the development of insurance loss curves and how to constrain our predictions within reasonable bounds. This is an apt problem for spline modelling in that the process is repeatable annually, but the domain of the support i.e. development in time is fixed year-on-year so the extrapolation problem is muted.

Insurance Loss Curves: Hierarchical Spline Models

Hierarchical models of spline fits allow us to characterise and sample from fits over a family of curves while improving out of sample predictive accuracy. We draw on car insurance losses data set discussed in Mick Cooney’s Stan case-study, but we simplify things for ourselves considerably by focusing on two types of loss and ensuring that each year under consideration has equal observations of the accruing losses.

loss_df = pd.read_csv('ppauto_pos.csv')
loss_df = loss_df[(loss_df['GRCODE'].isin([43, 353])) & (loss_df['DevelopmentYear'] < 1998)]

loss_df = loss_df[['GRCODE', 'AccidentYear', 'DevelopmentYear', 'DevelopmentLag', 'EarnedPremDIR_B', 'CumPaidLoss_B']]


loss_df.columns = ['grcode', 'acc_year', 'dev_year', 'dev_lag', 'premium', 'cum_loss']
loss_df['lr'] = loss_df['cum_loss'] / loss_df['premium']
loss_df = loss_df[(loss_df['acc_year'] <= 1992) & (loss_df['dev_lag'] <= 6)].reset_index(drop=True)

loss_df['year_code'] = loss_df['acc_year'].astype(str) + '_' + loss_df['grcode'].astype(str)
loss_df.sort_values(by=['year_code', 'acc_year', 'dev_lag'], inplace=True)
loss_df['standardised_premium'] = (loss_df.assign(mean_premium = np.mean(loss_df['premium']))
.assign(std_premium = np.std(loss_df['premium']))
.apply(lambda x: (x['mean_premium'] - x['premium']) /x['std_premium'], axis=1)
)

loss_df.head(12)

	grcode	acc_year	dev_year	dev_lag	premium	cum_loss	lr	year_code	standardised_premium
30	353	1988	1988	1	18793	4339	0.230884	1988_353	-0.347453
31	353	1988	1989	2	18793	9617	0.511733	1988_353	-0.347453
32	353	1988	1990	3	18793	11584	0.616400	1988_353	-0.347453
33	353	1988	1991	4	18793	12001	0.638589	1988_353	-0.347453
34	353	1988	1992	5	18793	12640	0.672591	1988_353	-0.347453
35	353	1988	1993	6	18793	12966	0.689938	1988_353	-0.347453
0	43	1988	1988	1	957	133	0.138976	1988_43	1.722344
1	43	1988	1989	2	957	333	0.347962	1988_43	1.722344
2	43	1988	1990	3	957	431	0.450366	1988_43	1.722344
3	43	1988	1991	4	957	570	0.595611	1988_43	1.722344
4	43	1988	1992	5	957	615	0.642633	1988_43	1.722344
5	43	1988	1993	6	957	615	0.642633	1988_43	1.722344

Plot the Loss Curves

Here we have plotted the developing loss curves from two different coded insurance products.

pivot = loss_df.pivot(index=['dev_lag'], columns=['grcode', 'acc_year'], values='lr')
fig, axs = plt.subplots(1, 2, figsize=(9, 7))
pivot.plot(figsize=(10, 6), ax=axs[0])
axs[0].set_title("Loss Ratios by Year");
for c in pivot.columns:
    if 43 in c:
        color='red'
    else: 
        color='grey'
    axs[1].plot(pivot[c], color=color, label=c)
axs[1].legend()
axs[1].set_title("Loss Ratio by Group");

We want to model these curves collectively as instances of draws from a distribution of loss curves. To do so we will specify a PyMC hierarchical (mixed) spline model. To do so we will have a spline basis for the global hyper parameters beta_g and the individual parameters for each curve. Here we define a convenience function to generate the basis splines.

def make_basis_splines(num_knots=3, max_dev=7):
    knot_list = np.linspace(0, max_dev, num_knots+2)[1:-1]
    dev_periods = np.arange(1, max_dev, 1)

    Bi = dmatrix(
        "bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True) - 1",
        {"dev_periods": dev_periods, "knots": knot_list},
    )

    Bg = dmatrix(
        "bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True) - 1",
        {"dev_periods": dev_periods, "knots": knot_list})


    return Bi, Bg

Bi, Bg = make_basis_splines()
Bg

DesignMatrix with shape (6, 7)
  Columns:
    ['bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[0]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[1]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[2]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[3]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[4]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[5]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[6]']
  Terms:
    'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)' (columns 0:7)
  (to view full data, use np.asarray(this_obj))

Next we specify a model maker function to create the various pooled, unpooled and hierarhical (mixed) models of the insurance curve data. Note that even though we’re specifying a hierarhical model we have not specified a hierarchy over the insurance codes, instead we have added this as a “fixed” effect feature into our regression model. The idea is that this fixed effect will capture the differences in expected baseline per group code of insurance product.

def make_model(loss_df, num_knots=3, max_dev=7, model_type='mixed'):
    Bi, Bg = make_basis_splines(num_knots, max_dev)
    observed = loss_df['lr'].values
    uniques, unique_codes = pd.factorize(loss_df['year_code'])
    coords= {'years': unique_codes, 
            'splines': list(range(Bi.shape[1])) ,
            'measurement': list(range(6)), 
            'obs': uniques
            }

    with pm.Model(coords=coords) as sp_insur:
        basis_g = pm.MutableData('Bg', np.asfortranarray(Bg))

        tau = pm.HalfCauchy('tau', 1)
        ## Global Hierarchical Spline Terms
        beta_g = pm.Normal("beta_g", mu=0, sigma=tau, 
        dims='splines')
        mu_g = pm.Deterministic("mu_g", pm.math.dot(basis_g, beta_g), dims='measurement')

        ## Individual or Year Specific Spline Modifications
        if model_type in ['mixed', 'unpooled']:
            sigma = pm.HalfCauchy('sigma_i', 1)
            basis_i = pm.MutableData('Bi', np.asfortranarray(Bi))
            beta = pm.Normal("beta", mu=0, sigma=sigma, dims=('splines', 'years'))
            mui = pm.Deterministic("mui", pm.math.dot(basis_i, beta), dims=('measurement', 'years'))
        
        ## Features
        prem = pm.MutableData('prem', loss_df['standardised_premium'].values)
        grcode = pm.MutableData('grcode', loss_df['grcode'] == 43)

        beta_prem = pm.Normal('beta_prem', 0, 1)
        beta_grcode = pm.Normal('beta_grcode', 0, 1)
        mu_prem = beta_prem*prem
        mu_grcode = beta_grcode*grcode

        ## Likelihood
        sigma = pm.TruncatedNormal("sigma", 1, lower=0.005)
        if model_type == 'mixed':
            mu = pm.Deterministic('mu',  mu_grcode + mu_prem + (mu_g.T + mui.T).ravel(), dims='obs')
            lr_likelihood = pm.Normal("lr", mu, sigma, observed=observed, dims=('obs'))
        elif model_type == 'pooled': 
             lr_likelihood = pm.Normal("lr",  np.repeat(mu_g, len(unique_codes)), sigma, observed=observed, dims='obs')
        elif model_type == 'unpooled':
            lr_likelihood = pm.Normal("lr",  mui.T.ravel(), sigma, observed=observed, dims=('obs'))


        ## Sampling
        idata_sp_insur = pm.sample(2000, return_inferencedata=True, target_accept=.99,
        idata_kwargs={"log_likelihood": True})
        idata_sp_insur = pm.sample_posterior_predictive(
            idata_sp_insur,extend_inferencedata=True)

    return idata_sp_insur, sp_insur


idata_sp_insur_unpooled, sp_insur_unpooled = make_model(loss_df, model_type='unpooled')
idata_sp_insur_pooled, sp_insur_pooled = make_model(loss_df, model_type='pooled')
idata_sp_insur_mixed, sp_insur_mixed = make_model(loss_df, model_type='mixed')

The model structure can be seen more clearly in this graph. The important features to note are the spline components and the manner of their combination. We are postulating a kind of hierarchy of development curves. We specify the “global” curve in 7 development steps and allow that each distinct year is a realisation of a curve that can be thought of as modifying the profile of “global” curve as mui is combined with mu_g. These structures are then combined into a more traditional regression and fed into our likelihood term.

pm.model_to_graphviz(sp_insur_mixed)

We can extract the effect of the differences grcodes and examine the baseline and annual spline related coefficients.

summary = az.summary(idata_sp_insur_mixed, var_names=['sigma', 'beta_grcode', 'beta_prem', 'beta_g', 'beta'])

summary

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
sigma	0.014	0.007	0.005	0.027	0.001	0.000	153.0	348.0	1.04
beta_grcode	0.231	0.057	0.123	0.337	0.001	0.001	4479.0	4712.0	1.00
beta_prem	-0.031	0.029	-0.086	0.023	0.000	0.000	4884.0	5002.0	1.00
beta_g[0]	0.132	0.069	0.006	0.265	0.001	0.001	3185.0	3812.0	1.00
beta_g[1]	0.070	0.395	-0.651	0.841	0.006	0.006	4464.0	3851.0	1.00
...	...	...	...	...	...	...	...	...	...
beta[6, 1990_43]	0.547	0.073	0.412	0.686	0.001	0.001	2884.0	4276.0	1.00
beta[6, 1991_353]	-0.098	0.070	-0.225	0.034	0.001	0.001	2873.0	4230.0	1.00
beta[6, 1991_43]	0.115	0.073	-0.016	0.257	0.001	0.001	3073.0	4496.0	1.00
beta[6, 1992_353]	-0.115	0.070	-0.250	0.012	0.001	0.001	2861.0	4547.0	1.00
beta[6, 1992_43]	-0.033	0.091	-0.212	0.135	0.002	0.001	3726.0	4461.0	1.00

80 rows × 9 columns

Again we can compare the performance metrics of the various models.

compare_df = az.compare({'unpooled': idata_sp_insur_unpooled, 
            'pooled': idata_sp_insur_pooled, 
            'mixed': idata_sp_insur_mixed})

az.plot_compare(compare_df)
compare_df

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
mixed	0	134.697068	61.393405	0.000000	1.000000e+00	1.382933	0.000000	True	log
unpooled	1	88.378689	73.735810	46.318379	0.000000e+00	1.267429	1.520694	True	log
pooled	2	-9.288780	5.995182	143.985848	6.050367e-10	5.414548	5.020442	False	log

The hierarchical splne model is by far the best fit to our data. What does it represent? Yes, splines are function approximation tools. Sure, extrapolation is dangerous, but in this domain where the range of the development curves is constrained they seem to offer an elegant way of articulating the structure of the risk profiles for various products.

Plot the Posterior Predictive Checks

We can check how well the model can recapture the observed data.

def plot_ppc_splines(idata):
    fig, axs = plt.subplots(2, 5, figsize=(20, 10), sharey=True)
    axs = axs.flatten()
    dev_periods = np.arange(1, 7, 1)
    uniques, unique_codes = pd.factorize(loss_df['year_code'])
    for y, c in zip(unique_codes, range(10)):
        az.plot_hdi(dev_periods, idata['posterior_predictive']['lr'].sel(obs=c), color='firebrick', ax=axs[c], fill_kwargs={'alpha': 0.2}, hdi_prob=.89)
        az.plot_hdi(dev_periods, idata['posterior_predictive']['lr'].sel(obs=c), color='firebrick', ax=axs[c], hdi_prob=0.5)
        axs[c].scatter(dev_periods, loss_df[(loss_df['year_code'] == y)]['lr'], color='k', label='Actual Loss Ratio')
        axs[c].plot(dev_periods, loss_df[(loss_df['year_code'] == y)]['lr'], color='k')
        axs[c].set_title(f"PPC 89% and 50% HDI: {y}")
        axs[c].set_ylabel("Loss Ratio")
        axs[c].legend();

plot_ppc_splines(idata_sp_insur_mixed);

Plot the Hierarchical Components

In the following plot we show similarly how to recapture the observed data, but additionally we can decompose the structure of the model in each case and extract baseline forecasts which would be our guide to future loss-ratio development curves in lieu of any other information.

mu_g = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))["mu_g"]

mu_i = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))["mui"]

mu  = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))['mu']

beta_grcode = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))['beta_grcode']

dev_periods = np.arange(1, 7, 1)
uniques, unique_codes = pd.factorize(loss_df['year_code'])

mosaic = """
         ABCDE
         FGHIJ
         KKKKK
"""
fig, axs = plt.subplot_mosaic(mosaic, sharey=True, 
figsize=(20, 15))
axs = [axs[k] for k in axs.keys()] 

mu_g_mean = mu_g.mean(dim='draws')
for y, c in zip(unique_codes, range(10)):
    group_effect = 0
    if '43' in y: 
        group_effect = beta_grcode.mean().item()
    mu_i_mean = mu_i.sel(years=y).mean(dim='draws')
    axs[c].plot(dev_periods, group_effect + mu_g_mean.values + mu_i_mean.values, label='Combined + E(grp_effect)', color='purple', linewidth=3.5)
    axs[c].plot(dev_periods, group_effect + mu_g_mean.values, label='E(Hierarchical Baseline)', color='red', linestyle='--')
    axs[c].plot(dev_periods,  group_effect + mu_i_mean.values, label='E(Year Specific Adjustment)', color='blue', linestyle='--')
    axs[c].scatter(dev_periods, loss_df[(loss_df['year_code'] == y)]['lr'], color='k', label='Actual Loss Ratio')
    az.plot_hdi(dev_periods,mu.sel(obs=c).T  , ax=axs[c], color='firebrick', fill_kwargs={'alpha': 0.2})
    az.plot_hdi(dev_periods, mu.sel(obs=c).T , ax=axs[c], color='firebrick', fill_kwargs={'alpha': 0.5}, hdi_prob=.50)
    axs[c].set_title(f"Components for Year {y}")
    axs[c].set_ylabel("Loss Ratio")
    if (c == 0):
         axs[c].legend()

axs[10].plot(dev_periods, mu_g_mean.values, label='E(Hierarchical Baseline)', color='black')
axs[10].plot(dev_periods, mu_g_mean.values + group_effect, label='E(Hierarchical Baseline) + E(grp_effect)', color='black', linestyle='--')
az.plot_hdi(dev_periods, mu_g.T.values, color='slateblue', ax=axs[10], fill_kwargs={'alpha': 0.2})
az.plot_hdi(dev_periods, mu_g.T.values + group_effect, color='magenta', ax=axs[10], fill_kwargs={'alpha': 0.2})
az.plot_hdi(dev_periods, mu_g.T.values, color='slateblue', ax=axs[10], hdi_prob=.5)
az.plot_hdi(dev_periods, mu_g.T.values  + group_effect, color='magenta', ax=axs[10], hdi_prob=.5)
axs[10].set_title("Baseline Forecast Loss Ratio")
axs[10].legend();

This plot is a bit clunky, because we’re mixing expectations and posterior distributions over the parameters. The point is just to highlight the “compositional” structure of our model.

Predicting Next Year’s Losses

A better way to interrogate the implications of the model is to “push” forward different data through the posterior predictive distribution and derive a kind of ceteris paribus rule for accrual of losses.

with sp_insur_mixed: 
    pm.set_data({'grcode': np.ones(len(loss_df)), 
    })
    idata_43 = pm.sample_posterior_predictive(idata_sp_insur_mixed, var_names=['lr'], extend_inferencedata =True)

with sp_insur_mixed: 
    pm.set_data({'grcode': np.zeros(len(loss_df))})
    idata_353 = pm.sample_posterior_predictive(idata_sp_insur_mixed, var_names=['lr'], extend_inferencedata=True)

idata_353

Sampling: [lr]

100.00% [8000/8000 00:00<00:00]

Sampling: [lr]

100.00% [8000/8000 00:00<00:00]

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:      (chain: 4, draw: 2000, splines: 7, years: 10, measurement: 6,
                    obs: 60)
  Coordinates:
    * chain        (chain) int64 0 1 2 3
    * draw         (draw) int64 0 1 2 3 4 5 6 ... 1994 1995 1996 1997 1998 1999
    * splines      (splines) int64 0 1 2 3 4 5 6
    * years        (years) <U8 '1988_353' '1988_43' ... '1992_353' '1992_43'
    * measurement  (measurement) int64 0 1 2 3 4 5
    * obs          (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 ... 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      beta_g       (chain, draw, splines) float64 0.08213 0.08078 ... 0.6575
      beta         (chain, draw, splines, years) float64 0.1395 ... -0.1751
      beta_prem    (chain, draw) float64 -0.0156 -0.03719 ... -0.0874 -0.06236
      beta_grcode  (chain, draw) float64 0.2377 0.2905 0.2288 ... 0.3751 0.3465
      tau          (chain, draw) float64 0.5886 0.5384 0.4783 ... 0.6427 0.4663
      sigma_i      (chain, draw) float64 0.2091 0.171 0.2062 ... 0.2064 0.1871
      sigma        (chain, draw) float64 0.01148 0.009831 ... 0.005943 0.005919
      mu_g         (chain, draw, measurement) float64 0.08213 0.3377 ... 0.6575
      mui          (chain, draw, measurement, years) float64 0.1395 ... -0.1751
      mu           (chain, draw, obs) float64 0.2271 0.509 ... 0.9015 0.9269
  Attributes:
      created_at:                 2024-05-27T06:09:05.927942
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3
      sampling_time:              227.7068202495575
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • splines: 7
    • years: 10
    • measurement: 6
    • obs: 60
  • Coordinates: (6)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

    • splines
      (splines)
      int64
      0 1 2 3 4 5 6

      array([0, 1, 2, 3, 4, 5, 6])

    • years
      (years)
      <U8
      '1988_353' '1988_43' ... '1992_43'

      array(['1988_353', '1988_43', '1989_353', '1989_43', '1990_353', '1990_43',
             '1991_353', '1991_43', '1992_353', '1992_43'], dtype='<U8')

    • measurement
      (measurement)
      int64
      0 1 2 3 4 5

      array([0, 1, 2, 3, 4, 5])

    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (10)
    • beta_g
      (chain, draw, splines)
      float64
      0.08213 0.08078 ... 0.1945 0.6575

      array([[[ 0.08213429,  0.08078411,  0.43921067, ...,  0.68826351,
                0.76152486,  0.57747165],
              [ 0.0729024 ,  0.34594954,  0.30222266, ...,  0.58894912,
                0.70131994,  0.56389405],
              [ 0.1296795 , -0.07301492,  0.59161604, ...,  0.89261377,
                0.32944684,  0.70257189],
              ...,
              [-0.04196952,  0.17650917,  0.4217089 , ...,  0.62258926,
                0.59316501,  0.79250849],
              [-0.00829219, -0.1684174 ,  0.68416978, ...,  0.81418464,
                0.38975666,  0.74581672],
              [ 0.0817568 ,  0.21642253,  0.4364298 , ...,  0.70190257,
                0.48913041,  0.62736766]],
      
             [[ 0.1372587 ,  0.23309973,  0.49454757, ...,  0.65137438,
                0.75772987,  0.61066381],
              [ 0.01118497, -0.501921  ,  0.70257164, ...,  0.86052648,
               -0.11635177,  0.5797538 ],
              [ 0.19565787,  0.27753886,  0.68404036, ...,  0.96958051,
                0.72825027,  0.73195019],
      ...
              [ 0.13938073, -0.02592998,  0.73699995, ...,  0.88935344,
                0.379609  ,  0.81252247],
              [ 0.11689267,  0.04623414,  0.67605048, ...,  0.90426575,
                0.4922117 ,  0.70663987],
              [ 0.12929436,  0.12568724,  0.61095241, ...,  0.75010102,
                0.58792975,  0.74250617]],
      
             [[ 0.1464027 ,  0.76694943, -0.06153482, ...,  0.15534104,
                1.32746341,  0.66072813],
              [ 0.14358004,  0.68940626, -0.11149118, ...,  0.17124399,
                1.18011474,  0.71314634],
              [ 0.05419155, -0.42883811,  1.14901708, ...,  1.26919344,
                0.01586078,  0.63687761],
              ...,
              [ 0.18723226,  0.48566063,  0.23265438, ...,  0.5794988 ,
                0.85015745,  0.86002011],
              [ 0.03751469,  0.0086789 ,  0.55346192, ...,  0.8412314 ,
                0.43493969,  0.65331407],
              [ 0.05265733, -0.19628186,  0.64580651, ...,  0.90007346,
                0.19448681,  0.65747667]]])

    • beta
      (chain, draw, splines, years)
      float64
      0.1395 -0.1399 ... -0.1164 -0.1751

      array([[[[ 1.39520480e-01, -1.39920512e-01,  1.59795716e-01, ...,
                -4.88250518e-02,  1.49757586e-01, -7.87841938e-02],
               [ 1.57822205e-01,  1.57221882e-01,  8.80230108e-02, ...,
                 1.84740484e-01,  1.40654109e-01,  1.05066007e-01],
               [ 1.71863836e-01, -3.89547722e-01,  2.10933769e-01, ...,
                -1.10459344e-02,  2.27594707e-02, -1.97138584e-01],
               ...,
               [-1.36108960e-01, -1.82625112e-01,  4.75130784e-02, ...,
                 4.45905407e-02, -3.94624785e-02, -1.40864962e-01],
               [ 6.71214512e-02, -3.95470614e-01, -1.96252840e-01, ...,
                -2.64734451e-02, -3.31353605e-01,  2.14460499e-02],
               [ 1.11116707e-01, -1.45940962e-01,  3.79732609e-02, ...,
                 1.98087329e-01, -8.13679795e-03,  8.28432092e-02]],
      
              [[ 1.46217147e-01, -1.45987241e-01,  1.33636350e-01, ...,
                -1.04141339e-01,  1.56999588e-01, -1.73189166e-01],
               [ 1.77012604e-01,  1.41028626e-02,  7.00230739e-02, ...,
                -6.63265264e-02,  5.91449369e-02, -2.69751242e-01],
               [ 1.41749001e-01, -3.79009870e-01,  1.25803026e-01, ...,
                 2.14529318e-02,  5.38979285e-02, -1.33607062e-01],
      ...
               [ 3.45479799e-02, -3.58625458e-01, -2.18827643e-01, ...,
                -1.90020093e-01, -2.01458145e-01, -3.92093686e-01],
               [-1.46373265e-01, -5.99317086e-02,  1.30108020e-01, ...,
                 9.35956673e-02, -9.62921164e-02, -1.04324354e-01],
               [ 1.15647935e-02, -2.22568566e-01, -4.62102677e-02, ...,
                -9.27173884e-03, -1.22732622e-01, -2.49213829e-01]],
      
              [[ 1.63621902e-01, -1.51455454e-01,  1.65518502e-01, ...,
                -1.56291791e-01,  1.56621714e-01, -2.49188309e-01],
               [-1.31496227e-02,  1.33838754e-01, -1.48588957e-01, ...,
                -8.37641665e-02,  9.22397066e-02, -3.70879910e-01],
               [ 2.55943473e-01, -4.43680242e-01,  2.87918409e-01, ...,
                -6.67540255e-02, -1.11330929e-03, -1.87938351e-01],
               ...,
               [ 1.05025544e-01, -2.90129799e-01,  1.05270857e-01, ...,
                -4.30757899e-02, -1.14967846e-01, -1.59326571e-01],
               [-3.26459908e-02, -8.45911318e-02, -1.55795510e-01, ...,
                 1.08681450e-01, -7.15263058e-02, -1.45134288e-01],
               [ 8.68094467e-03, -2.62496397e-01, -4.86758721e-02, ...,
                 9.04772502e-03, -1.16437070e-01, -1.75097569e-01]]]])

    • beta_prem
      (chain, draw)
      float64
      -0.0156 -0.03719 ... -0.06236

      array([[-0.01560439, -0.03719103, -0.04590153, ..., -0.04335828,
              -0.06037419,  0.00715623],
             [-0.08052324, -0.01663825, -0.07448925, ..., -0.04856727,
              -0.00828554, -0.03699675],
             [-0.01992275, -0.04479077, -0.05495835, ...,  0.03346062,
              -0.00432382, -0.03984384],
             [-0.06903017, -0.03561944, -0.06378486, ..., -0.04559814,
              -0.08739882, -0.06236288]])

    • beta_grcode
      (chain, draw)
      float64
      0.2377 0.2905 ... 0.3751 0.3465

      array([[0.23773592, 0.29052938, 0.22882446, ..., 0.22874581, 0.27521708,
              0.20095271],
             [0.31476531, 0.33120918, 0.10765162, ..., 0.31997823, 0.19133539,
              0.27806806],
             [0.21980155, 0.17918741, 0.25618125, ..., 0.16353416, 0.13860823,
              0.19789005],
             [0.26856768, 0.26101595, 0.17110103, ..., 0.2398666 , 0.37505705,
              0.34646743]])

    • tau
      (chain, draw)
      float64
      0.5886 0.5384 ... 0.6427 0.4663

      array([[0.58855469, 0.5383688 , 0.47833408, ..., 0.35456915, 1.00747773,
              0.63732983],
             [0.37739955, 0.53458212, 0.66920479, ..., 0.48638153, 0.59084796,
              0.76687933],
             [0.79058091, 0.78629887, 0.47156563, ..., 0.6155269 , 0.84141259,
              0.90084556],
             [0.74634296, 0.56718238, 0.63635355, ..., 0.53907184, 0.64268993,
              0.46630978]])

    • sigma_i
      (chain, draw)
      float64
      0.2091 0.171 ... 0.2064 0.1871

      array([[0.20909747, 0.17103069, 0.20618882, ..., 0.211108  , 0.19893185,
              0.21397901],
             [0.19959592, 0.21750556, 0.23931261, ..., 0.20211561, 0.18453507,
              0.20412235],
             [0.19160775, 0.19335863, 0.1767455 , ..., 0.19974999, 0.19592067,
              0.2032291 ],
             [0.20263706, 0.19317583, 0.18770624, ..., 0.19822884, 0.20642422,
              0.18713983]])

    • sigma
      (chain, draw)
      float64
      0.01148 0.009831 ... 0.005919

      array([[0.0114768 , 0.0098313 , 0.01093488, ..., 0.00679955, 0.00660215,
              0.00630935],
             [0.00763042, 0.00790286, 0.00785762, ..., 0.01620113, 0.01155728,
              0.01278793],
             [0.01091791, 0.01269563, 0.0112426 , ..., 0.0160019 , 0.00890353,
              0.00816588],
             [0.01126429, 0.01476696, 0.01066087, ..., 0.00525937, 0.00594278,
              0.00591885]])

    • mu_g
      (chain, draw, measurement)
      float64
      0.08213 0.3377 ... 0.6124 0.6575

      array([[[ 0.08213429,  0.33768236,  0.49079526,  0.58195553,
                0.6893857 ,  0.57747165],
              [ 0.0729024 ,  0.3523693 ,  0.47549742,  0.57555192,
                0.62509911,  0.56389405],
              [ 0.1296795 ,  0.36973404,  0.53621521,  0.64158577,
                0.66511154,  0.70257189],
              ...,
              [-0.04196952,  0.37081521,  0.53617268,  0.60818718,
                0.61363189,  0.79250849],
              [-0.00829219,  0.38796381,  0.54652576,  0.59643256,
                0.63414497,  0.74581672],
              [ 0.0817568 ,  0.38371482,  0.51608403,  0.6080783 ,
                0.61884953,  0.62736766]],
      
             [[ 0.1372587 ,  0.40331188,  0.4591502 ,  0.51765001,
                0.65439136,  0.61066381],
              [ 0.01118497,  0.26250702,  0.41669178,  0.47397997,
                0.47131783,  0.5797538 ],
              [ 0.19565787,  0.55795127,  0.69739504,  0.79814199,
                0.8594259 ,  0.73195019],
      ...
              [ 0.13938073,  0.47041618,  0.61011331,  0.66575445,
                0.68220337,  0.81252247],
              [ 0.11689267,  0.45871305,  0.59006843,  0.67156949,
                0.72612402,  0.70663987],
              [ 0.12929436,  0.45995572,  0.6102857 ,  0.66215402,
                0.68172986,  0.74250617]],
      
             [[ 0.1464027 ,  0.32512369,  0.50687362,  0.58588734,
                0.62142815,  0.66072813],
              [ 0.14358004,  0.26546548,  0.4564482 ,  0.55671318,
                0.57768676,  0.71314634],
              [ 0.05419155,  0.53789895,  0.62458472,  0.66990082,
                0.74346936,  0.63687761],
              ...,
              [ 0.18723226,  0.35784665,  0.45608484,  0.57637053,
                0.66758408,  0.86002011],
              [ 0.03751469,  0.38496915,  0.56650101,  0.66771488,
                0.68016691,  0.65331407],
              [ 0.05265733,  0.34718215,  0.49383633,  0.58342827,
                0.61237385,  0.65747667]]])

    • mui
      (chain, draw, measurement, years)
      float64
      0.1395 -0.1399 ... -0.1164 -0.1751

      array([[[[ 1.39520480e-01, -1.39920512e-01,  1.59795716e-01, ...,
                -4.88250518e-02,  1.49757586e-01, -7.87841938e-02],
               [ 1.65873605e-01, -1.93497148e-01,  1.47566232e-01, ...,
                 7.87630258e-02,  5.26650556e-02, -5.69859310e-02],
               [ 1.42490699e-01, -2.54958783e-01,  9.21905855e-02, ...,
                 1.17105881e-01, -7.24684266e-03,  4.86017675e-03],
               [ 3.83479319e-02, -1.91999665e-01,  3.42301192e-02, ...,
                 1.33318539e-01, -3.34724977e-02,  2.27443612e-02],
               [-3.72703490e-02, -2.46357200e-01, -3.30969747e-02, ...,
                 4.52311029e-02, -1.28371119e-01, -5.07836293e-02],
               [ 1.11116707e-01, -1.45940962e-01,  3.79732609e-02, ...,
                 1.98087329e-01, -8.13679795e-03,  8.28432092e-02]],
      
              [[ 1.46217147e-01, -1.45987241e-01,  1.33636350e-01, ...,
                -1.04141339e-01,  1.56999588e-01, -1.73189166e-01],
               [ 1.50806947e-01, -2.30291464e-01,  1.09203466e-01, ...,
                 9.20703581e-03,  4.54642291e-02, -1.56726472e-01],
               [ 1.22668583e-01, -2.49806973e-01,  1.13781261e-01, ...,
                 8.35633149e-02,  1.22118628e-03, -5.88513761e-02],
               [ 6.62855094e-02, -2.01179440e-01,  4.96226930e-02, ...,
      ...
                -9.33073690e-02, -1.27232209e-01, -3.13057795e-01],
               [-5.68682820e-02, -2.92825044e-01, -7.16499541e-02, ...,
                -9.83126153e-02, -1.67003142e-01, -3.19813169e-01],
               [-4.65304179e-02, -2.48490531e-01, -8.00679145e-02, ...,
                -8.15331960e-02, -1.62924280e-01, -2.87153932e-01],
               [ 1.15647935e-02, -2.22568566e-01, -4.62102677e-02, ...,
                -9.27173884e-03, -1.22732622e-01, -2.49213829e-01]],
      
              [[ 1.63621902e-01, -1.51455454e-01,  1.65518502e-01, ...,
                -1.56291791e-01,  1.56621714e-01, -2.49188309e-01],
               [ 1.36287416e-01, -2.29995142e-01,  1.06986574e-01, ...,
                -5.57859181e-02,  2.19595602e-02, -2.41298686e-01],
               [ 9.00781502e-02, -2.78835663e-01,  7.79656367e-02, ...,
                 2.00867494e-03, -3.20768389e-02, -1.75140835e-01],
               [ 3.93295341e-02, -2.29924865e-01,  1.67347175e-02, ...,
                 1.21352879e-02, -7.20658012e-02, -1.62981950e-01],
               [ 4.44384071e-02, -2.10066162e-01,  8.83961075e-04, ...,
                 1.70640485e-02, -9.33275239e-02, -1.55367387e-01],
               [ 8.68094467e-03, -2.62496397e-01, -4.86758721e-02, ...,
                 9.04772502e-03, -1.16437070e-01, -1.75097569e-01]]]])

    • mu
      (chain, draw, obs)
      float64
      0.2271 0.509 ... 0.9015 0.9269

      array([[[0.22707658, 0.50897777, 0.63870776, ..., 0.86695823,
               0.9008604 , 0.9225732 ],
              [0.2320417 , 0.5160984 , 0.61108815, ..., 0.85957629,
               0.88786053, 0.91962846],
              [0.22591164, 0.49834019, 0.63366749, ..., 0.87302862,
               0.8952962 , 0.91718248],
              ...,
              [0.21752702, 0.50863783, 0.62085732, ..., 0.85667073,
               0.90717039, 0.91667823],
              [0.24668158, 0.50946944, 0.61496311, ..., 0.85950557,
               0.91499729, 0.93394897],
              [0.22803856, 0.51225596, 0.60431136, ..., 0.85384508,
               0.91546043, 0.92367727]],
      
             [[0.20798215, 0.49653525, 0.61163146, ..., 0.84536501,
               0.90679466, 0.9304019 ],
              [0.21296999, 0.52098307, 0.6079544 , ..., 0.85788457,
               0.90174194, 0.91528221],
              [0.21005812, 0.50465123, 0.59585638, ..., 0.85754494,
               0.90510048, 0.92520922],
      ...
              [0.23575948, 0.52517243, 0.61730376, ..., 0.86419416,
               0.88562497, 0.92776999],
              [0.23289817, 0.51858843, 0.61101066, ..., 0.86935524,
               0.89006191, 0.92026702],
              [0.24012795, 0.50146092, 0.611026  , ..., 0.85900662,
               0.9094241 , 0.916668  ]],
      
             [[0.24430336, 0.49530459, 0.60862876, ..., 0.86711487,
               0.89899689, 0.89876914],
              [0.23968997, 0.50137434, 0.61292425, ..., 0.87145524,
               0.90762533, 0.91330743],
              [0.22874435, 0.5038266 , 0.63214802, ..., 0.86415242,
               0.900888  , 0.92221532],
              ...,
              [0.22737739, 0.50688812, 0.61785909, ..., 0.86241171,
               0.90291879, 0.93063882],
              [0.2253554 , 0.52076927, 0.61201075, ..., 0.86030662,
               0.90541789, 0.91650515],
              [0.23794743, 0.50513776, 0.60558268, ..., 0.86491745,
               0.9014776 , 0.92685024]]])

  • Indexes: (6)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

    • splines
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6], dtype='int64', name='splines'))

    • years
      PandasIndex

      PandasIndex(Index(['1988_353', '1988_43', '1989_353', '1989_43', '1990_353', '1990_43',
             '1991_353', '1991_43', '1992_353', '1992_43'],
            dtype='object', name='years'))

    • measurement
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5], dtype='int64', name='measurement'))

    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (6)

    created_at :

    2024-05-27T06:09:05.927942

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    227.7068202495575

    tuning_steps :

    1000

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 2000, obs: 60)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 1993 1994 1995 1996 1997 1998 1999
    * obs      (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 2 ... 7 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      lr       (chain, draw, obs) float64 0.233 0.5237 0.6346 ... 0.8981 0.9297
  Attributes:
      created_at:                 2024-05-27T06:09:08.031307
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • obs: 60
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (1)
    • lr
      (chain, draw, obs)
      float64
      0.233 0.5237 ... 0.8981 0.9297

      array([[[0.2329779 , 0.52369396, 0.63459651, ..., 0.88727935,
               0.89708185, 0.92159448],
              [0.22917617, 0.52664668, 0.59261549, ..., 0.86452905,
               0.90716367, 0.92848945],
              [0.21781397, 0.4993653 , 0.62570146, ..., 0.87851744,
               0.91024927, 0.90860314],
              ...,
              [0.21017767, 0.50642437, 0.61274687, ..., 0.86717561,
               0.90019066, 0.92583528],
              [0.25375665, 0.50843402, 0.60506852, ..., 0.85534182,
               0.91866865, 0.94337612],
              [0.22315311, 0.50562829, 0.60910443, ..., 0.85065712,
               0.91298724, 0.92813007]],
      
             [[0.21910589, 0.51068439, 0.60243355, ..., 0.82541821,
               0.91008109, 0.93325851],
              [0.20831853, 0.51436674, 0.61926284, ..., 0.85784094,
               0.91192378, 0.91751401],
              [0.211794  , 0.50545555, 0.59703336, ..., 0.84964635,
               0.90139699, 0.93754076],
      ...
              [0.20342139, 0.53073673, 0.58714448, ..., 0.85629361,
               0.86931632, 0.93170111],
              [0.22936355, 0.51296337, 0.62169933, ..., 0.86757595,
               0.87222132, 0.91036641],
              [0.23170526, 0.50702719, 0.615165  , ..., 0.86935216,
               0.89592565, 0.9289633 ]],
      
             [[0.23720474, 0.49880944, 0.60052934, ..., 0.86819616,
               0.89829526, 0.89245385],
              [0.21610787, 0.49419633, 0.59675283, ..., 0.87135881,
               0.91189684, 0.92522446],
              [0.21011217, 0.48862193, 0.62340207, ..., 0.85734297,
               0.90236241, 0.90941639],
              ...,
              [0.22319361, 0.50040966, 0.62049363, ..., 0.85878946,
               0.90237754, 0.93678927],
              [0.21677194, 0.51807386, 0.61236904, ..., 0.85692388,
               0.90624593, 0.91623522],
              [0.23988418, 0.5116527 , 0.59962956, ..., 0.86614767,
               0.89813962, 0.92968554]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (4)

    created_at :

    2024-05-27T06:09:08.031307

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 2000, obs: 60)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 1993 1994 1995 1996 1997 1998 1999
    * obs      (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 2 ... 7 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      lr       (chain, draw, obs) float64 3.493 3.52 1.659 ... 3.711 4.137 4.128
  Attributes:
      created_at:                 2024-05-27T06:09:06.204725
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • obs: 60
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (1)
    • lr
      (chain, draw, obs)
      float64
      3.493 3.52 1.659 ... 4.137 4.128

      array([[[ 3.49346525,  3.51967093,  1.65940483, ...,  3.3082412 ,
                3.51674037,  3.53525774],
              [ 3.69631071,  3.6046685 ,  3.55729917, ...,  3.70154431,
                2.3967742 ,  3.58347456],
              [ 3.4934784 ,  2.84680694,  2.35000913, ...,  2.77423593,
                3.29784217,  3.37659525],
              ...,
              [ 2.14258977,  3.96834899,  3.85707249, ...,  4.01314636,
                3.94561888,  3.42040059],
              [ 1.23863774,  4.04264235,  4.07774621, ...,  4.09852065,
                2.65095731,  3.06425839],
              [ 4.04510066,  4.14335   ,  2.31136224, ...,  3.81265832,
                2.42504945,  4.13947279]],
      
             [[-0.54742785,  1.97315188,  3.76142259, ...,  2.35947554,
                3.87719424,  3.6514542 ],
              [ 1.35250133,  3.23660373,  3.35059476, ...,  3.91158231,
                3.88923252,  3.25015959],
              [ 0.41505901,  3.52118601,  0.50965773, ...,  3.91012279,
                3.91261702,  3.9225444 ],
      ...
              [ 3.16969071,  2.86342788,  3.21451326, ...,  3.16348367,
                2.57445527,  3.19445904],
              [ 3.77677685,  3.50595189,  3.61919133, ...,  3.12638761,
                2.62018425,  3.69252388],
              [ 3.24809316,  3.09764755,  3.67232438, ...,  3.88885252,
                3.64764734,  3.43590068]],
      
             [[ 2.85754207,  2.50362844,  3.32921447, ...,  3.30786291,
                3.47806162,  0.97031193],
              [ 3.11861345,  3.05038605,  3.26872831, ...,  2.94087519,
                3.26203195,  3.01224386],
              [ 3.60209919,  3.34722452,  2.5311683 , ...,  3.50557154,
                3.58614045,  3.60045959],
              ...,
              [ 4.10655831,  3.90449352,  4.29030831, ...,  4.11874516,
                4.31624387,  3.63427395],
              [ 3.77393089,  3.05063185,  3.93392071, ...,  4.18257176,
                4.16737028,  3.31520253],
              [ 3.49856554,  3.58985088,  2.54068726, ...,  3.71137697,
                4.136818  ,  4.12777749]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (4)

    created_at :

    2024-05-27T06:09:06.204725

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:                (chain: 4, draw: 2000)
  Coordinates:
    * chain                  (chain) int64 0 1 2 3
    * draw                   (draw) int64 0 1 2 3 4 5 ... 1995 1996 1997 1998 1999
  Data variables: (12/17)
      perf_counter_diff      (chain, draw) float64 0.08216 0.0801 ... 0.08181
      lp                     (chain, draw) float64 184.4 196.5 ... 216.8 224.4
      step_size_bar          (chain, draw) float64 0.005908 0.005908 ... 0.005495
      smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan
      acceptance_rate        (chain, draw) float64 0.9998 0.9855 ... 0.944 0.9537
      reached_max_treedepth  (chain, draw) bool False False False ... True False
      ...                     ...
      largest_eigval         (chain, draw) float64 nan nan nan nan ... nan nan nan
      energy                 (chain, draw) float64 -141.4 -154.8 ... -177.6 -175.2
      tree_depth             (chain, draw) int64 10 10 10 10 10 ... 10 10 10 10 10
      diverging              (chain, draw) bool False False False ... False False
      index_in_trajectory    (chain, draw) int64 -352 387 -740 ... 790 521 -683
      process_time_diff      (chain, draw) float64 0.08216 0.08007 ... 0.0818
  Attributes:
      created_at:                 2024-05-27T06:09:05.940582
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3
      sampling_time:              227.7068202495575
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

  • Data variables: (17)
    • perf_counter_diff
      (chain, draw)
      float64
      0.08216 0.0801 ... 0.08166 0.08181

      array([[0.08215638, 0.08010133, 0.08184954, ..., 0.08156783, 0.08167758,
              0.08158583],
             [0.08136029, 0.08029138, 0.08061071, ..., 0.08166729, 0.08090975,
              0.08182967],
             [0.08155483, 0.08204283, 0.08018262, ..., 0.08148587, 0.08159787,
              0.08199237],
             [0.080971  , 0.08061054, 0.07909571, ..., 0.08166879, 0.08165904,
              0.08180558]])

    • lp
      (chain, draw)
      float64
      184.4 196.5 191.3 ... 216.8 224.4

      array([[184.36486414, 196.53933231, 191.27157196, ..., 200.00169387,
              204.82818005, 214.33525947],
             [192.7863688 , 187.12040812, 189.2459854 , ..., 155.91231505,
              178.6049379 , 187.65005752],
             [188.86778753, 183.74271399, 179.99640544, ..., 167.53100342,
              193.17434736, 201.62651467],
             [178.50910416, 177.79847515, 184.59234725, ..., 220.75603875,
              216.79913757, 224.43676879]])

    • step_size_bar
      (chain, draw)
      float64
      0.005908 0.005908 ... 0.005495

      array([[0.0059085 , 0.0059085 , 0.0059085 , ..., 0.0059085 , 0.0059085 ,
              0.0059085 ],
             [0.00804049, 0.00804049, 0.00804049, ..., 0.00804049, 0.00804049,
              0.00804049],
             [0.00597803, 0.00597803, 0.00597803, ..., 0.00597803, 0.00597803,
              0.00597803],
             [0.00549452, 0.00549452, 0.00549452, ..., 0.00549452, 0.00549452,
              0.00549452]])

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • acceptance_rate
      (chain, draw)
      float64
      0.9998 0.9855 ... 0.944 0.9537

      array([[0.99979105, 0.98550281, 0.97845882, ..., 0.98798177, 0.99585789,
              0.91300196],
             [0.88604944, 0.92966456, 0.9978866 , ..., 0.97677582, 0.99334874,
              0.9967715 ],
             [0.99347033, 0.99333916, 0.95811576, ..., 0.97819394, 0.99938068,
              0.99843286],
             [0.99895433, 0.99905734, 0.97817478, ..., 0.94058163, 0.94399245,
              0.95368363]])

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... True False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [ True, False, False, ...,  True, False,  True],
             [False, False, False, ..., False,  True, False]])

    • step_size
      (chain, draw)
      float64
      0.004065 0.004065 ... 0.007016

      array([[0.00406536, 0.00406536, 0.00406536, ..., 0.00406536, 0.00406536,
              0.00406536],
             [0.00551699, 0.00551699, 0.00551699, ..., 0.00551699, 0.00551699,
              0.00551699],
             [0.00537033, 0.00537033, 0.00537033, ..., 0.00537033, 0.00537033,
              0.00537033],
             [0.00701583, 0.00701583, 0.00701583, ..., 0.00701583, 0.00701583,
              0.00701583]])

    • perf_counter_start
      (chain, draw)
      float64
      2.361e+06 2.361e+06 ... 2.361e+06

      array([[2361318.79279567, 2361318.87508063, 2361318.95529054, ...,
              2361479.04572258, 2361479.12740954, 2361479.20918579],
             [2361312.089392  , 2361312.17086767, 2361312.25125425, ...,
              2361457.72276992, 2361457.80452825, 2361457.88554038],
             [2361312.71388937, 2361312.79556071, 2361312.87769787, ...,
              2361473.83453225, 2361473.91610921, 2361473.99781033],
             [2361315.01769779, 2361315.098796  , 2361315.17949475, ...,
              2361475.50722583, 2361475.58900379, 2361475.67075933]])

    • n_steps
      (chain, draw)
      float64
      1.023e+03 1.023e+03 ... 1.023e+03

      array([[1023., 1023., 1023., ..., 1023., 1023., 1023.],
             [1023., 1023., 1023., ..., 1023., 1023., 1023.],
             [1023., 1023., 1023., ..., 1023., 1023., 1023.],
             [1023., 1023., 1023., ..., 1023., 1023., 1023.]])

    • max_energy_error
      (chain, draw)
      float64
      -0.03173 -0.06478 ... 0.1507 0.15

      array([[-0.03173012, -0.06477837, -0.2075868 , ..., -0.24468182,
              -0.11603018,  0.29430899],
             [-0.84910791, -0.52014748, -0.76125186, ..., -0.32416101,
              -0.36622319, -0.5667363 ],
             [-0.15657836, -0.25783916, -0.18139475, ...,  0.06635242,
              -0.05819053, -0.06931609],
             [-0.13684815, -0.08643352, -0.15030015, ...,  0.17884379,
               0.15065315,  0.14999317]])

    • energy_error
      (chain, draw)
      float64
      -0.02417 0.02804 ... -0.06934

      array([[-0.02417412,  0.02804133, -0.01108127, ..., -0.03346141,
              -0.00129006,  0.24695244],
             [ 0.05038906,  0.22945603, -0.36579941, ...,  0.04978283,
               0.07012965, -0.27605762],
             [ 0.00575441, -0.04686598, -0.03784106, ...,  0.00566893,
               0.00287121, -0.0143361 ],
             [-0.00074252, -0.00241762, -0.06416166, ..., -0.02266332,
               0.07892408, -0.06933777]])

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • energy
      (chain, draw)
      float64
      -141.4 -154.8 ... -177.6 -175.2

      array([[-141.35781495, -154.82042184, -158.50830829, ..., -169.51286638,
              -161.92045136, -167.16629184],
             [-161.91301831, -144.81534601, -149.29314566, ..., -124.68284739,
              -126.4725391 , -140.03003424],
             [-142.80921945, -146.5387865 , -146.8058751 , ..., -129.37532231,
              -140.67959809, -156.69836175],
             [-139.25837416, -135.14965861, -136.52022668, ..., -188.92859374,
              -177.63265783, -175.21581418]])

    • tree_depth
      (chain, draw)
      int64
      10 10 10 10 10 ... 10 10 10 10 10

      array([[10, 10, 10, ..., 10, 10, 10],
             [10, 10, 10, ..., 10, 10, 10],
             [10, 10, 10, ..., 10, 10, 10],
             [10, 10, 10, ..., 10, 10, 10]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • index_in_trajectory
      (chain, draw)
      int64
      -352 387 -740 273 ... 790 521 -683

      array([[-352,  387, -740, ...,  643,  428, -650],
             [ 453,  487,  609, ...,  749, -425,  216],
             [-596,  472,  590, ..., -161,  452,  263],
             [-220,  272, -830, ...,  790,  521, -683]])

    • process_time_diff
      (chain, draw)
      float64
      0.08216 0.08007 ... 0.08164 0.0818

      array([[0.082157, 0.080075, 0.081789, ..., 0.081559, 0.081672, 0.08158 ],
             [0.08135 , 0.080263, 0.080593, ..., 0.081654, 0.080874, 0.081816],
             [0.081545, 0.082021, 0.080146, ..., 0.081464, 0.081571, 0.081977],
             [0.080926, 0.080607, 0.079093, ..., 0.081597, 0.081638, 0.0818  ]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

  • Attributes: (6)

    created_at :

    2024-05-27T06:09:05.940582

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    227.7068202495575

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (obs: 60)
  Coordinates:
    * obs      (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 2 ... 7 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      lr       (obs) float64 0.2309 0.5117 0.6164 0.6386 ... 0.859 0.9038 0.9244
  Attributes:
      created_at:                 2024-05-27T06:09:05.943015
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • obs: 60
  • Coordinates: (1)
    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (1)
    • lr
      (obs)
      float64
      0.2309 0.5117 ... 0.9038 0.9244

      array([0.23088384, 0.51173309, 0.61639972, 0.63858884, 0.67259086,
             0.68993774, 0.13897597, 0.34796238, 0.45036573, 0.59561129,
             0.64263323, 0.64263323, 0.2345894 , 0.47973401, 0.58713321,
             0.62354866, 0.63473718, 0.63742875, 0.25277402, 0.47253045,
             0.64005413, 0.69797023, 0.74776725, 0.80270636, 0.29829006,
             0.51463926, 0.59872363, 0.66151898, 0.68509768, 0.70209967,
             0.33072662, 0.79244053, 1.12088628, 1.31753014, 1.40029326,
             1.42440534, 0.25679105, 0.44820942, 0.53957139, 0.57970674,
             0.59681361, 0.59784754, 0.25876918, 0.65744596, 0.86254492,
             0.94998004, 0.98790852, 1.03097017, 0.24606574, 0.41036624,
             0.49781029, 0.54771644, 0.56508975, 0.57923865, 0.25779626,
             0.54739102, 0.76565216, 0.85900276, 0.90375243, 0.9244402 ])

  • Indexes: (1)
    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (4)

    created_at :

    2024-05-27T06:09:05.943015

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  
• constant_data

  <xarray.Dataset>
  Dimensions:       (Bg_dim_0: 6, Bg_dim_1: 7, Bi_dim_0: 6, Bi_dim_1: 7,
                     prem_dim_0: 60, grcode_dim_0: 60)
  Coordinates:
    * Bg_dim_0      (Bg_dim_0) int64 0 1 2 3 4 5
    * Bg_dim_1      (Bg_dim_1) int64 0 1 2 3 4 5 6
    * Bi_dim_0      (Bi_dim_0) int64 0 1 2 3 4 5
    * Bi_dim_1      (Bi_dim_1) int64 0 1 2 3 4 5 6
    * prem_dim_0    (prem_dim_0) int64 0 1 2 3 4 5 6 7 ... 52 53 54 55 56 57 58 59
    * grcode_dim_0  (grcode_dim_0) int64 0 1 2 3 4 5 6 7 ... 53 54 55 56 57 58 59
  Data variables:
      Bg            (Bg_dim_0, Bg_dim_1) float64 1.0 0.0 0.0 0.0 ... 0.0 0.0 1.0
      Bi            (Bi_dim_0, Bi_dim_1) float64 1.0 0.0 0.0 0.0 ... 0.0 0.0 1.0
      prem          (prem_dim_0) float64 -0.3475 -0.3475 -0.3475 ... -1.572 -1.572
      grcode        (grcode_dim_0) float64 0.0 0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0
  Attributes:
      created_at:                 2024-05-27T06:09:05.943907
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • Bg_dim_0: 6
    • Bg_dim_1: 7
    • Bi_dim_0: 6
    • Bi_dim_1: 7
    • prem_dim_0: 60
    • grcode_dim_0: 60
  • Coordinates: (6)
    • Bg_dim_0
      (Bg_dim_0)
      int64
      0 1 2 3 4 5

      array([0, 1, 2, 3, 4, 5])

    • Bg_dim_1
      (Bg_dim_1)
      int64
      0 1 2 3 4 5 6

      array([0, 1, 2, 3, 4, 5, 6])

    • Bi_dim_0
      (Bi_dim_0)
      int64
      0 1 2 3 4 5

      array([0, 1, 2, 3, 4, 5])

    • Bi_dim_1
      (Bi_dim_1)
      int64
      0 1 2 3 4 5 6

      array([0, 1, 2, 3, 4, 5, 6])

    • prem_dim_0
      (prem_dim_0)
      int64
      0 1 2 3 4 5 6 ... 54 55 56 57 58 59

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59])

    • grcode_dim_0
      (grcode_dim_0)
      int64
      0 1 2 3 4 5 6 ... 54 55 56 57 58 59

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59])

  • Data variables: (4)
    • Bg
      (Bg_dim_0, Bg_dim_1)
      float64
      1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 1.0

      array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 3.08571429e-01, 5.69339736e-01, 1.21488595e-01,
              6.00240096e-04, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 1.14285714e-02, 4.09819928e-01, 5.03721489e-01,
              7.50300120e-02, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 7.50300120e-02, 5.03721489e-01,
              4.09819928e-01, 1.14285714e-02, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 6.00240096e-04, 1.21488595e-01,
              5.69339736e-01, 3.08571429e-01, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

    • Bi
      (Bi_dim_0, Bi_dim_1)
      float64
      1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 1.0

      array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 3.08571429e-01, 5.69339736e-01, 1.21488595e-01,
              6.00240096e-04, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 1.14285714e-02, 4.09819928e-01, 5.03721489e-01,
              7.50300120e-02, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 7.50300120e-02, 5.03721489e-01,
              4.09819928e-01, 1.14285714e-02, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 6.00240096e-04, 1.21488595e-01,
              5.69339736e-01, 3.08571429e-01, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

    • prem
      (prem_dim_0)
      float64
      -0.3475 -0.3475 ... -1.572 -1.572

      array([-0.34745346, -0.34745346, -0.34745346, -0.34745346, -0.34745346,
             -0.34745346,  1.72234378,  1.72234378,  1.72234378,  1.72234378,
              1.72234378,  1.72234378, -0.3654406 , -0.3654406 , -0.3654406 ,
             -0.3654406 , -0.3654406 , -0.3654406 ,  1.40460971,  1.40460971,
              1.40460971,  1.40460971,  1.40460971,  1.40460971, -0.5486773 ,
             -0.5486773 , -0.5486773 , -0.5486773 , -0.5486773 , -0.5486773 ,
              1.12110923,  1.12110923,  1.12110923,  1.12110923,  1.12110923,
              1.12110923, -0.63582788, -0.63582788, -0.63582788, -0.63582788,
             -0.63582788, -0.63582788, -0.20123544, -0.20123544, -0.20123544,
             -0.20123544, -0.20123544, -0.20123544, -0.5779209 , -0.5779209 ,
             -0.5779209 , -0.5779209 , -0.5779209 , -0.5779209 , -1.57150713,
             -1.57150713, -1.57150713, -1.57150713, -1.57150713, -1.57150713])

    • grcode
      (grcode_dim_0)
      float64
      0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0

      array([0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0.,
             0., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1.,
             1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1., 0., 0., 0.,
             0., 0., 0., 1., 1., 1., 1., 1., 1.])

  • Indexes: (6)
    • Bg_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5], dtype='int64', name='Bg_dim_0'))

    • Bg_dim_1
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6], dtype='int64', name='Bg_dim_1'))

    • Bi_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5], dtype='int64', name='Bi_dim_0'))

    • Bi_dim_1
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6], dtype='int64', name='Bi_dim_1'))

    • prem_dim_0
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59],
            dtype='int64', name='prem_dim_0'))

    • grcode_dim_0
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59],
            dtype='int64', name='grcode_dim_0'))

  • Attributes: (4)

    created_at :

    2024-05-27T06:09:05.943907

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  

Even here though we want to average the curves over the specific years in the data and abstract a view of the model implications under different counterfactual settings. Or put another way, we want a view of the average development curve between years, not within a year. If we want to predict the novel insurance loss for next year we need some way to aggregate between the annual trajectories. Here we define a helper function to effect this step.

def get_posterior_predictive_curve(idata, prem=2, grcode=1):
    weighted_splines = [np.dot(np.asfortranarray(Bi), az.extract(idata['posterior']['beta'])['beta'].values[:, :, i]) for i in range(2000)]

    weighted_splines_1 = [np.dot(np.asfortranarray(Bg), az.extract(idata['posterior']['beta_g'])['beta_g'].values[:, i]) for i in range(2000)]

    beta_grcode = az.extract(idata['posterior']['beta_grcode'])['beta_grcode']

    beta_prem = az.extract(idata['posterior']['beta_prem'])['beta_prem']
    df1 = pd.DataFrame([beta_prem.values[i]*prem + beta_grcode.values[i]*grcode for i in range(2000)]).T


    ## How do we averaging over the years to get
    ## a view of a new development period?
    #df = pd.concat([pd.DataFrame(weighted_splines_1[i].T + weighted_splines[i].T).mean() for i in range(1000)], axis=1)

    ## Sample random group each draw from posterior
    df = pd.concat([pd.DataFrame((weighted_splines_1[i].T + weighted_splines[i].T)[np.random.choice(list(range(10))), :]) for i in range(2000)], axis=1)

    ## Average random subset of of groups
    #df = pd.concat([pd.DataFrame((weighted_splines_1[i].T + weighted_splines[i].T)[np.random.choice(list(range(10)), 5), :]).mean() for i in range(2000)], axis=1)

    df = df1.iloc[0].values + df

    return df

pred_df_1 = get_posterior_predictive_curve(idata_43, prem=1, grcode=1)
pred_df_0 = get_posterior_predictive_curve(idata_353, prem=1, grcode=0) 

fig, ax = plt.subplots(figsize=(9, 7), sharey=True)

az.plot_hdi(range(6), pred_df_0.values.T, ax=ax, color='slateblue', smooth=False, fill_kwargs={'alpha': 0.2})
az.plot_hdi(range(6), pred_df_1.values.T, ax=ax, color='firebrick', smooth=False, fill_kwargs={'alpha': 0.2})
az.plot_hdi(range(6), pred_df_0.values.T, ax=ax, color='slateblue', smooth=False, fill_kwargs={'alpha': 0.5}, hdi_prob=.5)
az.plot_hdi(range(6), pred_df_1.values.T, ax=ax, color='firebrick', smooth=False, fill_kwargs={'alpha': 0.5}, hdi_prob=.5)
ax.plot(pred_df_0.mean(axis=1), linestyle='-', color='k', linewidth=4, label='grcode 353 prem 1')

ax.plot(pred_df_1.mean(axis=1), linestyle='--', color='grey', linewidth=4, label='grcode 43 prem 1')

ax.set_title("Posterior Samples of the Trajectories \n Under different Counterfactual settings")
ax.set_ylabel("Loss Ratio")
ax.set_xlabel("Development Period")
ax.legend();

Note how we’ve varied the predictor inputs to showcase the expected realisation under different product codes. In this manner we can construct a hierarchical spline model with informed by the historic trajectories of the development curves but with predictions modulated by more recent information.

Multiple Smoother Regression Models

Non-linear spline relationships can be additively combined across multiple variables in simple and hierarchical regressions. To demonstrate how spline modelling can be further adapted to the multiple regression like cases we use the PISA data set discussed in this case study from Michael Clark. We’ll see how crucial it is to carefully assess the implications of these model fits.

The data set has been constructed using average Science scores by country from the Programme for International Student Assessment (PISA) 2006, along with GNI per capita (Purchasing Power Parity, 2005 dollars), Educational Index, Health Index, and Human Development Index from UN data. We want to model the overall outcome score as a function of these broad demographic features.

pisa_df = pd.read_csv("https://raw.githubusercontent.com/m-clark/generalized-additive-models/master/data/pisasci2006.csv")

pisa_df.head()

	Country	Overall	Issues	Explain	Evidence	Interest	Support	Income	Health	Edu	HDI
0	Albania	NaN	NaN	NaN	NaN	NaN	NaN	0.599	0.886	0.716	0.724
1	Argentina	391.0	395.0	386.0	385.0	567.0	506.0	0.678	0.868	0.786	0.773
2	Australia	527.0	535.0	520.0	531.0	465.0	487.0	0.826	0.965	0.978	0.920
3	Austria	511.0	505.0	516.0	505.0	507.0	515.0	0.835	0.944	0.824	0.866
4	Azerbaijan	382.0	353.0	412.0	344.0	612.0	542.0	0.566	0.780	NaN	NaN

The relationships displayed between each of these measures is not obviously linear, and as such could plausibly benefit from being modelled with splines.

g = sns.pairplot(pisa_df[['Overall', 'Income', 'Support', 'Health', 'Edu']],  kind="reg", height=1.5)
g.fig.suptitle("Pair Plot of Complex Relations", y=1.05);

We define three models for contrasting the implications. Note here how we have to define a seperate spline basis for each of the covariates. Here we create the knots for defining our basis on each covariate.

knots_income = np.linspace(np.min(pisa_df['Income']), np.max(pisa_df['Income']), 5+2)[1:-1]

knots_edu = np.linspace(np.min(pisa_df['Edu']), np.max(pisa_df['Edu']), 5+2)[1:-1]

knots_health = np.linspace(np.min(pisa_df['Health']), np.max(pisa_df['Health']), 5+2)[1:-1]

knots_income1 = np.linspace(np.min(pisa_df['Income']), np.max(pisa_df['Income']), 3+2)[1:-1]

knots_edu1 = np.linspace(np.min(pisa_df['Edu']), np.max(pisa_df['Edu']), 3+2)[1:-1]

knots_health1 = np.linspace(np.min(pisa_df['Health']), np.max(pisa_df['Health']), 3+2)[1:-1]

Now we initialise these models

formula = "Overall ~ Income + Edu + Health"
base_model = bmb.Model(formula, pisa_df, dropna=True)

formula_spline = """Overall ~ bs(Income, degree=3, knots=knots_income) + bs(Edu, degree=3, knots=knots_edu) + bs(Health, degree=3, knots=knots_health) """

formula_spline1 = """Overall ~ (1 | Country) + bs(Income, degree=3, knots=knots_income1) + bs(Edu, degree=3, knots=knots_edu1) + bs(Health, degree=3, knots=knots_health1) """

spline_model = bmb.Model(formula_spline, pisa_df, dropna=True)

spline_model1 = bmb.Model(formula_spline1, pisa_df, dropna=True)

spline_model1

Automatically removing 13/65 rows from the dataset.

Automatically removing 13/65 rows from the dataset.

Automatically removing 13/65 rows from the dataset.

       Formula: Overall ~ (1 | Country) + bs(Income, degree=3, knots=knots_income1) + bs(Edu, degree=3, knots=knots_edu1) + bs(Health, degree=3, knots=knots_health1) 
        Family: gaussian
          Link: mu = identity
  Observations: 52
        Priors: 
    target = mu
        Common-level effects
            Intercept ~ Normal(mu: 471.1538, sigma: 517.9979)
            bs(Income, degree=3, knots=knots_income1) ~ Normal(mu: [0. 0. 0. 0. 0. 0.], sigma: [1630.0781
                872.587   532.485   674.4864  570.5563  779.4582])
            bs(Edu, degree=3, knots=knots_edu1) ~ Normal(mu: [0. 0. 0. 0. 0. 0.], sigma: [965.2606 734.261
                559.679  590.0573 684.2515 649.6795])
            bs(Health, degree=3, knots=knots_health1) ~ Normal(mu: [0. 0. 0. 0. 0. 0.], sigma: [1842.3891
                815.7797  534.7546  671.1905  544.161   802.929 ])
        
        
        Group-level effects
            1|Country ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 517.9979))
        
        Auxiliary parameters
            sigma ~ HalfStudentT(nu: 4.0, sigma: 53.4654)

spline_model.build()
spline_model.graph()

base_idata = base_model.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})
spline_idata = spline_model.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True}, target_accept=.95)
spline_idata1 = spline_model1.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True}, target_accept=.99)

We can compare the simple regression approach to the spline based regression in the now usual way.

compare_df = az.compare({'spline': spline_idata, 'raw': base_idata, 'spline_hierarchy':spline_idata1})

compare_df

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
spline_hierarchy	0	-244.691895	38.487159	0.000000	1.000000e+00	3.901855	0.000000	True	log
spline	1	-254.101306	20.218465	9.409411	0.000000e+00	4.724572	2.840914	True	log
raw	2	-263.061455	9.617888	18.369559	1.054712e-13	10.083921	9.376156	True	log

The coefficients comparisons are harder

az.summary(base_idata)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Intercept	118.438	81.996	-41.166	266.033	1.522	1.086	2915.0	2986.0	1.0
Income	181.955	89.547	17.413	353.598	1.913	1.353	2196.0	2432.0	1.0
Edu	234.125	59.153	120.739	340.343	1.017	0.740	3414.0	2518.0	1.0
Health	30.395	140.966	-223.699	298.258	2.852	2.153	2450.0	2598.0	1.0
Overall_sigma	34.234	3.620	27.882	41.175	0.060	0.043	3757.0	2879.0	1.0

since it’s less clear what the spline coefficient terms mean.

az.summary(spline_idata)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Intercept	225.689	205.515	-133.222	630.021	6.036	4.292	1163.0	1653.0	1.0
bs(Income, degree=3, knots=knots_income)[0]	-355.159	1613.343	-3460.788	2642.535	56.357	39.865	819.0	1504.0	1.0
bs(Income, degree=3, knots=knots_income)[1]	194.905	122.719	-42.478	423.960	2.343	1.728	2748.0	2719.0	1.0
bs(Income, degree=3, knots=knots_income)[2]	48.748	123.638	-182.012	272.863	4.065	2.875	923.0	1431.0	1.0
bs(Income, degree=3, knots=knots_income)[3]	68.533	98.242	-118.858	250.496	3.369	2.383	853.0	1416.0	1.0
bs(Income, degree=3, knots=knots_income)[4]	210.290	103.892	17.347	406.753	3.626	2.565	825.0	1404.0	1.0
bs(Income, degree=3, knots=knots_income)[5]	222.546	105.553	31.163	427.878	3.589	2.539	866.0	1442.0	1.0
bs(Income, degree=3, knots=knots_income)[6]	202.395	105.402	0.453	394.110	3.629	2.567	846.0	1507.0	1.0
bs(Income, degree=3, knots=knots_income)[7]	81.592	103.313	-106.154	279.937	3.641	2.575	805.0	1390.0	1.0
bs(Edu, degree=3, knots=knots_edu)[0]	145.558	287.518	-410.399	654.146	8.525	6.029	1142.0	1509.0	1.0
bs(Edu, degree=3, knots=knots_edu)[1]	47.574	193.029	-323.504	392.406	5.595	3.957	1197.0	1656.0	1.0
bs(Edu, degree=3, knots=knots_edu)[2]	71.806	211.848	-312.889	469.050	6.246	4.417	1156.0	1647.0	1.0
bs(Edu, degree=3, knots=knots_edu)[3]	147.390	200.622	-224.182	518.214	5.972	4.224	1134.0	1603.0	1.0
bs(Edu, degree=3, knots=knots_edu)[4]	70.472	203.119	-314.815	427.662	6.018	4.256	1144.0	1535.0	1.0
bs(Edu, degree=3, knots=knots_edu)[5]	121.297	202.750	-264.164	485.547	5.998	4.243	1147.0	1575.0	1.0
bs(Edu, degree=3, knots=knots_edu)[6]	78.372	206.089	-296.582	457.669	6.012	4.252	1178.0	1713.0	1.0
bs(Edu, degree=3, knots=knots_edu)[7]	124.157	203.070	-258.062	490.727	6.008	4.249	1148.0	1532.0	1.0
bs(Health, degree=3, knots=knots_health)[0]	235.831	885.656	-1490.479	1862.011	31.104	22.002	813.0	1482.0	1.0
bs(Health, degree=3, knots=knots_health)[1]	162.805	389.636	-552.198	913.839	7.756	6.545	2529.0	2422.0	1.0
bs(Health, degree=3, knots=knots_health)[2]	14.967	120.574	-202.130	248.327	3.854	2.726	979.0	1657.0	1.0
bs(Health, degree=3, knots=knots_health)[3]	25.213	100.498	-170.867	204.686	3.411	2.413	871.0	1555.0	1.0
bs(Health, degree=3, knots=knots_health)[4]	-21.315	102.473	-223.005	158.665	3.577	2.530	822.0	1468.0	1.0
bs(Health, degree=3, knots=knots_health)[5]	-19.449	101.418	-213.704	165.809	3.541	2.505	823.0	1444.0	1.0
bs(Health, degree=3, knots=knots_health)[6]	-55.820	100.895	-252.216	126.193	3.547	2.509	811.0	1347.0	1.0
bs(Health, degree=3, knots=knots_health)[7]	0.840	101.892	-196.002	186.216	3.543	2.506	828.0	1469.0	1.0
Overall_sigma	24.766	3.303	18.816	30.884	0.085	0.060	1491.0	2068.0	1.0

We can also check that the model seems to recover the observed data well.

Plotting Posterior Predictive Checks

base_model.predict(base_idata, kind='pps')
spline_model.predict(spline_idata, kind='pps')
spline_model1.predict(spline_idata1, kind='pps')


fig, axs = plt.subplots(1, 3, figsize=(20, 6))
axs = axs.flatten()
az.plot_ppc(base_idata, ax=axs[0])
az.plot_ppc(spline_idata, ax=axs[1]);
az.plot_ppc(spline_idata1, ax=axs[2]);
axs[0].set_xlabel('')
axs[1].set_title("PPC: Spline Model");
axs[0].set_title("PPC: Regression Model");
axs[2].set_title("PPC: Hierarchical Spline Model");

Next we’ll dig into the spline basis features and decompose the predicted outcome and show how the outcome varies for levels in the inputs variables. The Bambi interpret package offers some functionality to assess the conditional predictions for varying values of the input variables.

Model Interpretability Plots in Bambi

We want to highlight the differences between the hierarchical and non-hierarchical multivariable spline models here. The shrinkage effects of hierarchically modelling the country intercepts is clearly evident. They helpfully constrain the poor extrapolation effects of the simpler spline model.

fig, axs = plt.subplots(3, 2, figsize=(9, 14))
axs = axs.flatten()
bmb.interpret.plot_predictions(
    spline_model, spline_idata, "Income", ax=axs[0]
)

axs[0].set_title("Non-Hierarchical Income")
axs[2].set_title("Non-Hierarchical Edu")
axs[4].set_title("Non-Hierarchical Health")

axs[1].set_title("Hierarchical Income")
axs[3].set_title("Hierarchical Edu")
axs[5].set_title("Hierarchical Health")


bmb.interpret.plot_predictions(
    spline_model, spline_idata, "Edu", ax=axs[2]
)

bmb.interpret.plot_predictions(
    spline_model, spline_idata, "Health", ax=axs[4]
);

bmb.interpret.plot_predictions(
    spline_model1, spline_idata1, "Income", 
    sample_new_groups=True, ax=axs[1]
)

bmb.interpret.plot_predictions(
    spline_model1, spline_idata1, "Edu", sample_new_groups=True, ax=axs[3]
)

bmb.interpret.plot_predictions(
    spline_model1, spline_idata1, "Health",sample_new_groups=True, ax=axs[5]
);

We can pull these types of values out into a table. Note here how we ask for predictions based on varying values of the Edu and Income variables where we keep the Health variable fixed at the mean value.

summary_df = bmb.interpret.predictions(
    spline_model,
    spline_idata,
    conditional={
        "Edu": list(np.linspace(0.5, .95, 10)),
        "Income": list(np.linspace(0.4, .95, 10)),
        },
)
summary_df


summary_df1 = bmb.interpret.predictions(
    spline_model1,
    spline_idata1,
    conditional={
        "Edu": list(np.linspace(0.5, .95, 10)),
        "Income": list(np.linspace(0.4, .95, 10)),
        },
    sample_new_groups=True
)
summary_df1

	Edu	Income	Country	Health	estimate	lower_3.0%	upper_97.0%
0	0.50	0.400000	Albania	0.885672	102.686176	-587.761967	750.373783
1	0.50	0.461111	Albania	0.885672	451.070276	-37.688251	905.442110
2	0.50	0.522222	Albania	0.885672	376.852378	-161.138067	945.914437
3	0.50	0.583333	Albania	0.885672	283.859330	-357.146695	933.659941
4	0.50	0.644444	Albania	0.885672	288.608418	-345.493652	920.643943
...	...	...	...	...	...	...	...
95	0.95	0.705556	Albania	0.885672	460.351293	413.197965	508.270287
96	0.95	0.766667	Albania	0.885672	516.885593	474.706704	561.845324
97	0.95	0.827778	Albania	0.885672	552.375825	506.285465	603.102537
98	0.95	0.888889	Albania	0.885672	499.921356	455.457677	550.896029
99	0.95	0.950000	Albania	0.885672	229.655207	49.386349	412.317948

100 rows × 7 columns

We can then plot these results based on the conditional realisations and see some interesting behaviour at the implausilble reaslisations. This should make us somewhat wary of using this model to extrapolate too far beyond the observable range of data

g = sns.relplot(data=summary_df, x="Income", y="estimate", hue="Edu")

g.fig.suptitle("Marginal Predictions of the Outcome Variable \n conditional on counterfactual values for Edu and Income", y=1.05);


g = sns.relplot(data=summary_df1, x="Income", y="estimate", hue="Edu")

g.fig.suptitle("Hierarchical Marginal Predictions of the Outcome Variable \n conditional on counterfactual values for Edu and Income", y=1.05);

However, we can see here how the hierarchical component shrinks the predicted values back towards a reasonable range! Constraining the poor extrapolation of the more naive spline model. This demonstrates something of the interplay between in-sample approximation and out of sample generalisation we saw above in the case of the insurance curve development. The additional structure insists on commensurate realisations under different counterfactual settings reflective of the fact that countries do not vary so radically.

The Spline Component Contributions

Finally, we’ll pull out the component contributions of each variable and see how they combine additively. This will also serve as a kind of posterior predictive check for each country as we can show the degree which posterior draws from each component sum to achieve a plausible mirror of the observed data.

Bincome = spline_model.response_component.design.common['bs(Income, degree=3, knots=knots_income)']

income_coefs = az.extract(spline_idata['posterior']['bs(Income, degree=3, knots=knots_income)'])['bs(Income, degree=3, knots=knots_income)']

Bedu = spline_model.response_component.design.common['bs(Edu, degree=3, knots=knots_edu)']

edu_coefs = az.extract(spline_idata['posterior']['bs(Edu, degree=3, knots=knots_edu)'])['bs(Edu, degree=3, knots=knots_edu)']


Bhealth = spline_model.response_component.design.common['bs(Health, degree=3, knots=knots_health)']

health_coefs = az.extract(spline_idata['posterior']['bs(Health, degree=3, knots=knots_health)'])['bs(Health, degree=3, knots=knots_health)']

income = np.dot(Bincome, income_coefs).T 
edu = np.dot(Bedu, edu_coefs).T
health = np.dot(Bhealth, health_coefs).T

intercept = az.extract(spline_idata['posterior']['Intercept'])['Intercept'].values

fig, ax = plt.subplots(figsize=(10, 7))
for i in range(100):
    if i == 1:
        ax.plot(income[i], label='Income Component', color='red')
        ax.plot(edu[i], label='Edu Component', color='blue')
        ax.plot(health[i], label='Health Component', color='darkgreen')
        ax.plot(intercept[i] + income[i] + edu[i] + health[i], label='Combined Components', color='purple')
    else: 
        ax.plot(income[i], alpha=0.1, color='red')
        ax.plot(edu[i], alpha=0.1, color='blue')
        ax.plot(health[i], alpha=0.1, color='darkgreen')
        ax.plot(intercept[i] + income[i] + edu[i] + health[i], color='purple', alpha=0.3)

ax.scatter(range(len(spline_idata['observed_data']['Overall'])), spline_idata['observed_data']['Overall'], label='Observed', color='grey', s=56, ec='black')
ax.set_title("PISA Outcomes \n Additive Spline Components", fontsize=20)
ax.legend();
ax.set_ylabel("Overall Score", fontsize=12)
ax.set_xticklabels(pisa_df.dropna(axis=0).reset_index()['Country'], fontsize=12);

We can see here how the combined components borrow the structure of the outcome variable primarily from the income variable component. The health measures have closer to zero additive contribution while the uncertainty in the educational component varies wildly. But the blue educational component here is used primarily as a scaling contribution which adds to the level of the outcome variable rather than distorting the shape. These initially opaque synthetic features of splines can offer deep insight into the structure of our data generating process when seen correctly.

Conclusion

We’ve seen the application of splines as univariate smoothers to approximate wiggly curves of arbitrary shape. We’ve also tried gaussian process approximations of the same univariate functions. The suggested flexibility of both methods is a strength, but we need to be careful where splines have a tendency to over-fit to individual curves. As such we have tried to show that we can incorporate spline basis modelling in a hierarchical bayesian model and recover more compelling posterior predictive checks and additionally derive predictiions from the mixed variant of the hierarhical model which helps us understand the implications of the data for generic forecasts of insurance loss curves. Finally we showed how splines can be used additively to model non-linear functions of multiple covariates. These are a powerful tool to interrogate complex non-linear relationships, but they offer interpolation of functions within a well understood domain. Applealing to these models for extrapolation needs to be done carefully.

We can, I think, draw a broad moral from this presentation. Observed phenomena admit description with a wide variety of mathematical abstraction. The “correctness” of the abstraction is never the right question. The manner in which one abstraction distils the structure of reality is not fundamentally any better than the manner in which an alternative scheme will do so. Predictive success achieved applying one abstraction can be traded off for interpretive tractability in another. What matters is the task at hand - solving the problem. Sometimes problems can be solved with new information, but mostly they’re resolved by re-arranging assumptions and phrasing the problem anew - seeing through a different lens, pulling the door handle where we previously pushed. Splines, hierarhical splines and gaussian process models all attempt to approximate an answer by rephrasing the question. The more of these tools we have, the better we are prepared to tackle novel problems and discern the pertinent structure in the world.

Citation

BibTeX citation:

@online{forde2024,
  author = {Nathaniel Forde},
  title = {GAMs and {GPs:} {Flexibility} and {Calibration}},
  date = {2024-04-07},
  langid = {en}
}

For attribution, please cite this work as:

Nathaniel Forde. 2024. “GAMs and GPs: Flexibility and Calibration.” April 7, 2024.