GitHub profile

What are we modeling?

Not settled yet. The open question isn't just where the Bayesian prior sits — it's what generative model of the data we're assuming in the first place. Two candidates are being developed side by side rather than picking one up front.

Option A — score-level uncertainty ("Bolt-On")

Keep the classical loadings/eigenvectors ($\mathbf{W}$) fixed, computed once by classical eigendecomposition.
$\mathbf{W}$ never enters the Bayesian model at all.

The likelihood target is the already-reduced 1-D component scores, one Normal per component, and the prior sits on that distribution's parameters, $(\mu_i, \sigma_i^2)$.

The Bayesian part is pure uncertainty quantification bolted onto a point estimate that classical PCA already produced — hence "Bolt-On". See Bayesian PCA (Bolt-On).

Option B — structure-level uncertainty (latent)

$\mathbf{W}$ is not fixed. The likelihood target is the raw demeaned yield matrix $\widetilde{\mathbf{Y}}_c$ itself, via $\widetilde{\mathbf{Y}}_c \sim \text{Normal}(\mathbf{0}, \mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I})$ — that is, a Normal whose covariance matrix is $\mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}$ rather than something computed separately — in the spirit of Tipping & Bishop's Probabilistic PCA and Bishop's Bayesian PCA. Here, the model is unsure about the components themselves — not just unsure about today's score on an already-agreed component.

Note there's no separate "compute the covariance matrix" step in this formulation. The covariance matrix an analyst would normally compute and eigendecompose never appears as an input here. Instead, $\mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}$ is worked out fresh from whatever $\mathbf{W}$ and $\sigma^2$ happen to be at each step of the sampler — it's an output of the model, not something fed into it. Doing the covariance/eigendecomposition step first and feeding the result in as fixed loadings is Option A, not this one.

The parameters this model actually infers are $\mathbf{W}$ and $\sigma^2$ — not the scores. Bishop's Bayesian PCA additionally puts a hierarchical prior on the columns of $\mathbf{W}$ (an ARD hyperparameter $\alpha_i$ per column, itself given a hyperprior), which lets the model automatically shrink away components the data doesn't need instead of the analyst fixing the number of components in advance. The 1-D scores can still be recovered afterward as a derived quantity — a closed-form conditional distribution given a posterior draw of $\mathbf{W}$ and $\sigma^2$ — but they are downstream of the model, not its target.

See Bayesian PCA (Latent).

Why this is more than a prior-placement choice

These are two different generative models of the data, not the same model with the prior moved around. That matters for the actual aim (systematic inclusion of expert judgement about future dynamics not in the historical data):

Both tracks are being kept alive in parallel until there's enough from each to compare and settle on one.