GitHub profile

Step 5 - Posterior Distribution of Parameters

Notation

Having confirmed in Step 4 that the chains converged, this step looks at what the posterior actually says about $\mu_i$ and $\sigma_i$ for each component — a full distribution (density plus a highest density interval), in place of the single fixed number classical PCA gives for each.

Each panel plots the marginal posterior for one parameter, with the classical point estimate overlaid as a reference line: $\bar{z}_i$ (the column mean of that component's scores, from Step 9) for $\mu_i$, and $\sqrt{\lambda_i}$ (the square root of that component's eigenvalue, from Step 8) for $\sigma_i$. Where the classical value falls within the bulk of the posterior, the uninformative-prior Bayesian model is reproducing classical PCA, as expected — this is the same check as Step 4's cross-check, just shown parameter-by-parameter rather than as a single number.

step 1 - plot the posterior distributions


%run ./bayesian/step5_posterior.py
posterior_plot_chart(trace, "output/snippets_png/bayesian_step5_posterior.png")

Python

filename: step5_posterior.py

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
from dataframes.step9_scores import scores
from dataframes.step7_eigendecomposition import eigenvalues

'''
posterior distribution plots for mu_i and sigma_i per component (Step 3's trace), with the
classical point estimate (column mean, sqrt(eigenvalue) from Step 8) overlaid as a reference
line - in place of classical's single fixed number per parameter.
'''

def classical_reference_values(components=("PC1", "PC2", "PC3")):
    df = scores()
    eigvals = eigenvalues().set_index("component")["eigenvalue"]
    ref = {}
    for component in components:
        ref[f"mu_{component}"] = df[component].mean()
        ref[f"sigma_{component}"] = np.sqrt(eigvals[component])
    return ref

def posterior_plot_chart(trace, filename, components=("PC1", "PC2", "PC3")):
    var_names = [f"mu_{c}" for c in components] + [f"sigma_{c}" for c in components]
    ref_val = classical_reference_values(components)
    az.plot_posterior(trace, var_names=var_names, ref_val=ref_val)
    plt.savefig(filename)
    return
posterior density and HDI for mu_i and sigma_i per component, with the classical point estimate overlaid as a reference line