The experiment: Counting red M&M’s

Let’s start this practical session by making inference about the proportion of red M&M’s in a very huge box. If you are following this practice in person, the instructor will give you a bag of M&M’s to make the practice more fun even if we have to be a bit flexible about the theoretical conditions of the experiment. If you are working on line it will be a bit more boring but we will comply better with the theoretical conditions of the experiment.

Representation of the sampling experiment.

For those of us who physically have the bag of M&M’s, the most appropriate experiment would be to take candies out of the bag one at a time and replace them, recording whether their colour is red or not, until the number n of candies defined in advance is reached. We consider n=20.

In case you are on line or you don’t want to carry out the experiment physically, you can use the following result of the experiment:

mm_sample <- data.frame(n = 20, r = 4) 

The sampling model is Binomial

Let Y be the random variable that describes the number of red M&M’s (successes) out of the n = 20 sampled. Given the proportion θ of red candies in the bag, Y is a binomial random variable with parameters n = 20 and probability θ associated with drawing a red candy.

Y ∣ θ ∼ Bi(n, θ)

whose conditional probability function, expectation and variance are:

• $P(Y=r \mid \theta)=\binom{n}{r}\,\theta^r\, (1-\theta)^{n-r}, \,\,r=0,1,\ldots, n$
• E(Y ∣ θ) = n θ,
• Var(Y ∣ θ) = n θ (1 − θ)

Binomial probability function for two values of θ and n = 20 observations.

A prior distribution for θ

Recall that the beta distribution is conjugate with respect to the binomial probability model. If we elicit a prior beta distribution Be(α₀, β₀) for θ, its density is

$$\pi(\theta \mid \alpha_0, \beta_0)=\frac{\Gamma(\alpha_0+\beta_0)}{\Gamma(\alpha_0)\,\Gamma(\beta_0)}\, \theta^{\alpha_0-1}\,(1-\theta)^{\beta_0-1}, \quad 0<\theta<1$$ with expectation and variance

• 𝔼(θ ∣ α₀, β₀) = α₀/(α₀ + β₀),
• Var(θ ∣ α₀, β₀) = α₀ β₀/[(α₀ + β₀)² (α₀ + β₀ + 1)]

We will start by working with the non-informative prior distribution Be(α₀ = 0.5, β₀ = 0.5) whose prior expectation and standard deviation is 𝔼(θ) = 0.5 and 𝕊𝔻(θ) = 0.3535, respectively. Its graphic is obtained via

curve(
  dbeta(x, 0.5, 0.5),
  col = "dodgerblue", lwd = 4,
  ylim = c(0, 4),
  xlab = expression(paste("Proportion of red candies ", theta)),
  ylab = 'prior'
)

The likelihood function of θ

The likelihood function of θ is a function of θ whose expression requires the observed data 𝒟 = {n, r}. It is defined as follows

$$L(\theta \mid \mathcal D) = P(Y=r \mid \theta, n) = \binom{n}{r}\, \theta^r\,(1- \theta)^{n-r},$$

where r is the observed number of successes in the sample (here, the number of red candies) of size n. In our first experiment, the likelihood function of θ for the data 𝒟₁ = {r₁ = 4 red candies in n₁ = 20 sampled} is

$$L(\theta \mid \mathcal D_1) = \binom{20}{4}\, \theta^4\,(1- \theta)^{16},$$

r <- mm_sample$r
n <- mm_sample$n
curve(
  choose(n, r) * x^r * (1 - x)^(n - r),
  col = "darkorange", lwd = 4,
  xlab = expression(paste("Proportion of red candies ", theta)),
  ylab = 'likelihood'
)

The posterior distribution of θ

The posterior distribution of θ is also a beta distribution with parameters

π(θ ∣ 𝒟) = Be(α = α₀ + r, β = β₀ + n − r).

In our case, the posterior distribution of the proportion of red candies is

π(θ ∣ 𝒟₁) = Be(α = 4.5, β = 16.5) whose graphic is

curve(
  dbeta(x, mm_sample$r + 0.5, mm_sample$n - mm_sample$r + 0.5),
  col = "darkgreen", lwd = 4,
  xlim = c(0, 1), ylim = c(0, 5), 
  xlab = expression(paste("Proportion of red candies ", theta)),
  ylab = 'posterior'
)

The posterior mean and posterior standard deviation of the proportion of red candies is

$$\mathbb{E}(\theta \mid \mathcal D)=\frac{4.5}{4.5+16.5}=0.2143, \quad\quad \mathbb{SD}(\theta \mid \mathcal D)=\sqrt{\frac{4.5 \cdot 16.5}{(4.5+16.5)^2 \cdot (4.5+16.5+1)}}=0.0875.$$ A 95% credible interval for θ is

qbeta(c(0.025, 0.975), 4.5, 16.5)
#> [1] 0.07152005 0.40822576

and the posterior probability that the proportion of red candies is between 0.1 and 0.3 is

pbeta(0.3, 4.5, 16.5) - pbeta(0.1, 4.5, 16.5)
#> [1] 0.7571733

You can see these posterior summaries in the Figure below

Posterior summaries for θ.

The posterior predictive distribution for the results of a new experiment

Recall that the posterior predictive distribution of the number of successes Y^′ in a new sample of size n^′ is a beta-binomial distribution Bb(α, β, n^′) with parameters α and β from the posterior distribution π(θ ∣ 𝒟) and the size n^′ of the new experiment. Its probability function is

$$P(Y^{\prime} = r^{\prime} \mid \mathcal{D}) = \binom{n^\prime}{r^\prime}\, \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \,\Gamma(\beta)} \,\frac{\Gamma(\alpha+r^{\prime})\,\Gamma(\beta+n^{\prime}-r^{\prime})}{\Gamma(\alpha+\beta+n^{\prime})}, \,\, r^{\prime}=0,1,\ldots, n^{\prime}$$

with expectation and variance,

• $\mathbb{E}(Y^{\prime}\mid \mathcal{D})=n^{\prime}\,\frac{\alpha}{(\alpha+\beta)},$
• $\mathbb{V}(Y^{\prime}\mid \mathcal{D})=n^{\prime}\,\frac{\alpha\,\beta\,(\alpha+\beta+n^{\prime})}{(\alpha+\beta)^{2}\,(\alpha+\beta+1)}.$

Suppose now that we are going to randomly draw n^′ = 10 new M&M’s from the bag. Before performing the experiment, the posterior predictive distributions quantifies the relative plausibility of each possible outcome. It is a beta-binomial distribution with parameters α = 4.5 and β = 16.5 from the posterior distribution π(θ ∣ 𝒟), and n^′ = 10, which is the size of the new experiment.

library(extraDistr)
curve(
  dbbinom(x, size = 10, mm_sample$r + 0.5, mm_sample$n - mm_sample$r + 0.5),
  col = "purple", lwd = 4, type = 'h', n = 11,
  xlim = c(0, 10), ylim = c(0, 0.4), 
  xlab = expression(paste("Number of red candies")),
  ylab = 'probability'
)

We observe that the most likely numbers of red candies in the new sample are 1 and 2. The posterior mean and posterior standard deviation of the number of red candies is this new sample is

𝔼(Y^′ ∣ 𝒟) = 2.14,    𝕊𝔻(Y^′ ∣ 𝒟) = 1.54.