Verifying that Monte Carlo posterior moments on scalar model parameters extend to the bivariate (and multivariate) parameter case. Obvious in retrospect but I am VERY VERY BAD at math.

bivariateMonteCarloPosterior.py

"""
We know how to use Monte Carlo to find moments of the posterior
```
P(parameter | data) ∝ P(data | parameter) * P(parameter)
```
That is, `posterior ∝ likelihood * prior`.

We generate parameter draws, and for each draw, assign a weight equal to
`P(data | parameters)`, since we have data. Then the weighted mean, weighted
variance, etc. give us the moments of the posterior.

This works for a scalar `parameter` but I want to know, how does it work
for multiple `parameters`?

Test via https://en.wikipedia.org/wiki/Normal-gamma_distribution which is conjugate
with a Normal with unknown mean and precision (tau = 1/sigma^2).
"""
from scipy.stats import norm, gamma  #type:ignore
import matplotlib.pylab as plt  # type:ignore
import numpy as np  #type:ignore


def normRv(mean, var, N):
  return norm.rvs(size=N, loc=mean, scale=np.sqrt(var))


(mu, lam, alpha, beta) = 1.0, 0.5, 10 * 1.4 + 1, 10.
N = 1_000_000
tau = gamma.rvs(alpha, scale=1 / beta, size=N)
mean = normRv(mu, 1 / (lam * tau), N)

Ndata = 5
data = normRv(2.0, 2.0, Ndata)

exampleFig, exampleAxs = plt.subplots(3, 1)
exampleAxs[0].hist(mean, bins=50, density=True)
exampleAxs[0].set_xlabel('mean')
exampleAxs[1].hist(tau, bins=50, density=True)
exampleAxs[1].set_xlabel('τ=1/σ^2')
exampleAxs[2].hist(data, bins=50, density=True)
exampleAxs[2].set_xlabel('data')

exampleFig.tight_layout()

sampleMean = np.mean(data)
sampleVar = np.var(data)
muNew = (lam * mu + Ndata * sampleMean) / (lam + Ndata)
lamNew = lam + Ndata
alphaNew = alpha + Ndata / 2
betaNew = beta + 0.5 * (Ndata * sampleVar + (lam * Ndata * (sampleMean - mu)**2) / (lam + Ndata))


def normalGammaMeans(mu, _lam, alpha, beta):
  return dict(mean=mu, tau=alpha / beta)


weights = np.exp(
    np.sum(np.array([norm.logpdf(x, loc=mean, scale=np.sqrt(1 / tau)) for x in data]), axis=0))
weightedMean = lambda w, x: np.sum(w * x) / np.sum(w)
print(normalGammaMeans(muNew, lamNew, alphaNew, betaNew))
print(dict(mcMean=weightedMean(weights, mean), mcTau=weightedMean(weights, tau)))
"""
Output:
```
{'mean': 2.822881709045586, 'tau': 1.4014492663032703}
{'mcMean': 2.8231642987861028, 'mcTau': 1.4018106759909776}
```
The analytical means for the mean and precision are very similar to the Monte Carlo ones.

Answer: it works exactly the same as with scalars. If you have data that's drawn from two
parameters,
- generate Monte Carlo draws from the two,
- compute the likelihood of the data for each pair of random samples: that's your weights,
- and compute weighted moments of the samples.
"""