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:
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))
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. =)