\documentclass[mathserif]{beamer}
\usetheme{Madrid}
\usepackage[spanish]{babel}\spanishdecimal{\raisebox{1.27ex}{,}}
\usepackage{dsfont}
\usepackage{pict2e} % por \line size unavailable en picture
\usepackage{mleftright}\mleftright
\usepackage{pgfplots} % para tikzpicture
\usepackage{etoolbox} % para \ifstrempty
\usetikzlibrary{decorations.pathreplacing,calligraphy} % para llaves
\def\qedsymbol{} % quita el símbolo final de las demos
\addtobeamertemplate{theorem begin}{\normalfont}{} % que no esté en cursiva
\newtheorem*{lema}{Lema}
\newtheorem*{teorema}{Teorema}
\newtheorem*{notacion}{Notación}
\newtheorem*{ejemplo}{Ejemplo}
\newcommand{\Dn}{\ensuremath{\Delta_n}} % estadígrafo de KyS
\renewcommand{\P}[1]{\Pr\left[#1\right]} % probabilidad
\newcommand{\casiseguro}%
{\ensuremath{\xrightarrow[n\rightarrow\infty]{\text{c.s.}}}}
\newcommand{\abs}[1]{\left|#1\right|} % valor absoluto
\newcommand{\F}[1]{F\left(#1\right)} % función de distribución F
\newcommand{\Fi}[1]{F^{-1}\left(#1\right)} % ojiva inversa
\newcommand{\f}[1]{f\left(#1\right)} % función de densidad f
\newcommand{\G}[2]{F_{#1}\ifstrempty{#2}{}{\left(#2\right)}} % ojiva Fx
\newcommand{\g}[2]{f_{#1}\ifstrempty{#2}{}{\left(#2\right)}} % densidad fx
\newcommand{\Fn}[1]{F_n\left(#1\right)} % función de distribución empírica
\newcommand{\N}{\mathbb{N}} % naturales
\newcommand{\R}{\mathbb{R}} % recta real
\renewcommand{\o}[2]{{#1}_{(#2)}} % estadígrafo de orden: \o x 1 = x_{(1)}
\newcommand{\uno}[2]{\mathds{1}_{#1}\left(#2\right)}
\newcommand{\sigue}{\mathbin{\sim}}
\newcommand{\binomial}[2]{\mathcal B \left(#1,#2\right)}
\newcommand{\gausiana}[2]{\mathcal N \left(#1,#2\right)}
\newcommand{\multigausiana}[3]{\mathcal N_{#1} \left(#2,#3\right)}
\newcommand{\uniforme}[2]{\mathcal U \left(#1,#2\right)}
\newcommand{\exponencial}[1]{\mathrm{Exp} \left(#1\right)}
\renewcommand{\d}{\mathrm{d}}
\newcommand{\integral}[4]{\int_{#1}^{#2}{#3}\,\d{#4}\,}
\newcommand{\derivada}[2]{\frac{\d #1}{\d #2}}
\newcommand{\conjuntotalque}[2]{\left\{#1\;\middle|\;#2\right\}}
\newcommand{\primitiva}[3]{\left|#1\right|_{#2}^{#3}}
\renewcommand{\le}{\leqslant}
\renewcommand{\ge}{\geqslant}
\begin{document}
\title{Simulación}
\author{Inferencia Estadística}
\maketitle

\begin{frame}
  \begin{block}{Simulación}
    \begin{itemize}
    \item Útil para aproximar la distribución exacta de un
      estadígrafo.
    \item Se obtienen muchos valores aleatorios del estadígrafo.
    \item Teorema GC garantiza que su ojiva empírica se aproxima
      uniformemente a la teórica.
%    \item Con la ojiva empírica se pueden aproximar probabilidades
%      relacionadas con el estadígrafo.
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}
  \begin{block}{Generación de seudoaleatorios \(\uniforme{0}{1}\)}
    \begin{itemize}
    \item Paso básico en procedimientos de simulación.
    \item Implementada en casi todo lenguaje de programación y
      calculadora.
    \item R: \texttt{runif(n)}; Octave: \texttt{rand}; Maxima:
      \texttt{random(1.0)}; Python: \texttt{import random;
        random.random()}; Common Lisp: \texttt{(random 1.0)}; Fortran:
      \texttt{call random\_number(r)}
    \item Algoritmo básico: generación por congruencia
      \begin{itemize}
      \item $m$ módulo, $a$ multiplicador, $c$ incremento, $X_0$ semilla\\
        (por ejemplo $m=2^{32}$, $a=1\,664\,525$, $c=1\,013\,904\,223$,\\
        {\scriptsize
          \url{https://es.wikipedia.org/wiki/Generador_lineal_congruencial}})
      \item secuencia seudoaleatoria
        $$X_{i+1}=\left[(a\cdot X_i+c)\bmod m\right]$$ 
      \item \emph{seudo} porque puede reproducirse a partir de la semilla
        (en R, mediante \texttt{set.seed})
      \end{itemize}
    \item \texttt{RNGkind} en R informa o establece el algoritmo.
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}
  \begin{block}{Método de la ruleta (Montecarlo original)\\
    para variables discretas finitas}
  \begin{itemize}
  \item Sea $X$ que toma valores $x_1,\dots,x_k$ con $$\P{X=x_i}=p_i\qquad
    P_i=\sum_{j=1}^i p_i\qquad P_0=0$$
  \item Sea $C_i=[P_{i-1},P_i)$ luego $C_1,\dots,C_k$ partición de $[0,1)$ con
    $$\text{amplitud de }C_i=p_i$$
  \item Algoritmo:
    \begin{enumerate}
    \item Generar $u$ según una $\uniforme{0}{1}$.
    \item Si $u\in C_i$ tomar $x_i$ como valor aleatorio de $X$.
    \end{enumerate}
  \item En R se usan habitualmente \texttt{rbinom}, \texttt{rpois},
    \texttt{rgeom}, $\dots$
  \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{Método de la ruleta (ejemplo)}
    \begin{itemize}
    \item Sea $X$ discreta que toma valores $\{1,1.5,3\}$ con
      probabilidades respectivas \(\bigl\{\frac12,\frac13,\frac16\bigr\}\).
    \item Objetivo: generar muestra de tamaño $n=5$ de $X$.
    \end{itemize}
    \begin{enumerate}
    \item Partir el intervalo $[0,1]$ en $C_1=\bigl[0,\frac12\bigr]$,
      $C_2=\bigl(\frac12,\frac56\bigr]$ y $C_3=\bigl(\frac56,1\bigr]$.
    \item Generar $5$ valores de una $\uniforme{0}{1}$, por ejemplo
      \((0.08417384,0.62186518,0.75045832,0.76283675,0.29542565)\).
    \item Clasificar los valores generados: $(C_1,C_2,C_2,C_2,C_1)$.
    \item Trasformar las clases en valores de $X$: $(1,1.5,1.5,1.5,1)$.
    \end{enumerate}
    \begin{itemize}
    \item En R:
\begin{verbatim}
  U <- runif (5)
  C <- cut (U, c(0, 1/2, 5/6, 1))
  X <- c(1,1.5,3)[C]
\end{verbatim}
      o simplemente (\texttt{TRUE} significa \emph{con reposición}):
\begin{verbatim}
  X <- sample (c(1,1.5,3), 5, TRUE, c(1/2,1/3,1/6))
\end{verbatim}
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}
   \begin{block}{Método de la transformación inversa}
   \begin{itemize}
   \item Sea $X$ continua con ojiva $F$.
   \item El método se basa en que $\F{x}\sigue\uniforme{0}{1}$.
   \item Algoritmo:
     \begin{enumerate}
     \item Generar $u$ según una $\uniforme{0}{1}$.
     \item Tomar $x=F^{-1}(u)$ como valor aleatorio de $X$.
     \end{enumerate}
   \end{itemize}
 \end{block}
 \begin{ejemplo}
   \[
   \begin{array}{cccl}
     X\sigue\exponencial{\lambda}&\implies&\F{x}=1-e^{-\lambda x}&
     \text{con $x>0$}\\
      &\implies& \Fi{u}=\frac{-\ln(1-u)}\lambda&\text{con $0<u<1$}
   \end{array}
   \]
 \end{ejemplo}
 \begin{block}{}
   \begin{itemize}
   \item En R, se usan habitualmente \texttt{rexp}, \texttt{rnorm},
     \texttt{rgamma}, \texttt{rbeta}, \texttt{rchisq}, \texttt{rt},
     \texttt{rf}.
   \end{itemize}
 \end{block}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{}
    \scriptsize
\begin{verbatim}
options (OutDec="'")
graf <- function (n)
{
    plot (pexp, 0, 5, xlab="Exp(1)", ylab="U(0,1)", main=paste("n =",n))
    X <- c ()
    for (i in 1:n)
    {
        u <- runif (1)
        x <- qexp (u)
        arrows (c(0,x), c(u,u), c(x,x), c(u,0), col=2)
        X <- c (X, x)
    }
    if (n>100)
    {
        d <- hist (X, breaks = sqrt(n), plot=FALSE)
        lines (d$mids, d$density, col=3, lwd=3)
    }
}
pdf ("graf/montecarlo%1d.pdf", onefile=FALSE, width=10)
graf(1); graf(1); graf(1); graf(1)
graf(10); graf(10); graf(10)
graf(100)
graf(1000)
dev.off ()
\end{verbatim}
    \end{block}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo1.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo2.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo3.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo4.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo5.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo6.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo7.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo8.pdf}
\end{frame}

\begin{frame}
  \includegraphics[width=\textwidth]{graf/montecarlo9.pdf}
\end{frame}

\begin{frame}
  \pgfplotsset{width=\textwidth,height=\textheight}
  \begin{tikzpicture}
    \begin{axis}[xmin=-1.5,ymin=-0.39]
      \addplot[blue,domain=0:5,samples=100]{1-exp(-x)};
      \draw[red,->](axis cs:0,0.5)--(axis cs:0.6931472,0.5);
      \draw[red,->](axis cs:0.6931472,0.5)--(axis cs:0.6931472,0);
      \draw[decorate, decoration = {calligraphic brace}]
      (axis cs:0,0)--(axis cs:0,0.5);
      \draw[decorate, decoration = {calligraphic brace}]
      (axis cs:0,0.5)--(axis cs:0,1);
      \fill(axis cs:0,0.25)circle(0pt)
      node [anchor=east]{\rotatebox{90}{$\Pr=\frac12$}};
      \fill(axis cs:0,0.75)circle(0pt)
      node [anchor=east]{\rotatebox{90}{$\Pr=\frac12$}};
      \draw[decorate, decoration = {calligraphic brace}]
      (axis cs:0.6931472,0)--(axis cs:0,0);
      \draw[decorate, decoration = {calligraphic brace}]
      (axis cs:5,0)--(axis cs:0.6931472,0);
      \fill(axis cs:0.35,0)circle(0pt)
      node [anchor=north]{$\Pr=\frac12$};
      \fill(axis cs:2.85,0)circle(0pt)
      node [anchor=north]{$\Pr=\frac12$};
      \draw[pen colour=blue,decorate, decoration = {calligraphic brace}]
      (axis cs:-.6,0)--(axis cs:-.6,1);
      \draw[pen colour=red,decorate, decoration = {calligraphic brace}]
      (axis cs:5,-.15)--(axis cs:0,-.15);
      \fill(axis cs:-.6,.5)circle(0pt)
      node [anchor=east,blue]{\rotatebox{90}{$\uniforme{0}{1}$}};
      \fill(axis cs:2.5,-.15)circle(0pt)
      node [anchor=north,red]{$\exponencial{\lambda=1}$};
      \fill(axis cs:3.5,.8)circle(0pt) node [blue]{$u=F(x)=1-e^{-\lambda x}$};
      \fill(axis cs:.8,0)circle(0pt)
      node [red,anchor=south west]{$x=F^{-1}(u)=\frac{-\ln(1-u)}\lambda$};
      \fill(axis cs:5,0)circle(0pt)
      node [anchor=south east]{\dots};
    \end{axis}
  \end{tikzpicture}
\end{frame}

 \begin{frame}
   \begin{block}{Método de Box y Muller}
   \begin{itemize}
   \item Sean $U_1$ y $U_2$ dos variables $\uniforme{0}{1}$ independientes.
   \item Sean
     \begin{eqnarray*}
       X&=&\sqrt{-2\ln U_1}\cos(2\pi U_2)\\
       Y&=&\sqrt{-2\ln U_1}\sen(2\pi U_2)\\
     \end{eqnarray*}
   \item Entonces $X$ y $Y$ son $\gausiana{0}{1}$ independientes.
   \item Disponible en R mediante \texttt{RNGkind(normal.kind="Box-Muller")}\\
     pero por omisión se usa la transformación inversa $\Phi^{-1}$ =
     \texttt{qnorm}.
   \end{itemize}
 \end{block}
\end{frame}

 \begin{frame}
   \begin{block}{Método de Box y Muller (ejemplo)}
   \begin{itemize}
   \item Objetivo: Generar una muestra de tamaño $n$ de valores aleatorios de
     una distribución gausiana multivariante
     $\vec X\sigue\multigausiana{4}{\vec\mu}{\Sigma}$.
   \item Iterar $n$ veces:
     \begin{enumerate}
     \item Sean $u_1$, $u_2$, $u_3$ y $u_4$ valores $\uniforme{0}{1}$
       independientes.
     \item Sean
     \begin{eqnarray*}
       y_1=\sqrt{-2\ln u_1}\cos(2\pi u_2)&\qquad&
       y_3=\sqrt{-2\ln u_3}\cos(2\pi u_4)\\
       y_2=\sqrt{-2\ln u_1}\sen(2\pi u_2)&\qquad&
       y_4=\sqrt{-2\ln u_3}\sen(2\pi u_4)
     \end{eqnarray*}
   \item Box-Muller $\implies$ 
     $(y_1,y_2)\sigue(y_3,y_4)
     \sigue\multigausiana{2}{\vec0}{I}$
   \item $u_1,u_2,u_3,u_4$ indep. $\implies$ $\vec y=(y_1,y_2,y_3,y_4)'$ indep.\ 
     y $\vec y\sigue\multigausiana{4}{\vec 0}{I}$
   \item $\Sigma$ matriz cov.\ $\implies$ $\Sigma$ def. positiva
     $\implies$ $\Sigma=U\Lambda U'$, $\Lambda$ diag.\, $U$ ortog.
   \item Sea $\vec x =\vec\mu+U\Lambda^{\frac12}\vec y
     \sigue\multigausiana{4}{\vec\mu}{\Sigma}$ pues
     $\mathrm{cov}(\vec x)=
     U\Lambda^{\frac12}\mathrm{cov}(\vec y)\Lambda^{\frac12}U'=\Sigma$
   \item Entonces $\vec x$ es un valor generado según la distribución
     de $\vec X$.
   \end{enumerate}
 \item En R puede usarse: \texttt{mvtnorm::rmvnorm(n,mu,Sigma)}
 \end{itemize}
\end{block}
\end{frame}

\newcommand{\integralMC}[1]{
\begin{frame}
  \begin{block}{Método de aceptación y rechazo (intro: integración
      de Montecarlo)}    
    \centering Objetivo: $\integral{0}{1}{\sen\abs{\frac1x}}{x}$
    
    \includegraphics[width=.9\textwidth]{graf/lisp/monte#1.png}
  \end{block}
\end{frame}}
\integralMC{0}
\integralMC{1}\integralMC{2}\integralMC{3}
\integralMC{4}\integralMC{5}\integralMC{6}
\newcommand{\rechazas}[1]{
\begin{frame}
  \begin{block}{Método de aceptación y rechazo (intro: simple)}
      \centering Objetivo: generar aleatorios
      según $\f{x}\propto\sen\abs{\frac1x}\uno{(0,1)}{x}$
    
    \includegraphics[width=.8\textwidth]{graf/lisp/recha#1.png}
  \end{block}
\end{frame}}
\rechazas{0}
\rechazas{1}\rechazas{2}\rechazas{3}
\rechazas{4}\rechazas{5}\rechazas{6}

% \begin{frame}
%    \begin{block}{Método de aceptación y rechazo (simple)}
%    \begin{itemize}
%    \item \label{simple:comienzo}
%      Sea $X$ continua con densidad $f$ con soporte acotado $[a,b]$.
%    \item Sea $C$ tal que $f\le C\cdot g$.
%    \item Algoritmo:
%      \begin{enumerate}
%      \item \label{comienzo} Generar $x$ a partir de $\uniforme{a}{b}$.
%      \item Sea $u$ generado según $\uniforme{0}{C}$.
%      \item Si $u\le f(x)$ tomar $x$ como valor de $X$. Si no, volver a
%        empezar.
%      \end{enumerate}
%    \end{itemize}
%  \end{block}
% \end{frame}
\begin{frame}
  \begin{columns}
    \column{0.5\textwidth}
    \begin{block}{Método de aceptación y rechazo (simple)}
      \begin{itemize}
      \item \label{simple:comienzo}
        Sea $X$ continua con densidad $f$ con soporte acotado $[a,b]$.
      \item<2-> Sea $C$ tal que $f\le C$.
      \item<3-> Algoritmo:
        \begin{enumerate}
        \item<3-> \label{comienzo}
          \uncover<4-8,10->{Generar $x$} a partir de $\uniforme{a}{b}$.
        \item<5-8,11-> 
          \uncover<6-8,12->{Sea $u$ generado} según $\uniforme{0}{C}$.
        \item<7-8,13> \alert<13>{Si $u\le f(x)$ tomar $x$ como valor de $X$.}
          \alert<8>{Si no, volver a empezar.}
        \end{enumerate}
      \end{itemize}
    \end{block}
    \column{0.45\textwidth}
    \begin{block}{}
      \pgfplotsset{width=1.1\textwidth,height=0.7\textheight}
    \begin{tikzpicture}
    \begin{axis}[xmin=-2.5, xmax=2.5, ymin=-.1, ymax=0.6]
      \addplot[blue,<-,thick] expression[domain=-2.5:-2,samples=10]{0};
      \addplot[blue,thick]expression[domain=-2:0,samples=10]{(x+2)/4};
      \addplot[blue,thick]expression[domain=0:2,samples=10]{(2-x)/4};
      \addplot[blue,->,thick]expression[domain=2:2.5,samples=10]{0};
      \fill(axis cs:-0.8,0.43) node[blue]{$f$};
      \only<2->{\addplot[green,thick]expression[domain=-2:2,samples=10]{.5};
        \fill(axis cs:-2.25,0.5) node[green]{$C$};}
      \only<3,9>{\draw [decorate, decoration = {calligraphic brace}]
        (axis cs:2,0) -- (axis cs:-2,0);
        \fill(axis cs:0,-0.04) node{\scriptsize$\uniforme{a}{b}$};}
      \only<4-8>{\fill[black](axis cs:1.1,0)circle(2pt) node[anchor=west]{$x$};}
      \only<10->{\fill[black](axis cs:-.4,0)circle(2pt) node[anchor=west]{$x$};}
      \only<5,11>{
        \draw [decorate, decoration = {brace}]
        (axis cs:-2.49,.5) -- (axis cs:-2.49,0);
        \fill(axis cs:-2.4,0.25) node[anchor=west]
        {\scriptsize$\uniforme{0}{C}$};}
    \only<6-8>{\fill[black](axis cs:-2.4,0.38)circle(2pt)
      node[anchor=north west]{$u$};}
    \only<12->{\fill[black](axis cs:-2.4,0.15)circle(2pt)
      node[anchor=south west]{$u$};}
    \only<7-8>{\draw(axis cs:-2.4,0.38)--(axis cs:1.1,0.38);
      \draw(axis cs:1.1,0)--(axis cs:1.1,0.38);}
    \only<8>{\draw(axis cs:-2.4,0.22)--(axis cs:1.1,0.22);}
    \only<8>{\fill[black](axis cs:-2.4,0.22)circle(2pt)
      node[anchor=south west]{$f(x)$};}
    \only<13>{\draw(axis cs:-2.4,0.15)--(axis cs:-.4,0.15);
      \draw(axis cs:-.4,0)--(axis cs:-.4,0.4);
      \draw(axis cs:-2.4,0.4)--(axis cs:-.4,0.4);
      \fill[black](axis cs:-2.4,0.4)circle(2pt)
      node[anchor=north west]{$f(x)$};}
  \end{axis}
\end{tikzpicture}
\end{block}
\end{columns}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{Método de aceptación y rechazo (simple)}\small
\begin{verbatim}
f.tri <- function (x, a)      # densidad triangular en [-a,+a]
    1/a^2 * ((a - x) * (x >= 0) * (x < a) +
             (a + x) * (x <  0) * (x > -a))

AR <- function (n, f, a, b, cota) { # método acepta/rechaza
    muestra <- numeric(n)
    i <- 1
    while (i <= n) {
        x <- runif (1, a, b); y <- runif (1, 0, cota)
        if (y <= f(x)) {
            muestra[i] <- x
            i <- i + 1
        }
    }
    muestra
}

muestra <- AR (1e5, function (x) f.tri (x, 2), -2, +2, 1/2)
\end{verbatim}
  \end{block}
\end{frame}

\newcommand{\rechazag}[1]{
\begin{frame}
  \begin{block}{Método de aceptación y rechazo (intro: simple)}
      \centering Objetivo: generar aleatorios
      según $\f{x}\propto\sen\abs{\frac1x}\uno{(0,1)}{x}$
    
    \includegraphics[width=.8\textwidth]{graf/lisp/acepta#1.png}
  \end{block}
\end{frame}}
\rechazag{0}
\rechazag{1}\rechazag{2}\rechazag{3}
\rechazag{4}\rechazag{5}

\begin{frame}
  \begin{columns}
    \column{0.5\textwidth}
    \begin{block}{Método de aceptación y rechazo (general)}
      \begin{itemize}
      \item Sea $X$ continua con densidad $f$.
      \item Sea $Y$ con densidad $g$, tal que es fácil generar seudoaleatorios
        de $Y$ y existe $C$ tal que $f\le C\cdot g$.
      \item Algoritmo:
        \begin{enumerate}
        \item Generar $y$ a partir de $g$.
        \item Generar $u$ a partir de $\uniforme{0}{C\cdot g(y)}$.
        \item Si $u\le f(y)$ tomar
          $y$ como valor de $X$. Si no, volver a empezar.
        \end{enumerate}
      \end{itemize}
    \end{block}
    \column{0.45\textwidth}
    \begin{block}{}
      \pgfplotsset{width=1.1\textwidth,height=0.7\textheight}
      \begin{tikzpicture}
        \begin{axis}[xmin=-2.5, xmax=2.5, ymin=-.1, ymax=0.6]
          \addplot[blue,<-,thick] expression[domain=-2.5:-2,samples=10]{0};
          \addplot[blue,thick]expression[domain=-2:0,samples=10]{(x+2)/4};
          \addplot[blue,thick]expression[domain=0:2,samples=10]{(2-x)/4};
          \addplot[blue,->,thick]expression[domain=2:2.5,samples=10]{0};
          \fill(axis cs:+0.6,0.43) node[blue]{$f$};
          \only<2->{\addplot[green,thick]
            expression[domain=-2:-1,samples=10]{.25};
            \addplot[green,thick]expression[domain=-1:1,samples=10]{.5};
            \addplot[green,thick]expression[domain=1:2,samples=10]{.25};
            \fill(axis cs:-1.8,0.53) node[green]{$C\cdot g$};}
          \only<4>{\fill[black](axis cs:1.1,0)circle(2pt)
            node[anchor=west]{$x$};
            \draw[thick,blue](axis cs:1.1,0)--(axis cs:1.1,0.225);
            \draw[thick, red](axis cs:1.1,0.225)--(axis cs:1.1,0.25);}
          \only<5>{\fill[black](axis cs:-.4,0)circle(2pt)
            node[anchor=west]{$x$};
            \draw[thick,blue](axis cs:-.4,0)--(axis cs:-.4,0.4);
            \draw[thick, red](axis cs:-.4,0.4)--(axis cs:-.4,0.5);}
        \end{axis}
      \end{tikzpicture}
    \end{block}
  \end{columns}
\end{frame}

\begin{frame}
  \begin{block}{Método de A/R: otro ejemplo}
    \begin{itemize}
    \item $X\sigue\gausiana{0}{1} \implies
      \g{X}{x}=\frac1{\sqrt{2\pi}}e^{\frac{-x^2}2}\quad \G{X}{}$ no invertible
    \item $\g{X}{}$ campaniforme con máximo $\frac1{2\pi}$ en $0$
    \item $x\notin [-1,1]\implies x^2>\abs{x}\implies
      \g{X}{x}<\frac1{\sqrt{2\pi}}e^{\frac{-\abs{x}}2}$
    \item $\forall\,x\in\R,\,\g{X}{x}\le h(x)$ con
      \[
        h(x)=
        \begin{cases}
          \frac1{\sqrt{2\pi}}&x\in[-1,1]\\
          \frac1{\sqrt{2\pi}}e^{\frac{-\abs{x}}2}&x\notin[-1,1]
        \end{cases}
      \]
      \item $\integral{-\infty}{+\infty}{h(x)}{x}=
        \sqrt{\frac2\pi}(2e^{-\frac12}+1) \approx 1.765\,767=C<\infty$
      \item Sea $Y$ la variable auxiliar con densidad $g=\frac hC$
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}
  \pgfplotsset{width=\textwidth,height=\textheight}
  \begin{tikzpicture}
    \begin{axis}
      \addplot [samples=100,domain=-4:4,thick] {1/sqrt(2*pi)*exp(-x^2/2)};
      \addplot [blue,samples=100,domain=-4:-1,thick]
      {1/sqrt(2*pi)*exp(-abs(x)/2)};
      \addplot [blue,samples=100,domain=1:4,thick]
      {1/sqrt(2*pi)*exp(-abs(x)/2)};
      \addplot [blue,samples=100,domain=-1:1,thick] {1/sqrt(2*pi)};
      \fill(axis cs:0,0.2)circle(0pt) node{$f\sigue\gausiana{0}{1}$};
      \fill(axis cs:-2,0.4)circle(0pt) node[blue]{$h=C\cdot g$};
      \fill(axis cs:-3,.15)circle(0pt) node[blue]{$h=C\cdot g$};
      \fill(axis cs:3,.15)circle(0pt) node[blue]{$h=C\cdot g$};
%      \draw[blue,dashed](-1,exp(-1/2)/sqrt(2*pi))--(-1,1/sqrt(2*pi));
      \draw[blue,dashed](axis cs:-1,0.242)--(axis cs:-1,0.399);
      \draw[blue,dashed](axis cs:1,0.242)--(axis cs:1,0.399);
    \end{axis}
  \end{tikzpicture}
\end{frame}

\begin{frame}
  \begin{block}{Método de A/R: otro ejemplo}
    \begin{itemize}
      \item existe inversa $H^{-1}$ de la acumulada de $h$,
        $H(x)=\int_{-\infty}^xh$, para aplicar el método de la
        transformación inversa con $0<y<C$
      \end{itemize}
      \[
          H^{-1}(y)=
          \begin{cases}
            \ln\bigl(\frac\pi2y^2\bigr)&\text{si }0<y<\sqrt{\frac2{\pi e}}\\
            y\sqrt{2\pi}-\frac2{\sqrt e}-1&\text{si }\sqrt{\frac2{\pi e}}\le y
            \le \sqrt{\frac{2}{\pi e}}+\sqrt{\frac2\pi}\\
            -2\ln\bigl(1+\frac2{\sqrt e}-y\sqrt{\frac\pi2}\bigr)&\text{si }
            \sqrt{\frac2{\pi e}}+\sqrt{\frac2\pi}<y<
            \bigl(\frac2{\sqrt e}+1\bigr)\sqrt{\frac2\pi}
          \end{cases}
        \]
    \begin{itemize}
    \item alternativamente, existe inversa $G^{-1}$ de la acumulada de
      $g$, $G(x)=\int_{-\infty}^xg$, para aplicar el mismo método con
      $0<y<1$
      \end{itemize}
      \[
          G^{-1}(y)=
          \begin{cases}
            2\ln[y(2+\sqrt e)]-1&\text{si }0<y<\frac1{2+\sqrt e}\\
            2y-1+\frac{4y-2}{\sqrt e}&\text{si }\frac1{2+\sqrt e}\le y
            \le \frac{\sqrt e+1}{\sqrt e+2}\\
            -2\ln\bigl(\frac{2-2y}{\sqrt e}-y+1\bigr)&\text{si }
            \frac{\sqrt e+1}{\sqrt e+2}<y<1
          \end{cases}
        \]
  \end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{Obtención de $G^{-1}$ con Maxima (\texttt{Ginversa.mac})}
\begin{verbatim}
f(x)  := 1/sqrt(2*%pi) * %e^(-x^2/2) $ /* f < h  */

h1(x) := 1/sqrt(2*%pi) * %e^(x/2)    $ /* x < -1 */
h2(x) := 1/sqrt(2*%pi)               $ /* -1<x<1 */
h3(x) := 1/sqrt(2*%pi) * %e^(-x/2)   $ /* x < -1 */

h(x)  := if x<-1 then g1(x) 
         elseif x<=1 then g2(x) else g3(x) $

/* plot2d (h(x), [x,-3,+3]) $ */
\end{verbatim}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{Obtención de $G^{-1}$ con Maxima (\texttt{Ginversa.mac})}
\begin{verbatim}
H1(x) := integrate (h1(t), t, -inf, x) $
H2(x) := H1(-1) + integrate (h2(t), t, -1, x) $
H3(x) := H2(+1) + integrate (h3(t), t, +1, x) $
C     :  H3(inf) $ /* aprox. 1,765767458879439 */
G1(x) := H1(x)/C $ G2(x) := H2(x)/C $ G3(x) := H3(x)/C
G(x)  := if x<-1 then G1(x) 
         elseif x<=1 then G2(x) else G3(x) $
plot2d (G(x), [x,-5,+5]) $

G_1(y)  := rhs (solve (G1(x)=y, x) [1]) $ /* G1 inversa */
G_2(y)  := rhs (solve (G2(x)=y, x) [1]) $
G_3(y)  := rhs (solve (G3(x)=y, x) [1]) $
G_(y)   := if y<G(-1) then G_1(y) 
           elseif y<=G(+1) then G_2(y) else G_3(y) $
Ginv(y) := G_(y) $ /* inversa definida en [0;1] */
plot2d (Ginv(y), [y,0.001,0.999]) $
\end{verbatim}
\end{block}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{Obtención de $G^{-1}$ con Maxima (\texttt{Ginversa.mac})}
    Sesión normal con \texttt{display2d:true}
    
    \scriptsize
\begin{verbatim}
 (%i32) G(-1), ratsimp;                                            
                                    sqrt(%e)                       
 (%o32)                          ---------------                   
                                 %e + 2 sqrt(%e)                   
 (%i33) G(+1), ratsimp;                                            
                                  %e + sqrt(%e)                    
 (%o33)                          ---------------                   
                                 %e + 2 sqrt(%e)                   
 (%i34) G_1(y), ratsimp;                                           
 (%o34)                    2 log(sqrt(%e) y + 2 y) - 1             
 (%i35) G_2(y), ratsimp;                                           
                          sqrt(%e) (2 y - 1) + 4 y - 2             
 (%o35)                   ----------------------------             
                                    sqrt(%e)                       
 (%i36) G_3(y), ratsimp;                                           
                                        sqrt(%e)                   
 (%o36)                2 log(- --------------------------)         
                               2 y + sqrt(%e) (y - 1) - 2          
\end{verbatim}
\end{block}
\end{frame}

\begin{frame}[fragile]
\begin{block}{Obtención de $G^{-1}$ con Maxima (\texttt{Ginversa.mac})}
  Sesión con \texttt{displa2d:false}
  para usar las expresiones con R
  {\scriptsize
\begin{verbatim}
 (%i55) display2d : false $
 (%i56) G(-1), ratsimp ;
 (%o56) sqrt(%e)/(%e+2*sqrt(%e))                             
 (%i57) G(1), ratsimp ;                                  
 (%o57) (%e+sqrt(%e))/(%e+2*sqrt(%e))                    
 (%i58) G_1(y), ratsimp ; 
 (%o58) 2*log(sqrt(%e)*y+2*y)-1            
 (%i59) G_2(y), ratsimp ;                  
 (%o59) (sqrt(%e)*(2*y-1)+4*y-2)/sqrt(%e)  
 (%i60) G_3(y), ratsimp ;                                 
 (%o60) 2*log(-sqrt(%e)/(2*y+sqrt(%e)*(y-1)-2))           
\end{verbatim}}
Para evitar sustituir las\ \texttt{ \%e} se puede
usar\texttt{ fortran} :
{\scriptsize
\begin{verbatim}
 (%i61) fortran (ratsimp (G_3(y))) $
2*log(-exp(1.0E+0/2.0E+0)/(2*y+exp(1.0E+0/2.0E+0)*(y-1)-2))
\end{verbatim}}
\end{block}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{Obtención de $G^{-1}$ con Maxima (\texttt{Ginversa.mac})}
    Sesión con \texttt{tex} para la presentación en PDF\LaTeX\scriptsize
\begin{verbatim}
 (%i63) tex(ratsimp(G(-1))) $                                          
 $${{\sqrt{e}}\over{e+2\,\sqrt{e}}}$$                                  
 (%i64) tex(ratsimp(G(1))) $                                           
 $${{e+\sqrt{e}}\over{e+2\,\sqrt{e}}}$$                                
 (%i65) tex(ratsimp(G_1(y))) $                                         
 $$2\,\log \left(\sqrt{e}\,y+2\,y\right)-1$$                           
 (%i66) tex(ratsimp(G_2(y))) $                                         
 $${{\sqrt{e}\,\left(2\,y-1\right)+4\,y-2}\over{\sqrt{e}}}$$           
 (%i67) tex(ratsimp(G_3(y))) $                                         
 $$2\,\log \left(-{{\sqrt{e}}\over{2\,y+\sqrt{e}\,\left(y-1\right)-2}} 
  \right)$$                                                            
\end{verbatim}
\end{block}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{Obtención de $G^{-1}$ con Maxima (\texttt{Ginversa.mac})}
    Sesión con \texttt{tex} para la presentación en PDF\LaTeX\scriptsize

 $$G(-1)={{\sqrt{e}}\over{e+2\,\sqrt{e}}}$$                                  
 $$G(1)={{e+\sqrt{e}}\over{e+2\,\sqrt{e}}}$$                                
 $$G_1^{-1}(y)=2\,\log \left(\sqrt{e}\,y+2\,y\right)-1$$                           
 $$G_2^{-1}(y)={{\sqrt{e}\,\left(2\,y-1\right)+4\,y-2}\over{\sqrt{e}}}$$           
 $$G_3^{-1}=2\,\log \left(-{{\sqrt{e}}\over{2\,y+\sqrt{e}\,\left(y-1\right)-2}} 
 \right)$$
\end{block}
\end{frame}

\begin{frame}[fragile]
  \begin{block}{Obtención de $H^{-1}$ con Maxima}
\begin{verbatim}
logcontract (solve (H1(x)=y, x)) ;
expand (solve (H2(x)=y, x)) ;
solve (rootscontract (ratsimp (H3(x))) = y, x);
\end{verbatim}
  \end{block}
\end{frame}

\begin{frame}
   \begin{block}{Método de aceptación y rechazo (demostración 1)}
   \begin{itemize}
   \item Sea $X'=\left(Y\mid\text{aceptación}\right)$ y su densidad
     \[
       f_{X'}(x)=f(Y=x\mid\text{aceptación})=
       \frac{\P{\text{aceptación}\mid Y=x}\cdot g(x)}{\P{\text{aceptación}}}
     \]
   \item Como
     \begin{eqnarray*}
       \P{\text{aceptación}}
       &=& \integral{}{}{\P{\text{aceptación}\mid Y=y}\,g(y)}{y}\\
       &=& \integral{}{}{\frac{f(y)}{C\cdot g(y)}\,g(y)}{y}=\frac1C
     \end{eqnarray*}
   \item se tiene
     \[
       f_{X'}(x) = f(Y=x\mid\text{aceptación})=
       \frac{\frac{f(x)}{C\cdot g(x)}\cdot g(x)}{\frac1C}=f(x)
     \]
   \item luego $X'$ tiene la misma distribución que $X$
   \end{itemize}
 \end{block}
\end{frame}

 \begin{frame}[t]
    \begin{block}{Método de aceptación y rechazo (demostración 2)}
    \begin{itemize}
    \item Sea $X'=\left(Y\mid\text{aceptación}\right)$ y su
      función de distribución
      \[
        \G{X'}{x}=\P{Y\le x\mid\text{aceptación}}=
        \frac{\P{\text{aceptación}\cap Y\le x}}{\P{\text{aceptación}}}
      \]
    \item Sean $(U\mid Y=y)\sigue\uniforme0{C\cdot g(y)}$
      y ``aceptación'' = ``$U\le f(y)$''
    \end{itemize}
      \[
        \P{\text{aceptación}\cap Y\le x} =
        \integral{-\infty}{x} {\integral{0}{f(y)}{\g{U,Y}{u,y}}{u}}{y}=
      \]
      \[=\integral{-\infty}{x}
        {\integral{0}{f(y)}{\g{U\mid Y=y}{u}\,g(y)}{u}}{y}
      =\integral{-\infty}{x}
      {\integral{0}{f(y)}{\frac1{C\cdot g(y)}}{u}\,g(y)}{y}=\]
    \[=\integral{-\infty}{x}
      {\frac{f(y)}{C\cdot g(y)}\,g(y)}{y}=\frac{\G{X}{x}}C\]
%    \item Como
%      \begin{eqnarray*}
%        \P{\text{aceptación}}
%        &=& \integral{}{}{\P{\text{aceptación}\mid Y=y}\,g(y)}{y}\\
%        &=& \integral{}{}{\frac{f(y)}{C\cdot g(y)}\,g(y)}{y}=\frac1C
%      \end{eqnarray*}
%    \item se tiene
%      \[
%        f_{X'}(x) = f(Y=x\mid\text{aceptación})=
%        \frac{\frac{f(x)}{C\cdot g(x)}\cdot g(x)}{\frac1C}=f(x)
%      \]
%    \item luego $X'$ tiene la misma distribución que $X$
  \end{block}
 \end{frame}

\end{document}
