\message{ !name(bondad.Rnw.tex)}
\message{ !name(bondad.Rnw) !offset(964) }
\subsection{Shapiro y Wilk}
Se trata de una prueba muy potente, por lo que se ha popularizado como
la prueba de gausianidad por antonomasia\footnote{Es la única
  disponible en R básico.}:
<<>>=
X <- 1:60
shapiro.test (X)
@ 

Su idea es similar al de Shapiro y Francia pero tiene en cuenta
también las covarianzas entre los estadígrafos de orden.
Considérese una muestra gausiana típica y sea
$\vec Z_{(\cdot)}=(Z_{(1)},\dots,Z_{(n)})^t$ la correspondiente
muestra ordenada.  Sean sus esperanzas y covarianzas
$$
\vec m=E\begin{pmatrix}Z_{(1)}\\\vdots\\Z_{(n)}\end{pmatrix}
\quad V=\text{Var}\begin{pmatrix}Z_{(1)}\\\vdots\\Z_{(n)}\end{pmatrix}
\quad Z_i\equiv\mathcal N(0,1)\;\forall\,i=1,\dots,n
$$
Sea $\vec X_{(\cdot)}=(X_{(1)},\dots,X_{(n)})^t$ una muestra ordenada proveniente
de una población $X$.
Si $X$ es gausiana $\mathcal N(\mu,\sigma)$, entonces se puede expresar
$$
X_{(i)}=\mu+\sigma Z_{(i)}
$$
Estimando $\mu$ y $\sigma$ por mínimos cuadrados
generalizados\footnote{Si $\vec y=X\vec\beta+\vec\epsilon$ con
  $E(\vec\epsilon)=\vec0$ y $\text{Var}(\epsilon)=V$, entonces
  $\hat\beta=(X^tV^{-1}X)^tX^tV^{-1}\vec y$.}  hay que minimizar los
residuos
``mahalanobizados''\footnote{\url{https://es.wikipedia.org/wiki/Distancia_de_Mahalanobis}}
$$
[\vec x_{(\cdot)}-\mu\vec1-\sigma\vec z_{(\cdot)}]^t
V^{-1}
[\vec x_{(\cdot)}-\mu\vec1-\sigma\vec z_{(\cdot)}]
=
$$
$$
=\left(\vec x_{(\cdot)}-
\begin{bmatrix}\vec1&\vec z_{(\cdot)}\end{bmatrix}
\begin{bmatrix}\mu\\\sigma\end{bmatrix}\right)^t V^{-1}
\left(\vec x_{(\cdot)}-
\begin{bmatrix}\vec1&\vec z_{(\cdot)}\end{bmatrix}
\begin{bmatrix}\mu\\\sigma\end{bmatrix}\right)=
$$
\newcommand{\Z}{{\mathcal Z}}
$$
=
\left(\vec x_{(\cdot)}-\Z\vec\beta\right)^t V^{-1}
\left(\vec x_{(\cdot)}-\Z\vec\beta\right)
$$
con
$$
\Z = \begin{bmatrix}\vec1&\vec z_{(\cdot)}\end{bmatrix}\qquad
\vec\beta=\begin{bmatrix}\mu\\\sigma\end{bmatrix}
$$
La estimación de $\vec\beta$ es
\begin{eqnarray*}
  \hat\beta
  &=&(\Z^tV^{-1}\Z)^{-1}\Z^tV^{-1}\vec x_{(\cdot)}\\
  \begin{bmatrix}
    \hat\mu\\\hat\sigma
  \end{bmatrix}
  &=&\begin{bmatrix}
    \vec1^tV^{-1}\vec1&\vec1^tV^{-1}\vec m\\
    \vec m^tV^{-1}\vec1&\vec m^tV^{-1}\vec m\end
  {bmatrix}^{-1}
  \begin{bmatrix}
    \vec1^t\\\vec m^t
  \end{bmatrix}
  V^{-1}\vec y\\
  % &=&\begin{bmatrix}
  %   \vec1^tV^{-1}\vec1&\vec m^tV^{-1}\vec1\\
  %   \vec m^tV^{-1}\vec1&\vec m^tV^{-1}\vec m
  % \end{bmatrix}^{-1}
  % \begin{bmatrix}
  %   \vec1^t\\\vec m^t
  % \end{bmatrix}
  % V^{-1}\vec y\\
  &=&\frac{\begin{bmatrix}
    \vec m^tV^{-1}\vec m&-\vec 1^tV^{-1}\vec m\\
    -\vec m^tV^{-1}\vec1&\vec 1^tV^{-1}\vec 1
  \end{bmatrix}
  \begin{bmatrix}
    \vec1^t\\\vec m^t
  \end{bmatrix}}
{\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m-(\vec1^tV^{-1}\vec m)^2}V^{-1}\vec y\\
\left<\text{son escalares}\right>
  &=&\frac{\begin{bmatrix}
    \vec m^tV^{-1}\vec m&-\vec m^tV^{-1}\vec 1\\
    -\vec 1^tV^{-1}\vec m&\vec 1^tV^{-1}\vec 1
  \end{bmatrix}
  \begin{bmatrix}
    \vec1^t\\\vec m^t
  \end{bmatrix}}
  {\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m-(\vec1^tV^{-1}\vec m)^2}V^{-1}\vec y\\
   &=&\frac{\begin{bmatrix}
     \vec m^tV^{-1}\vec m \vec1^t -\vec m^tV^{-1}\vec1\vec m^t\\
     -\vec 1^tV^{-1}\vec m\vec1^t+\vec 1^tV^{-1}\vec 1\vec m^t
   \end{bmatrix}}
   {\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m-(\vec1^tV^{-1}\vec m)^2}V^{-1}\vec y\\
   &=&\frac{\begin{bmatrix}
     \vec m^tV^{-1}(\vec m \vec1^t - \vec1\vec m^t)V^{-1}\vec y\\
     \vec 1^tV^{-1}(\vec 1\vec m^t-\vec m\vec1^t)V^{-1}\vec y
   \end{bmatrix}}
 {\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m-(\vec1^tV^{-1}\vec m)^2}\\
 \hat\mu
 &=&\frac{\vec m^tV^{-1}(\vec m \vec1^t - \vec1\vec m^t)V^{-1}\vec y}
 {\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m-(\vec1^tV^{-1}\vec m)^2}\\
 \hat\sigma
 &=&\frac{\vec 1^tV^{-1}(\vec 1\vec m^t-\vec m\vec1^t)V^{-1}\vec y}
 {\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m-(\vec1^tV^{-1}\vec m)^2}
\end{eqnarray*}
     Para distribuciones simétricas\footnote{En tal caso, $\vec m$ es
     simétrica, $m_i=-m_{n-i+1}$, y $V$ es bisimétrica (simétrica y persimétrica),
     luego
     también
     lo es $V^{-1}=[v_{i,j}]_{i,j}$
     y $v_{i,j}$ = $v_{j,i}$ = $v_{n-j+1,n-i+1}$ = $v_{n-i+1,n-j+1}$.  Entonces
     $\vec1^tV^{-1}\vec m$ = $\sum_{i=1}^n\sum_{j=1}^nv_{i,j}m_j$ =
     $\sum_{j=1}^{\lfloor n/2\rfloor}m_j
     \left(\sum_{i=1}^n v_{i,j}-\sum_{i=1}^n v_{i,n-j+1}\right)$
     + $[\text{$n$ impar}]\,m_{\frac {n+1}2}
     \left(\sum_{i=1}^n v_{i,\frac{n+1}2}-\sum_{i=1}^n v_{i,\frac{n+1}2}\right)$ =
     $\sum_{j=1}^{\lfloor n/2\rfloor}m_j
     \left(\sum_{i=1}^n v_{i,j}-\sum_{i=1}^n v_{n-i+1,n-j+1}\right)$
     + $0$ =
     $\sum_{j=1}^{\lfloor n/2\rfloor}m_j
     \left(\sum_{i=1}^n v_{i,j}-\sum_{i=1}^n v_{i,j}\right)$ = $0$
     (véase
     \url{https://math.stackexchange.com/questions/648856/how-to-prove-that-the-inverse-of-a-persymmetric-matrix-is-also-persymmetric}).}
     se tiene $\vec1^tV^{-1}\vec m$ =
     $0$ y entonces
     $$
     \hat\sigma=\frac{\vec 1^tV^{-1}\vec 1\vec m^t V^{-1}\vec y}
     {\vec1^tV^{-1}\vec1\vec m^tV^{-1}\vec m}
     =\frac {\vec m^t V^{-1}\vec y} {\vec m^tV^{-1}\vec m}
     $$

     El estadígrafo $W$ del contraste fue definido originalmente \cite{SW}
     como la estimación de la pendiente $\sigma$ de la regresión arriba planteada
     multiplicada por una constante normalizadora (tipificadora):
     $$
     W=\hat\sigma^2\frac{(\vec m^tV^{-1}\vec m)^2}{\vec m^tV^{-1}V^{-1}\vec m\,n\,S^2}
     =\frac{(\hat\sigma\vec m^tV^{-1}\vec m)^2}{\vec m^tV^{-1}V^{-1}\vec m}\frac1{n\,S^2}
     =\frac {(\vec m^t V^{-1}\vec y)^2} {\lVert\vec m^tV^{-1}\rVert^2}\frac1{n\,S^2}
     =\left(\frac {\vec m^t V^{-1}} {\lVert\vec m^tV^{-1}\rVert} \vec y \right)^2 \frac1{n\,S^2}
     $$
     siendo $S^2=\frac1n\sum_{i=1}^n(X_i-\bar X)^2$ la varianza muestral.
     Habitualmente se expresa de alguna de las siguientes maneras:
\begin{eqnarray*}
  W
  &=& \frac{\left[\sum_{i=1}^{\lfloor\frac n2\rfloor}a^{(n)}_{n-i+1}(X_{(n-i+1)}-X_{(i)})\right]^2}
      {nS^2}\\
  &=& \frac{\left[\sum_{i=1}^na^{(n)}_{i}(X_{(i)}-\bar X)\right]^2} {nS^2}\\
  &=& \frac{\left[\sum_{i=1}^na^{(n)}_{i}X_{(i)}\right]^2} {nS^2}
\end{eqnarray*}
donde
$$
\vec a^{(n)}=
\begin{pmatrix}a_1^{(n)}\\\vdots\\a_n^{(n)}\end{pmatrix}=
\frac{V^{-1}\vec m}{\sqrt{\vec m^tV^{-1}V^{-1}\vec m}}
$$
Se verifican las siguientes propiedades:
\begin{itemize}
\item $a^{n}_i = -a^{n}_{n-i+1}$
\item $0<\frac{na_1^2}{n-1}\le W\le1$
\end{itemize}

Para su cálculo se puede usar una función de R o aproximar mediante Montecarlo:
<<>>=
options(width=70)
n <- 50; x <- 1:n
shapiro.test (x)
## B <- 1e8 # repeticiones Montecarlo
### exige mucha memoria:
## Z <- t (replicate (1e8, sort(rnorm(n))))
## m <- colMeans (Z)
## V1 <- solve (var (Z))
### versión más ligera:
## m <- numeric (n)
## V <- matrix (0, n, n)
## for (i in 1:B)
## {
##     xi <- sort (rnorm (n))
##     m <- m + xi/B # medias
##     V <- V + tcrossprod (xi, xi) / B # medias de productos
## }
### para ahorrar tiempo
## save (m, V, file="shapiro.rda")
load ("shapiro.rda") # ojo: sirve sólo para n=50
V1 <- solve (V - tcrossprod (m, m)) # inversa de matriz de cov.
a <- drop (V1 %*% m) / sqrt (drop (m %*% V1 %*% V1 %*% m))
W <- function (x) sum (a * sort(x)) ^ 2 / sum ((x - mean(x)) ^ 2)
W (x)
mean (replicate (1e5, W (rnorm (n))) <= W (x)) # P-valor
@ 
\message{ !name(bondad.Rnw.tex) !offset(-200) }
