How can I include sample size information in my Bayesian inference model?

Question

I'm new to Bayesian inference, and have been following the book Bayesian Methods for Hackers. In this chapter, the author provides an example of how using small sample sizes can lead to extreme results, but doesn't say how to correct for this when modelling the problem.

This is how I'm currently modeling the problem with pymc. average_across_county is an array of the average person's height per county. There is also an array called population (unused here) which contains the population of each county.

mu  = pm.Uniform('height', 100, 200)
tau = pm.Uniform('precision', 0, 50)
hgt_dist = pm.Normal('hgt_dist', mu, tau, value=average_across_county, observed=True)

mcmc = pm.MCMC([mu, tau, hgt_dist])
mcmc.sample(50000, 2000)

mcmc.trace('height')[:].mean()

The error in above estimate of the true average height is about 10 times larger than a simple weighted average because I don't make use of the population data. How should I include the population variable in this model?

Answer 1

I think I should add some explanation. This example was to show estimates with small samples sizes can have extreme values. I think if I was to "model" this example, I would model it as follows:

For each county, there is an true $mu$, but a priori we don't know this. From other information, I can suggest that the mean is between 100 and 200. I have sample data from each county, and I want to use this to derive a posterior per county, and compare these posteriors together (perhaps against the global data).

county_samples = [
 [150., 160., 170., 150.],
 [140., 160., 170., 180., 150., 130.],
 [160., 140., 170.],
]

mus = [pm.Uniform('mu_%d' % i, 100, 200) for i in range(3)]
taus = [pm.Uniform('tau_%d' % i, 0, 50) for i in range(3)]
obs = [pm.Normal('obs_%d' % i, mus[i], taus[i], observed=True, value=county_samples[i]) for i in range(3)]

mcmc = pm.MCMC([mus, taus, obs])
mcmc.sample(50000,10000)


mu_0 = mcmc.trace('mu_0')[:]
mu_1 = mcmc.trace('mu_1')[:]
mu_2 = mcmc.trace('mu_2')[:]

plt.hist(mu_0, bins=50, histtype='stepfilled', label='mu_0')
plt.hist(mu_1, bins=50, histtype='stepfilled', label='mu_1', alpha=0.5)
plt.hist(mu_2, bins=50, histtype='stepfilled', label='mu_2', alpha=0.5)

So, even though the data has different sample sizes, we can't rule out that they are statistically different from each other. (Try playing around with the data to try to make that untrue.)

Your example uses the precomputed averages per country. This is not usually how I model in PyMC - I prefer to keep things granular, and as you can see, it fits the framework nicely.

Does this help?