So far in our chapters, we have created quite a few likelihood functions. When we introduced likelihoods at the very beginning1, we were careful to state that likelihoods are not probability functions. However, we know two things:
Because the likelihood function is non-negative and the area under this curve is finite, we can turn it into a distribution by a process called finding the marginal likelihood2 and dividing the likelihood function by this marginal likelihood.
More formally, consider a probability function \(f(x|\boldsymbol{\theta})\) with some finite but unknown parameter vector \(\boldsymbol{\theta} \in\mathbb{R}_p\). Further, consider an independent and identical sample, \(\textbf{x}\), of size \(n\) from this distribution, represented as \(x_i \overset{iid}{\sim} f(x|\boldsymbol\theta),\ i = 1, 2, \ldots, n\). The resulting likelihood is \[ \mathcal{L}(\boldsymbol\theta|\textbf{x}) = \prod_{i = 1}^n f(x_i|\boldsymbol\theta). \] In order to transform \(\mathcal{L}\) into a probability function, we must divide by the integral of \(\mathcal{L}\) over the support3 (all possible values) of \(\boldsymbol\theta\), where \(\mathcal{S}(\boldsymbol\theta)\) represents this support. For some distributions, especially those with a mixture of discrete and continuous parameters4, this support, \(\mathcal{S}(\boldsymbol\theta)\), may be quite complex.
Given this setup, the marginal likelihood is defined as \[ m(\boldsymbol\theta|\textbf{x}) = \int_{\mathcal{S}(\boldsymbol\theta)} \mathcal{L}(\boldsymbol\theta|\textbf{x}) \pi(\boldsymbol\theta) d\boldsymbol\theta. \] This \(\pi(\boldsymbol\theta)\) may come as a surprise, as we haven’t defined what it is. In Bayesian Statistics5, this \(\pi(\boldsymbol\theta)\) is known as a prior distribution6. In the context of Bayesian inference, this prior distribution represents all the expert knowledge about the unknown parameter vector \(\boldsymbol\theta\) that was known before the sample \(\textbf{x}\) was collected. However, for this class, we will make the statement that we don’t know much, if anything, about \(\boldsymbol\theta\), so we set \(\pi(\boldsymbol\theta) = 1\). Then, for our examples, \[ m(\boldsymbol\theta|\textbf{x}) = \int_{\mathcal{S}(\boldsymbol\theta)} \mathcal{L}(\boldsymbol\theta|\textbf{x}) d\boldsymbol\theta. \]
Finally, in order to transform the likelihood function of \(\boldsymbol\theta\) into a probability function of \(\boldsymbol\theta\), we will divide the likelihood by its integral. So, the probability function of \(\boldsymbol\theta\) given the observed data \(\textbf{x}\) is \[ f(\boldsymbol\theta|\textbf{x}) = \frac{\mathcal{L}(\boldsymbol\theta|\textbf{x})}{m(\boldsymbol\theta|\textbf{x})} = \frac{\mathcal{L}(\boldsymbol\theta|\textbf{x})}{\int_{\mathcal{S}(\boldsymbol\theta)} \mathcal{L}(\boldsymbol\theta|\textbf{x}) d\boldsymbol\theta}. \] We comment that distributions derived this way are true statistical distributions because they will always be non-negative, and an integral divided by itself equals 1 (so the total probability will be 1 by definition). In practice, this integral may be impossible to solve, so often numerical routines are used to estimate \(f\) directly. Such topics are beyond the scope of this course.
One small piece of theory we will need below is to swap the order of integration and summation (which we have done quite loosely up to this point). Basically, if we have two properties: 1) that \(f \ge 0\), and \(F < \infty\), then \[ \int \sum_n f_n(x) dx = \sum_n \int f_n(x) dx. \] This is an extension of Fubini’s Theorem7 which allows us to swap the order of summation and integration, as long as the function \(f\) is non-negative and the integral converges. Note that if \(f\) is a probability function of a statistical distribution, then we have both properties automatically.
We will begin this process by considering an independent and identical sample \(\textbf{x}\), with size \(n\), from an Binomial Distribution (with \(p\) unknown); that is \(x_i \overset{iid}{\sim} \text{Binom}(N, p),\ i = 1, 2, \ldots, n\). Further, we will take a play from the Negative Binomial distribution, and let \(k_i\) and \(r_i\) denote the number of successes and failures in Binomial sample \(i\), respectively. Thus, \[ \begin{align} \mathcal{L}(p|\textbf{x}) &= \prod_{i = 1}^n {N \choose x_i} p^{x_i} (1 - p)^{N - x_i} \\ \Longrightarrow \mathcal{L}(p|\textbf{r},\textbf{k}) &= \prod_{i = 1}^n {r_i + k_i \choose k_i} p^{k_i} (1 - p)^{r_i} \\ &= \left[ \prod_{i = 1}^n {r_i + k_i \choose k_i} \right] \times \left[ \prod_{i = 1}^n p^{k_i} \right] \times \left[ \prod_{i = 1}^n (1 - p)^{r_i} \right] \\ &= \left[ \prod_{i = 1}^n {r_i + k_i \choose k_i} \right] p^{S_k} (1 - p)^{S_r}, \end{align} \] where \[ S_k = \sum_{i = 1}^n k_i,\ \text{and}\ S_r = \sum_{i = 1}^n r_i. \]
As we discussed in our Formal Foundations section, we can marginalize this likelihood to create a probability function \(f\), of the parameter \(p\), given the sufficient statistics8 of the observed data \(S_k\) and \(S_r\). That is, \[ \begin{align} f(p|S_r, S_k) &= \frac{\mathcal{L}(p|\textbf{r},\textbf{k})}{m(p|\textbf{r},\textbf{k})} \\ &= \frac{\mathcal{L}(p|\textbf{r},\textbf{k})}{\int_{\mathcal{S}(p)} \mathcal{L}(p|\textbf{r},\textbf{k}) dp } \\ &= \frac{ \left[ \prod_{i = 1}^n {r_i + k_i \choose k_i} \right] p^{S_k} (1 - p)^{S_r} }{\int_0^1 \left[ \prod_{i = 1}^n {r_i + k_i \choose k_i} \right] p^{S_k} (1 - p)^{S_r} dp } \\ &= \frac{ \left[ \prod_{i = 1}^n {r_i + k_i \choose k_i} \right] p^{S_k} (1 - p)^{S_r} }{ \left[ \prod_{i = 1}^n {r_i + k_i \choose k_i} \right] \int_0^1 p^{S_k} (1 - p)^{S_r} dp } \\ &= \frac{ p^{S_k} (1 - p)^{S_r} }{ \int_0^1 p^{S_k} (1 - p)^{S_r} dp } \\ &= \frac{ p^{(S_k + 1) - 1} (1 - p)^{(S_r + 1) - 1} }{ \int_0^1 p^{(S_k + 1) - 1} (1 - p)^{(S_r + 1) - 1} dp } \\ &= \frac{ p^{\alpha - 1} (1 - p)^{\beta - 1} }{ \int_0^1 p^{\alpha - 1} (1 - p)^{\beta - 1} dp }, \end{align} \] where \(\alpha = S_k + 1\) (1 plus the total number of successes in the \(n\) Binomial trials) and \(\beta = S_r + 1\) (1 plus the total number of failures in the \(n\) Binomial trials).
This integral in the denominator should look familiar: it is the definition of the Complete Beta Function, which we covered in the “Formal Foundations” chapter on the Gamma and Beta functions. Thus, \[ \begin{align} f(p|S_r, S_k) &= \frac{ p^{\alpha - 1} (1 - p)^{\beta - 1} }{ \int_0^1 p^{\alpha - 1} (1 - p)^{\beta - 1} dp } \\ &= \left[ \int_0^1 p^{\alpha - 1} (1 - p)^{\beta - 1} dp \right]^{-1} p^{\alpha - 1} (1 - p)^{\beta - 1} \\ &= \left[ \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)} \right]^{-1} p^{\alpha - 1} (1 - p)^{\beta - 1} \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} p^{\alpha - 1} (1 - p)^{\beta - 1}, \end{align} \] which is the standard form of the Beta Distribution with \(p\in(0,1)\).
What are the allowed values for \(\alpha\) and \(\beta\)? Because we derived the Beta Distribution as the probability function of the Binomial Distribution’s parameter \(p\), we had an added restriction that \(S_k\) and \(S_r\) were non-negative integers (that is, \(S_k,\ S_r \in 0\cup\mathbb{N} = 0, 1, 2, 3, \ldots\)). This would impose the restriction that \(\alpha = S_k + 1\) and \(\beta = S_r + 1\) must be elements of \(\mathbb{N} = 1, 2, 3, \ldots\) (not including 0). However, notice that the form of \(f\) above uses Gamma functions, so it does not require that \(\alpha\) and \(\beta\) be restricted to the integers. What would it mean to have fractional/decimal counts of successes or failures? Well, some games allow for ties, which could be counted as half a success and half a failure. Other experiments could involve Likert-scale9 responses, where “strongly disagree” maps to 0, “strongly agree” maps to 1, but the values in between map to various fractions between 0 and 1.10 Thus, we state that \(\alpha\) and \(\beta\) simply need to be non-negative real numbers (\(\alpha,\beta\in\mathbb{R}^+\)).
set.seed(20150516)
xSymm <- rbeta(n = 500, shape1 = 10, shape2 = 10)
samplesSymm_ls <- list(
n10 = xSymm[1:10],
n30 = xSymm[1:30],
n60 = xSymm[1:60],
n500 = xSymm
)
xSkew <- rbeta(n = 500, shape1 = 5, shape2 = 1.5)
samplesSkew_ls <- list(
n10 = xSkew[1:10],
n30 = xSkew[1:30],
n60 = xSkew[1:60],
n500 = xSkew
)
range_num <- c(0, 1)
rm(xSymm, xSkew)PlotSharedDensity <- function(x, range_x, bandwidth = "nrd0") {
xDens_ls <- density(x, bw = bandwidth)
xHist_ls <- hist(x, plot = FALSE)
yLargest_num <- max(max(xDens_ls$y), max(xHist_ls$density))
hist(
x, prob = TRUE,
xlim = range_x, ylim = c(0, yLargest_num)
)
lines(xDens_ls, col = 4, lwd = 2)
}par(mfrow = c(2, 2))
PlotSharedDensity(
x = samplesSymm_ls$n10, range_x = range_num
)
PlotSharedDensity(
x = samplesSymm_ls$n30, range_x = range_num
)
PlotSharedDensity(
x = samplesSymm_ls$n60, range_x = range_num
)
PlotSharedDensity(
x = samplesSymm_ls$n500, range_x = range_num
)
par(mfrow = c(1, 1))
# , bandwidth = 0.005
par(mfrow = c(2, 2))
PlotSharedDensity(
x = samplesSkew_ls$n10, range_x = range_num
)
PlotSharedDensity(
x = samplesSkew_ls$n30, range_x = range_num
)
PlotSharedDensity(
x = samplesSkew_ls$n60, range_x = range_num
)
PlotSharedDensity(
x = samplesSkew_ls$n500, range_x = range_num
)
par(mfrow = c(1, 1))
Given all our work to derive this distribution, we can show that it is a distribution directly. First, note that because \(p\in(0,1)\) and \(\alpha,\beta > 0\), we have that \(f(p|\alpha,\beta) > 0\). Then, starting with the Riemann-Stieltjes integral and applying the definition of the Complete Beta Function11 we have \[ \begin{align} \int_{\mathcal{S}(p)} dF(p|\alpha,\beta) &= \int_0^1 \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} p^{\alpha - 1} (1 - p)^{\beta - 1} dp \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \int_0^1 p^{\alpha - 1} (1 - p)^{\beta - 1} dp \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \left[ \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)} \right] \\ &= 1. \end{align} \] Therefore, the Beta Distribution is a proper distribution.
When I was in grad school, some of my professors ignored deriving the Moment Generating Function of the Beta Distribution. When you look at its form in the back of a statistical theory textbook, you can probably see why (it’s unpleasant). However, we are going to derive it anyway. Many derivations you might find online involve using the Kummer’s Function of the First Kind12, which requires even more theoretical foundations to cover than I would ever want to write. Instead, we will opt for a longer derivation, but one that uses Formal Foundations that we have already covered, namely the Riemann-Stieljes Integral (that you should be comfortable with by now), the MacLaurin Series13 of \(e^x\), swapping the order of integration and summation via Fubini’s Theorem,14 the definition of the Complete Beta Function,15 and the Continued Recurrence Property of the Gamma Function.16
Let’s begin: \[ \begin{aligned} M_p(t) &= \int_{\mathcal{S}(p)} e^{tp}dF(p|\alpha,\beta) \\ &= \int_0^1 e^{tp} \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} p^{\alpha - 1} (1 - p)^{\beta - 1} dp \\ &\qquad\text{\emph{MacLaurin Series of }} e^{tp} \text{\emph{...}} \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \int_0^1 \left[ \sum_{k = 0}^{\infty} \frac{(tp)^k}{k!} \right] p^{\alpha - 1} (1 - p)^{\beta - 1} dp \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \int_0^1 \sum_{k = 0}^{\infty} \frac{t^k}{k!} p^k p^{\alpha - 1} (1 - p)^{\beta - 1} dp \\ &\qquad\text{\emph{Fubini's Theorem...}} \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \sum_{k = 0}^{\infty} \int_0^1 \frac{t^k}{k!} p^k p^{\alpha - 1} (1 - p)^{\beta - 1} dp \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \sum_{k = 0}^{\infty} \frac{t^k}{k!} \int_0^1 p^{\alpha + k - 1} (1 - p)^{\beta - 1} dp \\ &\qquad\text{\emph{Defn. of Beta Function...}} \\ &= \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \sum_{k = 0}^{\infty} \frac{t^k}{k!} \left[ \frac{\Gamma(\alpha + k)\Gamma(\beta)}{\Gamma(\alpha + \beta + k)} \right] \\ &= \sum_{k = 0}^{\infty} \frac{t^k}{k!} \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha + \beta + k)} \frac{\Gamma(\alpha + k)}{\Gamma(\alpha)} \right] \\ &\qquad\text{\emph{``Peel off'' the first summand...}} \\ &= \frac{t^0}{0!} \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha + \beta + 0)} \frac{\Gamma(\alpha + 0)}{\Gamma(\alpha)} \right] + \sum_{k = 1}^{\infty} \frac{t^k}{k!} \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha + \beta + k)} \frac{\Gamma(\alpha + k)}{\Gamma(\alpha)} \right] \\ &\qquad\text{\emph{Gamma Function Recurrence Property...}} \\ &= 1 + \sum_{k = 1}^{\infty} \frac{t^k}{k!} \left[ \frac{\Gamma(\alpha + \beta)}{ \prod_{j = 0}^{k - 1} (\alpha + \beta + j) \Gamma(\alpha + \beta) } \frac{ \prod_{j = 0}^{k - 1} (\alpha + j) \Gamma(\alpha) }{\Gamma(\alpha)} \right] \\ &= 1 + \sum_{k = 1}^{\infty} \frac{t^k}{k!} \left[ \frac{1}{ \prod_{j = 0}^{k - 1} (\alpha + \beta + j) } \frac{ \prod_{j = 0}^{k - 1} (\alpha + j) }{1} \right] \\ &= 1 + \sum_{k = 1}^{\infty} \frac{t^k}{k!} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right], \end{aligned} \] which is the standard form of the Beta Distribution’s MGF. We “peeled off” the first summand (incrementing from \(k = 0\) to \(k = 1\)) to make taking the first derivative easier. When we are ready to take the second derivative, we will “peel off” the second summand (incrementing from \(k = 1\) to \(k = 2\)) then.
Following our Binomial example with probability of success \(p_0 = 0.35\), let’s create a probability distribution from which this \(p_0\) is the most likely value. For the Beta Distribution, the formula for the mode17 is known, so we set it equal to 0.35, and we will have a frontier18 where one unknown variable can be written as a function of another: \[ \begin{aligned} 0.35 &= \frac{\alpha - 1}{\alpha + \beta - 2} \\ \Longrightarrow 0.35\alpha + 0.35\beta - 0.7 &= \alpha - 1 \\ \Longrightarrow 0.35\beta &= 0.65\alpha - 0.3 \\ \Longrightarrow \beta(\alpha) = \frac{13}{7}\alpha + \frac{6}{7}. \end{aligned} \] Notice that if we pick a value for one parameter of the Beta Distribution, while holding the mode constant, then the other value is determined.
Let’s pick a range of values for \(\alpha\), calculate the corresponding values of \(\beta\), and then look at the distribution of values as \(\alpha\) and \(\beta\) increase.
alpha_num <- 2^seq(from = 0, to = 6, length.out = 6)
beta_num <- (13/7) * alpha_num + (6/7)
set.seed(20150516)
par(mfrow = c(3, 2))
for (i in seq_along(alpha_num)) {
hist(
x = rbeta(n = 10000, shape1 = alpha_num[i], beta_num[i]),
xlim = c(0, 1),
main = paste0(
"alpha = ", round(alpha_num[i], 2),
"; beta = ", round(beta_num[i], 2)
),
xlab = NULL
)
abline(v = 0.35, col = "red", lwd = 2)
}
par(mfrow = c(1,1))Notice that for small values of \(\alpha\), the mode of the distribution is not the desired 0.35. However, as \(\alpha\) grows larger, the mode both gets closer to 0.35 and the distribution gets tighter around the desired value of 0.35.
The intuition for this goes back to the motivating derivation of the Beta Distribution. The parameters \(\alpha\) and \(\beta\) represent the total number of successes and failures, respectfully, in repeated independent Binomial experiments. The more total successes and failures we’ve observed, the more information we have about the true value of \(p\). Therefore, the distribution must shrink and become more narrow around the intended value of \(p_0 = 0.35\). In practice, larger parameters of the Beta Distribution represent more prior information known about the probability of success, \(p\).
Let’s generate some random data. Before we do this, we should think critically about this exercise. When, if ever, is it truly possible to “observe” a probability? We can observe successes and failures, but in the health science context, I can’t think of any case where we could actually take a random sample of probabilities themselves. However, this somewhat nonsensical step is necessary to move forward with our examples for the Method of Moments estimators and also to describe the forms of the solutions to the MLE exercises. So, we will press onward.
We will assume that we already know that in prior independent Binomial experiments, we observed a total of 64 successes and 120 failures. Let’s assume that we somehow “observed” \(n = 12\) probabilities from a Beta distribution with \(\alpha = 64 + 1\) and \(\beta = 120 + 1\). This would yield a theoretical mode value of \((65 - 1)/(65 + 121 - 2) \approx 0.348\). We can sample these values as shown below:
set.seed(20150516)
alphaParam_int <- 65
betaParam_int <- 121
(p_num <- rbeta(n = 12, shape1 = alphaParam_int, shape2 = betaParam_int))
[1] 0.396 0.355 0.337 0.341 0.290 0.334 0.344 0.376 0.365 0.345 0.299 0.324Thus, the 12 “observed” probabilities are 0.396, 0.355, 0.337, 0.341, 0.29, 0.334, 0.344, 0.376, 0.365, 0.345, 0.299, 0.324.
Let’s get to work: \[ \begin{aligned} M_p(t) &= 1 + \sum_{k = 1}^{\infty} \frac{t^k}{k!} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \\ \Longrightarrow M^{\prime}_p(t) &= 0 + \frac{\partial}{\partial t} \sum_{k = 1}^{\infty} \frac{t^k}{k!} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \\ &= \sum_{k = 1}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{\partial}{\partial t} \frac{t^k}{k!} \\ &= \sum_{k = 1}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{kt^{k - 1}}{k!} \\ &= \sum_{k = 1}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 1}}{(k - 1)!} \\ &\qquad\text{\emph{``Peel off''}}\ k = 1\ldots \\ &= \left\{ \left[ \prod_{j = 0}^{[1] - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{[1] - 1}}{([1] - 1)!} \right\} + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 1}}{(k - 1)!} \\ &= \left\{ \left[ \frac{\alpha + 0}{\alpha + \beta + 0} \right] \frac{t^0}{0!} \right\} + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 1}}{(k - 1)!} \\ &= \frac{\alpha}{\alpha + \beta} + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 1}}{(k - 1)!} \\ \Longrightarrow M^{\prime}_p(0) &= \frac{\alpha}{\alpha + \beta} + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{[0]^{k - 1}}{(k - 1)!} \\ &= \frac{\alpha}{\alpha + \beta} + 0 \\ &= \mathbb{E}[p]. \end{aligned} \]
Similarly, \[ \begin{aligned} M^{\prime}_p(t) &= \frac{\alpha}{\alpha + \beta} + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 1}}{(k - 1)!} \\ \Longrightarrow M^{\prime\prime}_p(t) &= \frac{\partial}{\partial t} \left\{ \frac{\alpha}{\alpha + \beta} + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 1}}{(k - 1)!} \right\} \\ &= 0 + \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{\partial}{\partial t} \frac{t^{k - 1}}{(k - 1)!} \\ &= \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{(k - 1)t^{k - 2}}{(k - 1)!} \\ &= \sum_{k = 2}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 2}}{(k - 2)!} \\ &\qquad\text{\emph{``Peel off''}}\ k = 2\ldots \\ &= \left\{ \left[ \prod_{j = 0}^{[2] - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{[2] - 2}}{([2] - 2)!} \right\} + \sum_{k = 3}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 2}}{(k - 2)!} \\ &= \left\{ \left[ \prod_{j = 0}^1 \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^0}{0!} \right\} + \sum_{k = 3}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 2}}{(k - 2)!} \\ &= \left\{ \left[ \frac{\alpha + 0}{\alpha + \beta + 0} \times \frac{\alpha + 1}{\alpha + \beta + 1} \right] [1] \right\} + \sum_{k = 3}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 2}}{(k - 2)!} \\ &= \left[ \frac{\alpha}{\alpha + \beta} \frac{\alpha + 1}{\alpha + \beta + 1} \right] + \sum_{k = 3}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{t^{k - 2}}{(k - 2)!} \\ \Longrightarrow M^{\prime\prime}_p(0) &= \left[ \frac{\alpha}{\alpha + \beta} \frac{\alpha + 1}{\alpha + \beta + 1} \right] + \sum_{k = 3}^{\infty} \left[ \prod_{j = 0}^{k - 1} \frac{\alpha + j}{\alpha + \beta + j} \right] \frac{[0]^{k - 2}}{(k - 2)!} \\ &= \left[ \frac{\alpha}{\alpha + \beta} \frac{\alpha + 1}{\alpha + \beta + 1} \right] + 0 \\ &= \mathbb{E}[p^2]. \end{aligned} \]
Thus, \[ \begin{aligned} \text{Var}[p] &= \mathbb{E}[p^2] - \left[\mathbb{E}[p]\right]^2 \\ &= \left[ \frac{\alpha}{\alpha + \beta} \frac{\alpha + 1}{\alpha + \beta + 1} \right] - \left[ \frac{\alpha}{\alpha + \beta} \right]^2 \\ &= \frac{\alpha(\alpha + 1)}{(\alpha + \beta)(\alpha + \beta + 1)}\frac{(\alpha + \beta)}{(\alpha + \beta)} - \frac{\alpha^2}{(\alpha + \beta)^2}\frac{(\alpha + \beta + 1)}{(\alpha + \beta + 1)} \\ &= \frac{\alpha(\alpha + 1)(\alpha + \beta) - \alpha^2(\alpha + \beta + 1)}{(\alpha + \beta)^2(\alpha + \beta + 1)} \\ &= \frac{\alpha(\alpha^2 + \alpha\beta + \alpha + \beta) - (\alpha^3 + \alpha^2\beta + \alpha^2)}{(\alpha + \beta)^2(\alpha + \beta + 1)} \\ &= \frac{(\alpha^3 + \alpha^2\beta + \alpha^2) + \alpha\beta - (\alpha^3 + \alpha^2\beta + \alpha^2)}{(\alpha + \beta)^2(\alpha + \beta + 1)} \\ &= \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)} \\ \end{aligned} \]
After these derivations, we have the following system of equations: \[ \bar{p} = \frac{\alpha}{\alpha + \beta};\ s^2 = \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)}. \] Let’s solve the first equation for \(\beta\), since it only appears once: \[ \begin{aligned} \bar{p} &= \frac{\alpha}{\alpha + \beta} \\ \Longrightarrow \alpha\bar{p} + \beta\bar{p} &= \alpha \\ \Longrightarrow \beta\bar{p} &= \alpha - \alpha\bar{p} \\ \Longrightarrow \beta &= \frac{\alpha}{\bar{p}} - \alpha. \end{aligned} \]
We will now substitute this value for \(\beta\), which depends on the known value of \(\bar{p}\) and the unknown value of \(\alpha\), into the second equation: \[ \begin{aligned} s^2 &= \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)} \\ &= \frac{\alpha\left[ \frac{\alpha}{\bar{p}} - \alpha \right]}{\left(\alpha + \left[ \frac{\alpha}{\bar{p}} - \alpha \right]\right)^2 \left(\alpha + \left[ \frac{\alpha}{\bar{p}} - \alpha \right] + 1\right)} \\ &= \frac{ \alpha^2\left[ \frac{1}{\bar{p}} - 1 \right] }{\left( \frac{\alpha}{\bar{p}} \right)^2 \left( \frac{\alpha}{\bar{p}} + 1 \right)} \\ &= \frac{ \alpha^2\left[ \frac{1}{\bar{p}} - 1 \right] }{ \alpha^2 \left( \frac{\alpha}{\bar{p}^3} + \frac{1}{\bar{p}^2} \right)} \times \frac{\bar{p}^3}{\bar{p}^3} \\ &= \frac{ \left[ 1 - \bar{p} \right] \bar{p}^2 }{ \frac{\alpha\bar{p}^3}{\bar{p}^3} + \frac{\bar{p}^3}{\bar{p}^2} } \\ &= \frac{ (1 - \bar{p})\bar{p}^2 }{\alpha + \bar{p}} \\ \Longrightarrow s^2(\alpha + \bar{p}) &= (1 - \bar{p})\bar{p}^2 \\ \Longrightarrow s^2\alpha &= (1 - \bar{p})\bar{p}^2 - s^2\bar{p} \\ \Longrightarrow \hat{\alpha} &= (1 - \bar{p})\frac{\bar{p}^2}{s^2} - \bar{p}. \end{aligned} \]
Finally, we can substitute this estimate for \(\alpha\), which is entirely in terms of the known quantities \(\bar{p}\) and \(s^2\), back into the first equation that we solved for \(\beta\). Thus, \[ \begin{aligned} \beta &= \frac{\alpha}{\bar{p}} - \alpha \\ &= \left[ \alpha \right]\left[ \frac{1}{\bar{p}} - 1 \right] \\ \Longrightarrow \hat{\beta} &= \left[ (1 - \bar{p})\frac{\bar{p}^2}{s^2} - \bar{p} \right] \left[ \frac{1}{\bar{p}} - 1 \right]. \end{aligned} \]
Now that we have these two equations in terms of the known quantities, we can find the Method of Moments estimates for \(\alpha\) and \(\beta\) given the \(n = 12\) “observed” data points:
pBar <- mean(p_num)
pS2 <- var(p_num)
(alphaHat_MoM <- (1 - pBar)*(pBar^2 / pS2) - pBar^2)
[1] 86.7
(betaHat_MoM <- alphaHat_MoM * (1 / pBar - 1))
[1] 167So, while the true parameter values used to generate this sample of \(n = 12\) probabilities were \(\alpha\) = 65 and \(\beta\) = 121, our Method of Moments estimates are \(\hat{\alpha}\) = 86.7 and \(\hat{\beta}\) = 166.7.
We have similar qualms here as in the Method of Moments section, in that it’s a strange thing to think about “observing” probabilities. However, we can still work through some of the MLE steps for the Beta Distribution. Let \(\textbf{p} = p_1, p_2, \ldots, p_n\) be an independent and identical sample from a Beta Distribution with unknown parameters \(\alpha\) and \(\beta\). We start with the Likelihood function: \[ \begin{aligned} \mathcal{L}(\alpha, \beta|\textbf{p}) &= \prod_{i = n}^n \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} p_i^{\alpha - 1} (1 - p_i)^{\beta - 1} \\ \Longrightarrow \ell(\alpha, \beta|\textbf{p}) &= \sum_{i = 1}^n \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} p_i^{\alpha - 1} (1 - p_i)^{\beta - 1} \right] \\ &= \sum_{i = 1}^n \left\{ \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \right] + \log \left[ p_i^{\alpha - 1} \right] + \log \left[ (1 - p_i)^{\beta - 1} \right] \right\} \\ &= n\log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \right] + (\alpha - 1) \sum_{i = 1}^n \log(p_i) + (\beta - 1)\sum_{i = 1}^n \log(1 - p_i). \end{aligned} \]
We immediately see two problems:
The partial derivative with respect to \(\alpha\) is found by: \[ \begin{aligned} \ell(\alpha, \beta|\textbf{p}) &= n\log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \right] + (\alpha - 1) \sum_{i = 1}^n \log(p_i) + (\beta - 1)\sum_{i = 1}^n \log(1 - p_i) \\ &= n\log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \right] - n\log \left[ \Gamma(\beta) \right] + (\alpha - 1) \sum_{i = 1}^n \log(p_i) + (\beta - 1)\sum_{i = 1}^n \log(1 - p_i) \\ \Longrightarrow \frac{\partial}{\partial \alpha} \ell(\alpha, \beta|\textbf{p}) &= n \frac{\partial}{\partial \alpha} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \right] - [0] + \sum_{i = 1}^n \log(p_i) + [0] \\ &= n \frac{\partial}{\partial \alpha} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \right] + n\overline{\log(p)} \\ 0 &\overset{\text{set}}{=} n \frac{\partial}{\partial \alpha} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \right] + n\overline{\log(p)} \\ \Longrightarrow -\overline{\log(p)} &= \frac{\partial}{\partial \alpha} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \right]. \end{aligned} \] We can’t get much further than this analytically. If either \(\alpha\) or \(beta\) are known to be an integer (for instance, if we know that these parameters truly represent integer counts of successes and failures in prior studies), then we can employ the Continued Recurrence Property of the Gamma Function as follows (but it doesn’t make the problem much easier): \[ \begin{aligned} -\overline{\log(p)} &= \frac{\partial}{\partial \alpha} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \right] \\ &= \frac{\partial}{\partial \alpha} \log \left[ \frac{\Gamma(\alpha) \prod_{j = 0}^{\beta - 1} (\alpha + j)}{\Gamma(\alpha)} \right] \\ &= \frac{\partial}{\partial \alpha} \log \left[ \prod_{j = 0}^{\beta - 1} (\alpha + j) \right] \\ &= \frac{\partial}{\partial \alpha} \sum_{j = 0}^{\beta - 1} \log (\alpha + j) \\ &= \sum_{j = 0}^{\beta - 1} \frac{\partial}{\partial \alpha} \log (\alpha + j) \\ &= \sum_{j = 0}^{\beta - 1} \frac{1}{\alpha + j} \\ &= \frac{1}{\alpha} + \frac{1}{\alpha + 1} + \ldots + \frac{1}{\alpha + \beta - 1}. \end{aligned} \] Now, this may look easier to work with (and it is easier computationally, because this is a smooth, positive, decreasing function for \(\alpha > 0\)), but notice that the unknown parameter \(\beta\) is now in the support of the computation for \(\alpha\). So, we are still in a situation where we need to estimate one parameter first in order to approximate a value for the second. If we have a guess for \(\beta\) then this becomes a manageable coding problem to estimate \(\alpha\). Also, note that because \(p\in(0,1)\), \(-\overline{\log(p)} \in (0,\infty)\), so there will be a solution to the problem.
The partial derivative with respect to \(\beta\) is almost identical, so I will skip some of the same steps we saw above: \[ \begin{aligned} \ell(\alpha, \beta|\textbf{p}) &= n\log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \right] + (\alpha - 1) \sum_{i = 1}^n \log(p_i) + (\beta - 1)\sum_{i = 1}^n \log(1 - p_i) \\ &= n\log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\beta)} \right] - n\log \left[ \Gamma(\alpha) \right] + (\alpha - 1) \sum_{i = 1}^n \log(p_i) + (\beta - 1)\sum_{i = 1}^n \log(1 - p_i) \\ \Longrightarrow \frac{\partial}{\partial \beta} \ell(\alpha, \beta|\textbf{p}) &= n \frac{\partial}{\partial \beta} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\beta)} \right] - [0] + [0] + \sum_{i = 1}^n \log(1 - p_i) \\ &= n \frac{\partial}{\partial \beta} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\beta)} \right] + n\overline{\log(1 - p)} \\ 0 &\overset{\text{set}}{=} n \frac{\partial}{\partial \beta} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\beta)} \right] + n\overline{\log(1 - p)} \\ \Longrightarrow -\overline{\log(1 - p)} &= \frac{\partial}{\partial \beta} \log \left[ \frac{\Gamma(\alpha + \beta)}{\Gamma(\beta)} \right] \\ &\qquad\text{\emph{Assuming one parameter can be an integer...}} \\ &= \frac{\partial}{\partial \beta} \log \left[ \frac{\Gamma(\beta) \prod_{j = 0}^{\alpha - 1} (\beta + j)}{\Gamma(\beta)} \right] \\ &= \frac{\partial}{\partial \beta} \sum_{j = 0}^{\alpha - 1} \log (\beta + j) \\ &= \sum_{j = 0}^{\alpha - 1} \frac{1}{\beta + j}. \end{aligned} \]
We now have a system of two equations and two unknowns (assuming that we can limit \(\alpha\) and \(\beta\) to the integers—if we cannot, then we have to use the Digamma functions19). That system is: \[ -\overline{\log(p)} = \sum_{j = 0}^{\beta - 1} \frac{1}{\alpha + j};\quad -\overline{\log(1 - p)} = \sum_{j = 0}^{\alpha - 1} \frac{1}{\beta + j}. \] Theoretically, we could solve this numerically. Often, we would use the Method of Moments estimates as the initial values for \(\alpha\) and \(\beta\).
The Continuous Uniform Distribution over \((0,1)\) is a special case of the Beta Distribution with \(\alpha = 1\) and \(\beta = 1\).
Try to write some code that will estimate \(\alpha\) and \(\beta\) from a vector of known “observed” probabilities, \(\textbf{p}\).
Other exercises to be determined.
https://gabriel.quarto.pub/stat-distributions-primer/chapters/bernoulli_20250310.html#the-likelihood-function↩︎
https://en.wikipedia.org/wiki/Support_(mathematics)#Support_of_a_distribution↩︎
For example, the Negative Binomial distribution has parameter vector \(\langle n,p \rangle\) when \(n\) is unknown, where \(p\in (0,1)\) and \(n\in\mathbb{N}\).↩︎
See this discussion for more details: https://math.stackexchange.com/questions/4244890/intuition-of-beta-distribution-with-less-than-one-parameters↩︎
Covered in our Formal Foundations chapter on the Gamma and Beta functions↩︎
https://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html↩︎
https://gabriel.quarto.pub/stat-distributions-primer/chapters/negative_binomial_20250310.html#taylormaclaurin-series↩︎
https://gabriel.quarto.pub/stat-distributions-primer/chapters/theory_gamma_function_20250707.html#the-complete-beta-function↩︎
https://gabriel.quarto.pub/stat-distributions-primer/chapters/theory_gamma_function_20250707.html#the-gamma-continued-recurrence-equation↩︎
For continuous distributions, the mode is the most likely value.↩︎
A constrained set of possible values, inspired by the concept of a production frontier in economics. See https://en.wikipedia.org/wiki/Production%E2%80%93possibility_frontier↩︎