This post relates to my original question here, but this time focusing on a more fundamental misunderstanding of how MCMC actually works.
When using Bayesian MCMC for parameter inference with a single model, it is usual to ignore the "evidence" integral and write:
$$P(\theta | D) \propto P(D | \theta) P(\theta)$$
For numerical reasons, it is also common to work with logarithms:
$$ln[P(\theta | D)] \propto ln[P(D | \theta)] + ln[P(\theta)]$$
I've used this successfully on a number of occasions, typically using Python and emcee. However, for reasons explained in my other question, I've recently been wondering about "normalising" my likelihood by dividing it by the number of observations (I realise this probably sounds odd, but please see my other question for why I thought it was worth a try).
My likelihood assumes i.i.d. Gaussian errors, so my log-likelihood is just the sum of the individual data likelihoods. My prior is uniform and I'm currently using an improper log-prior function, which simply returns zero if the parameter values are "within range" and '-inf'
otherwise.
I initially assumed that dividing by the number of data points wouldn't make any difference, because the RHS is still proportional to the log-posterior:
$$ln[P(\theta | D)] \propto \frac{1}{n}ln[P(D | \theta)] + ln[P(\theta)]$$
(My prior only ever returns 0 or '-inf'
, so multiplying that by $1/n$ makes no difference).
However, when I re-run my MCMC, I get very different results. As a sense check, I tried a simple OLS linear regression using the same code.
With my original (not normalised) likelihood, I get sensible results: the MCMC recovers the "true" parameters and the 95% CI broadly agrees with the Frequentist results from 'statsmodels'
:
However, with the "normalised" likelihood, I get very different (and obviously wrong) results:
I think I understand why this is, at least from the point of view of a simple Metropolis algorithm. The probability of accepting the proposed jump, $x_{i+1}$, compared to the current location, $x_i$, is:
$$\alpha = \frac{P(x_{i+1})}{P(x_i)}$$
Which, re-writing in terms of my log-posterior, $L$, becomes:
$$\alpha = \frac{exp(L(x_{i+1}))}{exp(L(x_i))} = exp(L(x_{i+1}) - L(x_i))$$
On the other hand, for my "normalised" version, all the likelihoods are divided by $n$, so this becomes:
$$\alpha = [exp(L(x_{i+1}) - L(x_i))]^\frac{1}{n}$$
which is different, and might explain why my uncertainty bounds suddenly become much larger. (I understand that emcee is more sophisticated than just basic Metropolis, but I guess a similar principle applies?).
My question:
What am I missing/misunderstanding to reconcile these two results, please? On the one hand, it's OK to ignore the "normalising constant" in Bayes' equation, because anything proportional to the log-posterior is sufficient. On the other hand, multiplying the log-posterior by a constant $\frac{1}{n}$ seems to have a dramatic effect on my results?
Thank you!