# Part 1: Probability distributions Question a: A contestant on a game show needs to answer 10 questions correctly to win the jackpot. However, if they get 4 incorrect answers, they are kicked off the show. Suppose one contestant consistently has a 80\% chance of correctly responding to any question. \begin{enumerate} \item What is the probability distribution? \item What is the probability that she will correctly answer 10 questions before 4 incorrect responses? \item Write out the R code to calculate (b). \end{enumerate} **Solution** 1. The experiment is modeled by the \textbf{Negative Binomial Distribution}, as it counts the number of failures ($k$) before a specified number of successes ($r$) occurs. \begin{itemize} \item Success: Correct answer, $p = 0.8$ \item Failure: Incorrect answer, $q = 1-p = 0.2$ \item Required successes: $r=10$ \item Losing condition: $k=4$ failures \end{itemize} 2. The contestant wins if she achieves $r=10$ successes before getting 4 failures. This means the total number of failures $k$ must be in the set $\{0, 1, 2, 3\}$. The probability mass function is $P(X=k) = \binom{k+r-1}{r-1} p^r q^k$. The total probability of winning is the sum over the allowed values of $k$: $$ P(\text{Win}) = \sum_{k=0}^{3} \binom{k+10-1}{10-1} (0.8)^{10} (0.2)^k $$ $$ P(\text{Win}) = \binom{9}{9}(0.8)^{10}(0.2)^0 + \binom{10}{9}(0.8)^{10}(0.2)^1 + \binom{11}{9}(0.8)^{10}(0.2)^2 + \binom{12}{9}(0.8)^{10}(0.2)^3 $$ $$ P(\text{Win}) \approx 0.1074 + 0.2147 + 0.2362 + 0.1890 \approx 0.7473 $$ 3. The calculation can be performed concisely using R's cumulative negative binomial function, \texttt{pnbinom}. \begin{verbatim} # Calculates the probability of 0, 1, 2, or 3 failures (q) # before 10 successes (size) are achieved. pnbinom(q = 3, size = 10, prob = 0.8) # Output: [1] 0.747323 \end{verbatim} Question b: A small town’s police department issues 5 speeding tickets per month on average. \begin{enumerate} \item Using a Poisson random variable, what is the likelihood that the police department issues 3 or fewer tickets in one month? \item What is the probability that 10 days or fewer elapse between two tickets being issued? \item Write out the R code to calculate (a), (b). \end{enumerate} **Solution** 1. First, we note that here $P(Y \leq 3)=P(Y=0)+P(Y=1)+\cdots+P(Y=3)$. Applying the probability mass function for a Poisson distribution with $\lambda=5$, we find that $$ \begin{aligned} P(Y \leq 3) &=P(Y=0)+P(Y=1)+P(Y=2)+P(Y=3) \\ &=\frac{e^{-5} 5^{0}}{0 !}+\frac{e^{-5} 5^{1}}{1 !}+\frac{e^{-5} 5^{2}}{2 !}+\frac{e^{-5} 5^{3}}{3 !} \\ &=0.27 . \end{aligned} $$ ```{r} # q is the number of events (3 or fewer). # lambda is the average rate of events per interval. ppois(q = 3, lambda = 5) ``` 2. We know the town's police issue 5 tickets per month. For simplicity's sake, assume each month has 30 days. Then, the town issues $\frac{1}{6}$ tickets per day. That is $\lambda=\frac{1}{6}$, and the average wait time between tickets is $\frac{1}{1 / 6}=6$ days. Therefore, $$ P(Y<10)=\int_{0}^{10} \frac{1}{6} e^{-\frac{1}{6} y} d y=0.81 $$ ```{r} pexp(10, rate = 1/6) ``` # Part 2: Statistical inference \begin{enumerate} \item (AoS 6.6.2) Let $X_{1}, \ldots, X_{n} \sim \operatorname{Uniform}(0, \theta)$ and let $\hat{\theta}=\max \left\{X_{1}, \ldots, X_{n}\right\}$. Find the bias, se and MSE of this estimator. \item (AoS 6.6.3) Let $X_{1}, \ldots, X_{n} \sim \operatorname{Uniform}(0, \theta)$ and let $\hat{\theta}=2 \bar{X}_{n}$. Find the bias, se and MSE of this estimator. \item Let $X_{1}, \ldots, X_{n} \sim \operatorname{Uniform}(0,1)$. Let $Y_{n}=\bar{X}_{n}^{2}$. Find the limiting distribution of $Y_{n}$. (Hint: CLT) \end{enumerate} **Solution** 1. The CDF $G$ of $\widehat{\theta}$ is $$ \begin{aligned} G(\widehat{\theta}) &=\mathbb{P}(\widehat{\Theta} \leq \widehat{\theta}) \\ &=\mathbb{P}\left(\max \left\{X_{1}, \ldots, X_{n}\right\} \leq \widehat{\theta}\right) \\ &=\prod_{i=1}^{n} \mathbb{P}\left(X_{i} \leq \widehat{\theta}\right) \\ &=F_{\theta}(\widehat{\theta})^{n} \\ &=\left(\frac{\widehat{\theta}}{\theta}\right)^{n} \end{aligned} $$ The density is therefore $$ g(\widehat{\theta})=\left(\frac{n}{\theta}\right)\left(\frac{\widehat{\theta}}{\theta}\right)^{n-1} $$ Thus, $$ \mathbb{E}_{\theta}(\widehat{\theta})=\int_{0}^{\theta} \widehat{\theta} g(\widehat{\theta}) d \widehat{\theta}=\frac{n \theta}{n+1} $$ and $$ \text { bias }=\frac{n \theta}{n+1}-\theta=-\frac{\theta}{n+1} . $$ Also, $$ \mathbb{E}_{\theta}\left(\widehat{\theta}^{2}\right)=\int_{0}^{\theta} \widehat{\theta}^{2} g(\widehat{\theta}) d \widehat{\theta}=\frac{n \theta^{2}}{n+2} $$ and so $$ \mathbb{V}_{\theta}(\widehat{\theta})=\frac{n \theta^{2}}{n+2}-\left(\frac{n \theta}{n+1}\right)^{2}=\frac{n \theta^{2}}{(n+2)(n+1)^{2}} $$ The mse is $$ \operatorname{bias}^{2}+\mathbb{V}=\left(\frac{\theta}{n+1}\right)^{2}+\frac{n \theta^{2}}{(n+2)(n+1)^{2}}=\frac{2 \theta^{2}}{(n+1)(n+2)} $$ 2. Recall that $\mathbb{E}\left(X_{i}\right)=\theta / 2, \mathbb{V}\left(X_{i}\right)=\theta^{2} / 12$. So $$ \mathbb{E}_{\theta}(2 \bar{X})=2 \mathbb{E}_{\theta}(\bar{X})=2 \frac{\theta}{2}=\theta $$ and hence bias $=0$. Now $$ \mathbb{V}_{\theta}(2 \bar{X})=4 \mathbb{V}_{\theta}(\bar{X})=\frac{4 \sigma^{2}}{n}=\frac{4 \theta^{2}}{12 n}=\frac{\theta^{2}}{3 n} . $$ Since this estimator is unbiased, $$ \mathrm{mse}=\mathbb{V}_{\theta}(\widehat{\theta})=\frac{\theta^{2}}{3 n} . $$ 3. $\mu=\mathbb{E}\left(X_{i}\right)=1 / 2$ and $\sigma^{2}=\mathbb{V}\left(X_{i}\right)=1 / 12$. By the CLT, $$ \frac{\sqrt{n}(\bar{X}-\mu)}{\sigma}=\sqrt{12 n}\left(\bar{X}-\frac{1}{2}\right) \rightsquigarrow N(0,1) . $$ Now $Y=g(\bar{X})$ where $g(s)=s^{2}$. And $g^{\prime}(s)=2 s$ and $g^{\prime}(\mu)=g^{\prime}(1 / 2)=$ $2(1 / 2)=1$. From the delta method, $$ \frac{\sqrt{n}(Y-g(\mu))}{\left|g^{\prime}(\mu)\right| \sigma}=\sqrt{12 n}\left(\bar{X}-\frac{1}{4}\right) \leadsto N(0,1) $$ # Part 3: Newton-Raphson implementation The $\mathrm{ABO}$-gene or $\mathrm{ABO}$-locus is on chromosome 9. It has 3 alleles (antigens) $(A, B, O$) and it determines 4 blood type $(A, B, A B, O)$. We have a large random sample obtained from Berlin (Bernstein 1925, Sham's book page 44): - $n_{A}=9123$ blood type $A$ - $n_{B}=2987$ blood type $B$ - $n_{A B}=1269$ blood type $A B$ - $n_{O}=7725$ blood type O For instance, $n_{A}=9123=n_{A A}+n_{A O}:$ Among 9123 blood type $A$ individuals, some have genotype $A A$ and the others have genotype $A O$. Our interest is to estimate the allele frequencies of alleles A, B, and O. i.e. $p=$ freq $($allele $A)$, $q=$ freq $($allele $B)$, $1-p-q=$ freq $($allele $O)$. 1. Write out the log-likelihood $L(p,q)$. 2. Is there a closed-form solution of this log-likelihood function? 3. Formulate the problem as a missing data problem and use the Newton-Raphson algorithm to find the MLEs, $\hat{p}$ and $\hat{q}$, that maximize the log-likelihood, $\ln L(p, q)$. \clearpage **Solution** 1. Let $X=\left(n_{A}, n_{B}, n_{A B}, n_{O}\right), \mathrm{X}$ follows the multinomial distribution, $$ L(p, q)=\left(\begin{array}{c} n \\ n_{A}, n_{B}, n_{A B}, n_{O} \end{array}\right)\left(p^{2}+2 p(1-p-q)\right)^{n_{A}}\left(q^{2}+2 q(1-p-q)\right)^{n_{B}}(2 p q)^{n_{A B}}\left((1-p-q)^{2}\right)^{n_{O}} $$ The general approach to estimate allele frequency is maximum likelihood estimation. To find MLE, take the derivatives of the log-likelihood function, and find the pair of $(p, q)$ that set the derivatives to 0. The log transformation is performed since finding a maximizer of $L(p,q)$ is equivalent to finding a maximizer of $\ln L(p,q)$. The log-likelihood is, \begin{equation} \ln L(p, q) \sim n_A \ln (p^2 + 2p(1- p-q )) + n_B \ln (q^2 + 2q(1 - p - q)) + n_{AB} \ln (2pq) + n_O \ln ((1- p -q)^2) \end{equation} Take the partial derivatives of the log-likelihood and set them to 0, \begin{equation} \frac{\partial \ln L(p,q)}{\partial p} = \frac{2(1 - p - q)}{p(2 - p - 2q)} n_A + \frac{2}{2p+q-2} n_B + \frac{1}{p} n_{AB} - \frac{2}{1-p-q} n_{O} =0 \end{equation} \begin{equation} \frac{\partial \ln L(p,q)}{\partial q} = \frac{2}{p+2q-2}n_A + \frac{2(1-p-q)}{q(2-2p-q)} n_B + \frac{1}{q}n_{AB} - \frac{2}{1-p-q} n_{O} =0 \end{equation} 2. It is hard to find the explicit form of the $(p, q)$ from equation (2) and (3) \textbf{directly}. We consider solve the problem \textbf{iteratively}. 3. Newton-Raphson algorithm is another method that can be applied to estimate the ABO allele frequency. The algorithm is initially designed to find the roots (or zeros) of a real-valued function iteratively. In the ABO blood type settings, maximizing the likelihood is equivalent to finding the roots of the derivative of log-likelihood function. Since it is hard to find the roots directly, one can approximate it iteratively. Here denote the parameters to be estimated as $\vec{\theta} = (p, q)$, denote the log-like hood shown in equation (1) as $f(\vec{\theta}) = \ln L(\theta)$, then the partial derivatives of the log-likelihood (score function) is $f'(\vec{\theta}) = f'(p,q)$, the second derivatives (Hessian matrix, or observed information $- I(\theta)$) as $f''(\vec{\theta}) = f''(p,q)$. The explicit forms of the score function and observed information are shown in the appendix. 4. Expectation-Maximum (EM) algorithm is a method for obtaining the Maximum likelihood estimates (MLE) of parameters iteratively. It usually contains two parts: \textbf{E-step (expectation)} computes an expected value of the log-likelihood using current estimate for the parameters; \textbf{M-step (Maximization)} calculates MLE based on the log likelihood in the E-step, then updates the estimates of the parameters. The EM algorithms are often used when the model contains unobserved data. In the ABO settings, one could observe the phenotype counts $(n_A, n_B, n_{AB}, n_O)$, while the genotype counts $n_{AA}$ or $n_{AO}$ in blood $A$ group and $n_{BB}$ or $n_{BO}$ in blood $A$ group and $B$ group are missing. This leads to a problem when applying direct counting (Sham, 1998) to estimate $p, q$, the allele frequency of $A$ and $B$ respectively. The following steps will show how this missing data problem could be solved by the EM algorithm. In each iteration $k$, firstly, the \textbf{E-step} computes the expected value of the log-likelihood $h(p,q)$ using the observed data $n_{\text{obs}} = (n_A, n_B, n_{AB}, n_O)$ and current parameter value $p^{(k)}, q^{(k)}$: \begin{equation} Q_k (p,q) = E[h_k (p,q) | n_{\text{obs}}, p^{(k)}, q^{(k)}] \end{equation} and $h_k (p,q) \sim 2 n_{AA} \log p^{(k)} + n_{AO} \log(2p^{(k)} (1-p^{(k)}-q^{(k)})) + 2 n_{BB} \log q^{(k)} + n_{BO} \log (2p^{(k)}(1-p^{(k)}-q^{(k)})) + n_{AB} \log(2p^{(k)}q^{(k)}) + 2 n_O \log(1-p^{(k)}-q^{(k)})$ Since the log-likelihood is linear w.r.t. the missing data, therefore, when taking the expectation for each component, the missing data can be imputed in this case. For example, since $E[2 n_{AA} \log p^{(k)} | n_{\text{obs}}, p^{(k)}, q^{(k)}] = 2 \log p^{(k)} E[n_{AA}| n_{\text{obs}}, p^{(k)}, q^{(k)}]$, and under the HWE assumption, $$E[n_{AA} | n_{\text{obs}}, p^{(k)}, q^{(k)}] = \frac{\text{freq}(A A)}{\text{freq}(A A)+\text{freq}(A O)} n_{A} =\frac{p^{(k)} p^{(k)}}{p^{(k)} p^{(k)}+2 p^{(k)}\left(1-p^{(k)}-q^{(k)}\right)} n_{A} := n^{(k)}_{AA}$$ By the similar calculation, we can get the imputed data $n_{AA}^{(k)}, n_{AO}^{(k)}, n_{BB}^{(k)}, n_{BO}^{(k)}$ for each iteration. Then the \textbf{M-step} computes the MLE based on the likelihood in equation (4). To be specific, $$\frac{\partial Q_k (p,q)}{\partial p} = \frac{2 n_{A A}^{(k)}+\ n_{A O}^{(k)}+n_{A B}}{p }-\frac{n_{A O}^{(k)}+ n_{B O}^{(k)}+n_{O}}{1-p-q} = 0$$ $$\frac{\partial Q_{k}(p,q)}{\partial q}=\frac{2 n_{B B}^{(k)}+ n_{B O}^{(k)}+n_{A B}}{q}-\frac{n_{A O}^{(k)}+ n_{B O}^{(k)}+n_{O}}{1-p-q} = 0$$ Thus, the updated values of parameters are $$p^{(k+1)}=\frac{2 n_{A A}^{(k)}+ n_{A O}^{(k)}+n_{A B}}{2 n} \qquad \quad q^{(k+1)}=\frac{2 n_{B B}^{(k)}+n_{B O}^{(k)}+n_{A B}}{2 n}$$ ```{r, echo=FALSE} nA <- 9123 nB <- 2987 nAB <- 1269 nO <- 7725 n <- nA + nB + nAB + nO ``` ```{r, echo=FALSE} log_like <- function(nA, nB, nAB, nO, p, q) { n <- nA + nB + nAB + nO nA*log(p*p + 2*p*(1-p-q)) + nB*log(q*q + 2*q*(1-p-q)) + nAB*log(2*p*q) + nO*log((1-p-q)^2) } ``` ```{r, echo=TRUE} df <- function(p, q) { dfp <- 2*(1 - p - q)*nA / (p*(2 - p - 2*q)) + 2*nB / (2*p + q - 2) + nAB/p - 2*nO / (1 - p - q) dfq <- 2*nA / (p + 2*q - 2) + 2*(1 - p - q)*nB / (q*(2 - 2*p - q)) + nAB/q - 2*nO / (1 - p - q) c(dfp, dfq) } d2f <- function(p, q) { pp <- -4*(1 - p - q)^2*nA / (p^2*(2 - p - 2*q)^2) - 2*nA / (p* (2 - p - 2*q)) - 4*nB / ((2 - 2*p - q)^2) - nAB / (p^2) - 2*nO / ((1 - p - q)^2) pq <- 4*(1 - p - q)*nA / (p*(2 - p - 2*q)^2) - 2*nA /(p*(2 - p - 2*q)) - 2*nB / ((2 - 2*p - q)^2) - 2*nO / ((1 - p - q)^2) qp <- -2*nA / ((2 - p - 2*q)^2) + 4*(1 - p - q)*nB / (q *(2 - 2*p - q)^2) - 2*nB / (q* (2 - 2*p - q)) - 2*nO / ((1 - p - q)^2) qq <- -4*nA / ((2 - p - 2*q)^2) - 4*(1 - p - q)^2*nB /((q^2*(2 - 2*p - q)^2)) - 2*nB / (q*(2 - 2*p - q)) - nAB / (q^2) - 2*nO / ((1 - p - q)^2) matrix(c(pp, pq, qp, qq), nrow = 2, ncol = 2, byrow = T) } max <- 10000 epsilon <- 1e-5 iter <- 0 p <- 0.3333333 q <- 0.3333333 diff1 <- 1 diff2 <- 1 nA <- 9123 nB <- 2987 nAB <- 1269 nO <- 7725 theta <- c(p, q) rp <- NULL rq <- NULL rlike <- NULL rp[1] <- p rq[1] <- q rlike[1] <- log_like(nA, nB, nAB, nO, p, q) while (diff1 > epsilon & diff2 > epsilon & iter < max) { p <- theta[1] q <- theta[2] theta.new <- theta - solve(d2f(p,q)) %*% df(p, q) diff1 <- abs(theta[1] - theta.new[1]) diff2 <- abs(theta[2] - theta.new[2]) theta <- theta.new log_lik <- log_like(nA, nB, nAB, nO, theta[1], theta[2]) iter <- iter + 1 rp[iter+1] <- theta[1] rq[iter+1] <- theta[2] rlike[iter+1] <- log_lik } knitr::kable( data.frame(c(0:(length(rp)-1)), rp, rq, rlike), col.names = c("iteration $k$", "$p^{(k)}$", "$q^{(k)}$", "log likelihood"), booktabs = TRUE, align = "cccr", caption = 'Results for Newton-Raphson algorithm' ) ```