5.1 Monte Carlo

(Metropolis and Ulam 1949) propuseram o chamado método de Monte Carlo, em alusão ao complexo de cassinos em Monte Carlo na Itália. A proposta “consiste em uma mistura de processos determinísticos e estocásticos” ao considerar extrações “de forma aleatória e independente” admitindo “duas classes: (1) produção de valores ‘aleatórios’ com sua distribuição de frequência igual àquelas que regem a mudança de cada parâmetro, (2) cálculo dos valores daqueles parâmetros que são determinísticos, ou seja, obtidos algebricamente dos demais”.

It may seem strange that the machine can simulate the production of a series of random numbers, but this is indeed possible. In fact, it suffices to produce a sequence of numbers between 0 and 1 which have a uniform distribution in this interval but otherwise are uncorrelated, i.e., pairs will have uniform distribution in the unit square, triplets uniformly distributed in the unit cube, etc., as far as practically feasible. (Metropolis and Ulam 1949, 339)

Exemplo 5.1 A função stats::runif é o gerador de uniformes padrão (entre 0 e 1) do R, anotado por \(\mathcal{U}(0,1)\) conforme Seção de Distribuição Uniforme Contínua.

# R base
par(mfrow=c(3,3))
for(i in 1:9){
  set.seed(i); u <- runif(1000)
  hist(u); rug(u)
}

Exercício 5.1 Veja a documentação de base::set.seed.

Exercício 5.2 Veja o Capítulo 6 de (Peng 2022).

5.1.1 Geração de variáveis aleatórias discretas

(Gamerman and Lopes 2006) sugerem procedimentos bastante simples para gerar valores das principais distribuições discretas a partir de valores de uma \(\mathcal{U}[0,1]\).

5.1.1.1 Distribuição Bernoulli

É possível gerar \(M\) valores de uma Bernoulli de parâmetro \(\theta\) da seguinte forma:

Obtenha \(u\) de uma \(\mathcal{U}(0,1)\).
Se \(u \le \theta\) faça \(x=1\), caso contrário \(x=0\).
Repita 1 e 2 \(M\) vezes.

Exemplo 5.2 Deseja-se gerar \(M=1000\) observações de uma Bernoulli de parâmetro \(\theta=0.75\).

# parâmetros
M <- 1000
theta <- 0.75
x <- vector(length = M)
# loop
for(i in 1:M){
  # passo 1
  u <- runif(1)
  # passo 2
  set.seed(3*i); x[i] <- as.numeric(u <= theta)
}
# verificando
tab <- table(x)
(ptab <- prop.table(tab))

## x
##     0     1 
## 0.227 0.773

barplot(ptab)

Exercício 5.3 Compare os resultados do Exemplo 5.2 com aqueles obtidos por stats::rbinom ou np.random.binomial.

5.1.1.2 Distribuição Binomial

É possível gerar \(M\) valores \(x\) de uma Binomial de parâmetros \(n\) e \(\theta\) da seguinte forma:

Obtenha \(u_1,u_2,\ldots,u_n\), todos de uma \(\mathcal{U}(0,1)\).
Se \(u_i \le \theta\) faça \(x_i=1\), caso contrário \(x_i=0\). Obtenha \(\sum_{i=1}^{n} x_i\).
Repita 1 e 2 \(M\) vezes.

Exemplo 5.3 Deseja-se gerar \(M=1000\) observações de uma Binomial de parâmetros \(n=10\) e \(\theta=0.75\).

# parâmetros
M <- 1000
n <- 10
theta <- 0.75
x <- vector(length = M)
# loop
for(i in 1:M){
  # passo 1
  u <- runif(n)
  # passo 2
  set.seed(i+314); x[i] <- sum(u < theta)
}
# verificando
tab <- table(x)
(ptab <- prop.table(tab))

## x
##     3     4     5     6     7     8     9    10 
## 0.001 0.018 0.066 0.136 0.259 0.283 0.177 0.060

barplot(ptab)

Exercício 5.4 Compare os resultados do Exemplo 5.3 com aqueles obtidos por stats::rbinom ou np.random.binomial.

5.1.2 Método da Transformação Inversa

Pode-se obter uma amostra de tamanho \(M\) de uma função com densidade \(f\) e função distribuição acumulada \(F(x) = \int_{-\infty}^{x} f(t) dt\) da seguinte maneira:

Obtenha \(u\) de uma \(\mathcal{U}(0,1)\).
Calcule \(F^{-1}(u)\).
Repita 1 e 2 \(M\) vezes.

Exercício 5.5 Veja o artigo Inverse transform sampling.

Exemplo 5.4 Considere a Distribuição Exponencial tal que \[ f(x|\lambda) = \lambda e^{-\lambda x}\] \[F(x|\lambda) = 1 - e^{-\lambda x}\] Desta forma a inversa de \(F\) é

\[\begin{equation} F^{-1}(u|\lambda) = -\frac{\log(1-u)}{\lambda} \tag{5.1} \end{equation}\]

# Passo 1: uniforme em [0,1]
set.seed(1); u <- runif(1000)
hist(u); rug(u)

# Passo 2: exponencial de taxa 3
lambda <- 3
x <- -log(1-u)/lambda
hist(x, 30, probability = TRUE); rug(x)
curve(dexp(x,3), add = TRUE, col = 'red')

mean(x); 1/lambda

## [1] 0.3373056

## [1] 0.3333333

var(x); 1/lambda^2

## [1] 0.1208113

## [1] 0.1111111

Exercício 5.6 Mostre que a inversa da FDA da exponencial pode ser dada pela Equação (5.1).

Exercício 5.7 Considere a distribuição triangular apresentada na Seção de Distribuição Triangular.

Simule a densidade de uma \(\mathcal{Tri}(0,.6,1)\) a partir do Método da Transformação Inversa.
Compare com os resultados de extraDistr::rtriang do R ou de numpy.random.triangular do Python.

Exercício 5.8 Sabe-se que se a variável \(U\) tem distribuição uniforme no intervalo (0,1), então \(X=\sigma \sqrt{-2 \log U}\) tem distribuição Rayleigh de parâmetro \(\sigma\).

Simule a densidade de uma \(\mathcal{Rayleigh}(3)\) a partir do Método da Transformação Inversa.
Compare com os resultados de extraDistr::rrayleigh do R ou de numpy.random.rayleigh do Pyhon.

5.1.3 Amostragem com rejeição

Pode-se obter uma amostra de uma função com densidade \(f\) a partir de uma função \(g\) de fácil obtenção da seguinte maneira:

Obtenha \(u\) de uma \(\mathcal{U}(0,1)\).
Obtenha um candidato \(x\) com distribuição \(g\).
Se \(u \le \frac{f(x)}{k\,g(x)}\) admite-se o candidato \(x\), caso contrário rejeita-se \(x\).
Repita 1 a 3 \(M\) vezes.

(Albert 2009, 117) aponta que “[a] amostragem por rejeição é um método geral para simular a partir de uma distribuição posteriori arbitrária, mas pode ser difícil de ser estabelecida, pois requer a construção de uma densidade de proposta [\(g(x)\)] adequada⁵.”

Exemplo 5.5 Adaptado da Seção 6.3 de (Peng 2022). Deseja-se amostrar de uma exponencial com \(\lambda=1\). Para isso considera-se \(U \sim Unif(0,10)\) e constante \(k=10\).

A seguir podem-se considerar os gráficos e as simulações.

# gráficos analíticos
curve(dexp(x, 1), 0, 10, col = 4, ylab = 'Densidade')
segments(0, 1, 10, 1, col = 1)

# elementos da simulação
fdist <- vector()
lambda <- 1
g <- 1/10
k <- 10
M <- 1000
j <- 0

# loop
ini <- Sys.time()
for(i in 1:M){
  
  # Passo 1: Simular U ~ Unif(0,1)
  set.seed(i); (U <- runif(1, 0, 1))
  
  # Passo 2: Simular um candidato X ~ g, i.e., X ~ Unif(0,10)
  set.seed(i+1); (X <- runif(1, 0, 10))
  
  # Apresentando X no gráfico
  rug(X)
  
  # Avaliando f no ponto candidato X
  f <- lambda*exp(-lambda*X)
  
  # Passo 3: decidindo se o candidato entra ou não na amostra
  # Se U > f/(k*g) indicar com um x vermelho, X não entra como ponto de f
  if(U > f/(k*g)){
    points(X, U, pch = 'x', col = 'red')
  }
  # Se U <= f/(k*g)), indicar com um x verde, X entra como ponto de f
  if(U <= f/(k*g)){
    points(X, U, pch = 'x', col = 'green')
    j <- j+1
    fdist[j] <- X
  }
  # print(i/M)
}
Sys.time()-ini

## Time difference of 0.09595108 secs

# legenda
legend('right', c('10*g(x) Unif(0,10)', 'f(x) Exp(1)'), lty = 1, 
       col = c(1, 4), box.col = 'white')

Note que o tamanho de fdist é 103, o que representa apenas 10.3% dos 1000 valores simulados, indicando neste caso um algoritmo pouco eficiente. Pode-se utilizar a simulação para calcular probabilidades, tal como \(Pr(X<1)\). Por se tratar de uma conhecida distribuição exponencial, pode-se comparar o valor estimado com o valor obtido analiticamente.

sum(fdist < 1)/length(fdist) # estimado

## [1] 0.6601942

1-exp(-1*1) # analítico via função distribuição acumulada

## [1] 0.6321206

pexp(1) # analítico via função stats::pexp

## [1] 0.6321206

Referências

Albert, Jim. 2009. Bayesian Computation with R. Springer.

Gamerman, Dani, and Hedibert F Lopes. 2006. Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. CRC press.

Metropolis, Nicholas, and Stanislaw Ulam. 1949. “The Monte Carlo Method.” Journal of the American Statistical Association 44 (247): 335–41. https://doi.org/10.2307/2280232.

Peng, Roger D. 2022. Advanced Statistical Computing. Work in Progress. https://bookdown.org/rdpeng/advstatcomp/.

“Rejection sampling is a general method for simulating from an arbitrary posterior distribution, but it can be difficult to set up since it requires the construction of a suitable proposal density.” (Albert 2009, 117)↩︎