Normalization Factor Wrong? (Bug?)

Question

I'm new to PyMC and Bayesian stuff in general, so I started off with what I thought was a very simple toy problem. I generated some normally-distributed noise with a given mean and standard deviation, and then played with PyMC and looked at the posterior distribution of possible mu/sigma.

It's possible that I just don't understand what I'm looking at, but shouldn't the standard deviation of the trace for mu converge on the standard error of the mean? It's currently off my a factor of sqrt(2).

Simply: Parent normal distribution with parameters mu0, sigma0. Child population sampled from parent distribution, resulting in mu1, sigma1.

SE = sigma1/sqrt(N) SD = sigma0/sqrt(N)

Shouldn't SE~SD~mcmc.trace['mu'].std() for sufficiently large N and # of iterations? What I actually see is:

SE~SD~mcmc.trace['mu'].std()/sqrt(2)

This is using PyMC 2.2 and here is the code:

%pylab inline
import pymc as mc

true_noise_mu = 50
true_noise_sigma = 5
true_noise_datapoints = 2000

noise_field = np.random.normal(true_noise_mu, true_noise_sigma, true_noise_datapoints)


data_mean = noise_field.mean()
data_std = noise_field.std()


mu = mc.Normal( 'mu', data_mean, 1/(data_std/sqrt(n_datapoints))**2)

# just guessed at tau here, but it shouldn't impact the result (nor does it), #only the time to converge, yes?
sigma = mc.Normal( 'sigma', 1/data_std**2, 0.01)


observation = mc.Normal( "obs", mu, sigma, value=noise_field, observed=True)

mcmc = mc.MCMC(model)
mcmc.sample( 400000, 20000, 1 )

mu_samples = mcmc.trace('mu')[:]
sigma_samples = mcmc.trace('sigma')[:]

print "SEM estimated from sample population: %.3f" % (noise_field.std()/sqrt(true_noise_datapoints))
print "SEM calculated from parent population: %.3f " % (true_noise_sigma / sqrt(true_noise_datapoints))
print "SEM calculated from posterior distribution: %.3f" % (mu_samples.std())
print "Ratio of posterior: sample estimate: %.3f" % ((mu_samples.std())/(noise_field.std()/sqrt(true_noise_datapoints)))
print "1/sqrt(2): %.3f" % (1/sqrt(2))

Last Print:

SEM estimated from sample population: 0.677
SEM calculated from parent population: 0.671 
SEM calculated from posterior distribution: 0.478
Ratio of posterior: sample estimate: 0.707
1/sqrt(2): 0.707

I'm probably just an idiot, but I would love to understand what I'm missing. =)

Answer 1

Try removing the sqrt(n_sata_points) from the variance of the mu prior. The results are that the ratio in question is (approximately) 1. This suggests the error is in that term.

Assume for the moment that the true std. deviation is known, call it $\sigma$. Also, define $\hat{\mu}$ as the empirical mean. You have assigned the prior

$$ \mu \sim \text{Normal}\left( \hat{\mu}, \frac{ \sigma^2}{N} \right)$$

to the unknown $\mu$.

The likelihood is Normal (Normally distributed data) as well. Hence, the posterior will also be Normal (by conjugancy, see here). The standard deviation of the posterior Normal, according to the wiki page, is:

$$\left( \frac{1}{ \frac{ \sigma^2 }{N} } + \frac{ N }{\sigma^2} \right)^{-1}$$

$$ = \frac{\sigma^2}{2N}$$

which I believe is where your sqrt 2 is coming from. (of course, this still assumes we know $\sigma$, but for large $N$ our estimate $\hat{\sigma} \approx \sigma$)

---

There is something I don't like about using $\frac{\sigma^2}{N}$ as the prior variance. It feels like you are double-counting your sample size (hence the additional factor of 2 cropping up).

I would suggest using $\hat{\sigma}$ or even better is a less subjective constant like 100 or something.