Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Questão 58, Problemas do Capítulo 8 do Livro "Mathematical Statistics and Data Analysis" de John A. Rice (3ª Edição)

0 votos
37 visitas
perguntada Abr 4 em Estatística por gustavobrangel (6 pontos)  
reclassificado Abr 4 por gustavobrangel

Se a frequência de genes estão em equilíbrio, os genótipos AA, Aa, aa ocorrem com probabilidades \((1-\theta)^2 \), \(2 \theta(1-\theta) \), e \( \theta^2 \), respectivamente. Plato et al. (1964) publicaram os seguintes dados a respeito do tipo de haptoglobina em uma amostra de 190 pessoas:

\begin{array}{|l|l|l|}
\hline
Hp1-1 & Hp1-2 & Hp2-2 \\ \hline
10 & 68 & 112 \\ \hline
\end{array}

a. Encontre a MLE de \(\theta \).
b. Encontre a variância assintótica do estimador de ML.
c. Encontre uma aproximação do intervalo de confiança de 99% para \(\theta \).
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).
e. Utilize o bootstrap para econtrar uma aproximação do intervalo de confiança de 99% para \(\theta \) e compare-o com a parte (c).

Compartilhe

1 Resposta

0 votos
respondida Abr 4 por gustavobrangel (6 pontos)  
republicada 5 dias atrás por gustavobrangel

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

...