Método da Transformação Integral da Probabilidade

Gerando números da distribuição exponencial

A função densidade de uma v.a. Exponencial é \[ f(x, \lambda) = \lambda\exp\{-\lambda x\}, \quad x>0, \lambda>0. \]

A função densidade de probabilidade acumulada é \[ F(x, \lambda) = 1-\exp\{-\lambda x\}. \]

A inversa da função \(F\), \(F^{-1}\), é portanto \[ F^{-1}(u) = -\frac{\log(1-u)}{\lambda}. \]

Pelo teorema da transformação integral da probabilidade, podemos gerar números aleatórios de variáveis aleatórias por meio da inversa da função de distribuição \(F(x)\).

##----------------------------------------------------------------------
## Simulando valores da Exponencial a partir de Uniformes.

Fx <- function(x, lambda){
    1-exp(-lambda*x)
}

iFx <- function(u, lambda){
    -log(1-u)/lambda
}

## Números uniformes. Poderia usar a rand0() da aula passada.
u <- runif(1000)

## Números de distribuição Exponencial.
x <- iFx(u, lambda=0.5)
head(x)

[1] 1.3846903 2.2494784 4.0682154 5.2400140 0.3723905 1.7382503

## Histograma? Não tem nada melhor não?
hist(x, prob=TRUE)
curve(0.5*exp(-0.5*x), add=TRUE, col=2)
rug(x)

## Hum! Pode ficar melhor?
plot(density(x, from=0))
curve(0.5*exp(-0.5*x), add=TRUE, col=2)
rug(x)

## Bem melhor.
plot(ecdf(x))
curve(Fx(x, lambda=0.5), add=TRUE, col=2, from=0)

## Pra ficar melhor ainda, fazer um PP-plot.
Pobs <- (1:length(x))/length(x) ## Freq. rel. acum. observadas.
Pteo <- Fx(sort(x), lambda=0.5) ## Prob acum. teóricas.

plot(Pobs~Pteo,
     xlab="Frequências relativas acumuladas teóricas",
     ylab="Probabilidades acumuladas observadas")
rug(x=Pteo, side=2)
rug(x=Pobs)

##----------------------------------------------------------------------
## O que acontece se usarmos um mau gerador de números uniformes, como o
## discutido na aula anterior?

u <- numeric(1001)
u[1] <- 0.33
for (i in 2:1000){
    u[i] <- 0.5+0.5*sin(2*pi*x[i-1])
}
u <- u[-1]

## Repita os gráficos da sessão anterior.

Gerando números aleatórios por meio da relação entre v.a.

Normal por Uniformes

##----------------------------------------------------------------------
## Simulando números da distribuição Normal tirando vantagem do Teorema
## do Limite Central.

## A soma de 6 Uniformes independentes retorna uma distribuição muito
## próxima da Normal.

## 6 mil uniformes.
u <- matrix(runif(6*1000), ncol=6)

x <- rowSums(u)
qqnorm(x)

Os valores gerados apresentam (satisfatória) distribuição Normal, porém não é uma Normal padrão. Como transformar em Normal padrão? Tente fazer isso. Dica: procure conhecer a média e variância de uma distribuição Uniforme.

Binomial por Bernoulli

##----------------------------------------------------------------------
## Simulando da binomial por meio da uniforme -> Bernoulli.

p <- 0.3 ## Prob de sucesso.
n <- 6   ## Número de ensaios.
u <- matrix(runif(n*1000), nrow=1000)

x <- rowSums(u < p)

plot(ecdf(x), cex=0)
curve(pbinom(x, size=n, prob=p),
      type="s", add=TRUE, col="2")

##----------------------------------------------------------------------
## Simulando também pela "inversa" da acumulada.

## Intervalos para classificar os números uniformes.
iFx <- c(0, pbinom(0:n, size=n, prob=p))

u <- runif(1000)

## Classificando.
x <- (findInterval(x=u, vec=iFx))-1

plot(ecdf(x), cex=0)
curve(pbinom(x, size=n, prob=p),
      type="s", add=TRUE, col="2")

Qual dos dois métodos para gerar valores de uma binomial é mais interessante? Em qual deles o custo é maior? Imagine a situação em que o \(n\) é grande como 100.

Poisson pela Exponencial

##----------------------------------------------------------------------
## Simulando a Poisson por meio da Exponencial.

## O número de eventos dentro de um intervalo unitário é uma Poisson se
## o tempo para ocorrência de um evento tiver distribuição Exponencial.

## Gerando números da Exponencial.
iFx <- function(u, lambda){
    -log(1-u)/lambda
}

lambda <- 2
y <- iFx(runif(5000), lambda=lambda)

## Acumular.
y <- cumsum(y)

## Truncar para o inteiro acima.
y <- ceiling(y)

## Contar ocorrências em cada intervalo unitário.
x <- tabulate(y)
length(x)

[1] 2589

plot(ecdf(x), cex=0)
curve(ppois(x, lambda=lambda),
      type="s", add=TRUE, col="2")

Qual o problema desse gerador de Poisson? Teremos o mesmo número de números Poisson se um valor de \(\lambda\) maior for usado? Como deve ser a implementação para gerar exatos 1000 números?

Aproximação da inversa por finitos segmentos de reta

##----------------------------------------------------------------------
## Usando recursos númericos para simular da normal via método da
## transformação integral da probabilidade.

## Alguns pontos no domínio da normal padrão.
x <- seq(-3,3,l=6)

## A função acumulada da normal.
Fx <- pnorm(x)

## Aproximação e real.
plot(Fx~x, type="o")
curve(pnorm(x), add=TRUE, col=2)

## Aproximação da F(x) por segmentos de reta.
aFx <- approxfun(x=x, y=Fx)

xi <- seq(-3,3, by=0.1)
Fxi <- aFx(xi)
plot(Fxi~xi)

## Aproximação da inversa de F(x) por segmentos de reta.
aiFx <- approxfun(x=Fx, y=x)

u <- runif(1000)
x <- aiFx(u)

plot(ecdf(x), cex=0.1)
curve(pnorm(x), add=TRUE, col=2)

## Valores NA. De onde vieram?
c(sum(is.na(u)), sum(is.na(x)))

[1] 0 5

Encontre uma explicação para a quantidade de valores NA obtidos. Tente aumentar o número de intervalos e verifique se existe melhoria quanto à distribuição de probabilidades. Amplie também os limites usados.

options(width=90)
devtools::session_info()
Sys.time()

 setting  value                       
 version  R version 3.2.2 (2015-08-14)
 system   x86_64, linux-gnu           
 ui       X11                         
 language en_US                       
 collate  en_US.UTF-8                 
 tz       Brazil/East                 
 date     2015-10-15                  

 package      * version date       source        
 devtools       1.9.1   2015-09-11 CRAN (R 3.2.2)
 digest         0.6.8   2014-12-31 CRAN (R 3.2.2)
 evaluate       0.8     2015-09-18 CRAN (R 3.2.2)
 formatR        1.2.1   2015-09-18 CRAN (R 3.2.2)
 htmltools      0.2.6   2014-09-08 CRAN (R 3.2.0)
 knitr        * 1.11    2015-08-14 CRAN (R 3.2.2)
 lattice      * 0.20-33 2015-07-14 CRAN (R 3.2.1)
 latticeExtra * 0.6-26  2013-08-15 CRAN (R 3.1.0)
 magrittr       1.5     2014-11-22 CRAN (R 3.2.2)
 memoise        0.2.1   2014-04-22 CRAN (R 3.1.0)
 RColorBrewer * 1.1-2   2014-12-07 CRAN (R 3.2.2)
 rmarkdown    * 0.7     2015-06-13 CRAN (R 3.2.1)
 stringi        0.5-5   2015-06-29 CRAN (R 3.2.2)
 stringr        1.0.0   2015-04-30 CRAN (R 3.2.2)
 yaml           2.1.13  2014-06-12 CRAN (R 3.1.0)
[1] "2015-10-15 18:55:40 BRT"

CE 089 - Estatística Computacional II