BUGS-type calculations in Mathematica

Question

Recently I have been teaching myself how to Bayesian calculations with the BUGS language (JAGS, in particular). However, I find myself wondering how one might use Mathematica to do similar calculations (perhaps based around NIntegrate). I would be most interested in building up models using the awesome combined computing/graphing facilities that Mathematica provides. I have a copy of "Bayesian Logical Data Analysis for the Physical Sciences" and while it is an excellent and enjoyable book, I think that it more focused on "analytic" solutions than "numerical"--it seems that more complicated models would quickly outpace the ability to obtain analytic solutions.

For example, consider the problem of inferring the difference between two rates, taken from Lee and Wagenmaker's book:

model {
# Prior on Rates
theta1 ~ dbeta(1,1)
theta2 ~ dbeta(1,1)

# Observed Counts
k1 ~ dbin(theta1,n1)
k2 ~ dbin(theta2,n2)

# Difference between Rates
delta <- theta1 - theta2
}

Given k1, n1, k2, and n2, how might I translate this into Mathematica to get, e.g. the posterior density of theta, the difference in rates, I suppose p(delta | k1,k2,n1,n2) ? (For example, k1=5, k2=7, n1=n2=10). Knowing how this is done would go a long way toward me being able to build and solve other models!

Answer 1

The posterior for theta1 is beta(6,6) and the posterior for theta2 is beta(8,4). Thus the posterior for delta can be approximated by taking draws from these two distributions and subtracting these.

theta1Draws = RandomVariate[BetaDistribution[6, 6], 10000];

theta2Draws = RandomVariate[BetaDistribution[8, 4], 10000];

deltaDraws = theta1Draws - theta2Draws;

You can calculate the mean and the 95% posterior intervals for delta from the sample of draws stored in deltaDraws

In[7]:= Mean[deltaDraws]

Out[7]= -0.168601

In[10]:= Quantile[deltaDraws, 0.025]

Out[10]= -0.532054

In[9]:= Quantile[deltaDraws, 0.975]

Out[9]= 0.213398