Bayesian Principal Component Analysis (Bolt-On)

This is one of two candidate approaches to adding a Bayesian layer, kept side by side rather than settled on — see What are we modeling? for why, and for the other option.

Contents

• Steps
  • Data prep (reused from Classical PCA)
  • Bayesian-specific steps
• Notes

Steps

We go through the following steps to produce a Bayesian principal component analysis.

Data prep (reused from Classical PCA)

Data prep through deriving a stress (loading the raw data, truncating, logging, differencing, de-meaning, covariance, eigendecomposition, variance explained, projecting to scores, fitting/testing the score distribution, and reversing out the stress) is unchanged from Classical PCA — the loadings/eigenvectors are kept fixed, so the Bayesian model consumes the same component scores as its data rather than recomputing them. That pipeline is still growing on the classical side, so rather than continuing its step numbers here (and having to renumber every time a classical step is added), the Bayesian-specific steps below use their own numbering, starting at 1.

Bayesian-specific steps

1. Specify the prior — a non-informative prior on $(\mu_i, \sigma_i^2)$ for each component's score distribution: a flat prior on $\mu_i$, a Jeffreys prior $p(\sigma_i^2) \propto 1/\sigma_i^2$ on the variance (see Priors)
2. Build the probabilistic model — a Normal likelihood for each component's scores combined with the Step 1 priors, defined in PyMC
3. Sample the posterior (MCMC) — run the sampler to draw from the joint posterior of $(\mu_i, \sigma_i^2)$ for each component
4. Convergence diagnostics — check the sampler actually converged before trusting the draws: trace plots, $\hat{R}$, effective sample size (arviz)
5. Posterior distribution of parameters — the fitted posterior for $\mu_i, \sigma_i^2$ per component, in place of classical's single point estimate for each
6. Posterior predictive distribution — for each posterior draw of $(\mu_i, \sigma_i)$ from Step 3, simulate a score $\sim \text{Normal}(\mu_i, \sigma_i)$; the resulting spread of simulated scores, across all the draws, is the score-level distribution that folds in parameter uncertainty (the closed-form conjugate route would instead give this analytically as a Student-t)
7. Compare to classical output — stresses/percentiles derived from the posterior predictive vs. the classical normal fit. This is aim 3 from the Aims page: with an uninformative prior, the two should come out close, which acts as a sanity check on the Bayesian machinery before any expert judgement is introduced in a later, informative-prior pass

Notes

!!! in pca the model is a model of 3 principal components.... and we fit this to a distributions.... so the prior must be a parameter of that distribution.

!! the challenge as i see it !!

• we need to be able to formulate priors that indicate our measure of uncertainty
• a simple adjustment might be to say we expect rates to be lower at all terms by a constant amount (a downward shift)
• each principal component however is difficult to interpret because it is a vector in many directions
• so interpreting what each parameter of principal component model is non intuitive

• and yet bayesian framework IS setting priors for the parameters of the principal components

• it is possible to reverse engineer component scores that give our expected shift

• which in turn we could perhaps use as a mean value to set parameters
• but do we actually want to produce the set of possible parameters... in some ways no !!
• perhaps we can use it to answer question how likely expected outcome given historical data...
• but even then the different parameter values that give something close to expected outcome -- difficult to anticpate combinations
• and we havent even discussed combinations... does this need to be a joint probability distribution.

i think we should keep eigenvectors set at what they are...
so the challenge is reverse engineering parameter adjustments ....