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.
Exercício 5.2 Veja a documentação de base::set.seed
.
Exercício 5.3 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
set.seed(3*i); u <- runif(1)
# passo 2
x[i] <- as.numeric(u <= theta)
}
# verificando
tab <- table(x)
(ptab <- prop.table(tab))
## x
## 0 1
## 0.228 0.772
Exercício 5.4 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
set.seed(i+314); u <- runif(n)
# passo 2
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.258 0.284 0.177 0.060
## [1] 9.536743e-07 2.861023e-05 3.862381e-04 3.089905e-03 1.622200e-02 5.839920e-02 1.459980e-01 2.502823e-01
## [9] 2.815676e-01 1.877117e-01 5.631351e-02
Exercício 5.5 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.6 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 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')
## [1] 0.3333333
## [1] 0.3373056
## [1] 0.1111111
## [1] 0.1208113
## [1] 0.7768698
## [1] 0.7768698
## [1] 0.777
Exercício 5.7 Mostre que a inversa da FDA da exponencial pode ser dada pela Equação (5.1).
Exercício 5.8 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 denumpy.random.triangular
do Python.
Exercício 5.9 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\). Realizando a transformação de variáveis conforme Seção 2.1, obtém-se \(f(x|\sigma)=\frac{x}{\sigma^2}e^{-x^2/(2\sigma^2)}\) e \(F(x|\sigma)=1-e^{-x^2/(2\sigma^2)}\).
- Obtenha \(F^{-1}(u)\).
- Simule a densidade de uma \(\mathcal{Rayleigh}(1)\) a partir do Método da Transformação Inversa.
- Compare com os resultados de
extraDistr::drayleigh
do R ounumpy.random.rayleigh
do Pyhon. - Calcule \(P(X<0.5)\) via simulação e compare com a probabilidade teórica.
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)\)] adequada5.”
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.09337902 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.
## [1] 0.6601942
## [1] 0.6321206
## [1] 0.6321206
Referências
“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)↩︎