# Trouble replicating simple example of Bayesian inference

On pages 20-21 of John Kruschke's Doing Bayesian Data Analysis book (2nd ed.), there is an introductory illustration of Bayesian inference. We know that balls can have four sizes: 1, 2, 3 and 4, but the manufacturing process is not perfect, so the empirically observed balls are of sizes 1.77, 2.23 and 2.7.

Now, the assumptions are:

• the three observed balls were produced as balls of the same nominal size,
• the priors in this example are 0.25 for each ball type,
• the distribution of ball size variability is normal and centered on ball size for each ball type.

We would like to measure which nominal size is most probable for this sample: 1, 2, 3 or maybe 4.

At this stage of the book, this problem is only used to show how posterior changes when data is observed, and no calculations are provided. The Author does offer precisely stated result, though:

...there is 56% probability that the balls are size 2, 31% probability that the balls are size 3, 11% probability that the balls are size 1, and only 2% probability that the balls are size 4.

I tried to replicate this result and I failed. My logic was as follows:

I have: $$\theta \in \{1, 2, 3, 4\}$$ and: $$data \in \{1.77, 2.23, 2.7\}$$.

Given Bayes' formula:

$$p(\theta|data) = \frac{p(data|\theta) \times p(\theta)}{p(data)} = \frac{p(data|\theta) \times p(\theta)}{\displaystyle\sum_{i=1}^{4} p(data|\theta_i) \times p(\theta_i)}$$

Since $$p(\theta) = 0.25$$ for all values, I can factor it out in the denominator and then cancel, so I am left with:

$$p(\theta|data) = \frac{p(data|\theta)}{\displaystyle\sum_{i=1}^{4} p(data|\theta_i) }$$

For numerator, I have to count the following:

$$p(data|\theta) = \mathcal{N}(data_1|\theta,1) \times \mathcal{N}(data_2|\theta,1) \times \mathcal{N}(data_3|\theta,1)$$

For example:

$$p(1.77, 2.23, 2.7|2) = \mathcal{N}(1.77|2,1) \times \mathcal{N}(2.23|2,1) \times \mathcal{N}(2.7|2,1)$$

This translates to Python as:

import numpy as np
from scipy.stats import norm
mean_2 = np.prod(norm.pdf([1.77, 2.23, 2.7], loc=2, scale=1))


So, the full calculation of $$p(\theta|data)$$ should look like:

mean_1 = np.prod(norm.pdf([1.77, 2.23, 2.7], loc=1, scale=1))
mean_2 = np.prod(norm.pdf([1.77, 2.23, 2.7], loc=2, scale=1))
mean_3 = np.prod(norm.pdf([1.77, 2.23, 2.7], loc=3, scale=1))
mean_4 = np.prod(norm.pdf([1.77, 2.23, 2.7], loc=4, scale=1))
posterior_2 = mean_2 / (mean_1 + mean_2 + mean_3 + mean_4)


My assumption is that posterior_2 should have a value close to 56% mentioned by John Kruschke, but it's actually close to 64% (after multiplying by 100).

I see two reasons for this discrepancy:

• there is some fault in my logic at this early stage of learning Bayesian data analysis (do you see any?),
• there's not enough information in the book to replicate the results.

The latter is possible because the Author does not state the value of standard deviation for these normal distributions. I am assuming it's 1.

• I think, you are doing everything correctly. Not knowing the variance seems to be the problem. Assuming $1$ is quite big from a modelling point of view; $1$ is also the distance between two nominal sizes. – Jonas Jan 12 at 22:39
• Found by trial-and-error: if you change scale=1.16, the numbers agree. Direct link to the example: books.google.pl/… – Tim Jan 12 at 22:56