# Bayesian analysis from summary statistics in PyMC

I'm trying to reproduce this example from Confidence intervals vs Bayesian intervals by ET Jaynes: Assuming dispersion is standard deviation, how does one model that in PyMC? If I had the 31 + 61 samples, I'd just set them to the observed values of a generative model such as the one explained by Kruschke on his BEST procedure.

One idea I had was to create a deterministic function that would calculate the "dispersion" (standard deviation) of each group and set the two values to be observed. However, I don't know how to use the information about the 31 and 61 samples in this case. Ignoring sample size seems wrong to me, so I'm stuck.

Is it possible to model something like that in PyMC?

• Can someone explain why I was downvoted? – Pedro Tabacof Feb 2 '15 at 16:55

I think that it is possible to model this in PyMC, but in my skimming of the Jaynes article, it does not seem to have all of the detail necessary. It does include the answer

Maybe you can reverse engineer the model. Here is my guess at a start for you:

import pymc as pm

@pm.stochastic
def R1(value=1500):
return 1./value

@pm.stochastic
def R2(value=1500):
return 1./value

@pm.observed
def d1(value=2237, R=R1, nu=31):
x = nu * value/R
return pm.chi2_like(x, nu)

@pm.observed
def d2(value=1347, R=R2, nu=61):
x = nu * value/R
return pm.chi2_like(x, nu)

m = pm.MCMC([R1, R2, d1, d2])
m.sample(iter=200000, burn=100000, thin=10)

• Thank you very much for the answer and the code, I took some time to understand it but now it makes sense, there are just two things I don't understand: Why are you using custom priors for R1 and R2? Also, shouldn't the degrees of freedom used be subtracted by one since they are from the sample size? – Pedro Tabacof Feb 3 '15 at 16:55
• I thought the custom priors for R1 and R2 might help if you want to figure out what Jaynes would do. Jeffries, perhaps? That is not implemented in PyMC, AFAIK. Regarding the degrees of freedom, you are probably right, and that will make my answer a little closer to Jaynes value. Here is a notebook that makes that change and gets a results. – Abraham D Flaxman Feb 3 '15 at 18:25