# 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'
)
```
