Estatı́stica Básica
10.2 Resolução
-
10.2.1 Vamos dividir nossos resultados nos diferentes testes ou procedimentos gráficos. Vamos iniciar com os testes.
-
a) Considerando a amostra de tamanho \(n\) \(=\) \(15\), temos os resultados apresentados na tabela seguinte, para o teste de Shapiro-Wilk. Aplicamos as fórmulas de estimação dos \(a_i\)’s apresentadas no Livro, sendo escolhidas as de maior acurácia. Também apresentamos os resultados das estimativas usando os métodos extremamente acurados apresentados por Royston (1993) e detalhadamente descritos em Ferreira (2018). Fizemos duas funções em R para realizarmos estes cálculos com mais facilidade. A função bSW.test, aplica o teste de Shapiro-Wilk com a versão aproximada, em que boa parte das expressões foram obtidas para este Livro apenas. A segunda função SW.test, utiliza uma versão mais acurada do teste, pois se baseia na proposta de Royston (1993).
A duas funções para aplicação do teste de SW, juntamente com a função já existente no próprio R, e os dados do exemplo são:
# Function to apply the Shapiro-Wilk's # normality test - version of my book bSW.test <- function(x) { n <- length(x) pis <- ((1:n) - 0.375) / (n + 0.25) mi <- qnorm(pis) if (n == 3) ai <- c(-2^0.5/2, 0, 2^0.5/2) else { ais <- 2 * mi[2:(n-1)] if (n<=20) a1 <- gamma(n/2)/(2*sqrt(2)*gamma((n+1)/2))+ gamma((n+1)/2)/(2*sqrt(2)*gamma((n+2)/2)) else if (n<=50) a1 <- gamma((n+1)/2)/(sqrt(2)*gamma((n+2)/2)) else { # n > 50 a12 <- mi[1]^2/sum(mi^2) a1 <- 30*gamma((n+1)/2)/(sqrt(2)*gamma((n+2)/2))/n+ (n-30)/n*a12 } a1s <- a1 / (1-2*a1)*sum(ais^2) ais <- c(-(a1s^0.5),ais,(a1s^0.5)) ai <- ais/(sqrt(sum(ais^2))) } w <- sum(ai * sort(x))^2 / ((n-1)*var(x)) if (n == 3) { valor.p <- 6 / pi * (asin(sqrt(w)) - asin(sqrt(0.75))) } else { if (n <= 11) { lamb <- -2.273+0.459*n U <- log(1-w) if (lamb>U) Y <- -log(lamb-U) else Y<- Inf muy <- 0.544 - 0.39978*n + 0.025054*n^2 - 0.0006714*n^3 sigY <- exp(1.3822 - 0.77857*n + 0.062767*n^2- 0.0020322*n^3) Zc <- (Y - muy) / sigY valor.p <- 1 - pnorm(Zc) } else { u <- log(n) Y <- log(1-w) muy <- -1.5861 - 0.31082*u - 0.083751*u^2 + 0.0038915*u^3 sigY <- exp(-0.4803-0.082676*u+0.0030302*u^2) Zc <- (Y - muy) / sigY valor.p <- 1 - pnorm(Zc) } } plot(ai,sort(x),xlab="a's",ylab="x's") abline(lm(sort(x)~ai),col="red") return(list(w=w,p.value=valor.p)) } # Function to apply the Shapiro-Wilk's # normality test - version of Royston(1993) SW.test <- function(x) { n <- length(x) if (n == 3) a <- c(-sqrt(2)/2,0,sqrt(2)/2) else { p <- ((1:n) - 0.375) / (n + 0.25) m <- qnorm(p) sm2 <- sum(m^2) c <- m / sqrt(sm2) u <- 1 / sqrt(n) if (n <= 5) { an <- c[n]+0.221157*u-0.14798*u^2- 2.071190*u^3+4.434685*u^4-2.706056*u^5 phi <- (sm2 - 2 * m[n]^2) / (1 - 2 * an^2) a <- c(-an,m[2:(n-1)]/phi^0.5,an) } else { an <- c[n]+0.221157*u-0.14798*u^2- 2.071190*u^3+4.434685*u^4-2.706056*u^5 an1 <- c[n-1]+0.042981*u-0.293762*u^2- 1.752461*u^3+5.682633*u^4-3.582633*u^5 phi <- (sm2 - 2 * m[n]^2 - 2 * m[n-1]^2) / (1 - 2 * an^2 - 2 * an1^2) a <- c(-an,-an1,m[3:(n-2)]/phi^0.5,an1,an) } } w <- (sum(a*sort(x)))^2 / ((n - 1) * var(x)) if (n == 3) p.value <- 6/pi*(asin(sqrt(w))-asin(sqrt(0.75))) else { # n>3 if (n <= 11) { gam <- -2.273+0.459 * n Y <- -log(gam - log(1 - w)) muY <- 0.5440-0.39978*n+ 0.025054*n^2-0.0006714*n^3 sigY<-exp(1.3822-0.77857*n+ 0.062767*n^2-0.0020322*n^3) Z <- (Y - muY) / sigY p.value <- 1 - pnorm(Z) } else { Y <- log(1- w) u <- log(n) muY <- -1.5861-0.31082*u- 0.083751*u^2+0.0038915*u^3 sigY<-exp(-0.4803-0.082676*u+0.0030302*u^2) Z <- (Y - muY) / sigY p.value <- 1 - pnorm(Z) } } plot(a,sort(x),xlab="a's",ylab="x's") abline(lm(sort(x)~a),col="red") return(list(w=w,p.value=p.value)) } # Example: Exercise 10.5.1 x <- c(8.52,9.07,11.90,9.43,4.32, 5.77,8.52,5.95,5.56,9.02, 12.08,8.06,9.96,6.02,8.69) bSW.test(x) SW.test(x) shapiro.test(x)#R function for SW's testOs valores da estatística do teste, usando a versão do Livro e de Royston (1993), são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} W=& \dfrac {\left (\displaystyle \sum _{i=1}^n \hat {a}_i X_{(i)}\right )^2}{\displaystyle \sum _{i=1}^{n} X^2_i-\dfrac {\left (\displaystyle \sum _{i=1}^{n} X_i\right )^2}{n}} =\dfrac {8,306^2}{72,99}\\ =& 0,9452 \end{align*} e
\(\seteqnumber{0}{10.}{0}\)\begin{align*} W=& \dfrac {\left (\displaystyle \sum _{i=1}^n \hat {a}_i X_{(i)}\right )^2}{\displaystyle \sum _{i=1}^{n} X^2_i-\dfrac {\left (\displaystyle \sum _{i=1}^{n} X_i\right )^2}{n}} =\dfrac {8,296^2}{72,99}\\ =& 0,9430, \end{align*} respectivamente.
Uma vez obtido o valor da estatística, usando a versão do Livro, devemos obter determinadas quantidades necessárias, que para este caso, \(n>11\), são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} x =& \ln (n) = \ln (15) = 2,7081\\ Y=& \ln (1-W)=-2,9036. \end{align*} sendo a média e desvio padrão de \(Y\), são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} \mu _Y=& -2,9647 & \textrm { e }&& \sigma _Y&= 0,5056. \end{align*}
Portanto, o valor padronizado de \(Y\) é
\(\seteqnumber{0}{10.}{0}\)\begin{align*} Z_c=& \dfrac {Y-\mu _Y}{\sigma _Y}=\dfrac {-2,9036+2,9647}{0,5056}= 0,1210. \end{align*}
Logo,
\(\seteqnumber{0}{10.}{0}\)\begin{align*} \textrm {valor-}p =& P(Z>0,1210) = 0,4519. \end{align*}
Portanto, não se deve rejeitar a hipótese nula de normalidade no nível nominal de significância de \(5\%\).
Para a segunda opção, o teste mais acurado usando os resultados de Royston (1993), temos que \(Z_c\) \(=\) 0,1975 e o valor-\(p\) \(=\) 0,4217. Assim, a conclusão é a mesma anterior. Assim, podemos concluir que nossa aproximação é bem razoável. É claro que fizemos simulações Monte Carlo para validar e não apenas um simples exemplo como este, antes de publicar o Livro. Sugerimos aos leitores que utilizem a versão mais acurada, seja a segunda função do programa acima ou a função shapiro.test do programa R, também apresentada no programa acima. A vantagem de nossa função em relação a do R é que ele faz um gráfico tipo quantil-quantil, em que um dos quantis esta relacionado com a esperança das estatísticas de ordem normalizadas (quantidades \(a_i\)’s) e o outro, são as estatísticas de ordem observadas. Este gráfico, usando o cálculo usando as expressões de Royston (1993) está apresentado a seguir:
Verificamos que a linha de mínimo quadrados se ajusta razoavelmente bem aos pontos observados \((\hat {a}_i\), \(x_{(i)})\), indicando a normalidade dos dados. Fizemos uma simulação com \(n\) \(=\) \(1000\) dados de uma distribuição gama(0,5) e a reta de mínimos quadrados ajustada foge totalmente da dispersão dos pontos, indicando, como esperado, a não normalidade dos dados.
-
b) O segundo caso a ser considerado são os testes de desvios de assimetria e curtose de D’Agostino et al. (1990). Após o exemplo, apresentamos novamente três funções para aplicação automática destes testes, que são trabalhosos de serem aplicados, mas não difíceis.
A média e o desvio padrão amostrais são \(\bar {X}\) \(=\) 8,1913 dias e \(S\) \(=\) 2,2833 t/ha. As estimativas da assimetria e curtose são \(\sqrt {b_1}\) \(=\) 0,041948, \(b_2\) \(=\) 2,2116 e \(g_2\) \(=\) \(-\)0,59356.
-
i) Teste de assimetria
As hipóteses nula e alternativa para o teste de assimetria são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} H_0: \sqrt {\beta _1}=0 \qquad \qquad \textrm { versus } \qquad \qquad H_1: \sqrt {\beta _1}\ne 0. \end{align*}
-
1) Primeira alternativa
A estatística do teste é
\(\seteqnumber{0}{10.}{0}\)\begin{align*} Z_1=& \sqrt {b_1}\sqrt {\dfrac {(n+1)(n+3)}{6(n-2)}}\\ =& 0,041948 \sqrt {\dfrac {(15+1)(15+3)}{6(15-2)}} =& 0,080605. \end{align*}
Como \(Z_1\) \(\le \) \(Z_{0,025}\) \(=\) 1,96, não se rejeita \(H_0\) e conclui-se que a distribuição é simétrica. Obs. O programa só vai aplicar esta aproximação se \(n<7\), que não é o caso.
-
2) Segunda alternativa
Esta é uma aproximação melhorada em relação à primeira alternativa. Para isso, é necessário calcular os seguintes valores
\(\seteqnumber{0}{10.}{0}\)\begin{align*} B=& \dfrac {3(n^2+27n-70)(n+1)(n+3)}{(n-2)(n+5)(n+7)(n+9)} \\ =& \dfrac {3(15^2+27\times 15-70)(15+1)(15+3)}{(15-2)(15+5)(15+7)(15+9)} = 3,5245, \\ C=& \sqrt {2(B-1)}-1= \sqrt {2(3,5245-1)}-1=1,2470,\\ D=& \sqrt {C} =\sqrt {1,2470} = 1,1167,\\ E=& \dfrac {1}{\sqrt {\ln (D)}} = \dfrac {1}{\sqrt {\ln (1,1167)}} = 3,0101,\\ F=& \dfrac {Z_1}{\sqrt {2/(C-1)}} = \dfrac {0,080605}{\sqrt {2/(1,1167-1)}} =0,028326. \end{align*}
Logo,
\(\seteqnumber{0}{10.}{0}\)\begin{align*} Z_{1c} =& E \ln \left (F+\sqrt {F^2+1}\right )\\ =& 3,0101\times \ln \left (0,028326+\sqrt {0,028326^2+1}\right )\\ =& 0,0853. \end{align*}
Como o valor \(|Z_{1c}|\) \(\le \) 1,96, a conclusão é, portanto, a mesma do item anterior, ou seja, não há desvios de assimetria na distribuição, considerando o nível nominal de significância de \(5\%\). Neste segundo caso (como poderia ter sido feito no primeiro) computamos o valor-\(p\) que foi 0,93206. Logo, como o valor-\(p\)>0,05, a hipótese nula não deve ser rejeitada. Esta aproximação é mais adequada e o valor da estatística difere do valor anterior, que é 0,0806.
-
-
ii) Hipótese sobre curtose
As hipóteses nula e alternativa são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} H_0: {\beta _2}=3 \qquad \qquad \textrm { versus } \qquad \qquad H_1: {\beta _2}\ne 3. \end{align*}
-
1) Primeira alternativa
A estatística do teste é
\(\seteqnumber{0}{10.}{0}\)\begin{align*} Z_2 =& \left (b_2-3+ \dfrac {6}{n+1}\right ) \sqrt {\dfrac {(n+1)^2(n+3)(n+5)}{24 n(n-2)(n-3)}} \\ =& \left (2,2116-3+ \dfrac {6}{15+1}\right ) \sqrt {\dfrac {(15+1)^2(15+3)(15+5)}{24\times 15(15-2)(15-3)}}\\ =& -0,52954. \end{align*}
Como \(|Z_2|\) \(<\) 1,96, então \(H_0\) não deve ser rejeitada, ou seja, a distribuição não pode ser considerada diferente da mesocúrtica, no nível nominal de significância de \(5\%\).
-
2) Segunda alternativa
Esta é uma opção melhorada do teste e a estimativa gama da curtose é
\(\seteqnumber{0}{10.}{0}\)\begin{align*} g_2 =& \dfrac {(n-1)\left [(n+1)b_2-3(n-1)\right ]}{(n-2)(n-3)} \\ =& \dfrac {(15-1)\left [(15+1)2,2116-3(15-1)\right ]}{(15-2)(15-3)}\\ =& -0,59356. \end{align*}
As demais quantidades necessárias para obter a melhor aproximação são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} H=& \dfrac {(n-2)(n-3)|g_2|}{(n+1)(n-1)} \sqrt {\dfrac {(n+1)^2(n+3)(n+5)}{24n (n-2)(n-3)}} \\ =& \dfrac {(15-2)(15-3)|-0,59356|}{(15+1)(15-1)} \sqrt {\dfrac {(15+1)^2(15+3)(15+5)}{24\times 15(15-2)(15-3)}}\\ =& 0,52954,\\ J=& \dfrac {6(n^2-5n+2)}{(n+7)(n+9)}\sqrt {\dfrac {6(n+3)(n+5)}{n (n-2)(n-3)}} \\ =& \dfrac {6(15^2-5\times 15+2)}{(15+7)(15+9)} \sqrt {\dfrac {6(15+3)(15+5)}{15 (15-2)(15-3)}}\\ =& 1,6595,\\ K=& 6 +\dfrac {8}{J}\left [\dfrac {2}{J}+\sqrt {1+\dfrac {4}{J^2}}\right ]\\ =& 6 +\dfrac {8}{1,6595} \left [\dfrac {2}{1,6595}+\sqrt {1+\dfrac {4}{1,6595^2}}\right ]\\ =& 19,3591,\\ L=& \dfrac {1-2K^{-1}}{1+H\sqrt {2(K-4)^{-1}}}\\ =& \dfrac {1-2\times 19,3591^{-1}}{1+0,52954 \sqrt {2(19,3591-4)^{-1}}}\\ =& 0,75283. \end{align*}
Portanto, a estatística do teste é
\(\seteqnumber{0}{10.}{0}\)\begin{align*} Z_{2c} =& \dfrac {1-2(9K)^{-1}-\sqrt [3]{L}}{\sqrt {2(9K)^{-1}}}\\ =& \dfrac {1-2(9\times 19,3591)^{-1}-\sqrt [3]{0,75283}}{\sqrt {2(9 \times 19,3591)^{-1}}}\\ =& 0,7357. \end{align*}
Como \(|Z_{2c}|\) \(<\) 1,96, então \(H_0\) não deve ser rejeitada, ou seja, a distribuição não pode ser considerada diferente da mesocúrtica.
-
-
iii) Teste conjunto de assimetria e curtose:
As hipóteses nula e alternativa para o teste conjunto de assimetria e curtose são
\(\seteqnumber{0}{10.}{0}\)\begin{align*} H_0: \sqrt {\beta _1}=0 \textrm { e } \beta _2=3 \textrm { versus } H_1: \sqrt {\beta _1}\ne 0 \textrm { ou } \beta _2\ne 3. \end{align*}
A estatística do teste é
\(\seteqnumber{0}{10.}{0}\)\begin{align*} \chi ^2_c =& 0,08525^2+0,73566^2\\ =& 0,5485. \end{align*} Considerando a distribuição qui-quadrado com \(\nu =2\) graus de liberdade, temos o seguinte valor-\(p\) \(=\) \(P(\chi ^2_{2}\ge \)0,5485) \(=\) 0,7602. Como o valor é inferior a 0,05, não devemos rejeitar \(H_0\) de distribuição simétrica e mesocúrtica. Assim, temos que não há evidências de desvios de normalidade considerando o nível nominal de \(5\%\). Este resultado, embora seja um teste indireto de normalidade, pois simetria e mesocurtose são condições necessárias, mas não suficientes, concordam com o teste de Shapiro-Wilk, anteriormente apresentado. Um protocolo que recomendamos é aplicar o teste conjunto. Se não houver significância, encerramos o processo e concluímos a favor da normalidade. Caso contrário, aplicamos individualmente cada teste (assimetria e curtose) para identificarmos o desvio ou os desvios de normalidade que ocorreu ou ocorreram.
O programa R para estes casos é:
# Function to apply the D'Agostino's test # for departure from normality, based on # skewness departure: improved DagSkewNormTest <- function(x) { n <- length(x) srb1 <- (sum((x-mean(x))^3)/n) / ((n-1)*var(x)/n)^1.5 z1 <- srb1 * ((n+1)*(n+3)/(6*(n-2)))^0.5 if (n < 7) z1c <- z1 else { B <- 3*(n^2+27*n-70)*(n+1)*(n+3) / ((n-2)*(n+5)*(n+7)*(n+9)) C <- sqrt(2*(B-1))-1 D <- sqrt(C) E <- 1 / sqrt(log(D)) FF <- z1 / sqrt(2 / (C - 1)) z1c <- E*log(FF +sqrt(FF^2 + 1)) } p.value <- 2 * (1 - pnorm(abs(z1c))) return(list(zc = z1c, p.value = p.value, Skewness=srb1)) } # Function to apply the D'Agostino's test # for departure from normality, based on # kurtosis departure: improved DagKurtNormTest <- function(x) { n <- length(x) b2 <- (sum((x-mean(x))^4)/n) / ((n-1)*var(x)/n)^2 g2 <- (n-1)*((n+1)*b2-3*(n-1))/((n-2)*(n-3)) z2 <- (b2 - 3 + 6 / (n + 1)) * ((n+1)^2*(n+3)*(n+5)/(24*n*(n-2)*(n-3)))^0.5 if (n < 4) z2c <- z2 else { H <- (n-2)*(n-3)*abs(g2)/((n+1)*(n-1)) * ((n+1)^2*(n+3)*(n+5)/(24*n*(n-2)*(n-3)))^0.5 J <- 6*(n^2-5*n+2)/((n+7)*(n+9)) * (6*(n+3)*(n+5)/ (n*(n-2)*(n-3)))^0.5 K <- 6 + 8/J*(2/J+sqrt(1+4/J^2)) L <- (1-2*1/K)/(1+H*sqrt(2/(K-4))) z2c <- (1-2/(9*K)-L^(1/3))/ sqrt(2/(9*K)) } p.value <- 2 * (1 - pnorm(abs(z2c))) return(list(zc = z2c, p.value = p.value, Kurtosis=b2, gKurt=g2)) } # Function to apply the D'Agostino's # ominibus test for Skewnes and Kurtosis # departure: indirect normality test OminbusNormTest <- function(x) { chi2c <- DagSkewNormTest(x)$zc^2 + DagKurtNormTest(x)$zc^2 nu <- 2 p.value <- 1 - pchisq(chi2c, nu) return(list(chi2c = chi2c, p.value = p.value)) } # Example: Exercise 10.5.1 x <- c(8.52,9.07,11.90,9.43,4.32, 5.77,8.52,5.95,5.56,9.02, 12.08,8.06,9.96,6.02,8.69) OminbusNormTest(x) DagSkewNormTest(x) DagKurtNormTest(x) -
-
c) Computando a correlação entre as colunas \(\hat {m}_i\) e \(X_{(i)}\) da tabela do item (i), temos \(r_Q\) \(=\) 0,97477. O valor crítico, obtido na Tabela 10.2 (Livro), com \(n=15\) e \(\alpha =\)0,05 é 0,9389. Como o valor calculado é maior que o tabelado, não há razão para rejeitar a hipótese nula de normalidade.
Para os procedimentos gráficos temos os seguintes resultados:
\(j\) \(P_j\) \(\hat {m}_j\) \(\hat {P}_j\) \(x_{(j)}\) 1 0,0410 -1,7394 0,0396 4,32 2 0,1066 -1,2450 0,1165 5,56 3 0,1721 -0,9458 0,1362 5,77 4 0,2377 -0,7137 0,1548 5,95 5 0,3033 -0,5150 0,1625 6,02 6 0,3689 -0,3349 0,4763 8,06 7 0,4344 -0,1651 0,5592 8,52 8 0,5000 0,0000 0,5592 8,52 9 0,5656 0,1651 0,5894 8,69 10 0,6311 0,3349 0,6464 9,02 11 0,6967 0,5150 0,6548 9,07 12 0,7623 0,7137 0,7128 9,43 13 0,8279 0,9458 0,7887 9,96 14 0,8934 1,2450 0,9536 11,90 15 0,9590 1,7394 0,9610 12,08 -
i) Para o Q-Q plot, basta que plotemos \(\hat {m}_j\) versus \(x_{(j)}\). O gráfico obtido está apresentado abaixo:
-
ii) Para o Q-Q plot alternativo, basta que plotemos \(P_j\) versus \(x_{(j)}\). O gráfico obtido está apresentado abaixo:
-
iii) Finalmente, o P-P plot é obtido plotando-se \(\hat {P}_j\) versus \(x_{(j)}\). O gráfico obtido está apresentado abaixo:
Os três procedimentos gráficos apresentaram gráficos com um bom ajuste da reta de mínimos quadrados à dispersão de pontos em cada caso, indicando, visualmente (daí tem subjetividade) que não há desvios de normalidade. As funções em R estão apresentadas na sequência:
# Function to obtain the normal Q-Q plots # and rQ, normal quantil-quantil correlation normQQPlot <- function(x) { n <- length(x) Pj <- ((1:n)-0.375) / (n + 0.25) mj <- qnorm(Pj) Xj <- sort(x) par(mfrow = c(1, 2)) plot(mj, Xj, xlab="m's",ylab="x's", main="Normal Q-Q Plot") abline(lm(Xj~mj),col="red") ######## Normal Probability Plot ##### plot(Pj, Xj, xlab="p's",ylab="x's", main="Normal Probability Plot") abline(lm(Xj~Pj),col="red") rQ <- cor(Xj,mj) return(list(rQ=rQ)) } # Function to obtain the normal P-P plots # and rQ, normal quantil-quantil correlation normPPPlot <- function(x) { n <- length(x) Pj <- ((1:n)-0.375) / (n + 0.25) Xj <- sort(x) Sn <- sqrt(var(x) * (n - 1) / n) Phj <- pnorm((Xj -mean(x))/Sn) plot(Phj, Pj, xlab="Est. P's",ylab="P's", main="Normal P-P Plot") abline(lm(Pj~Phj),col="red") } # Example: Exercise 10.5.1 x <- c(8.52,9.07,11.90,9.43,4.32, 5.77,8.52,5.95,5.56,9.02, 12.08,8.06,9.96,6.02,8.69) normQQPlot(x) normPPPlot(x) -
-
10.2.2 Vamos repetir o processo anterior para o novo conjunto de dados. Vamos escolher os procedimentos mais adequados e fornecer menos detalhes dos cálculos efetuados, uma vez que a maior parte dos detalhes necessários já foi apresentada no exercício anterior.
-
a) Teste de qui-quadrado:
A distribuição de frequências está apresentada a seguir, incluindo as frequências esperadas pelo modelo normal considerando uma média de 3,8397 e desvio padrão de 3,3154 (valores estimados na amostra):
----------------------------------------------------------------------- Classes (i) Ponto médio X(i) FO(i) FE(i) ----------------------------------------------------------------------- [-0,9268, 1,5286) 0,3009 12 9,7140 [ 1,5286, 3,9841) 2,7564 13 10,9798 [ 3,9841, 6,4396) 5,2119 7 10,6467 [ 6,4396, 8,8951) 7,6673 6 6,1125 [ 8,8951,11,3506) 10,1228 1 2,0763 [11,3506,13,8060) 12,5783 1 0,4697 -----------------------------------------------------------------------
O valor da estatística do teste é:
\(\seteqnumber{0}{10.}{0}\)\begin{align*} \chi ^2_c =& \sum _{i=1}^{k} \dfrac {\left (F_i-FE_i\right )^2}{FE_i}= \sum _{i=1}^{6} \dfrac {\left (F_i-FE_i\right )^2}{FE_i}\\ =& \dfrac {(12-9,7140)^2}{9,7140}+\dfrac {(13-10,9798)^2}{10,9798} + \dfrac {(7-10,6467)^2}{10,6467}+ \\ & \,\, +\dfrac {(6-6,1125)^2}{6,1125}+ \dfrac {(1-2,0763)^2}{2,0763}+\dfrac {(1-0,4697)^2}{0,4697}\\ =& 3,3175. \end{align*}
O valor tabelado é \(\chi ^2_{0,975;\nu =4}\) \(=\) 11,1433. Assim, como o valor calculado de qui-quadrado é menor que o tabelado, não há evidências significativas para rejeitar a hipótese nula. Este teste tem, em geral, um desempenho muito fraco, apresentando altas taxas de erros do tipo I e II. Neste caso, os dados foram simulados de uma população exponencial e a hipótese nula de normalidade, que é falsa, deveria ser rejeitada. Não foi, indicando que houve um erro do tipo II.
-
b) Usando a opção de Royston (1993), obtivemos \(W\) \(=\) 0,8656 e valor-\(p\) \(=\) 0,00022. Como o valor-\(p\) é menor que o nível nominal de \(5\%\), devemos rejeitar a hipótese nula de normalidade. O resultado levou a uma decisão correta e foi contrária a decisão do teste qui-quadrado, que como dissemos anteriormente tem baixo desempenho.
-
c) A correlação Q-Q plot foi \(r_Q\) \(=\) 0,9335. O valor tabelado para \(\alpha \) \(=\) \(5\%\) e \(n\) \(=\) \(40\) é 0,9726. Como o valor calculado é menor que o valor tabelado, devemos rejeitar a hipótese de normalidade no nível nominal de significância de \(5\%\). Este resultado concorda com o do teste de Shapiro-Wilk e discorda do de qui-quadrado.
Os procedimentos gráficos forneceram os seguintes resultados:
-
i) Gráficos Q-Q plot e probabilidade normal:
Há um indicativo de não normalidade dos dados.
-
ii) P-P plot:
Há um indicativo de não normalidade dos dados.
A comparação com os resultados dos testes de assimetria e curtose do exemplo 10.5 do Livro, indicam que os diversos testes concordam entre si. A exceção se deu com o teste qui-quadrado. O conjunto de dados e a chamada das funções apropriadas está apresentado a seguir para que o leitor possa reproduzir estes resultados ou, até mesmo, utilizar outra base de dados de seu interesse.
# Example 10.5 from the book x <-c(8.5209, 4.1871, 2.5163, 1.9133, 8.7796, 5.9117, 0.7608, 12.0372, 5.6255, 6.3604, 5.0679, 3.0310, 1.1281, 1.3850, 12.5783, 2.0292, 3.6014, 7.8288, 1.3829, 1.9344, 0.8642, 8.5144, 4.9774, 0.5759, 1.0414, 0.3009, 1.7809, 2.5638, 5.3587, 2.3066, 1.5298, 8.1052, 2.6035, 0.5953, 1.5033, 3.1507, 1.6889, 0.4451, 0.4750, 8.6276) normQQPlot(x) normPPPlot(x) -
-
10.2.3 O estimador de máxima verossimilhança do parâmetro \(\lambda \) do modelo exponencial é \(1/\bar {X}\). Embora sugerido o uso da estimativa de máxima verossimilhança, que no caso é \(1/\bar {X}\) \(=\) 1/3,83971 \(=\) 0,2604363, sabemos que a correlação é uma quantidade que é invariante a mudanças de escala e, portanto, não precisamos utilizar esta quantidade. Se, por exemplo, \(Y\sim \exp (1)\), então \(X\) \(=\) \(Y/\lambda \), \(\lambda >0\), tem distribuição exponencial com parâmetro \(\lambda \), ou seja, \(X\) \(\sim \) exponencial\((\lambda )\), que podemos provar facilmente usando o método do Jacobiano da transformação. Assim, a mudança de escala não influi na distribuição e, portanto, poderemos simular de uma distribuição exponencial padrão \((Y\sim \textrm {exp}(\lambda ))\), que a correlação \(r_Q\) não será afetada pela simples mudança de escala, quando consideramos uma exponencial com parâmetro \(\lambda \), de uma variável \(X\) exponencial qualquer.
-
(a) A proposta do teste é:
-
i) computamos \(p_j\) \(=\) \((j-0,375)/n(+0,25)\) para todo \(j=1\), \(2\), \(\cdots \), \(n\);
-
ii) ordenamos os valores de \(x\), \(x_{(j)}\), \(j=1\), \(2\), \(\cdots \), \(n\);
-
iii) estimamos a média das estatísticas de ordem da exponencial por \(m_j\) \(=\) \(F_Y^{-1}(p_j)\), em que \(F_Y^{-1}\) é a função de distribuição inversa de uma exponencial com \(\lambda =1\), considerando que \(Y\) \(\sim \) exponencial\((1)\), para todo \(j=1\), \(2\), \(\cdots \), \(n\);
-
iv) computamos a correlação amostral \(r_Q\) entre \(p_j\) e \(m_j\);
-
v) plotamos \(m_j\) (na abscissa) versus \(x_{(j)}\) (na ordenada) para obtermos o Q-Q plot;
-
vi) simulamos dados da distribuição exponencial com \(\lambda =1\) e mesmo tamanho amostral da amostra original. Repetimos, para cada amostra simulada os passos (ii), (iii) e (iv), com os mesmos valores de \(p_j\)’s do passo (i) e armazenamos os valores de \(r_Q\)’s obtidos. Denotamos estes valores simulados por \(r_{Q_i}\), para todo \(i=1\), \(2\), \(\cdots \), \(N\);
-
vii) calculamos o valor-\(p\) por
\(\seteqnumber{0}{10.}{0}\)\begin{align*} \textrm {valor-}p =& \dfrac {\displaystyle \sum _{i=1}^{N} I(r_{Q_j}\le r_Q)}{N}, \end{align*} em que \(I()\) é a função indicadora. Escolhemos \(N\) como sendo ao menos igual a \(1000\) simulações Monte Carlo;
-
viii) retornamos \(r_Q\) e o valor-\(p\).
Como o procedimento é essencialmente computacional, desenvolvemos uma função R para executar esta tarefa.
# Function to obtain the Q-Q plot dor exponential # distribution and to apply the Q-Q plot # correlation test expQQrQTest <- function(x, N = 10000) { n <- length(x) pj <- ((1:n)-0.375) / (n + 0.25) mj <- qexp(pj) xj <- sort(x) plot(mj, xj, xlab="m's",ylab="x's", main="Exponential Q-Q Plot") abline(lm(xj~mj),col="red") rQ <- cor(xj,mj) randomrQ <- function(n, pj) { xs <- rexp(n) msj <- qexp(pj) xsj <- sort(xs) rQs <- cor(xsj,msj) return(rQs) } nullRqDist <- matrix(c(n), N, 1) nullRqDist <- apply(nullRqDist,1,randomrQ,pj) p.value <- length(nullRqDist[nullRqDist<=rQ])/N return(list(rQ = rQ, p.value=p.value)) } # Example 10.5 from the book x <-c(8.5209, 4.1871, 2.5163, 1.9133, 8.7796, 5.9117, 0.7608, 12.0372, 5.6255, 6.3604, 5.0679, 3.0310, 1.1281, 1.3850, 12.5783, 2.0292, 3.6014, 7.8288, 1.3829, 1.9344, 0.8642, 8.5144, 4.9774, 0.5759, 1.0414, 0.3009, 1.7809, 2.5638, 5.3587, 2.3066, 1.5298, 8.1052, 2.6035, 0.5953, 1.5033, 3.1507, 1.6889, 0.4451, 0.4750, 8.6276) N <- 10000 expQQrQTest(x, N)Considerando \(N\) \(=\) \(10000\) simulações para este caso, obtivemos \(r_Q\) \(=\) 0,9967996 e valor-\(p\) \(=\) 0,5283. Assim, considerando o nível nominal de significância de \(5\%\) não devemos rejeitar a hipótese de distribuição exponencial para o conjunto de dados em questão. Sabemos que estes dados foram simulados de uma distribuição exponencial e, portanto, o resultado confirma o que já sabíamos previamente. O gráfico Q-Q plot deste caso é:
Podemos observar que há um bom ajuste da reta de mínimos quadrados pela dispersão de pontos, o que indica que a distribuição exponencial é adequada para este caso. Veja que fizemos este mesmo Q-Q plot para este conjunto de dados no exercício anterior, considerando como hipótese nula que os dados eram normais. O resultado foi que a hipótese de normalidade (visual e subjetivamente) devia ser rejeitada.
-
-
(b) Fizemos simulações Monte Carlo de duas formas: i) considerando dados de uma exponencial com parâmetro \(\lambda \) e amostra de tamanho \(n\). Como a correlação é invariante em relação a mudança de \(\lambda \), apenas o valor de \(n\) foi alterado. Assim, simulamos dados com \(\lambda =2\) e \(n\) \(=\) \(5\), \(50\), \(100\) e \(1000\); b) consideramos amostras da normal padrão, com os mesmos valores de \(n\) anteriores. Em ambos os caso \(N=2000\) valores foram usados para obter os valores-\(p\).
A proporção de rejeições de \(H_0\) foram computadas em ambos os casos em \(M=2000\) simulações Monte Carlo. No primeiro caso, situação (a), as rejeições indicam as taxas de erro tipo I e no segundo caso, situação (b), o poder do teste. Os resultados estão apresentados na tabela a seguir apenas para o nível \(\alpha \) \(=\) 0,05 de significância:
\(n\) Exponencial\((\lambda =2)\) Normal Padrão 5 0,0515 0,1895 50 0,0450 0,9620 100 0,0500 1,0000 1000 0,0515 1,0000 Os resultados mostram que o teste controla adequadamente a taxa de erro tipo I e possui poder elevado. As pequenas diferenças para o valor de \(5\%\) são atribuídas ao erro Monte Carlo. Claro que restringimos nosso estudo de poder a apenas uma distribuição, a normal. Um estudo mais amplo envolveria muitos níveis de significância, mas valores de \(n\) e outras distribuições sob a hipótese alternativa.
-