a. Encontre o MLE de \(\theta \)
Sejam \(X_1, X_2 \dots, X_{190}\) variáveis aleatórias que representam o tipo de haptoglobina de cada uma das 190 pessoas da amostra.
Para encontrarmos a estimativa por máxima verossimilhança (MLE) de \(\theta \), precisamos encontrar o ponto \(\theta = \theta_0 \) que maximiza a função de verossimilhança. Sejam \(n_1, n_2, n_3 \) a quantidade de cada um dos tipos de haptoglobina na amostra, a função de verossimilhança é definida da seguinte forma:
\[
\begin{align}
ver(\theta) &= f(x_1, \dots, x_{190}|\theta) \\
&= f(x_1|\theta) \dots f(x_{190}|\theta) && (x_1, \dots, x_{190} \; \text{são amostras iid})\\
&= P(X=x_1 | \theta) \dots P(X=x_{190}|\theta) \\
&= (1-\theta)^{2n_1} [2\theta (1-\theta)]^{n_2} \theta^{2n_3} \\
&= 2^{n_2} \theta^{n_2 + 2n_3} (1-\theta)^{2n_1 + n_2}
\end{align}
\]
Como o logaritmo natural é uma função monotônica, podemos aplicá-lo na função de verossimilhança sem alterar seu ponto de máximo. Em seguida, podemos encontrar o ponto crítico dessa função e verificar se ele é um ponto de máximo.
\[
\mathcal{L}(\theta) = \ln[ver(\theta)] = n_2\ln(2) + (n_2 + 2n_3) \ln(\theta) + (2n_1 + n_2)\ln(1-\theta) \\
\frac{\partial \mathcal{L(\theta)}}{\partial \theta} = \frac{n_2 + 2n_3}{\theta} - \frac{2n_1 + n_2}{1-\theta} = 0 \iff \\
\boxed{\theta^{MLE} = \frac{n_2 + 2n_3}{2*(n_1 + n_2 + n_3)} = \frac{292}{380} = 0.7684}
\]
Vamos checar se esse ponto crítico é realmente um ponto de máximo global tomando a segunda derivada da função e verificando seu sinal:
\[
\frac{\partial^2 \mathcal{L\theta)}}{\partial \theta^2} = -\frac{n_2 + 2n_3}{\theta^2} - \frac{2n_1 + n_2}{(1-\theta)^2} < 0 \qquad \forall \theta \in (0, 1).
\]
Assim, \(\mathcal{L}(\theta)\) é estritamente côncava e o ponto \(\theta^{MLE} = 0.7684\) é a estimativa de máxima verossimilhança para o parâmetro.
b. Encontre a variância assintótica do estimador de ML.
A variância assintótica de um estimador \(\theta^{MLE}\) de máxima verossimilhança é dada por:
\[
\begin{align}
Var(\theta^{MLE}) &= -\frac{1}{E[\mathcal{L}^{\prime\prime}(\theta)]} \\
&= -\frac{1}{-\frac{n_2 + n_3}{\theta^2} -\frac{2n_1 + n_2}{(1-\theta)^2}} \\
&= \frac{1}{\frac{292}{0.7684^2} + \frac{88}{(1-0.7684)^2}} \\
&= \boxed{0.0004683}
\end{align}
\]
c. Encontre uma aproximação do intervalo de confiança de 99% para \(\theta\).
Assintoticamente, o estimador de máxima verossimilhança tem distribuição normal, que é simétrica. Assim, o intervalo de confiança de 99% para \(\theta^{MLE}\) é:
\[
[\theta^{MLE} \pm z(0.995) * s_{ \theta^{MLE}}]
\]
em que \(s_{ \theta^{MLE}}\) é a estimativa do erro padrão de \(\theta^{MLE}\) e z(0.995) é o percentil 99,5% da distribuição normal padrão.
\[
s^2_{\theta^{MLE}} = \sqrt{0.0004683} = 0.02164
\]
Podemos utilizar o Python para consultar a tabela da distribuição normal
import scipy.stats
print(f'Z-score for the 99.5th percentile is {scipy.stats.norm.ppf(0.995):0.4f}')
Z-score for the 99.5th percentile is 2.5758
Assim, o intervalo de confiança de interesse é dado por
\[
[0.7684 - 2.5758*0.002164, 0.7684 + 2.5758*0.002164] = \boxed{ [0.7127, 0.8242]}
\]
d. Utilize o bootstrap para encontrar uma aproximação do desvio padrão do estimador de MLE e compare-o com o resultado da parte (b).
Da parte (a), já sabemos que \(\theta^{MLE} = \frac{n_2 + 2n_3}{2*190}\)
import numpy
# given a sample of size 190, return its MLE of theta
def get_theta_mle(x):
n2, n3 = sum(x == 'Hp1-2'), sum(x == 'Hp2-2')
return (n2 + 2*n3) / 380
# ======= ======= ======= parameters setup
N_bootstrap = 100_000
sample_size = 190
population = ['Hp1-1', 'Hp1-2', 'Hp2-2']
theta = 292 / 380 # = 0.7684
probs = [(1 - theta) ** 2, 2 * theta * (1 - theta), theta ** 2]
# ======= ======= ======= calculate theta's standard deviation via bootstrap
# get a random sample of size 190x100000 from ['Hp1-1', 'Hp1-2', 'Hp2-2'] according to the given probabilities
# each one of the 100000 column is a different sample of haptglobins with size 190
# matrix shape: 190x100000
bootstrap_samples = numpy.random.choice(population, size=(sample_size, N_bootstrap), p=probs)
# evaluate MLE of theta for each column (i.e. each sample)
# matrix shape: 100000x1
bootstrap_thetas = numpy.apply_along_axis(func1d=get_theta_mle, axis=0, arr=bootstrap_samples)
print(f'Theta\'s standard deviation estimation via bootstrap method is {bootstrap_thetas.std()}')
Theta's standard deviation estimation via bootstrap method is 0.02169
As estimativas do desvio-padrão de \(\theta^{MLE}\)$ pelos dois métodos foram bem próximas, divergindo apenas a partir a 5ª casa decimal.
e. Utilize o bootstrap para econtrar uma aproximação do intervalo de confiança de 99% para θ e compare-o com a parte (c).
Aproveitando o código da parte (d), o intervalo de confiança é encontrado por meio dos quantis referentes a 0,5% e 99,5% da amostra de theta calculados em cada 1 das 100000 iterações de bootstrap
# sort array in ascending order
bootstrap_thetas_sorted = numpy.sort(bootstrap_thetas)
# get the 0.5% quantile estimate
bootstrap_ci_estimated_inf = bootstrap_thetas_sorted[int(0.005*N_bootstrap)]
# get the 99.5% quantile estimate
bootstrap_ci_estimated_sup = bootstrap_thetas_sorted[int(0.995*N_bootstrap)]
print(f'Theta\'s 99% confidence interval is [{bootstrap_ci_estimated_inf:.4f},{bootstrap_ci_estimated_sup:.4f}]')
Theta's 99% confidence interval is [0.7105,0.8211]
As estimativas do intervalo de confiança para \(\theta\) por ambos os métodos também foram bem próximas, divergindo apenas a partir da 3ª casa decimal