\documentclass[a4paper]{article}

\usepackage[spanish]{babel}
\usepackage[utf8]{inputenc}
\usepackage{fontenc}
\newcommand{\sigue}{\;\hookrightarrow\;}
\renewcommand{\L}{\mathcal L}
\renewcommand{\(}{\left(}
\renewcommand{\)}{\right)}
\renewcommand{\[}{\left[}
\renewcommand{\]}{\right]}
\newcommand{\RC}[1]{\left\{\vec x\;\middle|\;#1\right\}}
\usepackage{Sweave}
\begin{document}
\section{Enunciado}
$$X\sigue N(\mu,\sigma)$$
$$
H_0\colon\;\mu=\mu_0\,\cap\,\sigma=\sigma_0\qquad
H_1\colon\;\mu=\mu_1\,\cap\,\sigma=\sigma_1
$$
\section{Resolución}
$$
f=\frac1{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}
$$
$$
\L=\frac1{\sigma^n\sqrt{2\pi}^n}e^{-\frac{\sum(x_i-\mu)^2}{2\sigma^2}}
$$
\begin{eqnarray*}
  \mbox{R.C.}
  &=&\RC{\frac{\L_1}{\L_0}>k}=
      \RC{\(\frac{\sigma_0}{\sigma_1}\)^n
      e^{-\[\frac{\sum(x_i-\mu_1)^2}{2\sigma_1^2}-\frac{\sum(x_i-\mu_0)^2}{2\sigma_0^2}\]}>k}\\
  &=&\RC{{\frac{\sum(x_i-\mu_1)^2}{2\sigma_1^2}-
      \frac{\sum(x_i-\mu_0)^2}{2\sigma_0^2}}<k'}\\
\left\langle*\right\rangle
  &=&\RC{{\frac{\sum(x_i-\mu_1)^2}{\sigma_1^2}-
      \frac{\sum(x_i-\mu_0)^2}{\sigma_0^2}}<k''}\\
  &=&\RC{{{\sum(x_i-\mu_1)^2}-
      \frac{\sigma_1^2}{\sigma_0^2}{\sum(x_i-\mu_0)^2}}<k'''}\\
  \left<c=\frac{\sigma_1^2}{\sigma_0^2}\right>
  &=&\RC{{{\sum(x_i-\mu_1)^2}- c{\sum(x_i-\mu_0)^2}}<k'''}\\
  &=&\RC{{{\sum(x_i^2+\mu_1^2-2\mu_1x_i)}- c{\sum(x_i^2+\mu_0-2\mu_0x_i)}}<k'''}\\
  &=&\RC{\sum x_i^2+n\mu_1^2-2\mu_1\sum x_i-
      \(c\sum x_i^2+cn\mu_0-2c\mu_0\sum x_i\)<k'''}\\
  &=&\RC{\sum x_i^2+n\mu_1^2-2\mu_1\sum x_i- c\sum x_i^2-cn\mu_0+2c\mu_0\sum x_i<k'''}\\
  &=&\RC{(1-c)\sum x_i^2+n\mu_1^2+2(c\mu_0-\mu_1)\sum x_i-cn\mu_0<k'''}\\
  &=&\RC{(1-c)\sum x_i^2+2(c\mu_0-\mu_1)\sum x_i<k''''}
\end{eqnarray*}
Dados unos parámetros concretos, por ejemplo
\mathcode`.="0027 % apOstrofo decimal
$$
\mu_0=10\qquad\mu_1=12
\qquad\sigma_0=1\qquad\sigma_1=2\qquad n=25
\qquad\alpha=0.05
$$
la frontera de la R.C. puede aproximarse mediante montecarlo:
\begin{Schunk}
\begin{Sinput}
> mu0 <- 10
> mu1 <- 12
> sigma0 <- 1
> sigma1 <- 2
> c <- sigma1^2/sigma0^2
> alfa <- 0.05
> n <- 25
> nr <- 1e6 # número de repeticiones montecarlo
> d1 <- replicate (nr, # distribución montecarlo
+ {
+     x <- rnorm(n,mu0,sigma0)
+     (1-c)*sum(x^2)+2*(c*mu0-mu1)*sum(x)
+ })
> (k.... <- quantile(d1,alfa)) # k''''
\end{Sinput}
\begin{Soutput}
      5% 
6373.366 
\end{Soutput}
\end{Schunk}
%% Alternativamente, puede usarse la expresión $\left\langle*\right>$
%% donde
%% $$
%% \frac{\sum(x_i-\mu_0)^2}{\sigma_0^2}
%% =\frac{n S_0^2}{\sigma_0^2}\sigue\chi^2_n
%% \qquad
%% \frac{\sum(x_i-\mu_1)^2}{\sigma_1^2}
%% =\frac{n S_1^2}{\sigma_1^2}\sigue\chi^2_{n,d}
%% $$
%% con $S_j^2=\frac1n\sum (X_i-\mu_j)^2$ y siendo
%% $$
%% d=\frac{|\mu_0-\mu_1|}{\sigma_0}
%% $$
%% el parámetro de descentralización de una $\chi^2$:
%% <<>>=
%% d2 <- sigma0^2/sigma1^2 * rchisq(nr, n, abs(mu0-mu1)/sigma0) -
%%       rchisq(nr,n) # distribución montecarlo de <*>
%% (k.. <- quantile(d2,alfa)) # k''
%% @ 
%% Tipificando las dos distribuciones se comprueba que tienen la misma forma:
%% <<fig=TRUE>>=
%% s1 <- scale(d1)[,1]
%% s2 <- scale(d2)[,1]
%% plot(ecdf(s1), lwd=5)
%% plot(ecdf(s2), col=2, add=TRUE)
%% @ 
%% <<fig=TRUE>>=
%% plot(density(s1), lwd=5)
%% lines(density(s2), col=2)
%% @ 
\section{Casos particulares}
\subsection{$\sigma_0=\sigma_1\iff c=1$}
\begin{eqnarray*}
  \mbox{R.C.}
  &=&\RC{(1-c)\sum x_i^2+2(c\mu_0-\mu_1)\sum x_i<k''''}\\
  &=&\RC{(1-1)\sum x_i^2+2(\mu_0-\mu_1)\sum x_i<k''''}\\
  &=&\RC{2(\mu_0-\mu_1)\sum x_i<k''''}\\
  &=&\RC{(\mu_0-\mu_1)\sum x_i<k''''}\\
  \left<\mbox{supóngase }\mu_1>\mu_0\right>
  &=&\RC{\sum x_i>k^{(5)}}\\
  &=&\RC{\bar x>k^{(6)}}
\end{eqnarray*}
Se obtiene la misma región crítica que intuitivamente.
\subsection{$\mu_0=\mu_1$}
Supóngase $\sigma_1>\sigma_0$, luego $c>1$.
Intuitivamente se podría pensar en una R.C. basada en un estimador de $\sigma$:
$$
\RC{S^2>k_1}=\RC{\hat S^2>k_2}
$$
Aplicando Neyman-Pearson:
\begin{eqnarray*}
  \mbox{R.C.}
  &=&\RC{(1-c)\sum x_i^2+2(c\mu_0-\mu_1)\sum x_i<k''''}\\
  &=&\RC{(1-c)\sum x_i^2+2(c\mu_0-\mu_0)\sum x_i<k_4}\\
  &=&\RC{(1-c)\sum x_i^2+2(c-1)\mu_0\sum x_i<k_4}\\
  &=&\RC{(1-c)\sum x_i^2-2(1-c)\mu_0\sum x_i<k_4}\\
  \left<1-c<0\right>
  &=&\RC{\sum x_i^2-2\mu_0\sum x_i>k_5}\\
  &=&\RC{\sum x_i^2-2\mu_0\sum x_i+n\mu_0^2>k_6}\\
  &=&\RC{\sum (x_i^2-2\mu_0x_i+\mu_0^2)>k_6}\\
  &=&\RC{\sum (x_i-\mu_0)^2>k_6}\\
  &=&\RC{\frac{\sum (x_i-\mu_0)}n^2>k_7}\\
  &=&\RC{S_0^2>k_7}
\end{eqnarray*}
$S_0^2=\frac1n\sum (X_i-\mu_0)^2$ sería un estimador de la varianza cuando
la esperanza fuese conocida.  El estadístico asociado tendría distribución:
$$
\frac{n S_0^2}{\sigma^2}\sigue\chi^2_n\qquad\iff\qquad
S_0^2\sigue\frac{\sigma^2}n\chi^2_n
$$
Se gana un grado de libertad con respecto a la distribución asociada a la cuasivarianza:
$$
\frac{(n-1) \hat S^2}{\sigma^2}\sigue\chi^2_{n-1}\qquad\iff\qquad
\hat S^2\sigue\frac{\sigma^2}{n-1}\chi^2_{n-1}
$$
Eso supone obtener un contraste con más potencia; de hecho, por Neyman-Pearson,
se trata del contraste más potente posible.  Por ejemplo:
\begin{Schunk}
\begin{Sinput}
> mu0 <- 0; mu1 <- 0; sigma0 <- 5; sigma1 <- 10
> n <- 10; alfa <- 0.05
> ## RC1 = [cuasivarianza > k1]   RC2 = [S0^2 > k2]
> (k1 <- sigma0^2 * qchisq (1-alfa, n-1) / (n-1))
\end{Sinput}
\begin{Soutput}
[1] 46.99716
\end{Soutput}
\begin{Sinput}
> (k2 <- sigma0^2 * qchisq (1-alfa, n) / n)
\end{Sinput}
\begin{Soutput}
[1] 45.7676
\end{Soutput}
\begin{Sinput}
> 1 - pchisq ((n-1) * k1 / sigma1^2, n-1) # potencia RC1
\end{Sinput}
\begin{Soutput}
[1] 0.8956509
\end{Soutput}
\begin{Sinput}
> 1 - pchisq (n * k2 / sigma1^2, n) # potencia RC2
\end{Sinput}
\begin{Soutput}
[1] 0.9176019
\end{Soutput}
\end{Schunk}
\end{document}
