# Sampling from predictive distribution with unknown mean and variance

I am following Pattern Recognition and Machine Learning and trying to implement Bayesian linear regression with unknown mean and variance, so that the posterior is given by a Normal-Gamma distribution $NG(w,\beta\mid w_0,\Lambda^{-1},a,b)$

I update the precision matrix by $\Lambda = \Lambda_0 + X^TX$, the weights by $w = \Lambda^{-1}(\Lambda_0w_0 + X^Ty)$, the shape parameter by $a = a_0 + \frac{n}{2}$, and the rate parameter by $b = b_0 + \frac{1}{2}\sum_i{(y_i - w^Tx_i)^2}$

After fitting the parameters to a simple polynomial curve, I then use the parameters to sample from the predictive distribution given by a Student's t-distribution in Python like so, where the variance according to the book should be equal to $(\beta\Lambda)^{-1} = (\frac{a}{b}\Lambda)^{-1}$

x1 = np.linspace(left, right, 250)
chol = np.linalg.cholesky(np.linalg.inv(posterior.L * (posterior.a / posterior.b)))
for k in range(100):
w_s = np.dot(chol, np.random.standard_t(2. * posterior.a, size=len(posterior.w))) + posterior.w
y1 = w_s[3] * x1**3 + w_s[2] * x1**2 + w_s[1] * x1 + w_s[0]
plt.plot(x1, y1, color='red', alpha=0.2)


However, the sampled lines do not seem right as they are far too narrow given the data as you can see here (blue dots are the data that is being fit to, the green line is the maximum a posteriori, and the red lines are the samples)

For some reason I had the intuition to remove the $a$ from the variance term, so that it is instead given by $(\frac{1}{b}\Lambda)^{-1}$, like so

chol = np.linalg.cholesky(np.linalg.inv(posterior.L / posterior.b))


After doing that, the samples actually seem to make sense

But I'm at a loss as to why. I'm sure it's something obvious, but my brain doesn't seem to want to help me. Can anyone point it out to me? Thanks.

I reproduced the results, using the wikipedia page at https://en.wikipedia.org/wiki/Bayesian_linear_regression as a guide, instead of Pattern Recognition, and I got the exact same results, eg for 25 samples:

The full jupyter notebook (mostly using torch, rather than numpy, just a bit of numpy for the gaussian/IG sampling), is at https://github.com/hughperkins/pub-prototyping/blob/f1398e2244aa74da914534e26a3c164823d52f5c/maths/bayesian_linear_regression.ipynb

I hypothesized that it's not that the results are wrong/incorrect, but that given these data, the most probable lines are actually fairly well centred between the points. It would seem more probable that the mean line is fairly central, with just highish variance, than that the line passes through some of the outliers, perhaps?

To test this hypothesis, I tried with just 7 results, and the lines were more varied, and looked fairly plausible to me:

So, I'm going to say that, given we both independently worked through the formulae, from independent sources of formulae, and writing the code entirely from scratch, and get the same results, that:

• the results are in fact correct
• given the number of points we have in your example, and my first example, the most likely mean lines are in fact fairly well centred between the points
• if you want to see lines that vary more, then you might need to reduce the number of points perhaps? As per my second example

tldr: your code and results seem correct, as far as I can see

• Thanks for taking the time to reproduce the results from scratch. I think it makes sense to me now - for some reason I had it in my head that the predictive distribution for w would be influenced more heavily by the variance of y, but after reading your answer and thinking about it more that wasn't a correct assumption. – rojmor Oct 8 '17 at 22:14