Complementa slides de GNA.

_site.yml

@@ -31,7 +31,7 @@ navbar:
         href: slides/02-docum-codigo.pdf
       - text: "GNA Uniformes"
         href: slides/03-gna-uniforme.pdf
-      - text: "GNA de v.a. Discretas"
+      - text: "GNA de v.a. não uniformes"
         href: slides/04-gna-nao-uniforme.pdf
       - text: "----------"
       - text: "Arquivos complementares"

index.Rmd

@@ -86,6 +86,7 @@ Estatística, está disponível em [Ementas DEST
     Karl Sigman.
   * <http://www.cs.cmu.edu/~tcortina/15-105sp09/Unit02PtC.pdf>;
   * [Bootstrap inference - Francisco Cribari-Neto](http://conteudo.icmc.usp.br/CMS/Arquivos/arquivos_enviados/SECAO-POSGRAD_87_bootstrap-slides.pdf).
+  * [Can You Behave Randomly?](http://faculty.rhodes.edu/wetzel/random/intro.html) - Robert R. Coveyou.
 
 # Atividades e avaliações
 

slides/04-gna-nao-uniforme.Rnw

@@ -41,7 +41,7 @@ source("config/setup.R")
 
   \begin{itemize}
   \item Mostrar a GNA de v.a. discretas.
-  \item Introduzir abordagens para GNA de v.a. contínuas.
+  \item Apresentar o método da transformação integral de probabilidades.
   \end{itemize}
 \end{frame}
 
@@ -213,6 +213,189 @@ curve(pgeom(x - 1, p = 0.5), add = TRUE, type = "s", col = 2)
 
 \end{frame}
 
+\begin{frame}[fragile]{Uma implementação que visa eficiência}
+<<eval = FALSE>>=
+random_geo <- function(n, p) {
+    q <- 1 - p
+    pq <- p * q
+    x <- 1
+    Fx <- p
+    u_vec <- sort(runif(n))
+    x_vec <- integer(n)
+    for (i in 1:n) {
+        while (u_vec[i] > Fx) {
+            x <- x + 1
+            Fx <- Fx + pq
+            pq <- pq * q
+        }
+        x_vec[i] <- x
+    }
+    x_vec <- sample(x_vec)
+    return(x_vec)
+}
+@
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+
+\begin{frame}{GNA de v.a. contínuas}
+
+  \begin{itemize}
+  \item A GNA de v.a. discretas é feito por busca direta no intervalo ao
+    qual o valor uniforme pertence na distribuição acumulada.
+  \item Em v.a. contínuas o mesmo raciocínio aplica-se.
+  \item Se a função de distribuição $F(x)$ tiver inversa $F^{-1}(u)$,
+    então a geração de números da v.a. $X$ com distribuição $F(x)$ é
+    muito simples.
+  \item Essa abordagem é conhecida como \textbf{método da transformação
+      integral da probabilidade} (MTIP).
+  \end{itemize}
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{O método da transformação integral da probabilidade}
+
+  \begin{itemize}
+  \item Considere que $X$ é uma v.a. contínua com função de distribuição
+    $F(x) = \Pr(X \leq x)$ para todo $x$.
+
+  \item Gera-se $U \sim \text{U}(0, 1)$.
+
+  \item O método da transformação integral da probabilidade estabelece
+    que a variável aleatória $W = F^{-1}(U)$ terá distribuição
+    $F(x)$. Ou seja, $W$ é na realidade $X$.
+
+  \item Para verificar isso, note que
+    \begin{equation}
+      \Pr(X \leq x) = \Pr(F^{-1}(U) \leq x) = \Pr(U \leq F(x)) = F(x)
+    \end{equation}
+    considerando que $F^{-1}(u) \leq x$ e $u \leq F(x)$ coincidem para
+    todo $u$ e $x$.
+  \end{itemize}
+
+  \vfill
+
+  \myurl{http://www.portalaction.com.br/simulacao-monte-carlo/metodo-da-transforma-inversa}\\
+  \myurl{https://en.wikipedia.org/wiki/Probability_integral_transform}
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{Distribuição exponencial}
+
+  Uma v.a. $X$ segue o modelo exponencial com parâmetro $\lambda > 0$ se sua densidade é dada por
+
+  \begin{equation}
+    f(x) = \lambda e^{-\lambda x} \cdot I(x > 0),
+  \end{equation}
+  sendo $\lambda > 0$. Denota-se por $X \sim \text{Exp}(\lambda)$.
+
+  A média e variância são
+  \begin{equation}
+    \text{E}(X) = \frac{1}{\lambda} \quad \text{e} \quad \text{V}(X) = \frac{1}{\lambda^2}.
+  \end{equation}
+
+\end{frame}
+
+\begin{frame}{Geração de números da Exponencial}
+  A função de distribuição da Exponencial é
+  \begin{equation}
+    F(x) = \int_{0}^{x} f(v)\, \text{d}v = 1 - \exp\{-\lambda x\}.
+  \end{equation}
+
+  Dessa forma, a inversa de $F(x)$ é
+  \begin{align*}
+    u &= 1 - \exp\{-\lambda x\}\\
+    1 - u &= \exp\{-\lambda x\}\\
+    \log(1 - u) &= -\lambda x\\
+    -\frac{\log(1 - u)}{\lambda} &= x\\
+    \therefore \quad x &= F^{-1}(u) = -\frac{\log(1 - u)}{\lambda}.
+  \end{align*}
+
+\end{frame}
+
+\begin{frame}[fragile, allowframebreaks]{Implementação em R}
+
+<<>>=
+random_exp <- function(n, lambda) {
+    u <- runif(n)
+    x <- -log(1 - u)/lambda
+    return(x)
+}
+@
+
+<<ecdf_exp, eval = FALSE>>=
+x <- random_exp(n = 1000, lambda = 0.5)
+
+Fx <- function(x, lambda) {
+    1 - exp(-lambda * x)
+}
+
+plot(ecdf(x))
+curve(Fx(x, lambda = 0.5), add = TRUE, col = 2, from = 0)
+legend("right", legend = c("Teórica", "Empírica"),
+       lty = 1, col = 1:2, bty = "n")
+@
+
+  \framebreak
+
+<<echo = FALSE, eval = TRUE, fig.dim = c(4, 4), out.width = "0.7\\linewidth">>=
+<<ecdf_exp>>
+@
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{Exercícios}
+
+  Obter o GNA das distribuições especificadas abaixo usando o MTIP.
+
+  \begin{itemize}
+  \item $F(x) = 1-(x-1)^{2}, \quad 0 < x < 1.$
+  \item $F(x) = x^{\theta}, \quad 0 < x < 1, \theta > 1.$
+  \item $F(x) = 1-\exp\{-x^2/2\tau^2\}, \quad x > 0, \tau > 0.$
+  \item $f(x) = \frac{1}{2} \sin(x), \quad 0 < x < \pi.$
+  \item $f(x) = \frac{\sec^2(x)}{\sqrt{3}} , \quad 0 < x < \pi/3.$
+  \end{itemize}
+
+  Dica: use programas de matemática simbólica para verificar os resultados.
+
+  Wolfram Alpha: \url{http://www.wolframalpha.com/}.
+\end{frame}
+
+\begin{frame}{Solução}
+
+    \begin{itemize}
+    \item $F^{-1}(u) = 1 - \sqrt{1 - u}$ ou apenas
+      $F^{-1}(u) = 1 - \sqrt{u}$ por simetria.
+  \item $F^{-1}(u) = u^{1/\theta}$
+  \item $F^{-1}(u) = \sqrt{-2\tau\cdot\log(1-u)}.$
+  \item $F(x) = \frac{1-\cos(x)}{2}$ e $F^{-1}(u) = \arccos(1 - 2u)$.
+  \item $F(x) = \frac{\tan(x)}{\sqrt{3}}$ e $F^{-1}(u) = \arctan(u\sqrt{3})$.
+  \end{itemize}
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{Quando usar o MTIP para GNA?}
+
+  \hi{Sempre que possível.}
+
+  \begin{itemize}
+  \item O método funciona quando $F^{-1}(u)$ existe.
+  \item Existem distribuições onde $F(x)$ é não inversão analítica.
+  \item Existem distribuições onde $F(x)$ não existe (apenas $f(x)$).
+  \end{itemize}
+
+  Alternativas
+
+  \begin{itemize}
+  \item Pode-se aproximar numericamente $F(x)$ e/ou $F^{-1}(u)$ para GNA.
+  \item Pode-se considerar métodos não baseados em $F^{-1}(u)$.
+  \end{itemize}
+\end{frame}
+
 % \begin{frame}
 %   Formalização dos conceitos
 %   \url{http://www.portalaction.com.br/simulacao-monte-carlo/metodo-da-transforma-inversa}
@@ -233,15 +416,15 @@ curve(pgeom(x - 1, p = 0.5), add = TRUE, type = "s", col = 2)
 
       \hi{Próxima aula}
       \begin{itemize}
-      \item GNA de variáveis aleatórias contínuas.
-      \item Método da transformação integral da probabilidade.
+      \item O método da aceitação-rejeição.
       \end{itemize}
 
       \hi{Avisos}
       \begin{itemize}
-      \item Sabatina 02 já está no Moodle!
+      \item Sabatina 03 já está no Moodle!
       \end{itemize}
 
+      \vspace{3em}
     \end{minipage}
 
 \end{frame}

slides/config/setup.R

 library(knitr)
-opts_chunk$set(dev.args = list(family = "Palatino"))
+opts_chunk$set(dev.args = list(family = "Palatino"),
+               fig.align = "center")
 # thm <- knit_theme$get("dusk")
 thm <- knit_theme$get("solarized-light")
 knit_theme$set(thm)