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

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

0 votos
49 visitas
perguntada Abr 11 em Estatística por gustavobrangel (6 pontos)  

Um naturalista inglês coletou dados sobre o comprimento dos ovos de cuco, medindo-os com precisão de 0,5 mm. Examine a normalidade dessa distribuição (a) construindo um histograma e superpondo a função de densidade normal, (b) plotando a probabilidade normal e (c) construindo um hanging rootgram.

\begin{array}{|c|c|}
\hline
Comprimento & Frequência \\ \hline
18.5 & 0 \\
19.0 & 1 \\
19.5 & 3 \\
20.0 & 33 \\
20.5 & 39 \\
21.0 & 156 \\
21.5 & 152 \\
22.0 & 392 \\
22.5 & 288 \\
23.0 & 286 \\
23.5 & 100 \\
24.0 & 86 \\
24.5 & 21 \\
25.0 & 12 \\
25.5 & 2 \\
26.0 & 0 \\
26.5 & 1 \\ \hline
\end{array}

Compartilhe

1 Resposta

0 votos
respondida Abr 11 por gustavobrangel (6 pontos)  
editado Abr 11 por gustavobrangel

Vamos utilizar o Python para construir todos os gráficos.

(a) Construir um histograma e superpor a função de densidade normal

A construção do histograma é imediata, mas precisamos ajustá-lo para que ele e a densidade da normal fiquem na mesma escala. Assim, vamos examinar a frequência relativa de cada comprimento observado.

A função densidade teórica da normal é construída assumindo que a média e a variância sejam iguais às da amostra

# generating data
df = pandas.DataFrame({
'length': numpy.arange(start=18.5, stop=27, step=0.5),
'frequency': [0, 1, 3, 33, 39, 156, 152, 392, 288, 286, 100, 86,
              21, 12, 2, 0, 1]
})

# number of obs is used to calculate relative frequency
n_obs = df.frequency.sum()

fig, ax = pyplot.subplots()

# histogram of eggs' length according to its relative frequency
ax.bar(x=df.length, height=df.frequency / n_obs / 0.5,
       width=0.5, color='C0')

sample_mean = (df.length * df.frequency).sum() / df.frequency.sum()

sample_var = sum([row.frequency * (row.length - sample_mean)**2 
                  for row in df.itertuples()]) / (n_obs - 1)
sample_std = sample_var ** 0.5    

# build a normal pdf from length=18.5 to length=26.5
x_length = numpy.linspace(start=df.length.values[0], 
               stop=df.length.values[-1], num=10000)
normal_pdf = stats.norm.pdf(x, sample_mean, sample_std)

ax.plot(x_length, normal_pdf, color='C1')

# ---------- ---------- ---------- graph details
ax.set_xlabel('Length')
ax.set_ylabel('Relative frequency')
ax.set_title("Eggs' length histogram vs normal density plot")
ax.xaxis.set_ticks(numpy.arange(19, 27, 1))

A imagem será apresentada aqui.

Comparando o histograma com a função densidade teórica, parece haver mais assimetria nos dados do que o esperado, sendo um indicativo de que o comprimento dos ovos não segue uma distribuição gaussiana.

(b) Plotar a probabilidade normal

O gráfico de probabilidade normal pode ser interpretado como uma comparação de percentis entre a distribuição empírica e a gaussiana. Assim, se estivermos certos quanto ao processo gerador dos dados empíricos, esse gráfico será uma linha reta com inclinação de 45º.

Como nossos dados estão com precisão de 0,5 mm, ou seja, estão discretizados, o gráfico de probabilidade normal fica um pouco mais difícil de interpretar. Assim, para facilitar a análise, resolvi mostrá-lo juntamente com esse mesmo tipo de gráfico, mas dessa vez usando dados simulados oriundos de uma distribuição normal (n = 100_000).

# sorting observation by ascending order
# ordered_obs = [19, 19.5, 19.5, 19.5, 20, ..., 25.5, 26.5]
ordered_obs = sorted(df.length.repeat(df.frequency))

# theoretical z-score distribution
normal_quantiles = [stats.norm.ppf(k / (n_obs + 1))
                    for k in range(1, n_obs + 1)]

fig, ax = pyplot.subplots()

ax.scatter(normal_quantiles, ordered_obs, color='C0', 
           label="Empirical normal probability")

# ---------- simulating a normal probability plot for comparison
n_sim = 100_000

# sampling from a normal distribuition and 
# discretizing it to a 0.5 precision mm as the exercise data
random_normal = (2*numpy.random.normal(sample_mean, 
                 sample_var**0.5, size=n_sim)).round() / 2

normal_quantiles_sim = [stats.norm.ppf(k / (n_sim + 1)) 
                        for k in range(1, n_sim + 1)]

ax.plot(normal_quantiles_sim, sorted(random_normal), color='C1',
        label='Simulated normal probability')

# ---------- ---------- graph details
ax.set_xlabel('Z-score')
ax.set_ylabel('Length')
ax.set_title("Normal probability plot: acual vs expected \ 
                (simulated)")
ax.legend()

A imagem será apresentada aqui.

Assim como o anterior, também é possível ver por esse gráfico que esse os dados são mais assimétricos do que o esperado para uma distribuição normal, já que os percentis de ambas as distribuiçẽs estão divergindo. Essa divergência, porém, pode ter ocorrido por mera chance, principalmente porque estamos fazendo a comparação com uma distribuição simulada.

(c) Construa um hanging rootgram

Outra forma de compararmos a distribuição dos dados empírica com uma distribuição teórica é por meio da diferença entre a quantidade de amostras observadas \(n_j\) e a quantidade esperada \(\hat n_j\) em cada intervalo \(j \). O gráfico dessa diferença é conhecido como hanging histogram.

Vamos, então, definir os seguintes intervalos para cada um dos nossos comprimentos:
\({[18.25, 18.75), [18.75, 19.25), \cdots, [25.75, 26.25), [26.25, 26.75]}\).

Sejam \(x_{j-1}\) o limite inferior e \(x_j\) o limite superior de um intervalo, a probabilidade de uma amostra cair nesse intervalo, dado que os dados são normalmente distribuiídos é:

\[ \hat p_j = \Phi \Bigg ( \frac{ x_j - \bar x }{\hat \sigma} \Bigg ) - \Phi \Bigg (\frac{ x_{j-1} - \bar x }{\hat \sigma} \Bigg) \quad j \in \{1, 2, \cdots, 17\} \]

O número de amostras esperado dentro de cada intervalo é \(\hat n_j = n \hat p_j\), em que \(n\) é o tamanho total da amostra. O gráfico abaixo mostra o hanging histogram, ou seja, a diferença \(n_j - \hat n_j\) para cada um dos 17 intervalos. Esse gráfico, entretanto, é difícil de interpretar porque a variância nos intervalos não é constante. Como a contagem no centro da distribuição é mais alta, assim também é sua variância.

A variância em cada um dos intervalos é dada por:

\[ \begin{align} Var(n_j - \hat n_j) &= Var(n_j) = np_j(1 - p_j) \\ &= np_j - np^2_j \\ &\approx np_j \qquad (\text{para } p \text{ pequeno}) \end{align} \]

O que precisamos fazer é aplicar uma função \(g\) nos nossos dados \(X\) tal que \(Y = g(X)\) tenha variância \(\sigma^2_Y\) constante.

Fazendo a expansão de Taylor de primeira ordem de \(g(X)\) em torno de \(\mu_X\), ficamos com:
\[ Y \approx g(\mu_X) + (X - \mu_X)g^{\prime}(\mu_X) \]

Assim:

\[ \begin{align} \sigma^2_Y &\approx [g^{\prime}(\mu_X)]^2\sigma^2_X \end{align} \]

Logo, seja \(\bar x_j = np_j\) a média em cada um dos intervalos, precisamos escolher uma função \(g\) tal que \(Var(n_j) \approx [g^{\prime}(\bar x_j)]^2 np_j = [g^{\prime}(\bar x_j)]^2 \bar x_j\) não dependa mais de \(\bar x_j\). Se a função for \(f(x) = \sqrt{x}\), ficamos com:

\[ \begin{align} Var(\sqrt{n_j}) \approx \frac{1}{4\bar x_j} \bar x_j = \frac{1}{4} \end{align} \]

# expected frequency of an interval 
# is the difference between upper and lower bound z score
# times the number of observation
df['upper_bound_zscore'] = stats.norm.cdf(df.length + .25, loc=sample_mean, scale=sample_std)
df['expected_frequency'] = df.upper_bound_zscore.diff().fillna(0) * n_obs

fig, axes = pyplot.subplots(2, 1, sharex=True)

# axes[0] is the hanging histogram
axes[0].bar(x=df.length, height=df.frequency - df.expected_frequency, width=0.5, color='C0')

# axes[1] is the hanging rootogram
axes[1].bar(x=df.length, height=df.frequency**0.5 - df.expected_frequency**0.5, width=0.5, color='C0')

# ---------- ---------- graph details
axes[1].set_xlabel('Length')
axes[0].set_ylabel('Difference between actual and expected values')
axes[1].set_ylabel('Difference between actual and expected values')
axes[0].set_title('Hanging histogram')
axes[1].set_title('Hanging rootogram')
axes[0].xaxis.set_ticks(numpy.arange(19, 27, 1))

A imagem será apresentada aqui.

Analisando esse último gráfico, vemos que realmente existe uma assimetria entre os dados do exercício e a distribuição normal já que a diferença em três invervalos supera 2 (em módulo), enquanto o desvio-padrão é 0.5, ou seja, é bastante improvável que essa difereça tenha ocorrido por mero acaso.

...