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))

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

| cite | improve this question | | | | |
  • $\begingroup$ An quick detail is that you are using a Normal to model the standard deviation. This us a poor choice as the Normal can return negative values (which a standard deviation forbids, and running your code throws a ZeroProbability: Stochastic obs's value is outside its support error.) A better choice is a Uniform, i.e.sigma = mc.Uniform("sigma", 0, 50)**(-2). $\endgroup$ – Cam.Davidson.Pilon Apr 19 '13 at 19:35
  • 1
    $\begingroup$ EDIT: actually, you are interested in the evaluation of Uniform(..), so I would code it as sigma= Uniform("sigma", 0, 50); precision = sigma**(-2) and use precision in the observation variable. $\endgroup$ – Cam.Davidson.Pilon Apr 19 '13 at 19:40

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.

| cite | improve this answer | | | | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.