\documentclass{article}
\usepackage[utf8]{inputenc}
%% em en spanishdecimal porque atkinson tiene el ojo muy alto
\usepackage[spanish]{babel}\spanishdecimal{\text{\raisebox{.6em}{$,$}}}
\usepackage{amsmath,amssymb}
\usepackage[margin=2cm]{geometry}
\usepackage{atkinson}
\begin{document}
\title{Ejercicios de contrastes mediante NP, KR, RV}
\maketitle
\renewcommand{\L}{\mathcal{L}}
\renewcommand{\le}{\leqslant}
\renewcommand{\leq}{\leqslant}
\renewcommand{\ge}{\geqslant}
\renewcommand{\geq}{\geqslant}
\newcommand{\sigue}{\;\hookrightarrow\;}
\newenvironment{resoluciOn}{\begin{quotation}\sffamily}{\end{quotation}}
\newcommand{\uno}
  {Sea $X$ una variable aleatoria con distribución $B(5,p)$, de la
  que se extrae una muestra aleatoria simple de tamaño $30$ con los
  siguientes resultados:
  
  (3,2,3,4,1,4,3,4,4,3,4,2,1,2,4,1,1,0,3,2,2,4,2,2,4,3,2,3,1,1)
  \begin{enumerate}
  \item Aplica el método de Neyman-Pearson para contrastar las
    hipótesis $H_0:p=0.50$ y $H_1:p=0.70$ al nivel de significación
    $\alpha=0.05$ y calcula el p-valor del contraste.
  \item Contrasta $H_0: p \le 0.50$ y $H_1: p > 0.50$ al nivel de
    significación $\alpha=0.05$. Calcula el p-valor del contraste y
    representa la función potencia.
\end{enumerate}}
\newcommand{\dos}
{Sea $X$ una variable aleatoria con distribución
  geométrica, $G(p)$, que contabiliza el número de pruebas hasta que
  aparece el primer éxito y de la que se extrajo una muestra de tamaño
  $40$ con los siguientes resultados:
  \begin{center}
    1  2  4  5  6  7  8 10 12 13 14 16 17 19 24 29 30 31 39 
    7  5  4  1  1  1  3  3  4  1  1  1  1  2  1  1  1  1  1
    \end{center}
    \begin{enumerate}
\item Aplica el método de Neyman-Pearson para contrastar las hipótesis
  $H_0:p=0.50$ y $H_1:p=0.30$ al nivel de significación $\alpha=0.05$
  y calcula el p-valor del contraste.
\item Contrasta $H_0:p \ge 0.50$ y $H_1:p < 0.35$ al nivel de
  significación $\alpha=0.05$. Calcula el p-valor del contraste y
  representa la función potencia.
\end{enumerate}}
\newcommand{\tres}
{
  Sea $X$ una variable aleatoria con distribución
  beta $\beta(a,1)$ de la que se extrajo una muestra
  de tamaño $20$
  con los siguientes resultados:

  (0.96,0.83,0.49,0.35,0.76,0.76,0.92,0.98,0.43,0.91,%
   0.77,0.97,0.51,0.53,0.59,0.66,0.96,0.94,0.78,0.97)
\begin{enumerate}
\item Contrasta las hipótesis $H_0:a=2$ y $H_1:a=2.3$ al nivel de
  significación $\alpha=0.01$ y calcula el p-valor del contraste.
\item Contrasta $H_0:a \ge 2$ y $H_1: a <2$ al nivel de significación
  $\alpha=0.01$. Calcula el p-valor del contraste y representa la
  función potencia.
\end{enumerate}
}
\newcommand{\cuatro}
{
Considera una población $X\equiv N(\mu,\sigma)$ con
  $\mu$ y $\sigma$ desconocidas.
  Usa el test de la razón de verosimilitudes para contrastar:
  \begin{enumerate}
  \item $H_0: \mu=\mu_0$ frente a $H_1:\mu\neq\mu_0$.
  \item $H_0: \mu\le\mu_0$ frente a $H_1:\mu>\mu_0$.
  \item $H_0: \sigma\le\sigma_0$ frente a $H_1:\sigma>\sigma_0$.
  \end{enumerate}
}
\section{Enunciados}
\begin{enumerate}
\item \uno
\item \dos
\item \tres
\item \cuatro
\end{enumerate}
\section{Resoluciones}
\begin{enumerate}
\item \uno
  \begin{resoluciOn}
    \begin{displaymath}
      X\sigue B(5,p) \implies f={5\choose x}p^x(1-p)^{5-x}
    \end{displaymath}
    \begin{displaymath}
      \implies \L=\prod{5\choose x_i}p^{\sum x_i}(1-p)^{5n-\sum x_i}
      \propto p^{\sum x_i}(1-p)^{5n-\sum x_i}
    \end{displaymath}
    \begin{displaymath}
      \implies \L
      \propto \left(\frac p{1-p}\right)^{\sum x_i}(1-p)^{5n}
    \end{displaymath}
    Para el contraste $H_0:p=0.5$, $H_1:p=0.7$ se puede usar el lema
    de Neyman y Pearson:
    \[
      {\rm RC} = \left[\frac{\L_1}{\L_0}>k\right]
      = \left[\left(\frac{0.7/0.3}{0.5/0.5}\right)^{\sum x_i}>k'\right]
      = \left[\sum x_i>k''\right]
    \]
    \[ H_0\implies\sum X_i \sigue B(5n,0.5) \]
\begin{verbatim}
> x <- c(3,2,3,4,1,4,3,4,4,3,4,2,1,2,4,
         1,1,0,3,2,2,4,2,2,4,3,2,3,1,1)
> (n <- length(x))
[1] 30
> (pvalor <- 1 - pbinom(sum(x)-.5, 5*n, .5))
[1] 0.5325193
\end{verbatim}

    Para el contraste $H_0:p\leq0.5$, $H_1:p>0.5$ se puede usar el
    teorema de Karlin y Rubin considerando la razón de verosimilitud
    monótona en $T=\sum X_i$.  Sea $p_1<p_2$. Entonces
    \[
      \frac{\L(p_2)}{\L(p_1)}\propto
      \left(\frac{p_2/(1-p_2)}{p_1/(1-p_1)}\right)^t\quad\mbox{creciente en }t
    \]
    porque $\frac{p_2/(1-p_2)}{p_1/(1-p_1)}=\frac{p_2}{p_1}\frac{1-p_1}{1-p_2}
    > \frac{p_2}{p_1}>1$.
    Por tanto, la región crítica tiene la forma $[T>k]$ y el p-valor es
    el mismo que en el apartado anterior.

    Para representar la función potencia:
\begin{verbatim}
> alfa <- .05
> k <- qbinom(1-alfa, 5*n, .5)
> potencia <- function (p) 1-pbinom(k-.5, 5*n, p)
> library(txtplot)
> txtplot(pes, potencia(pes))
\end{verbatim}
    \begin{samepage}
\begin{verbatim}
    +--+---------+----------+----------+---------+----------+--+
  1 +                                   *********************  +
    |                                  **                      |
0.8 +                                  *                       +
    |                                 *                        |
    |                                 *                        |
0.6 +                                *                         +
    |                                *                         |
0.4 +                                                          +
    |                               *                          |
    |                               *                          |
0.2 +                              *                           +
    |                             **                           |
  0 +  ***************************                             +
    +--+---------+----------+----------+---------+----------+--+
       0        0.2        0.4        0.6       0.8         1
\end{verbatim}
  \end{samepage}
\begin{verbatim}
> pes <- seq(.45,.65,.01)
> txtplot(pes, potencia(pes))
\end{verbatim}
  \begin{samepage}
\begin{verbatim}
  1 +--+------------+------------+------------+------------+---+
    |                                              *  *  * *   |
    |                                            *             |
0.8 +                                         *                +
    |                                      *                   |
    |                                    *                     |
0.6 +                                 *                        +
    |                               *                          |
0.4 +                                                          +
    |                            *                             |
    |                         *                                |
0.2 +                       *                                  +
    |                 *  *                                     |
  0 +  * *  *  * *  *                                          +
    +--+------------+------------+------------+------------+---+
     0.45          0.5         0.55          0.6         0.65
\end{verbatim}
  \end{samepage}
\begin{verbatim}
> pes <- seq(.45,.65,.001)
> txtplot(pes, potencia(pes))
\end{verbatim}
  \begin{samepage}
\begin{verbatim}
  1 +--+------------+------------+------------+------------+---+
    |                                              *********   |
    |                                          *****           |
0.8 +                                       ****               +
    |                                     ***                  |
    |                                   ***                    |
0.6 +                                ***                       +
    |                              ***                         |
0.4 +                            ***                           +
    |                          ***                             |
    |                       ****                               |
0.2 +                    ****                                  +
    |                *****                                     |
  0 +  ***************                                         +
    +--+------------+------------+------------+------------+---+
     0.45          0.5         0.55          0.6         0.65   
\end{verbatim}
  \end{samepage}
\end{resoluciOn}
\item \dos
  \begin{resoluciOn}
    \noindent
\begin{eqnarray*}
  X\sigue G(p)&\implies& f(x\mid p) = (1-p)^{x-1}p
                         \cdot I(x\in\mathbb N=\{1,2,3,...\})\\
              &\implies& \L(\vec x\mid p) = (1-p)^{\sum x_i}(1-p)^{-n}p^n
                          \cdot I(\vec x\in\mathbb N^n)
\end{eqnarray*}
Omito a partir de ahora la indicatriz $I(\vec x\in\mathbb N^n)$.

Para $H_0:p=0.50$ y $H_1:p=0.30$ se tiene la región crítica
\begin{eqnarray*}
  \rm{R.C.}
  &=&\left\{\vec x\;\middle|\;
      \frac{\L(\vec x\mid0.3)}{\L(\vec x\mid0.5)}>c\right\}\quad\exists\,c\\
  &=&\left\{\vec x\;\middle|\;
      \frac{0.7^{\sum x_i}0.7^{-n}0.3^n}{0.5^{\sum x_i}0.5^{-n}0.5^n}>c\right\}
  \quad\exists\,c\\
  &=&\left\{\vec x\;\middle|\;
      \left(\frac{0.7}{0.3}\right)^{\sum x_i}>c'\right\}
  \quad\exists\,c'\\
  &=&\left\{\vec x\;\middle|\;
      \sum x_i>c''\right\}
      \quad\exists\,c''\\
  H_0&\implies& \sum_{i=1}^n X_i-n \sigue{\rm BN}(n,0.5)\\
  &\implies&c''=n+\left<\text{cuantil $1-\alpha$ de }{\rm BN}(n,0.5)\right>
             {}=n+53\\
  &\implies&\text{p-valor}=\Pr\left[{\rm BN}(n,0.5)\ge\sum x_i-n\right]
             \approx0<\alpha
\end{eqnarray*}
Se rechaza por tanto la hipótesis nula.
La forma de la región crítica tiene sentido porque $E(X)=1/p$ y, a menor
$p$, se esperan valores mayores de $X$.
\begin{verbatim}
> x <- scan()
1:     1  2  4  5  6  7  8 10 12 13 14 16 17 19 24 29 30 31 39 
20:    7  5  4  1  1  1  3  3  4  1  1  1  1  2  1  1  1  1  1
39: 
Read 38 items
> n <- length(x) # 38
> alfa <- 0.05
> qnbinom(1-alfa,n,0.5)
[1] 53
> sum(x)
[1] 327
> 1-pnbinom(327-n-0.5, n, 0.5)
[1] 0
\end{verbatim}
Para $H_0:p \ge 0.50$ y $H_1:p < 0.35$, consideramos primeramente el contraste
$H_0:p \ge 0.50$ y $H_1:p < 0.50$ para aplicar Karlin-Rubin. Sea $p_1<p_2$.
Entonces
$$
\frac{\mathcal L(p_2)}{\mathcal L(p_1)}=
\frac{(1-p_2)^{\sum x_i}(1-p_2)^{-n}p_2^n}{(1-p_1)^{\sum x_i}(1-p_1)^{-n}p_1^n}
\propto\left(\frac{1-p_2}{1-p_1}\right)^{\sum x_i}
\text{decreciente en }\sum x_i
$$
porque 
$\frac{1-p_2}{1-p_1}<1$, luego hay razón de verosimilitud
monótona\footnote{Podría considerarse también que es
creciente en $-\sum x_i$.} y la
región crítica del contraste uniformemente más potente tiene\footnote{
  Considérese:
\begin{itemize}
\item $H_1:p>p_0$
  \begin{itemize}
  \item $X\sigue B\implies\frac{\L_1}{\L_0}$ creciente en $\sum x_i$
    $\implies\text{R.C.}=[\sum x_i>c]$.
  \item $X\sigue G\implies\frac{\L_1}{\L_0}$ decreciente en $\sum x_i$
    $\implies\text{R.C.}=[\sum x_i<c]$.
  \end{itemize}
\item $H_1:p<p_0$, $X\sigue G$ $\implies$ $\text{R.C.}=[\sum x_i>c]$.
\end{itemize}
}
la forma $[\sum X_i>c]$.  La hipótesi alternativa del enunciado implica
un subconjunto de la del contraste K-R, por lo que nos sirve la misma región
crítica para disponer del contraste uniformemente más potente.
Por tanto, el p-valor sería
\begin{eqnarray*}
  \Pr\left[\sum X_i\ge\sum x_i\;\middle|\; p=0.5\right]
  &=&\Pr\left[\sum X_i-n\ge\sum x_i-n\;\middle|\; p=0.5\right]\\
  &=&\Pr\left[{\rm BN}(n,0.5)\ge\sum x_i-n\right]
  \approx0
\end{eqnarray*}
como en el apartado anterior.

La representación de la función potencia en el espacio paramétrico
$(0,0.35)\cup[0.5,1)$:
\begin{verbatim}
> pot <- function (p) 1-pnbinom(53.5,n,p)
> pes <- c(seq(0,.35,.01),seq(.5,1,.01))
> txtplot(pes, pot(pes))

    +-+----------+----------+----------+---------+----------+--+
  1 +  ******************                                      +
    |                    *                                     |
0.8 +                                                          +
    |                                                          |
    |                                                          |
0.6 +                                                          +
    |                                                          |
0.4 +                                                          +
    |                                                          |
    |                                                          |
0.2 +                                                          +
    |                                                          |
  0 +                            ****************************  +
    +-+----------+----------+----------+---------+----------+--+
      0         0.2        0.4        0.6       0.8         1
\end{verbatim}
\end{resoluciOn}
\item \tres
  \begin{resoluciOn}
    Para $H_0:a=2$, $H_1:a=2.3$ se usa Neyman-Pearson:
    $$X\sigue\beta(a,1)\implies f\propto x^{a-1}
    \implies\L\propto\prod x_i^{a-1}$$
    $$\frac{\L_1}{\L_0}\propto\prod x_i^{1.3}>c\iff\prod x_i>c'$$
    Tiene sentido porque, a mayor $a$, más se acerca la distribución a $1$.

    Cálculo aproximado del p-valor mediante Montecarlo:
\begin{verbatim}
> x <- c(0.96,0.83,0.49,0.35,0.76,0.76,0.92,0.98,0.43,0.91,
         0.77,0.97,0.51,0.53,0.59,0.66,0.96,0.94,0.78,0.97)
> n <- length(x)
> mean (replicate(1e6,prod(rbeta(n,2,1))) >= prod(x))
[1] 0.042949
> mean (replicate(1e6,prod(rbeta(n,2,1))) >= prod(x))
[1] 0.042787
\end{verbatim}
    por lo que p-valor${}\approx0.043$.

    Cálculo exacto del p-valor:
    $$\prod x_i>c'\iff\sum\ln x_i>c''\iff-\sum\ln x_i<c'''$$
    $$X\sigue\beta(a,1)\implies-\ln X\sigue{\rm Exp}(a)
    \implies-\sum\ln X_i\sigue\gamma(n,a)$$
\begin{verbatim}
> pgamma(-sum(log(x)),n,2)
[1] 0.04298876
\end{verbatim}
    luego a nivel $\alpha=0.01$ no se rechaza $H_0$.

    Para $H_0:a\geq2$, $H_1:a<2$ se usa Karlin-Rubin.
    Sea $a_1<a_2$. Entonces:
    $$
    \frac{\L(a_2)}{\L(a_1)}\propto\prod x_i^{a_2-a_1}
    \text{ es creciente en $\prod x_i$ porque }a_2-a_1>0
    $$
    luego R.C.${}=\bigl[\prod x_i<c\bigr]=
    \bigl[-\sum\ln x_i>c'\bigr]$ y el p-valor es:
\begin{verbatim}
> 1 - pgamma(-sum(log(x)),n,2)
[1] 0.9570112
\end{verbatim}
    mucho mayor que $\alpha$, luego no se rechaza $H_0$.
  \end{resoluciOn}
\item \cuatro
  \begin{resoluciOn}
    Para $H_0: \mu=\mu_0$ frente a $H_1:\mu\neq\mu_0$:
    $$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}$$
    $$\hat\mu_{\rm MV}=\bar x\qquad
    \hat\sigma_{\rm MV,H_0}^2=\frac1n\sum(x_i-\mu_0)^2=:s_0^2\qquad
    \hat\sigma_{\rm MV,H_1}^2=\frac1n\sum(x_i-\bar x)^2=s^2$$
    $$\Lambda=\frac{\L(\mu_0,s_0)}{\L(\bar x,s)}
    =\frac{\frac1{s_0^n\sqrt{2\pi}^n}e^\frac{-\sum(x_i-\mu_0)^2}{2s_0^2}}
    {\frac1{s^n\sqrt{2\pi}^n}e^\frac{-\sum(x_i-\bar x)^2}{2s^2}}
    =\left(\frac{s^2}{s_0^2}\right)^\frac n2
    =\left(\frac{s_0^2}{s^2}\right)^{-\frac n2}
    =\left(\frac{\sum(x_i-\mu_0)^2}{\sum(x_i-\bar x)^2}\right)^{-\frac n2}={}$$
    $${}=\left(\frac{\sum(x_i-\bar x+\bar x-\mu_0)^2}
      {\sum(x_i-\bar x)^2}\right)^{-\frac n2}
    ={\underbrace{\left(1+\frac{\sum(\bar x-\mu_0)^2} {\sum(x_i-\bar x)^2}
      \right)}_{>1}}^{-\frac n2}
    < c {}\iff1+\frac{\sum(\bar x-\mu_0)^2} {\sum(x_i-\bar x)^2}
    > c'{}\iff{}$$
    $$\iff\frac{\sum(\bar x-\mu_0)^2} {\sum(x_i-\bar x)^2}=
    \frac{n(\bar x-\mu_0)^2} {\sum(x_i-\bar x)^2}=
    \frac{\left(\frac{\bar x-\mu_0}{\sigma_0/\sqrt n}\right)^2}
    {\frac{\sum(x_i-\bar x)^2}{\sigma_0^2}}\stackrel{H_0}=
    \frac{N(0,1)^2}{\chi^2_{n-1}{\rm\; indep.}}
    > c''\iff$$
    $$\iff \left(\frac{N(0,1)}{\sqrt{\frac{\chi^2_{n-1}}{n-1}}}
    {\rm\;independientes}\right)^2 = t_{n-1}^2 = F_{1,n-1}
  > c'''$$
  por lo que equivale al habitual contraste $t$ bilateral
  para una muestra basado
  en $$\frac{\bar x - \mu_0}{s/\sqrt{n-1}}\sigue t_{n-1}$$

  Para $H_0: \mu\le\mu_0$ frente a $H_1:\mu>\mu_0$:
  $$
  \Lambda=
  \begin{cases}
    1&\bar x\leq\mu_0\\
    \left(\frac{s^2}{s_0^2}\right)^{\frac n2}&\bar x>\mu_0\\
  \end{cases}
  $$
  luego $$\Lambda<c\iff\frac{\bar x-\mu_0}{\hat s/\sqrt n}>t_{n-1,1-\alpha}$$
  De nuevo, el caso habitual: contraste $t$ unilateral.

  Para $H_0: \sigma\le\sigma_0$ frente a $H_1:\sigma>\sigma_0$:
  $$
  \Lambda=
  \begin{cases}
    1&s\leq\sigma_0\\
    \frac{\L(\bar x,\sigma_0)}{\L(\bar x,s)}&s>\sigma_0\\
  \end{cases}
  $$
  Ahora,
  $$
  \frac{\L(\bar x,\sigma_0)}{\L(\bar x,s)}=
 \frac{\frac1{\sigma_0^n\sqrt{2\pi}^n}e^\frac{-\sum(x_i-\bar x)^2}{2\sigma_0^2}}
 {\frac1{s^n\sqrt{2\pi}^n}e^\frac{-\sum(x_i-\bar x)^2}{2s^2}}
 =\left(\frac{s}{\sigma_0}\right)^n
 e^{-\frac12\sum(x_i-\bar x)^2\left(\frac1{\sigma_0^2}+\frac1{s^2}\right)}=
 \left(\frac{s^2}{\sigma_0^2}\right)^{\frac n2}
 e^{-\frac n2\frac{s^2}{\sigma_0^2}}e^{-\frac n2}
 =:g\left(\frac{s^2}{\sigma_0^2}\right)
 $$
 con $g(x)=x^{n/2}e^{-nx/2}e^{-n/2}$.
\begin{verbatim}
> n <- 10
> g <- function (x) x^(n/2) * exp(-n*x/2) * exp(-n/2)
> library(txtplot)
> x <- seq (0, 10, .1)
> txtplot (x, g(x))
      +--+---------+----------+---------+----------+---------+--+
      |      **                                                 |
4e-05 +      * *                                                +
      |        *                                                |
      |     *   *                                               |
3e-05 +          *                                              +
      |     *                                                   |
      |          *                                              |
2e-05 +           *                                             +
      |    *      *                                             |
      |            *                                            |
1e-05 +    *       **                                           +
      |   *         **                                          |
    0 +  **           ****************************************  +
      +--+---------+----------+---------+----------+---------+--+
         0         2          4         6          8        10   
> x <- seq (0, 3, .1)
> txtplot (x, g(x))
      +--+-------+--------+--------+--------+-------+--------+--+
      |                 * * *                                   |
4e-05 +                *      *                                 +
      |                        *                                |
      |              *           *                              |
3e-05 +                            *                            +
      |            *                                            |
      |                              *                          |
2e-05 +                               *                         +
      |          *                      *                       |
      |                                   *                     |
1e-05 +         *                           **                  +
      |       *                                * * *            |
    0 +  ** *                                       * * * ** *  +
      +--+-------+--------+--------+--------+-------+--------+--+
         0      0.5       1       1.5       2      2.5       3   
> 
\end{verbatim}
Comprobemos que tiene un único máximo en $x=1$:
\begin{verbatim}
Maxima 5.46.0 https://maxima.sourceforge.io
using Lisp ECL 21.2.1
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) g(x) := x^(n/2) * exp(-n*x/2) * exp(-n/2) $
(%i2) diff(g(x),x) ;
                              n x                      n x
                           (- ---) - n/2            (- ---) - n/2
                 n/2 - 1       2              n/2       2
              n x        %e                n x    %e
(%o2)         -------------------------- - ----------------------
                          2                          2
(%i3) solve (%, x) ;
                                         n - 2
                                         -----
                                  n/2      2
(%o3)                           [x    = x     ]

(%i4) factor(%o2);
                                               n x
                                            (- ---) - n/2
                                  n/2 - 1       2
                       n (x - 1) x        %e
(%o4)                - ----------------------------------
                                       2
(%i5) solve(%,x) ;
(%o5)                           [x = 0, x = 1]
\end{verbatim}
Por tanto, la función es creciente hasta $x=1$ y decreciente a partir de ahí,
por lo que
$$
\Lambda<c\iff \frac{s^2}{\sigma_0^2}<c_1\cup\frac{s^2}{\sigma_0^2}>c_2\iff
{s^2}<c_1\;\cup\;{s^2}>c_2
$$
y la región crítica coincide con la del contraste intuitivo basado en
$(n-1)\hat S/\sigma_0^2\sigue\chi^2_{n-1}$.
  \end{resoluciOn}
\end{enumerate}
\end{document}
