Muitas distribuições podem ser geradas a partir de uma variável uniforme em \((0,1)\). O método da transformação inversa é um dos procedimentos mais importantes da simulação estatística baseado neste princípio.
Suponha \(X\) uma v.a. discreta com f.p.
\[ P(X = x_j) = p_X(x_j) = p_j, \qquad j = 0,1,\ldots, \qquad \mbox{com} \qquad \sum_{j} p_j = 1. \]
Neste caso, se \(U \sim \text{U}(0,1)\), note que
\[ P(a < U < b) = b - a, \qquad \forall\, 0< a < b <1. \]
Logo,
\[ \begin{aligned} P(X = x_j) &= p_j = \sum_{i=0}^{j} p_i - \sum_{i=0}^{j-1} p_i = P\left(\sum_{i=0}^{j-1} p_i < U < \sum_{i=0}^{j} p_i\right). \end{aligned} \]
O que significa que, se sabemos gerar \(U\), podemos criar um procedimento para gerar \(X\) a partir dos valores gerados para \(U\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < p_0\), tome \(X = x_0\); caso contrário,
se \(u < p_0 + p_1\), tome \(X = x_1\); caso contrário,
se \(u < \sum_{i=0}^{2} p_i\), tome \(X = x_2\); caso contrário,
…
se \(u < \sum_{i=0}^{j} p_i\), tome \(X = x_j\); caso contrário,
…
Repita os passos 1 e 2 quantas vezes forem necessárias até obter a amostra do tamanho desejado.
Simular valores para a v.a. \(X\) com f.p. dada por
| \(x\) | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| \(p_X(x)\) | 0,20 | 0,15 | 0,25 | 0,40 |
Solução:
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < 0,20\), tome \(X = 1\); caso contrário,
se \(u < 0,35\), tome \(X = 2\); caso contrário,
se \(u < 0,60\), tome \(X = 3\); caso contrário,
tome \(X = 4\).
n = 10000
x = NULL
for (i in 1:n){
u = runif(1)
if(u<0.20){x[i] = 1}else{
if(u<0.35){x[i] = 2}else{
if(u<0.60){x[i] = 3}else{
x[i] = 4
}
}
}
}
table(x)/n## x
## 1 2 3 4
## 0.1986 0.1533 0.2469 0.4012Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < 0,40\), tome \(X = 4\); caso contrário,
se \(u < 0,65\), tome \(X = 3\); caso contrário,
se \(u < 0,85\), tome \(X = 1\); caso contrário,
tome \(X = 2\).
O Algoritmo 1.2 é mais eficiente, uma vez que a chance de verificar menos condições é maior.
Se os \(x_j\), \(j \geq 0\), são ordenados, isto é,
\[ x_0 < x_1 < x_2 < \cdots, \]
e \(F_X\) representa a f.d.a. de \(X\), então
\[ F_X(x_j) = \sum_{i=0}^{j} p_i, \]
e, portanto,
\(X\) será igual a \(x_0\) se \(U<F_X(x_0)\) e
\(X\) será igual a \(x_j\) se \(F_X(x_{j-1}) < U \leq F_X(x_j)\), \(j\geq 1\)
\(\qquad\qquad\qquad\qquad \Leftrightarrow F_X^{-1}(U) \in [x_{j-1},\,x_j)\)
Daí o nome do método.
Suponha \(X\) uma v.a. com distribuição uniforme discreta no conjunto \(\{1,2,3,...,n\}\). Implemente o algoritmo da transformação inversa para gerar \(n=100\) valores de \(X\).
Simular valores de uma v.a. \(X \sim \text{Ber}(p)\), para \(p\) fixo.
Solução:
Sejam \(X\) e \(U\) duas v.a. discretas tais que \(X \sim \text{Ber}(p)\) e \(U \sim \text{U}(0,1)\).
Neste caso, \(\{X=1\}\) e \(\{U<p\}\) são eventos equiprováveis, uma vez que \[P(X=1) = P(U<p) = p.\]
Dessa forma, se queremos gerar uma amostra de tamanho \(n\) para \(X\), basta gerar valores \(u_i \sim \text{U}(0,1)\) e tomar \(x_i=1\) sempre que \(u_i<p\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < p\), tome \(X = 1\); caso contrário,
tome \(X = 0\).
n = 10000
p = 0.3
u = runif(n)
x = ifelse(u < p,1,0)
table(x)/n## x
## 0 1
## 0.7058 0.2942De forma equivalente, podemos utilizar a f.d.a. de \(X\): \[ F_X(x) = \left\{ \begin{array}{ll} 0, & x<0\\ 1-p, & 0\leq x<1\\ 1, & x \geq 1; \end{array} \right. \]
tomando \(\quad\) \(x_i=0\) sempre que \(0<u_i\leq 1-p\) \(\quad\) e \(\quad\) \(x_i=1\) se \(1-p<u_i\leq 1\) (ou, equivalentemente, \(0<1-u_i\leq p\)).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < 1-p\), tome \(X = 0\); caso contrário,
tome \(X = 1\).
n = 10000
p = 0.3
u = runif(n)
x = ifelse(u < 1-p,0,1)
table(x)/n## x
## 0 1
## 0.7026 0.2974# Escrevendo como uma função
rbern.inv = function(n,p){
u = runif(n)
x = ifelse(u < 1-p,0,1)
x
}Seja \(X\) a v.a. que conta o número de realizações independentes de um ensaio de Bernoulli com probabilidade de sucesso \(p\) necessárias até que ocorra o primeiro sucesso.
Simular valores da v.a. \(X\), para \(p\) fixo.
Solução:
Neste caso, vamos denotar \(X \sim \text{Geo}(p)\), cuja função de probabilidade é dada por \[p_X(x) = p(1-p)^{x-1} = p q^{x-1}, \quad x = 1, 2, 3, \dots,\] onde denotamos \(q=1-p\).
Aplicando a definição de função de distribuição acumulada,
\[F_X(x) = \sum_{j=1}^{x} P_X(j) = \sum_{j=1}^{x} p q^{j-1} = \frac{p}{q} \cdot \sum_{j=1}^{x} q^j\]
Para simplificar o resultado anterior, vamos utilizar a expressão da soma de uma quantidade finita de termos da PG de razão \(q\).
Seja \(S_n = \sum_{j=1}^{n} q^j\).
\[\begin{eqnarray*} S_n &=& q + q^2 + q^3 + \dots + q^n \\ q S_n &=& \qquad q^2 + q^3 + \dots + q^{n+1} \end{eqnarray*}\]
Diminuindo uma equação da outra,
\[S_n - q S_n = q - q^{n+1} \implies (1-q)S_n = q - q^{n+1} \implies S_n = \frac{q - q^{n+1}}{1-q} = \frac{q(1-q^n)}{1-q}\]
Donde segue que
\[F_X(x) = \frac{p}{q} \cdot \frac{q(1-q^x)}{1-q} = 1-q^x\]
Assim, utilizando o método da transformação inversa, gerando \(U \sim \text{U}(0,1)\), tomamos
\[X = j \qquad \text{se} \qquad F_X(j-1)=1-q^{j-1} < u \le F_X(j)=1-q^j\] ou, equivalentemente,
\[X = j \qquad \text{se} \qquad q^j \le 1-u < q^{j-1}.\]
Na prática, a expressão acima indica que
\[X = \min\{j: 1-u > q^j\} = \min\{j: \log(1-u) > \log(q^j)=j\log(q)\} \underbrace{=}_{0<q<1} \min\left\{j: j > \dfrac{\log(1-u)}{\log{q}}\right\}\]
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Tome \(x = \left\lfloor \dfrac{\log(1-u)}{\log{q}} \right\rfloor\).
rgeom.inv = function(n,p){
u = runif(n)
x = floor(log(1-u)/log(1-p))
x
}
x = rgeom.inv(1000,0.3)
table(x)[1:10]## x
## 0 1 2 3 4 5 6 7 8 9
## 268 227 162 116 70 38 42 20 21 10Simular valores de uma v.a. \(X \sim \text{Poi}(\lambda)\), para \(\lambda\) fixo.
Solução:
Seja \(X\) uma v.a. discreta tal que \(X \sim \text{Poi}(\lambda)\), cuja função de probabilidade é dada por
\[p_X(x) =
\dfrac{e^{-\lambda}\lambda^x}{x!}, \quad x = 1, 2, 3,
\dots,\]
Neste caso não há forma simplificada para a f.d.a.
\[F_X(x) = \sum_{j=0}^x p_X(j) =
\sum_{j=0}^x\dfrac{e^{-\lambda}\lambda^j}{j!},\]
no entanto, note que \[p_X(j+1)
= \dfrac{e^{-\lambda}\lambda^{j+1}}{(j+1)!}
= \dfrac{e^{-\lambda}\lambda\lambda^j}{(j+1)j!}
=\dfrac{\lambda}{j+1} \dfrac{e^{-\lambda}\lambda^j}{j!}
=\dfrac{\lambda}{j+1}p_X(j);\] o que induz uma forma recursiva
para implementação do Algoritmo original.
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\) e tome \(j=0\), \(p = p_0 = e^{-\lambda}\) e \(F = p\).
Se \(u < F\), tome \(X = j\); caso contrário,
tome \(j = j+1\), \(\qquad p = \dfrac{\lambda}{j}p \qquad\) e \(\qquad F = F+p\)
se \(u < F\), tome \(X = j\); caso contrário, retorne ao passo 2.
rpois.inv = function(n,lambda){
x = NULL
for (i in 1:n){
u = runif(1)
j = 0
p = exp(-lambda)
F = p
while(u>=F){
j = j+1
p = (lambda/j)*p
F = F + p
}
x[i] = j
}
x
}
x = rpois.inv(100000,5)
table(x)[1:10]## x
## 0 1 2 3 4 5 6 7 8 9
## 732 3380 8556 13928 17559 17555 14578 10364 6530 3629Crie uma função que simule valores de uma v.a. discreta \(X \sim \text{Bin}(n,p)\).
Dica: Utilize uma estratégia semelhante à do caso Poisson, explorando a relação entre \(p_X(j+1)\) e \(p_X(j)\).
Compare os geradores implementados com as funções nativas do R. Para isto utilize gráficos e medidas resumo sobre as amostras obtidas através das diferentes funções.
ROSS, Sheldon M. Simulation. 6. ed. London: Academic Press, 2022.