Simulação de Números Aleatórios
Mauro Campos. Departamento de Estatística, UFES
01 de agosto de 2018
Referência: S. Ross. Simulation, 2nd edition. New York: Academic Press, 1997.
1 Motivação
Resultado. Seja \(U\sim U(0,1)\) e \(F\) uma função de distribuição acumulada. Se \(X\) é definida por \(X=F^{-1}(U)\), então \(X\sim F\), onde \(F^{-1}(u)=\inf\{x:F(x)\geq u\}\) é a inversa generalizada de \(F\). Além disso, segue que \(F(X)\sim U(0,1)\).
Este resultado, bem conhecido da Teoria de Probabilidades, informa que se uma amostra aleatória da distribuição \(U(0,1)\) está disponível, então esta amostra pode ser usada para se obter uma amostra aleatória de uma distribuição \(F\). É importante lembrar que uma amostra aleatória simples da distribuição \(U(0,1)\) é simplesmente uma sequência de variáveis aleatórias independentes e uniformemente distribuídas no intervalo \((0,1)\).
É simples verificar o resultado acima. Em geral, observe que \[ \Pr(X\leq x)=\Pr(F^{-1}(U)\leq x)=\Pr(U\leq F(x))=F(x). \] Logo \(X\sim F\).
2 Números aleatórios \(\times\) números pseudo-aleatórios
Números aleatórios no intervalo \((0,1)\) formam uma amostra aleatória da distribuição \(U(0,1)\).
Números pseudo-aleatórios no intervalo \((0,1)\) formam uma sequência de valores, normalmente gerados determinísticamente por um computador, que para todos os fins práticos imitam números aleatórios no intervalo \((0,1)\) (ou seja, imitam uma amostra aleatória da distribuição \(U(0,1)\)). Tecnicamente, isto significa que se um método é usado para computacionalmente gerar números pseudo-aleatórios no intervalo \((0,1)\), os valores gerados devem necessarimente passar por testes de aderência e de independência para serem utilizados como uma amostra aleatória da distribuição \(U(0,1)\) em simulações computacionais de processos ou sistemas.
3 O método congruencial linear
O método congruencial linear é usado para gerar números pseudo-aleatórios \(u_{1},u_{2},u_{3},\ldots\) no intervalo \((0,1)\). O método consiste de duas etapas:
- Obtenha uma sequência \(x_{1},x_{2},x_{3},\ldots\) de valores usando a recorrência \(x_{k+1}=(a\cdot x_{k} + b) \mod m\) para \(k=1,2,3,\ldots\), onde \(a,b\) e \(m\) são inteiros positivos e \(x_{1}\) é escolhido tal que \(0\leq x_{1}<m\).
- Obtenha agora a sequência \(u_{1},u_{2},u_{3},\ldots\), onde \(u_{k}=x_{k}/m\) para \(k=1,2,3,\ldots\).
A recorrência na primeira etapa do método significa que \(x_{k+1}\) é igual ao resto após a divisão de \((a\cdot x_{k}+b)\) por \(m\). Os números gerados por esta recorrência estão entre \(0\) e \(m-1\), inclusive. Portanto, a sequência \(u_{1},u_{2},u_{3},\ldots\) são valores no intervalo \((0,1)\).
4 Exemplos
Exemplo. Quando \(a=6\), \(b=7\), \(m=5\) e \(x_{1}=2\).
kmax <- 15
a <- 6;b <- 7;m <- 5
x <- rep(NA,kmax);x[1] <- 2
for(k in 1:(kmax-1)){x[k+1] <- (a*x[k]+b)%%m}
u <- x/m;u
## [1] 0.4 0.8 0.2 0.6 0.0 0.4 0.8 0.2 0.6 0.0 0.4 0.8 0.2 0.6 0.0
Exemplo. Quando \(a=30\), \(b=60\), \(m=100\) e \(x_{1}=64\).
kmax <- 25
a <- 30;b <- 60;m <- 100
x <- rep(NA,kmax);x[1] <- 64
for(k in 1:(kmax-1)){x[k+1] <- (a*x[k]+b)%%m}
u <- x/m;u
## [1] 0.64 0.80 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60
## [15] 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60 0.60
Exemplo. Quando \(a=1\), \(b=1\), \(m=10\) e \(x_{1}=3\).
kmax <- 30
a <- 1;b <- 1;m <- 10
x <- rep(NA,kmax);x[1] <- 3
for(k in 1:(kmax-1)){x[k+1] <- (a*x[k]+b)%%m}
u <- x/m;u
## [1] 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## [18] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.0 0.1 0.2Estes exemplos mostram que o método congruencial linear é dependente dos valores escolhidos para \(a,b,m\) e \(x_{1}\).
Uma boa receita. Quando \(a=7^{5}\), \(b=0\), \(m=2^{31}-1\) e \(x_{1}=231\).
kmax <- 700
a <- 7^(5);b <- 0;m <- 2^(31)-1
x <- rep(NA,kmax);x[1] <- 231
for(k in 1:(kmax-1)){x[k+1] <- (a*x[k]+b)%%m}
u <- x/m
classes <- seq(0,1,0.1)
hist(u,breaks=classes,freq=F,right=F,col="gray",ylim=c(0,1.4),xlab="u",ylab="Densidade",main="")
points(seq(0,1,0.1),dunif(seq(0,1,0.1),0,1),type="l",col="blue")
Existem muitas outras receitas para \(a,b,m\) e \(x_{1}\) que proporcionam bons geradores de números pseudo-aleatórios.
Daqui pra frente, pra simplificar a linguagem, números pseudo-aleatórios no intervalo \((0,1)\) serão chamados simplesmente de números aleatórios no intervalo \((0,1)\).
5 Números aleatórios no R
Vimos nas seções anteriores que números aleatórios no intervalo \((0,1)\) são os blocos construtivos básicos para gerar valores de uma distribuição \(F\) não-uniforme. No R números aleatórios no intervalo \((0,1)\) são gerados usando a função runif().
runif(1)
## [1] 0.2143544
runif(5)
## [1] 0.33518252 0.41547661 0.29175395 0.66492202 0.09129638
runif(1,min=0,max=1)
## [1] 0.9211488
runif(5,min=0,max=1)
## [1] 0.5261263 0.4449835 0.5647774 0.6036906 0.3647829
set.seed(123)
u <- runif(700)
classes <- seq(0,1,0.1)
hist(u,breaks=classes,freq=F,right=F,col="gray",ylim=c(0,1.4),xlab="u",ylab="Densidade",main="")
points(seq(0,1,0.1),dunif(seq(0,1,0.1),0,1),type="l",col="blue")
6 Introdução à integração Monte Carlo
Uma aplicação imediata para números aleatórios no intervalo \((0,1)\) é a estimação de integrais. Vamos começar com a situação mais simples que é estimar a seguinte integral: \[ I=\int_{0}^{1}g(x)dx. \] É simples ver que \[ I=E(g(U)) \] onde \(U\sim U(0,1)\). Portanto, pela lei dos grandes números segue que \[ \bar{g}_{N}=\frac{1}{N}\sum_{k=1}^{N}g(U_{k})\stackrel{p}{\rightarrow}E(g(U))=I \hspace{1.5cm}N\rightarrow\infty \] onde \(U_{1},\ldots,U_{N}\stackrel{iid}{\sim}U(0,1)\). Este resultado garante que se \(N\) for suficientemente grande, então com alta probabilidade o valor de \(I\) é bem aproximado por \(\bar{g}_{k}\), ou seja, \[ I\approx\bar{g}_{N}=\frac{1}{N}\sum_{k=1}^{N}g(u_{k}) \] onde \(u_{1},\ldots,u_{N}\) são \(N\) números aleatórios no intervalo \((0,1)\). Este processo de estimação da integral \(I\) é chamado integração Monte Carlo.
O seguinte algoritmo pode ser aplicado:
Algoritmo: estimação de I
Entrada: tol, Nmax
Saída: I.hat
1 soma <- 0; N <- 0; I.hat <- 0
2 N <- N+1
3 I.old <- I.hat
4 gere um número aleatório u
5 soma <- soma + g(u)
6 I.hat <- soma/N
7 dif <- I.hat-I.old
8 se (dif < tol) ou (N>Nmax) pare. Caso contrário volte para o passo 2
9 retorna I.hat
Vejamos um exemplo ilustrativo simples: estimar \(I=\int_{0}^{1}3x^{2}dx=1\).
tol <- 0.00001; Nmax <- 1000
soma <- 0; N <- 0; I.hat <- 0
repeat{
N <- N+1
I.old <- I.hat
soma <- soma + 3*(runif(1)^(2))
I.hat <- soma/N
dif <- abs(I.hat-I.old)
if((dif<tol)|(N>Nmax)){break}
}
N
## [1] 420
I.hat
## [1] 0.9669233
c(I.hat,I.old,dif)
## [1] 9.669233e-01 9.669164e-01 6.892665e-06De forma alternativa, poderíamos fazer também:
N <- 1000
u <- runif(N)
g <- 3*(u^(2))
I.hat <- mean(g)
I.hat
## [1] 0.9974833Cada uma das seguintes integrais \[ \int_{a}^{b}h(x)dx\hspace{1cm}\int_{0}^{+\infty}h(x)dx\hspace{1cm}\int_{-\infty}^{+\infty}h(x)dx \] pode ser estimada por integração Monte Carlo após uma transformação de variável apropriada.
Por fim, a ideia geral da integração Monte Carlo pode ser extendida para o caso de integrais multidimensionais. De fato, basta ver que \[ I=\int_{0}^{1}\cdots\int_{0}^{1}g(x_{1},\ldots,x_{n})dx_{1}\cdots dx_{n}=E[g(U_{1},\ldots,U_{n})] \] e que \[ I\approx\frac{1}{N}\sum_{k=1}^{N}g(u_{1}^{k},\ldots,u_{n}^{k}) \] onde \(u_{1}^{k},\ldots,u_{n}^{k}\) são \(n\) números aleatórios independentes para cada \(k\) variando de \(1\) até \(N\).
7 Estimando \(\pi\)
Basta ver que \[ \frac{\pi}{4}=I=\Pr(A)=\int_{0}^{1}\int_{0}^{1}I_{A}(x,y)dxdy\hspace{2cm} A=\{(x,y)\in[0,1]^{2}:(x-0.5)^{2}+(y-0.5)^{2}<0.25\}. \] Portanto \[ \hat{\pi}\approx 4\cdot\hat{I}. \] Assim
N <- 20000
ux <- runif(N)
uy <- runif(N)
d2 <- (ux-0.5)^(2)+(uy-0.5)^(2)
aux1 <- which(d2<0.25)
aux2 <- length(aux1); aux2
## [1] 15712
I.hat <- aux2/N
I.hat
## [1] 0.7856
pi.hat <- 4*I.hat
pi.hat
## [1] 3.1424