Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
You've seen this? –  J. M. Oct 3 '12 at 19:05
Great resource! It also has several integration strategies (Integrate, MCMC, Gibbs, Interpolation) –  Eric Brown Oct 4 '12 at 1:32
This is a notebook from the 2007 Wolfram Technology Conference demonstrating many Bayesian statistical procedures. –  gwr Nov 13 '12 at 15:14
add comment

1 Answer

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
share|improve this answer
Thanks for posting, Asim. I think I understand that a Beta Distribution would be the posterior since Beta is conjugate with Binomial. I think I'm looking for some more pedagogical, symbolic expressions, e.g. showing Bayes' rule in action. –  Eric Brown Oct 4 '12 at 4:26
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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