Referência: S. Ross. Simulation, 2nd edition. New York: Academic Press, 1997.
Neste material estamos interessados em estudar métodos para resolver o seguinte problema:
Problema. Seja \(X\) uma variável aleatória (v.a.) discreta cuja função de probabilidade é dada por \[ f_{X}(x)=\Pr(X=x)=\left\{ \begin{array}{ll} p_{i}, & x=x_{i} \\ 0, & \mbox{caso contrário} \end{array} \right. \] para \(i=0,1,2,\ldots\). Desejamos então simular valores de \(X\). Em outras palavras, desejamos instruir um computador para gerar automaticamente valores de \(X\) de acordo com a “lei” \(f_{X}\).
A notação descrita acima é geral e será utilizada ao longo de toda exposição. Dois métodos serão estudados:
No R o bloco construtivo básico para simular valores de v.a.’s é a função runif(), que gera números aleatórios no intervalo \((0,1)\). Tais números são realizações de v.a.’s uniformemente distribuídas no intervalo \((0,1)\).
runif(1)
## [1] 0.865903
runif(5)
## [1] 0.09602833 0.59096798 0.30033937 0.86098779 0.54884286
runif(1,min=0,max=1)
## [1] 0.4817448
runif(5,min=0,max=1)
## [1] 0.58637598 0.02505666 0.92893268 0.22404233 0.61434796O método da transformação inversa é um método geral para simular valores de uma v.a.. Este método é baseado no seguinte resultado.
Resultado. Seja \(U\sim U(0,1)\). Para qualquer função de distribuição \(F\), a v.a. \(X\) definida por \[ X=F^{-1}(U) \] possui distribuição \(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 afirma que podemos simular um valor de uma v.a. \(X\) (\(X\) com distribuição \(F\)) aplicando o seguinte algoritmo:
Algoritmo: método da transformação inversa
Entrada: invGF() # inversa generalizada de F
1 Gere um número aleatório u
2 Calcule x <- invGF(u)
Saída: xÉ 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\). Em particular, o método também é válido no caso de v.a.’s discretas. Seja \(X\) uma v.a. discreta tal como descrita na Introdução. Quando aplicando o método da transformação inversa para simular valores de \(X\), obtemos o seguinte resultado: \[ X=F^{-1}(U)=\left\{ \begin{array}{ll} x_{0}, & 0\leq U<p_{0} \\ x_{1}, & p_{0}\leq U<p_{0}+p_{1} \\ x_{2}, & p_{0}+p_{1}\leq U<p_{0}+p_{1}+p_{2} \\ \vdots & \\ x_{i}, & \sum_{j=0}^{i-1}p_{j}\leq U<\sum_{j=0}^{i}p_{j} \\ \vdots & \end{array} \right. \] Por fim, note que \[ \Pr(X=x_{i})=\Pr\left(\sum_{j=0}^{i-1}p_{j}\leq U<\sum_{j=0}^{i}p_{j}\right)=f(x_{i})=p_{i}. \]
Exemplo. Simule um valor da v.a. discreta \(X\) cuja função de probabilidade é dada por
| X | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| P(X=x) | 0.20 | 0.15 | 0.25 | 0.40 |
Em R o método pode ser implementado tal como segue:
u <- runif(1)
if(u<0.20){
x <- 1
} else {
if(u<0.35){
x <- 2
} else {
if(u<0.60){
x <- 3
} else {
x <- 4
}
}
}
x
## [1] 3
u
## [1] 0.445029A implementação será certamente mais eficiente se o encadeamento condicional for realizado dos eventos mais prováveis para os menos prováveis. Neste caso, temos que \[ X=\left\{ \begin{array}{ll} 4, & 0\leq U<0.40 \\ 3, & 0.40\leq U<0.65 \\ 1, & 0.65\leq U<0.85 \\ 2, & \mbox{caso contrário.} \end{array} \right. \]
u <- runif(1)
if(u<0.40){
x <- 4
} else {
if(u<0.65){
x <- 3
} else {
if(u<0.85){
x <- 1
} else {
x <- 2
}
}
}
x
## [1] 1
u
## [1] 0.7239768Por fim, se desejar simular \(n\) valores de \(X\) (\(n>1\)), utilize um loop.
fsimul.vad <- function(u){
if(u<0.40){
out <- 4
} else {
if(u<0.65){
out <- 3
} else {
if(u<0.85){
out <- 1
} else {
out <- 2
}
}
}
return(out)
}
n <- 1000
amostra <- rep(NA,n)
for(k in 1:n){u <- runif(1); amostra[k] <- fsimul.vad(u)}
table(amostra) # freq. absoluta de cada possível valor de X na amostra
## amostra
## 1 2 3 4
## 213 161 264 362
freq <- (as.vector(table(amostra)))/n
freq # freq. relativa de cada possível valor de X na amostra
## [1] 0.213 0.161 0.264 0.362
c(0.20,0.15,0.25,0.40) # probabilidades teóricas para comparação
## [1] 0.20 0.15 0.25 0.40Um caso particular importante é a simulação da distribuição uniforme discreta com \(m\) valores possíveis. Ou seja, suponha que \(X\) é uma v.a. discreta assumindo qualquer um dos seguintes valores \(1,2,\ldots,m\) e que \[ p_{i}=\Pr(X=i)=\frac{1}{m} \] para todo \(i=1,2,\ldots,m\). Aplicando o método da transformação inversa temos que \[ X=i\Leftrightarrow\frac{i-1}{m}\leq U < \frac{i}{m}. \] Portanto, \[ X=i\Leftrightarrow i-1\leq mU < i. \] Em outras palavras, \[ X=\lfloor mU\rfloor + 1 \] onde \(\lfloor a\rfloor\) representa o piso do número real \(a\).
Exemplo. Simule um valor da distribuição uniforme discreta assumindo valores no conjunto \(\{1,2,\ldots,10\}\). Basta fazer \[ X=\lfloor 10U\rfloor + 1. \]
Em R o método pode ser implementado tal como segue:
m <- 10
x <- floor(m*runif(1))+1
x
## [1] 7Por fim, se desejar simular \(n\) valores desta uniforme discreta (\(n>1\)), basta fazer as seguintes operações:
n <- 1000
m <- 10
amostra <- floor(m*runif(n))+1
table(amostra) # freq. absoluta de cada possível valor de X na amostra
## amostra
## 1 2 3 4 5 6 7 8 9 10
## 106 105 91 99 95 116 105 86 103 94
freq <- (as.vector(table(amostra)))/n
freq # freq. relativa de cada possível valor de X na amostra
## [1] 0.106 0.105 0.091 0.099 0.095 0.116 0.105 0.086 0.103 0.094Uma v.a. discreta \(X\) possui distribuição geométrica com probabilidade de sucesso \(p\) se \[ p_{i}=\Pr(X=i)=pq^{i-1} \] onde \(0\leq p\leq 1\), \(q=1-p\) e \(i=1,2,\ldots\). Alguns autores definem a distribuição geométrica com probabilidade de sucesso \(p\) como \[ p_{j}=\Pr(Y=j)=pq^{j} \] onde \(0\leq p\leq 1\), \(q=1-p\) e \(j=0,1,2,\ldots\). A relação entre estas duas definições é bem simples. De fato, temos que \(X=Y+1\) ou \(Y=X-1\) onde \(X\) é o número de repetições de um ensaio de Bernoulli até o primeiro sucesso e \(Y\) é o número de repetições de um ensaio de Bernoulli antes do primeiro sucesso. A simulação de valores da distribuição geométrica com probabilidade de sucesso \(p\) segue da seguinte equação: \[ X=\lfloor log(U)/log(1-p)\rfloor + 1 \] onde \(\lfloor a\rfloor\) representa o piso do número real \(a\).
Exemplo. Simule valores da distribuição geométrica com probabilidade de sucesso \(p=0.7\).
n <- 10000
p <- 0.7
amostra <- floor(log(runif(n))/log(1-p))+1
tabela.freq <- table(amostra) # freq. absoluta
freq.rel <- (as.vector(tabela.freq))/n
y <- 0:9;x <- y+1;prob.real <- round(dgeom(y,0.7),4)
# probabilidades teóricas
data.frame(x,y,prob.real)
## x y prob.real
## 1 1 0 0.7000
## 2 2 1 0.2100
## 3 3 2 0.0630
## 4 4 3 0.0189
## 5 5 4 0.0057
## 6 6 5 0.0017
## 7 7 6 0.0005
## 8 8 7 0.0002
## 9 9 8 0.0000
## 10 10 9 0.0000
# comparação com as probabilidades teóricas
tabela.freq
## amostra
## 1 2 3 4 5 6 7 8
## 6958 2118 653 182 65 18 5 1
freq.rel
## [1] 0.6958 0.2118 0.0653 0.0182 0.0065 0.0018 0.0005 0.0001Uma v.a. discreta \(X\) possui distribuição de Poisson com parâmetro \(\lambda>0\) se \[ p_{i}=\Pr(X=i)=\frac{\lambda^{i}}{i!}e^{-\lambda} \] onde \(i=0,1,2,\ldots\) e \(\lambda\) é a taxa média de ocorrências por unidade de tempo. O símbolo \(e\) é definido por \(e=\lim_{n\rightarrow\infty}=[1+(1/n)]^{n}\approx 2.7183\). Um ponto importante para aplicar o método da transformação inversa para simular valores da distribuição de poisson é saber que as probabilidades \(p_{i}\) podem ser calculadas recursivamente sabendo que \(p_{0}=e^{-\lambda}\) e que \[ p_{i+1}=\frac{\lambda}{i+1}p_{i} \] para todo \(i=0,1,2,\ldots\). Assim, aplicando o método da transformação inversa temos que \[ X=i\Leftrightarrow\sum_{j=0}^{i-1}p_{j}\leq U<\sum_{j=0}^{i}p_{j}. \]
Exemplo. Simule um valor da distribuição de Poisson com \(\lambda=2\). Em R o método acima pode ser implementado tal como segue:
lamb <- 2
u <- runif(1)
flag <- 0
i <- 0; F <- p <- exp(-lamb)
repeat{
if(u<F){
x <- i; flag <- 1
} else {
p <- (lamb/(i+1))*p; F <- F+p
i <- i+1
}
if(flag==1){break}
}
x
## [1] 0
u
## [1] 0.08775602
round(dpois(0:10,2),4) # probabilidades teóricas para comparação
## [1] 0.1353 0.2707 0.2707 0.1804 0.0902 0.0361 0.0120 0.0034 0.0009 0.0002
## [11] 0.0000
round(ppois(0:10,2),4) # função de distribuição acumulada
## [1] 0.1353 0.4060 0.6767 0.8571 0.9473 0.9834 0.9955 0.9989 0.9998 1.0000
## [11] 1.0000Por fim, se desejar simular \(n\) valores da distribuição de Poisson (\(n>1\)), utilize um loop.
fsimul.pois <- function(u,lamb){
flag <- 0
i <- 0; F <- p <- exp(-lamb)
repeat{
if(u<F){
out <- i; flag <- 1
} else {
p <- (lamb/(i+1))*p; F <- F+p
i <- i+1
}
if(flag==1){break}
}
return(out)
}
n <- 10000; lamb <- 2
amostra <- rep(NA,n)
for(k in 1:n){u <- runif(1); amostra[k] <- fsimul.pois(u,lamb)}
table(amostra) # freq. absoluta de cada possível valor de X na amostra
## amostra
## 0 1 2 3 4 5 6 7 8 9 10 12
## 1298 2738 2784 1786 876 353 121 33 5 4 1 1
freq <- (as.vector(table(amostra)))/n
freq # freq. relativa de cada possível valor de X na amostra
## [1] 0.1298 0.2738 0.2784 0.1786 0.0876 0.0353 0.0121 0.0033 0.0005 0.0004
## [11] 0.0001 0.0001
round(dpois(0:10,2),4) # probabilidades teóricas para comparação
## [1] 0.1353 0.2707 0.2707 0.1804 0.0902 0.0361 0.0120 0.0034 0.0009 0.0002
## [11] 0.0000Esta seção fica como exercício. Apresente um script em R para simular valores da distribuição binomial com parâmetros \(n\) e \(p\). Neste caso \[ p_{i}=\Pr(X=i)=\binom{n}{i}p^{i}q^{n-i} \] para todo \(i=1,2,\ldots,n\), onde \(n\) é o número total de ensaios dicotômicos, \(p\) é a probabilidade de sucesso e \(q=1-p\).
Suponha que temos um método para simular valores de uma v.a. com distribuição \(g\). É possível usar este método como uma base para gerar valores de uma outra v.a. com distribuição \(f\). Basta gerar \(y\) de \(g\) e então aceitar este valor gerado como se fosse de \(f\), com uma probabilidade proporcional a \(f(y)/g(y)\). Este método geral para simular valores de uma v.a. é conhecido como o método da aceitação e rejeição.
Determine uma constante \(c\) tal que \[ f(x)\leq cg(x) \] para todo \(x\) onde \(f(x)>0\). O método da aceitação e rejeição para simular valores da distribuição \(f\) pode ser descrito pelo seguinte algoritmo:
Algoritmo
Entrada: f e g # distribuições
1 Gere um número aleatório y da distribuição g
2 Gere um número aleatório u
3 Se u < f(y)/cg(y), então faça x<-y. Caso contrário volte para 1.
Saída: xResultado. O método da aceitação e rejeição simula uma v.a. \(X\) tal que \(X\sim f\). Além disso, o número de iterações do algorítmo que são necessárias para gerar um valor simulado de \(X\) é uma variável aleatória com distribuição geométrica com média \(c\).
Novamente é bem simples verificar o resultado acima. Seja \(A\) o evento “o valor \(y\) gerado no passo \(1\) é aceito”. Segue da regra de Bayes que: \[ \Pr(X=y)=\Pr(Y=y|A)=\frac{\Pr(Y=y)\Pr(A|Y=y)}{\Pr(A)}. \] Como \[ \Pr(A|Y=y)=\Pr\left(U\leq\frac{f(y)}{cg(y)}\right)=\frac{f(y)}{cg(y)} \] segue do teorema da probabilidade total que \(\Pr(A)\) é dada por \[ \Pr(A)=\sum_{y}\Pr(Y=y)\Pr(A|Y=y)=\sum_{y}g(y)\cdot\frac{f(y)}{cg(y)}=\frac{1}{c}. \] Logo, \[ \Pr(X=y)=\frac{g(y)[f(y)/cg(y)]}{1/c}=f(y). \] Portanto, \(X\sim f\). Além disso, a cada iteração do algoritmo, o valor gerado no passo \(1\) pode ser aceito ou não (ou seja, a cada iteração o evento \(A\) ocorre ou não). As iterações são independentes uma da outra, visto que \(Y\) (no passo \(1\)) é independente \(U\) (no passo \(2\)). Se \(N\) denota o número de iterações do algorítmo que são necessárias para gerar um valor simulado da variável \(X\), segue que \(N\sim\mbox{Geo}(p)\), com \(p=\Pr(A)=1/c\), ou seja, \[ \Pr(N=n)=\left(1-\frac{1}{c}\right)^{n-1}\frac{1}{c},\hspace{1cm}n=1,2,3,\ldots \] e \(E(N)=1/\Pr(A)=c\). Concluímos que o número de iterações que são necessárias para gerar um valor simulado da distribuição de \(X\) é uma v.a. com distribuição geométrica com média \(c\). Assim, quanto menor for o valor de \(c\) para que \(c\cdot g\) envelope \(f\), mais eficiente será a simulação.
Exemplo. Deseja-se simular valores de uma v.a. discreta \(X\) cuja função de probabilidade é dada por:
| X | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| P(X=x) | 0.12 | 0.09 | 0.08 | 0.12 | 0.10 | 0.09 | 0.09 | 0.10 | 0.10 | 0.11 |
f.simul <- function(p.vec){
flag <- 0
repeat{
px <- p.vec; m <- length(p.vec); py <- rep(1/m,m)
c <- max(px/py)
u1 <- runif(1); y <- floor(m*u1)+1
u2 <- runif(1); aux <- (m*px[y])/c
if(u2<=aux){out <- y;flag <- 1}
if(flag==1){break}
}
return(out)
}
p.vec <- c(12,9,8,12,10,9,9,10,10,11)/100
f.simul(p.vec)
## [1] 3
#
# Se desejar n valores simulados, use um loop
#
p.vec <- c(12,9,8,12,10,9,9,10,10,11)/100
n <- 10000
amostra <- rep(NA,n)
for(k in 1:n){amostra[k] <- f.simul(p.vec)}
table(amostra) # freq. absoluta de cada possível valor de X na amostra
## amostra
## 1 2 3 4 5 6 7 8 9 10
## 1141 885 801 1191 1035 912 897 1048 986 1104
freq <- (as.vector(table(amostra)))/n
freq # freq. relativa de cada possível valor de X na amostra
## [1] 0.1141 0.0885 0.0801 0.1191 0.1035 0.0912 0.0897 0.1048 0.0986 0.1104
round(freq,2)
## [1] 0.11 0.09 0.08 0.12 0.10 0.09 0.09 0.10 0.10 0.11
p.vec # probabilidades teóricas para comparação
## [1] 0.12 0.09 0.08 0.12 0.10 0.09 0.09 0.10 0.10 0.11permuta <- function(k,p){
num.i <- floor(k*runif(1))+1
aux <- p[k];p[k] <- p[num.i];p[num.i] <- aux
return(p)
}
n <- 10
p <- 1:n
k <- n
repeat{
if(k==n){print(p)}
p <- permuta(k,p);print(p)
k <- k-1
if(k==1){break}
}
## [1] 1 2 3 4 5 6 7 8 9 10
## [1] 1 2 3 4 5 6 7 10 9 8
## [1] 1 2 9 4 5 6 7 10 3 8
## [1] 1 2 9 10 5 6 7 4 3 8
## [1] 1 2 7 10 5 6 9 4 3 8
## [1] 1 2 7 10 6 5 9 4 3 8
## [1] 1 2 7 6 10 5 9 4 3 8
## [1] 6 2 7 1 10 5 9 4 3 8
## [1] 6 7 2 1 10 5 9 4 3 8
## [1] 7 6 2 1 10 5 9 4 3 8