2 Geração de variáveis aleatórias

Neste capítulo, apresentamos técnicas práticas que podem produzir variáveis aleatórias de distribuições padrão e não padrão usando um programa de computador.

Dada a disponibilidade de um gerador uniforme em R, conforme explicado na Seção 2.1.1, não lidamos com a produção específica de variáveis aleatórias uniformes. As técnicas mais básicas relacionam a distribuição a ser simulada de uma variável uniforme por uma transformada ou uma propriedade probabilística particular, como na Seção 2.2, enquanto a mais genérica é uma versão de simulação do método de tentativa e erro, descrito na Seção 2.3 sob o nome do método Accept-Reject.

Em todos os casos, os métodos contam com a disponibilidade de sequências de gerações uniformes independentes fornecidas pelo gerador R residente, runif.


2.1 Introdução

Os métodos desenvolvidos neste texto e resumidos sob a denominação de métodos de Monte Carlo baseiam-se principalmente na possibilidade de produzir, com um computador, um fluxo supostamente infinito de variáveis aleatórias para distribuições conhecidas ou novas. Tal simulação é, por sua vez, baseada na produção de variáveis aleatórias uniformes no intervalo \((0, 1)\).

Embora não estejamos diretamente preocupados com a mecânica de produzir essas variáveis aleatórias uniformes, porque os geradores uniformes existentes podem ser considerados “perfeitos”, confiaremos totalmente nesses geradores para produzir outras variáveis aleatórias. Em certo sentido, a distribuição uniforme \(U[0,1]\) fornece a representação probabilística básica da aleatoriedade em um computador e os geradores para todas as outras distribuições requerem uma sequência de variáveis uniformes a serem simuladas.

Como já apontado na Seção 1.4 do Capítulo 1, o R possui um grande número de funções internas que irão gerar as variáveis aleatórias padrão listadas na Tabela 1.1. Por exemplo,

set.seed(8645)
rgamma(n=3, shape=2.5, rate=4.5)
## [1] 0.4823078 0.7784044 0.8130656
produz três gerações independentes de uma distribuição \(Gama(5/2,9/2)\) com todas as devidas garantias de representar esta distribuição.

Portanto, é contraproducente, ineficiente e até mesmo perigoso gerar a partir dessas distribuições padrão usando qualquer coisa que não seja os geradores R residentes. Os princípios desenvolvidos nas seções a seguir são, no entanto, essenciais para lidar com distribuições menos padrão que não são incorporadas ao R.


2.1.1 Simulação Uniforme

O gerador uniforme básico em R é a função runif, cuja única entrada obrigatória é o número de valores a serem gerados. Os outros parâmetros opcionais são min e max, que caracterizam os limites do intervalo que suporta o uniforme. O padrão é min=0 e max=1.

Por exemplo,

runif(100, min=2, max=5)
##   [1] 4.489150 2.033353 3.152120 3.752938 2.523232 2.260375 2.806874 4.321926
##   [9] 3.971142 2.532192 3.282657 3.753058 2.235750 4.334596 4.710294 4.725567
##  [17] 2.045297 2.156042 4.172345 4.569682 4.874380 2.938462 3.423675 2.204636
##  [25] 4.048736 3.418593 4.725319 3.750835 4.713688 3.659483 3.398340 4.941957
##  [33] 2.870286 4.514118 4.807172 2.393737 2.853960 4.628588 3.712513 3.036174
##  [41] 3.913611 2.071040 3.446080 2.695356 4.176157 4.057998 3.812711 2.329913
##  [49] 2.098854 2.425561 3.759616 3.653614 2.812323 2.832957 4.295760 2.871476
##  [57] 4.030831 3.201702 4.127869 2.381369 3.900356 3.470281 3.773113 4.982890
##  [65] 2.664752 2.550644 3.080258 4.654706 3.155671 4.337826 4.163387 3.767221
##  [73] 2.097771 2.464136 2.678777 2.833337 3.221942 3.651476 2.660196 2.176000
##  [81] 3.848358 2.241907 4.565161 4.667912 4.045428 3.937666 3.227411 2.653083
##  [89] 4.819891 4.493923 3.267842 2.977529 2.365820 2.979139 3.674945 4.348825
##  [97] 4.615357 3.977588 3.961450 2.146076
produzirá 100 variáveis aleatórias distribuídas uniformemente entre 2 e 5.

A rigor, todos os métodos que veremos e isso inclui runif, produzem números pseudo-aleatórios em que não há aleatoriedade envolvida; baseado em um valor inicial \(u_0\) de uma sequência uniforme \(U(0, 1)\) e uma transformação \(D\), o gerador uniforme produz uma sequência \(\{u_i\} = \{D^i(u_0)\}\) de valores em \((0, 1)\); mas o resultado tem as mesmas propriedades estatísticas de uma sequência i.i.d..

Embora testes extensivos dessa função tenham sido realizados para garantir que ela produza variáveis uniformes para todos os propósitos (consulte, por exemplo, Robert and Casella, 2004), uma verificação rápida das propriedades desse gerador uniforme é observar um histograma dos \(X_i\)’s, um gráfico dos pares \((X_i,X_{i+1})\) e a função de autocorrelação estimada, pois qualquer gerador de variável aleatória sofre de uma autocorrelação residual e bons algoritmos reduzirão isso a um valor insignificante.

O código R usado para produzir a saída na Figura 2.1 é apresentado aseguir e mostra que runif é aparentemente aceitável para esta avaliação casual. 

Nsim=10^4 #number of random numbers
x=runif(Nsim)
x1=x[-Nsim] #vectors to plot
x2=x[-1] #adjacent pairs
par(mfrow=c(1,3),mar=c(4,4,1,1))
hist(x); box()
plot(x1,x2, cex=0.6, xlab=expression(x[1]),ylab=expression(x[2]),
     main="Pairwise plot")
y=acf(x, plot=F)
plot(y[1:40]);grid()

Figura. 2.1. Histograma (à esquerda), gráfico de pares (centro) e função de autocorrelação estimada (à direita) de uma sequência de \(10^4\) números aleatórios uniformes gerados por runif.

Conforme apontado no comentário anterior, runif não envolve aleatoriedade per se. A produção de runif(Nsim) é melhor descrita como uma sequência determinística baseada em um ponto inicial aleatório. Uma ilustração extrema desse fato é obtida por meio da função R set.seed, que usa seu único argumento inteiro para definir quantas sementes forem necessárias.

Por exemplo,

set.seed(1)
runif(5)
## [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
set.seed(1)
runif(5)
## [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
set.seed(2)
runif(5)
## [1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393
mostra que definir a semente determina todos os valores subsequentes produzidos pelo gerador aleatório.

Na esmagadora maioria dos casos, não definimos a semente, que é então escolhida de acordo com o tempo atual. Mas em configurações onde precisamos reproduzir exatamente a mesma sequência de simulações aleatórias, por exemplo, para comparar dois procedimentos ou duas velocidades, é razoável definir um valor fixo da semente.


2.1.2 A transformada inversa

Existe uma transformação simples, às vezes útil, conhecida como transformação integral de probabilidade, que nos permite transformar qualquer variável aleatória em uma variável aleatória uniforme e, mais importante, vice-versa.

Por exemplo, se \(X\) tem densidade \(f\) e função de distribuição \(F\), então temos a relação \[ F(x) = \int_{-\infty}^x f(t) \mbox{d}t, \] e se definirmos \(U = F(X)\), então \(U\) é uma variável aleatória distribuída de uma uniforme \(U(0,1)\).

Isto é porque \[ P(U\leq u) = P\big(F(X)\leq F(x)\big) = P\big(F^{-1}(F(X))\leq F^{-1}(F(x))\big) = P(X\leq x), \] onde assumimos que \(F\) tem inversa. Essa suposição pode ser relaxada (ver Robert and Casella, 2004), mas vale para a maioria das distribuições contínuas.


Exemplo 2.1

Sabemos que se \(X\sim Exp(1)\), então \(F(x)=1-e^{-1}\). Resolvendo para \(x\) em \(u=1-e^{-x}\), temos \(x=-\log(1-x)\).

Portanto, se \(U\sim U(0,1)\), \[ X=-\log(U)\sim Exp(1), \] ambos, \(U\) e \(1-U\) são uniformes.

O código R a seguir compara a saída da transformada inversa de probabilidade com a saída do rexp. Os ajustes de ambos os histogramas ao seu limite exponencial não são distinguíveis na Figura 2.2. 

set.seed(546)
Nsim=10^4 # number of random variables
U=runif(Nsim)
X=-log(U) # transforms of uniforms
Y=rexp(Nsim) # exponentials from R
par(mfrow=c(1,2), mar=c(4,4,1,1)) # plots
hist(X,freq=F,main="Exp from Uniform",breaks=10,ylim=c(0,1))
grid();box()
curve(dexp(x), col="brown",add=T,lwd=3)
hist(Y,freq=F,main="Exp from R",ylim=c(0,1),breaks=10);grid();box()
curve(dexp(x), col="brown",add=T,lwd=3)

Figura. 2.2. Histogramas de variáveis aleatórias exponenciais usando a transformada inversa (à esquerda) e usando o comando R rexp (à direita), com a densidade \(Exp(1)\).


A geração de variáveis aleatórias uniformes é, portanto, um determinante chave do comportamento dos métodos de simulação para outras distribuições de probabilidade, uma vez que essas distribuições podem ser representadas como uma transformação determinística de variáveis aleatórias uniformes.


2.2 Métodos gerais de transformação

Quando uma distribuição com densidade \(f\) é vinculada de maneira relativamente simples a outra distribuição fácil de simular, essa relação pode frequentemente ser explorada para construir um algoritmo para simular variáveis de \(f\).


Exemplo 2.2

No Exemplo 2.1, vimos como gerar uma variável aleatória exponencial a partir de uma uniforme. Agora ilustramos algumas das variáveis aleatórias que podem ser geradas a partir de uma distribuição exponencial.

Se os \(X_i\) são variáveis aleatórias iid \(Exp(1)\), então três distribuições padrão podem ser derivadas como \[\begin{equation} \tag{2.1} \begin{array}{rcl} Y & = & \displaystyle 2\sum_{j=1}^\nu X_j \, \sim \, \chi^2_{2\nu}, \qquad \nu\in\mathbb{N}^* \\ Y & = & \displaystyle \beta\sum_{j=1}^a X_j \, \sim \, Gama(a,\beta), \qquad a\in\mathbb{N}^* \\ Y & = & \displaystyle \dfrac{\displaystyle \sum_{j=1}^a X_j}{\displaystyle \sum_{j=1}^{a+b} X_j} \, \sim \, Beta(a,b), \qquad a,b\in\mathbb{N}^*, \end{array} \end{equation}\] onde \(\mathbb{N}^*=\{1,2,\cdots\}\).

Por exemplo, para gerar variáveis aleatórias \(\chi^2_6\), poderíamos usar o código R: 

test1=function(){
  U=runif(3*10^4)
  U=matrix(data=U,nrow=3) # matrix for sums
  X=-log(U) # uniform to exponential
  X=2*apply(X,2,sum) # sum up to get chi squares
}
test2=function(){X=rchisq(10^4,df=6)}

Obviamente, isso não é tão eficiente quanto chamar rchisq, como pode ser verificado pelo código R: 

system.time(test1());system.time(test2())
##   usuário   sistema decorrido 
##      0.01      0.00      0.02
##   usuário   sistema decorrido 
##         0         0         0


Muitas outras derivações de distribuições padrão são possíveis ao aproveitar as propriedades probabilísticas existentes, conforme mostrado no Exercício 2.12.

Essas transformações são bastante simples de usar e, portanto, costumam ser as favoritas em nossas ilustrações. No entanto, há limites para sua utilidade, tanto no escopo das variáveis que podem ser geradas dessa maneira pense, por exemplo, em uma distribuição qui-quadrada com um número não par de graus de liberdade, quanto na eficiência da geração. Para qualquer distribuição específica, algoritmos eficientes foram desenvolvidos. Assim, se R tiver uma distribuição embutida, quase sempre vale a pena usá-la, como mostra o Exemplo 2.2. Além disso, o método de transformação descrito acima não pode atingir todas as distribuições; por exemplo, não podemos obter uma normal padrão.


2.2.1 Um gerador normal

Uma maneira de obter simulações da variável aleatória normal usando uma transformada é com o algoritmo Box-Muller, desenvolvido para a geração de variáveis \(N(0, 1)\).


Exemplo 2.3

Se \(U_1\) e \(U_2\) forem iid \(U(0,1)\), as variáveis \(X_1\) e \(X_2\), definidas por \[ X_1= \sqrt{-2\log(U_1)} \cos(2\pi U_2) \quad \mbox{e} \quad X_2=\sqrt{-2\log(U_1)} \sin(2\pi U_2), \] são então iid \(N(0,1)\) em virtude de uma simples mudança de argumento variável.

Observe que este não é o gerador implementado em R qnorm, que usa por padrão a transformada inversa de probabilidade, com base em uma representação muito precisa da inversa da distribuição acumulad normal, até 16 dígitos!. É, no entanto, possível, se não recomendado, mudar o gerador normal para a versão Box–Muller ou mesmo para o Kinderman-Ramage.


Em comparação com algoritmos grosseiramente aproximativos baseados no Teorema do Limite Central (CLT), o algoritmo Box-Muller é exato, produzindo duas variáveis aleatórias normais a partir de duas variáveis aleatórias uniformes, sendo a única desvantagem, em velocidade a necessidade de calcular funções transcendentais como log, cos e sin.

A simulação de uma distribuição normal multivariada \(N_p(\mu,\Sigma)\), onde \(\Sigma\) é uma matriz \(p\times p\) simétrica e definida positiva, pode ser derivada do gerador rnorm genérico usando uma decomposição de Cholesky de \(\Sigma\), ou seja, \(\Sigma = AA^T\) e tomando a transformada por \(A\) de um vetor normal iid de dimensão \(p\) leva a um vetor normal \(N_p(0,\Sigma)\).

Existe, no entanto, uma função R que replica essas etapas, chamada rmnormt e disponível na biblioteca mnormt (Genz and Azsalini, 2009). Esta biblioteca também permite o cálculo da probabilidade de hipercubos através da função sadmvn, como em

library(mnormt)
Sigma=matrix(c(1,2,0,2,5,0.5,0,0.5,3), 3, 3)
sadmvn(low=c(1,2,3),upp=c(10,11,12),mean=rep(0,3),varcov=Sigma)
## [1] 0.006014523
## attr(,"error")
## [1] 3.58177e-07
## attr(,"status")
## [1] "normal completion"
onde \(B\) é uma matriz definida positiva. Isso é bastante útil, pois a derivação analítica dessa probabilidade é quase sempre impossível.


2.2.2 Distribuições discretas

Em seguida, nos voltamos para a geração de variáveis aleatórias discretas, onde temos um algoritmo “para todos os fins”. Novamente, usando o princípio da transformada inversa da Seção 2.1.2, podemos de fato construir um algoritmo genérico que funcionará formalmente para qualquer distribuição discreta.

Para gerar \(X\sim P_\theta\), onde \(P_\theta\) é suportada pelos números inteiros, podemos calcular, de uma vez, assumindo que podemos armazená-los, as probabilidades \[ p_0\sim P_\theta(X\leq 0), \quad p_1=P_\theta(X\leq 1),\quad p_2=P(X\leq 2), \quad \cdots \] e então, gerar \(U\sim U(0,1)\) e selecionar \[ X=k \quad \mbox{se} \quad p_{k-1}<U<p_k\cdot \]


Exemplo 2.4

Para gerar \[ X\sim Binomial(10, 0.3), \] os valores de probabilidade são obtidos por pbinom(k,10,0.3) como:

pbinom(q=0:10,size=10,prob=0.3)
##  [1] 0.02824752 0.14930835 0.38278279 0.64961072 0.84973167 0.95265101
##  [7] 0.98940792 0.99840961 0.99985631 0.99999410 1.00000000

e para gerar \(X\sim Poisson(7)\),

ppois(q=0:20,lambda=7)
##  [1] 0.000911882 0.007295056 0.029636164 0.081765416 0.172991608 0.300708276
##  [7] 0.449711056 0.598713836 0.729091268 0.830495937 0.901479206 0.946650377
## [13] 0.973000227 0.987188607 0.994282798 0.997593420 0.999041817 0.999638216
## [19] 0.999870149 0.999955598 0.999985505
a sequência é interrompida quando atinge 1 com um determinado número de decimais, por exemplo, \(p_{20} = 0.999985\).


Algoritmos específicos são geralmente mais eficientes, como mostrado no Exemplo 2.5, mas principalmente por causa do problema de armazenamento. Muitas vezes podemos melhorar o algoritmo acima por uma escolha criteriosa de quais probabilidades calculamos primeiro. Por exemplo, se quisermos gerar variáveis aleatórias a partir de uma distribuição de Poisson com média \(\lambda = 100\), o algoritmo acima é lamentavelmente ineficiente.

Isso ocorre porque esperamos que a maioria de nossas observações esteja no intervalo \(\lambda\pm \sqrt{\lambda}\), lembre-se de que \(\lambda\) é tanto a média quanto a variância da distribuição Poisson, e para \(\lambda = 100\) esse intervalo é \((70,130)\). Assim, começar em 0 quase sempre produzirá 70 testes de \(p_{k-1} < U < p_k\) que são inúteis porque quase certamente serão rejeitados. Um primeiro remédio é “ignorar” o que está fora de um intervalo altamente provável, como \((70,130)\) no exemplo atual, como \[ P(X<70)+P(X>130)=0.00268\cdot \]

Formalmente, devemos encontrar um limite inferior e superior para tornar essa probabilidade pequena o suficiente, mas informalmente \(\pm 3\sigma\) funciona bem.


Exemplo 2.5

Mostramos um código R que pode ser usado para gerar variáveis aleatórias Poisson para grandes valores de lambda. A sequência \(t\) contém os valores inteiros no intervalo em torno da média. 

set.seed(832)
Nsim=10^4; lambda=100
spread=3*sqrt(lambda)
t=round(seq(max(0,lambda-spread),lambda+spread,1))
prob=ppois(t, lambda)
X=rep(0,Nsim)
for (i in 1:Nsim){
  u=runif(1)
  X[i]=t[1]+sum(prob<u) }
par(mfrow=c(1,1),mar=c(4,4,1,1))
hist(X);grid();box()

mean(X);var(X)
## [1] 99.8771
## [1] 98.0866

A última linha do programa verifica em qual intervalo a variável aleatória uniforme caiu e atribui o valor de Poisson correto a \(X\). Veja o Exercício 2.14 para outras distribuições.


Uma solução mais formal para a ineficiência de iniciar as probabilidades cumulativas em \(p_0\) é começar a partir da moda da distribuição discreta \(P_\theta\) e explorar os valores vizinhos até que a probabilidade cumulativa seja 1 até um erro de aproximação. Os \(p_k\)’s são então indexados pelos valores visitados em vez dos inteiros, mas a validade do método permanece completa.

Existem algoritmos específicos para quase todas as distribuições e geralmente são bastante rápidos. Assim, reforçamos mais uma vez que, se R tiver a distribuição que lhe interessa, o mais sensato é utilizá-la. Mais uma vez, a comparação do código do Exemplo 2.5 com o rpois residente mostra como essa simples implementação pode ser ineficiente. Se test3 corresponder ao acima e test4 a rpois, os tempos de execução são dados por 

test3=function(){
  Nsim=10^4; lambda=100
  spread=3*sqrt(lambda)
  t=round(seq(max(0,lambda-spread),lambda+spread,1))
  prob=ppois(t, lambda)
  X=rep(0,Nsim)
  for (i in 1:Nsim){
    u=runif(1)
    X[i]=t[1]+sum(prob<u) }
}
test4=function(){rpois(10^4,100)}
system.time(test3()); system.time(test4())
##   usuário   sistema decorrido 
##      0.02      0.02      0.03
##   usuário   sistema decorrido 
##         0         0         0

No entanto, o R não lida com todas as distribuições de que precisaremos; portanto, abordagens como a acima podem ser úteis. Veja o Exercício 2.15 para alguns algoritmos específicos.


2.2.3 Representações de mistura

Às vezes, uma distribuição de probabilidade pode ser representada naturalmente como uma distribuição de mistura; ou seja, podemos escrevê-lo na forma \[\begin{equation} \tag{2.2} f(x)=\int_\mathcal{Y} g(x|y)p(y) \mbox{d}y \qquad \mbox{ou} \qquad f(x)=\sum_{y\in\mathcal{Y}} p_i f_i(x), \end{equation}\] dependendo se o espaço auxiliar \(\mathcal{Y}\) é contínuo ou discreto, onde \(g\) e \(p\) são distribuições padrão que podem ser facilmente simuladas.

Para gerar uma variável aleatória \(X\) usando tal representação, podemos primeiro gerar uma variável \(Y\) a partir da distribuição mista e depois gerar \(X\) a partir da distribuição condicional selecionada.

Significa que, \[ \mbox{se} \; y\sim p(y) \;\mbox{e} \; X\sim f(x|y), \; \mbox{então} \; X\sim f(x), \; \mbox{caso seja contínua}, \\ \mbox{se} \; \gamma\sim P(\gamma=i)=p_i \;\mbox{e} \; X\sim f_\gamma(x), \; \mbox{então} \; X\sim f(x), \; \mbox{caso seja discreta}\cdot \]

Por exemplo, podemos escrever a densidade \(t\)-Student com \(\nu\) graus de liberdade \(T_\nu\), como uma mistura, onde \[ X|y N(0,\nu/y) \quad \mbox{e} \quad Y\sim \chi^2_\nu\cdot \]

Gerar a partir de uma distribuição \(T_\nu\) poderia equivaler a gerar a partir de uma distribuição \(\chi^2_\nu\) e, em seguida, a partir da distribuição normal correspondente. Obviamente, usar rt é um pouco mais eficiente, como você pode verificar via system.time.


Exemplo 2.6

Se \(X\) é uma variável aleatória binomial negativa, \(X\sim Neg(n,p)\), então \(X\) tem a representação da mistura \[ X \, | \,y\sim P(y) \] e \[ Y\sim Gama(n,\beta), \] onde \(\beta= (1-p)/p\).

O código R a seguir é gerado a partir dessa mistura e produz a Figura 2.3, onde o ajuste para a função de probabilidade binomial negativa também é mostrado. 

Nsim=10^4
n=6;p=.3
y=rgamma(Nsim,n,rate=p/(1-p))
x=rpois(Nsim,y)
hist(x,main="",freq=F,col="grey",breaks=40)
lines(1:50,dnbinom(1:50,n,p),lwd=2,col="sienna")
grid(); box()
Figura. 2.3. Histograma de \(10^4\) variáveis aleatórias binomiais negativas \(Neg(6,0.3)\) geradas a partir da representação da mistura juntamente com a função de probabilidade.



2.3 O método Accept-Reject

Existem muitas distribuições para as quais o método de transformação inversa e até mesmo transformações gerais falharão em gerar as variáveis aleatórias necessárias. Para esses casos, devemos recorrer a métodos indiretos; isto é, métodos nos quais geramos uma variável aleatória candidata e só a aceitamos sujeita a passar em um teste. Como veremos, esta classe de métodos é extremamente poderosa e nos permitirá simular virtualmente qualquer distribuição.

Esses chamados métodos de Accept-Reject ou Aceitar-Rejeitar exigem apenas que conheçamos a forma funcional da densidade \(f\) de interesse, chamada densidade alvo, até uma constante multiplicativa. Usamos uma densidade \(g\) mais simples para simular, chamada densidade instrumental ou candidata, para gerar a variável aleatória para a qual a simulação é realmente feita.

As únicas restrições que impomos a essa densidade candidata \(g\) são que

  1. \(f\) e \(g\) têm suportes compatíveis, isto é, \(g(x) > 0\) quando \(f(x) > 0\).

  2. Existe uma constante \(M\) com \(f(x)/g(x)\leq M\) para todo \(x\).

Neste caso, \(X\) pode ser simulado da seguinte forma.

Primeiro, geramos \(Y\sim g\) e, independentemente, geramos \(U\sim U(0,1)\). Se \[ U\leq \dfrac{1}{M}\dfrac{f(Y)}{g(Y)}, \] então definimos \(X = Y\). Se a desigualdade não for satisfeita, descartamos \(Y\) e \(U\) e começamos de novo. Sucintamente, a representação algorítmica do método Accept-Reject é a seguinte:


Algoritmo 1 Método Accept-Reject
1. Gerar \(Y\sim g\), \(U\sim U(0,1)\);

  1. Aceitar \(X = Y\) se \(U\leq f(Y)/M g(Y)\);

  2. Retornar para 1 caso contrário.


A implementação no R deste algoritmo é direta: se randg é uma função que fornece gerações a partir da densidade \(g\), no mesmo espírito de rnorm ou rt, uma versão R simples do Algoritmo 1 é


u=runif(1)*M
y=randg(1)
while (u>f(y)/g(y)){
 u=runif(1)*M
 y=randg(1)}
que produz uma única geração \(y\) de \(f\).

Por que esse método funciona? Um cálculo de probabilidade direto mostra que a função de distribuição da variável aleatória aceita \(P\left(Y\leq x \, | \, U\leq f(Y)/\big(Mg(Y)\big)\right)\), é exatamente a função de distribuição de \(X\). Ou seja, \[ \begin{array}{rcl} P\left(Y\leq x \, | \, U\leq f(Y)/\big(Mg(Y)\big)\right) & = & \dfrac{P\big(Y\leq x, U\leq f(Y)/Mg(Y)\big) \big)}{P\big( U\leq f(Y)/Mg(Y)\big)} \\ & = & \dfrac{\displaystyle \int_{-\infty}^x \int_0^{f(y)/Mg(y)} \mbox{d}u \, g(y)\mbox{d}y}{\displaystyle \int_{-\infty}^\infty \int_0^{f(y)/Mg(y)} \mbox{d}u \, g(y)\mbox{d}y} \\ & = & \dfrac{\displaystyle \int_{-\infty}^x \Big(f(y)/Mg(y)\big) g(y)\mbox{d}y}{\displaystyle \int_{-\infty}^\infty \Big(f(y)/Mg(y)\Big) g(y)\mbox{d}y} \\ & = & \dfrac{\displaystyle \int_{-\infty}^x f(y)\mbox{d}y}{\displaystyle \int_{-\infty}^\infty f(y)\mbox{d}y} \, = \, P(X\leq x), \end{array} \] onde usamos o fato de que a integral uniforme é igual ao seu limite superior.

Apesar de simular apenas a partir de \(g\), a saída desse algoritmo é exatamente distribuída a partir de \(f\). O método Accepts-Reject é aplicável em qualquer dimensão, desde que \(g\) seja uma densidade no mesmo espaço que \(f\).

Observe o cancelamento dos \(g(y)\)’s e dos \(M\)’s nas integrais acima. Também decorre dessa representação que não precisamos nos preocupar em normalizar constantes. Contanto que saibamos \(f/g\) até uma constante, \(f/g\approx \widetilde{f}/\widetilde{g}\), o algoritmo pode ser implementado se um limite superior \(\widetilde{M}\) puder ser encontrado em \(\widetilde{f}/\widetilde{g}\). As constantes que faltam realmente são absorvidas em \(M\).


Exemplo 2.7

O Exemplo 2.2 não forneceu um algoritmo geral para simular variáveis aleatórias \(Beta(\alpha,\beta)\). Podemos, no entanto, construir um algoritmo baseado no método Accept-Reject, usando como distribuição instrumental a distribuição \(U(0,1)\), quando ambos \(\alpha\) e \(\beta\) são maiores que 1.

A função rbeta genérica não impõe essa restrição. O limite superior \(M\) é então o máximo da densidade beta, obtido, por exemplo, por optimize ou seu alias optimise: 

optimize(f=function(x){dbeta(x,2.7,6.3)}, interval=c(0,1),maximum=T)$objective
## [1] 2.669744

Como a densidade candidata \(g\) é igual a um, o valor proposto \(Y\) é aceito se \[ M\times U < f(Y), \] ou seja, se \(M\times U\) estiver abaixo da densidade beta \(f\) naquela realização.

Observe que gerar \[ U\sim U(0,1) \] e multiplicar por \(M\) é equivalente a gerar \(U\sim U(0,M)\). Para \(\alpha= 2.7\) e \(\beta = 6.3\), uma implementação R alternativa do algoritmo Accept-Reject é:

Nsim=2500
a=2.7;b=6.3
M=2.67
u=runif(Nsim,max=M)   # uniform over (0,M)
y=runif(Nsim)         # generation from g
x=y[u<dbeta(y,a,b)]   #accepted subsample
e o painel esquerdo na Figura 2.4 mostra os resultados da geração de 2500 pares \((Y,U)\) de \(U(0,1)\times U(0,M)\). 

Nsim = 10^3
a = 2.7
b = 6.3
M = 2.67
y = runif(Nsim)
u = runif(Nsim, max = M)
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(y, u, col = "grey", pch = 19, cex = 0.4, ylab = expression(u.g(y)))
points(y[u < dbeta(y, a, b)], u[u < dbeta(y, a, b)], pch = 19, cex = 0.4)
curve(dbeta(x, a, b), col = "sienna", lwd = 2, add = T)
abline(h = M, col = "gold4", lwd = 2)
grid()
M = 1.68
y = rbeta(Nsim, 2, 6)
u = runif(Nsim, max = M)
labels = u < dbeta(y, a, b)/dbeta(y, 2, 6)
plot(y, u * dbeta(y, 2, 6), col = "grey", pch = 19, cex = 0.4, ylab = expression(u.g(y)))
points(y[labels], u[labels] * dbeta(y[labels], 2, 6), pch = 19, cex = 0.4)
curve(dbeta(x, a, b), col = "sienna", lwd = 2, add = T)
curve(M * dbeta(x, 2, 6), col = "gold4", lwd = 2, add = T)
grid()

Figura. 2.4. Geração de variáveis aleatórias \(X\sim Beta(2.7, 6.3)\): Usando o algoritmo Accept-Reject, 2500 \((Y,U)\) propostas foram geradas a partir de \(g\) e \(U(0,M)\), respectivamente, e os pontos \((Y, Ug(Y))\) foram representados com pontos cinza. No painel esquerdo, \(Y\sim U(0,1)\) e 36% das variáveis aleatórias candidatas foram aceitas e representadas com pontos pretos. No painel direito, \(Y\sim Beta(2, 6)\) e 58% dos valores simulados foram aceitos e representados de forma semelhante com pontos pretos. Em ambos os painéis, \(f\) e \(Mg\) também são plotados.

Os pontos pretos \((Y, Ug(Y))\) que caem sob a densidade \(f\) são aqueles para os quais aceitamos \(X = Y\) e rejeitamos os pontos cinzas \((Y, Ug(Y))\) que ficam fora. É novamente claro a partir desta representação gráfica que os pontos pretos são uniformemente distribuídos sobre a área sob a densidade \(f\). Como a probabilidade de aceitação de uma dada simulação é \(1/M\) (Exercício 2.5), com \(M = 2.67\) aceitamos aproximadamente \(1/2.67\) = 37% dos valores.


Na implementação do algoritmo Aceitar-Rejeitar acima, o número total de tentativas Nsim é fixo, o que significa que o número de valores aceitos é uma variável aleatória binomial com probabilidade \(1/M\). Em vez disso, na maioria dos casos, o número de valores aceitos é fixo, mas essa implementação pode ser explorada como em


x=NULL
while (length(x)<Nsim){
  y=runif(Nsim*M)
  x=c(x,y[runif(Nsim*M)*M<dbeta(y,a,b)])}
x=x[1:Nsim]

Observe que usar y=u=runif(Nsim* M) no programa produziria um viés, já que \(y\) e \(u\) teriam os mesmos valores.

Algumas propriedades-chave do algoritmo Accept-Reject, que sempre devem ser consideradas ao usá-lo, são as seguintes:

  1. Apenas a razão \(f/M\) é necessária, então o algoritmo não depende da constante de normalização.

  2. O limite \(f\leq M g\) não precisa ser rígido; o algoritmo permanece válido, se menos eficiente, quando \(M\) é substituído por qualquer constante maior.

  3. A probabilidade de aceitação é \(1/M\), então \(M\) deve ser o menor possível para um dado esforço computacional.

A eficiência de um determinado algoritmo Accept-Reject pode ser medida em termos de sua probabilidade de aceitação, pois quanto maior essa probabilidade, menos simulações desperdiçadas de \(g\). Em termos absolutos, isso deve ser ponderado pelo custo computacional de produzir um valor de \(g\), pois, caso contrário, a melhor escolha para \(g\) seria \(f\).


Exemplo 2.8 (Continuação do Exemplo 2.7)

Em vez disso, considere simular \(Y\sim Beta(2,6)\) como uma distribuição proposta. Esta escolha de \(g\) é aceitável, pois

optimize(f=function(x){dbeta(x,2.7,6.3)/dbeta(x,2,6)},maximum=T,interval=c(0,1))$objective
## [1] 1.671808

Esta modificação da proposta leva assim a um valor menor de \(M\) e uma taxa de aceitação correspondentemente maior de 58% do que com a proposta uniforme.

O painel direito da Figura 2.4 mostra o resultado do algoritmo Accept-Reject correspondente e ilustra o ganho de eficiência trazido pela simulação de pontos em um conjunto menor.


Como mostrado no Exemplo 2.8, alguma otimização do algoritmo Accept-Reject é possível escolhendo a densidade candidata \(g\) em uma família paramétrica e determinando o valor do parâmetro que minimiza o limite \(M\).

Às vezes pode acontecer que a complexidade da otimização seja muito cara em termos de análise ou tempo de computação. No primeiro caso, a construção do algoritmo ótimo ainda deve ser realizada quando o algoritmo for submetido a uso intensivo.

Este é, por exemplo, o caso da maioria dos geradores aleatórios em R, como pode ser verificado em help. No segundo caso, na maioria das vezes é preferível explorar o uso de outra família de distribuições instrumentais \(g\), ver Exercício 2.22. Uma aplicação particular do algoritmo Accept-Reject encontrou um nicho na genética de populações e é chamada de ABC, seguindo a denominação proposta por Pritchard et al. (1999).

A versão principal deste algoritmo é um algoritmo Acceptr-Reject adequado para problemas Bayesianos, onde uma distribuição a posteriori \(\pi(\theta|x_0) \approx \pi(\theta)f(x_0|\theta)\) deve ser simulada para uma função de verossimilhança \(f(x|\theta)\) que não está disponível, mas pode ser simulada. O algoritmo ABC então gera valores a partir da a priori e da verossimilhança até que a observação simulada seja igual à observação original \(x_0\):


Repetir
Gerar \(\theta\sim \pi(\theta)\) e \(X\sim f(x|\theta)\)
até \(X = x_0\).


Este algoritmo é, portanto, válido, mas só se aplica em configurações onde \(\pi(\theta)\) é uma priori adequado e onde \(P_\theta (X = x_0)\) tem uma probabilidade positiva de ocorrer. Mesmo em genética de população onde \(X\) é uma variável aleatória discreta, o tamanho do espaço de estado é frequentemente tal que este algoritmo não pode ser implementado.

A proposta de Pritchard et al. (1999) é então substituir a condição exata de aceitação \(X = x_0\) por uma condição aproximada \(d(X, x_0)<\epsilon\), onde \(d\) é uma distância e \(\epsilon\) um nível de tolerância. Embora inevitável, essa etapa de aproximação torna o método ABC difícil de recomendar de maneira geral, embora trabalhos mais recentes o reformulem em uma estrutura não paramétrica que visa aproximar a função de verossimilhança \(f(x|\theta)\) (Beaumont et al., 2002).

Uma crítica ao algoritmo Accept-Reject é que ele gera simulações “inúteis” da proposta \(g\) ao rejeitar, mesmo aquelas necessárias para validar a saída como sendo gerada a partir do alvo \(f\). Veremos no Capítulo 3 como o método de amostragem por importância (Seção 3.3) pode ser usado para ignorar este problema.


2.4 Exercícios