Introdução:

Método da Transformação Inversa:

Suponha que deseja-se construir um gerador da vari?vel aleat?ria discreta \(X\), com fun??o de probabilidade \(f_X(x_j)=\operatorname{P}\{X = x_j\}=p_j \text{, } j = 0, 1, \dots\), com \(\sum_{j}p_j = 1\). A expressão de sua função de distribuição acumulada é \(F_X(x)=\operatorname{P}\{X \leq x\}\), ou seja:

\[\begin{equation} F_X(x) = \begin{cases} 0, & \text{se } x < x_0 \\ p_0, & \text{se } x_0 \leq x < x_1 \\ p_0 + p_1, & \text{se } x_1 \leq x < x_2 \\ \vdots & \\ \sum_{i=1}^{j-1}p_i & \text{se } x_{j-1} \leq x < x_j \\ \vdots \end{cases} \end{equation}\]

Embora não se possa estabelecer a inversa da função de sitribuição acumulada no caso das variáveis aleatórias discretas, é possível aplicar o Método da Transformação Inversa. Para tanto, gera-se um número \(U\) ao acaso no intervalo \((0, 1)\). Como \(\operatorname{P}\{a \leq U \leq b \} = b -a\), para \(0 < a < b < 1\), tem-se que \[ \operatorname{P}\{X = x_j \} = F_X(x_j) - F_X(x_{j-1}) = \operatorname{P}\left\{\sum_{i=0}^{j-1}p_i\leq U < \sum_{i=0}^j p_i\right\}=p_j \] ou seja, toma-se:

\[\begin{equation} X = \begin{cases} x_0, & \text{se } U < p_0 \\ x_1, & \text{se } p_0 \leq U < p_0 + p_1 \\ \vdots & \\ x_j & \text{se } \sum_{i=0}^{j-1}p_i \leq U < \sum_{i=0}^{j}p_i \\ \vdots \end{cases} \end{equation}\]

Método de Transformação Inversa (Discretas) - Procedimento A

  1. Gerar um número ao acaso \(U \sim \mathcal{Unif.}(0,1)\).

  2. Se \(U < p_0\), fazer \(X=x_0\) e parar.

  3. Se \(U < p_0 + p_1\), fazer \(X=x_1\) e parar.

  4. Se \(U < p_0 + p_1 + p_2\), fazer \(X=x_2\) e parar

  5. \(\vdots\)

Em outras palavras, depois da geração do número aleatório \(U\), o valor de \(X\) é determinado buscando-se o intervalo \((F_X(x_{j-1}), F_X(x_j))\) que contém \(U\). Ou seja, o algoritmo acima pode ser resumido como:

  1. Gerar um número ao acaso \(U \sim \mathcal{Unif.}(0,1)\).

  2. Fazer \(X=x_j\) se \(F_X(x_{j-1}) \leq U < F_X(x_j)\).

É importante salientar que o tempo de processamento para gerar números aleatórios pelo método descrito acima será proporcional à quantidade de intervalos que devem ser verificados. Por essa razão, pode ser razoável considerar os possíveis valores de \(x_j\) em ordem decrescente de \(p_j\).

Exemplo:

Simular a variável aleatória \(X\), com função de probabilidade \(p_j = \operatorname{P}\{ X = j\}\), tal que \(p_1 = 0,20\), \(p_2 = 0,15\), \(p_3 = 0,25\) e \(p_4 = 0,40\).

Procedimento A

  1. Gerar um número aleatório ao acaso, \(U\).

  2. Se \(U < 0.20\), então \(X = 1\) e pare,

  3. Se \(U < 0.35\), então \(X = 2\) e pare.

  4. Se \(U < 0.60\), então \(X = 3\) e pare

  5. Caso contrário, então \(X = 4\).

## [1] 0.6732237 4.0000000

Observe o gráfico abaixo e verifique que para \(U = 0.6732237\), temos que \(X = 4\).

Inversa de função de distribuição acumulada discreta

Inversa de função de distribuição acumulada discreta

Um procedimento mais eficiente, com ordenação decrescente das massas de cada ponto, seria o seguinte:

Procedimento B

  1. Gerar um número aleatório ao acaso, \(U\).

  2. Se \(U < 0.40\), então \(X = 4\) e pare,

  3. Se \(U < 0.65\), então \(X = 3\) e pare.

  4. Se \(U < 0.85\), então \(X = 1\) e pare

  5. Caso contrário, então \(X = 2\).

Gerador de Variável Aletória Uniforme Discreta:

Na construção de gerador de números aleatórios de uma distribuição uniforme discreta não é necessário proceder a uma busca por intervalos. Seja \(X\) uma variável aleatória em que todos seus valores possíveis, \(1, 2, \dots, n\), são equiprováveis, ou seja, \(\operatorname{P}\{X = j \}= \frac{1}{n}\), para \(j = 1, 2, \dots, n\). Aplicando o procedimento de inversão de função de distribuição acumulada discreta, temos que: \[ X = j \text{, se } \frac{j-1}{n} \leq U < \frac{j}{n} \]

Ou seja,

\[ X = j \text{, se } j-1 \leq nU < j \]

Ou seja,

\[ X = \lceil nU \rceil \]

em que, \(\lceil y \rceil\) é o menor inteiro maior ou igual a \(y\). De maneira mais formal, \(\lceil y \rceil = \max\{n \in \mathbb{N}; \, n \geq y\}\).

Gerador de Variável Aleatória Geométrica:

Seja \(X\) uma variável aleatória geométrica com parâmetro \(p\). \(X\) pode ser associada à contagem de tentativas independentes, com mesma probabilidade de sucesso \(p\), até a ocorrência do primeiro sucesso. Sua função de probabilidade é: \[ p_i = \operatorname{P}\{X = i \} = pq^{i-1} \text{, }i \geq 1, \text{em que } q = 1-p \] Sua função de distribuição acumulada é:

\[\begin{align} F_X(j-1) = \operatorname{P}\{X \leq j-1 \} & = \sum_{i=1}^{j-1} \operatorname{P}\{X = i \} = 1 - \operatorname{P}\{X > j-1 \} \\ & = 1 - \operatorname{P}\{\text{as primeiras $j-1$ tentantivas fracassam}\} \\ & = 1 - q^{j-1}\text{, } j \geq 1. \end{align}\]

Um procedimento para gerar o valor de \(X\) é:

Gerador de Variável Aleatória Geométrica - Procedimento

  1. Gerar um número aleatório ao acaso, \(U\).

  2. Tornar \(X = j\) em que \(1-q^{j-1} \leq U < 1 - q^j\).

O passo 2 é equivalente a:

\[ q^j < 1 - U \leq q^{j-1}, \]

ou seja, podemos definir \(X\) por:

\[ X = \min\{j; \, q^j < 1 - U \} \]

Considerando que a função logarítmo é monótona, ou seja, se \(a < b\), então \(\log(a) < \log(b)\), a expressão para obter \(X\) pode ser dada por:

\[\begin{align} X & = \min\{j; \, j\log(q) < \log(1 - U) \} \\ & = \min\left\{j; \, j \frac{\log(1-U)}{\log(q)} \right\}. \end{align}\]

Perceba que o sinal da desigualdade muda, pois \(\log(q) <0\) já que \(0 < q < 1\). Podemos expressar \(X\) como:

\[ X = \left\lceil \frac{\log(1-U)}{\log(q)} \right\rceil \]

Finalmente, como \(1-U\) também é distribuída uniformemente no intervalo \((0,1)\), então:

\[ X \equiv \left\lceil \frac{\log(U)}{\log(q)} \right\rceil \]

tem distribuição geométrica com parâmetro \(p\).

Gerador de Sequência de Variáveis Aleatórias Independentes de Bernoulli:

Seja \(X_1, X_2, \dots, X_n\) uma sequência de variáveis aleatórias independentes, todas com distribuição de Bernoulli com parâmetro \(p\). A sequência é facilmente obtida gerando-se uma sequência de \(n\) números aleatórios \(U_1, U_2, \dots, U_n\), tomando-se:

\[ X_i = \begin{cases} 1, & \text{ se } U_i \leq p \\ 0, & \text{ se } U_i < p \end{cases}. \]

Há uma abordagem mais eficiente. Imagine que essa variáveis aleatórias de Bernoulli representam o resultado de tentativas em que \(X_i=1\), caso a tentativa \(i\) seja um sucesso ou, \(X_i = 0\) em caso contrário. Para gerar essas tentativas quando \(p\leq 1/2\), use o resultado obtido com o gerador da subseção anterior, para gerar a variável aleatória \(N\), igual à quantidade de tentativas até o primeiro sucesso, quando todas as tentatovas têm probabilidade de sucesso \(p\). Suponha que o valor simulado é \(N=j\). Se \(j>n\), faça \(X_i = 0\), para \(i = 1, 2, \dots, n\); se \(j\leq n\), faça \(X_i = \dots = X_{j-1}= 0, X_j = 1\); e, se \(j<n\), repita o procedimento para obter as \(n-j\) variáveis aleatórias de Bernoulli remanescentes. como desejamos gerar simultaneamente tantas variáveis de Bernoulli quanto possível, para \(p > 1/2\) geramos o número de tenativas até a primeira falha, ao invés de até o primeiro sucesso.

Essa mesma ideia pode ser aplicada quando \(X_i\) são variáveis aleatórias de Bernoulli independentes, mas não identicamente distribuídas. Para cada \(i = 1, \dots, n\), seja \(U_i\) o menos provável dos dois valores possíveis de \(X_i\). Ou seja, \(U_i = 1\) se \(\operatorname{P}\{X_i = 1\} \leq 1/2\), e \(U_i = 0\), caso contrário. Além disso, seja $p_i = {X_i = u_i } e seja \(q_i = 1 − p_i\) . Simula-se a sequência de Bernoullis gerando primeiro o valor de \(X\), em que \(X = j\) para \(j = 1, \dots , n\), quando a tentativa \(j\) for a primeira tentativa que resulta em um valor improvável, e \(X = n+1\) se nenhuma das \(n\) tentativas resultar em seu valor improvável. Para gerar X, faça \(q_{n+1}=0\) e note que:

\[ \operatorname{P}\{X >j\} = \prod_{i=1}^j q_i \text{, } j= 1, \dots, n+1. \]

Portanto,

\[ \operatorname{P}\{X \leq j\} = 1 - \prod_{i=1}^j q_i \text{, } j= 1, \dots, n+1. \]

Consequentemente, pode-se simular \(X\) por meio da geração de um número aleatório \(U\), fazendo-se:

\[ X = \min \left\{j; \, U \leq 1 - \prod_{i=1}^j q_i \right\} \]

Se \(X=n+1\), a sequência de variáveis aleatórias de Bernoulli simulada será \(X_i = 1 - U_i\), para \(i = 1, \dots, n\). Se \(X=j\) j n$, faça \(X_i = 1 - u_i\), para \(i = 1, \dots, j-1\) e \(X_j = u_j\). Se \(j<n\), então gere de maneira similar os valores remanescentes de \(X_{j+1}, \dots, X_n\).

Reuso de Números aleatórios:

O procedimento discutido anteriormente seja mais eficiente para gerar os resultados de \(n\) tentativas independentes do que gerar uma variável aleatória uniforme para cada tentativa. Entretanto, em teoria, pode-se usar um único número aleatório para gerar todos os \(n\) resultados. Para fazer isso, Gere um número \(U\) ao acso e faça:

\[ X_1 = \begin{cases} 1, & \text{ se } U \leq p_1 \\ 0, & \text{ se } U < p_1 \end{cases}. \]

Sabe-se que a distribuição condicional de \(U\) dado que \(U \leq p\) é uma uniforme no intervalo \((0, p)\). Conseqüentemente, dado que \(U \leq p_1\), a razão \(\frac{U}{p_1}\) é uma uniforme no intervalo \((0, 1)\). De maneira similar, sabe-se que a distribuição condicional de \(U\) dado que \(U > p\) é uma uniforme no intervalo \((p, 1)\). Assim, \(U\), condicionada a \(U > p1\), tem-se que a razão \(\frac{U-p_1}{1-p_1}\) é uniforme em (0,1). Dessa maneira, em teoria, podemos usar um único número aleatório \(U\) para gerar os resultados das \(n\) tentativas da seguinte forma:

Geração com Reuso de Números Aleatórios - Procedimento

  1. \(I\) = 1.

  2. Gerar um número aleatório ao acaso, \(U\).

  3. Se \(U \leq p_I\), faça \(X_1 = I\), caso contrário \(X_I=0\).

  4. Se \(I=n\), pare.

  5. Se \(U\leq p_I\), faça \(U= \frac{U}{p_I}\), caso contrário \(U= \frac{U-p_I}{1-p_I}\).

  6. \(I = I+1\).

  7. Retorne o passo 3.

Há, no entanto, um problema prático ao reutilizarmos um único número aleatório. Os computadores fornecem os números aleatórios apenas até uma certa quantidade de casas decimais. Dessa maneira, os erros de arredondamento podem levar a variáveis transformadas menos uniformes. Por exemplo, suponha no procedimento indicado que todo \(p_i = 0,5\). Então \(U\) é transformado em \(2U\) se \(U \leq 0,5\), ou \(2U − 1\) se \(U > 0,5\). Consequentemente, se o último dígito de \(U\) for \(0\), ele permanecerá \(0\) na próxima transformação. Além disso, se o penúltimo dígito se tornar \(5\), ele será transformado em \(0\) na próxima iteração e, portanto, a partir de então os últimos 2 dígitos sempre serão \(0\), e assim por diante. Assim, após um grande números de iterações, todos os números aleatórios gerados podem terminar iguais a \(1\) ou \(0\). (Uma possível solução seria usar \(2U − 0,999...9\) em vez de \(2U − 1\).)

Eficiência da transformação Inversa para Variáveis Discretas:

Em vez de usar uma tabela com os pontos de massa da distribuição, podemos utilizar outros métodos eficientes para buscar o valor da função inversa da distribuição. Muitas vezes, a busca pode ser melhorada a partir do conhecimento das probabilidades dos pontos. A idéia básica é começar a partir de um ponto com alta probabilidade de satisfazer a relação de inversão de distribuição discreta (\(F_X(x_{j-1}) \leq U < F_X(x_j)\)). Dessa maneira, a moda da distribuição é uma escolha óbvia para iniciar o procedimento de busca, especialmente se a probabilidade de ocorrência do valor modal for bastante alta.

Para muitas distribuições discretas de interesse, pode haver uma relação recursiva simples entre as probabilidades de pontos de massa adjacentes:

\[ F_X(x_j) = f(F_X(x_{j-1})) \text{, para } j = 1, 2, \dots \]

onde \(f\) é uma função simples (assumindo-se \(x_0\) é o menor valor com massa positiva). Por exemplo, na distribuição de Poisson temos que:

\[ F_X(i) = \frac{\lambda F_X(i-1)}{i} \]

Para este caso, Kemp (1981) propõe duas abordagens: build-up search (busca acumulada) e chp-down. Note que qualquer um desses métodos (detalhados abaixo) pode ser modificado para iniciar em qualquer ponto, como, por exemplo, a moda.

Algoritmo de busca build-up para distribuições discretas

No método de busca build-up, função de distribuição acumulada da variável aleatória discreta é construída pelo cálculo recursivo da função de probabilidade de cada ponto.

Procedimento de busca Build-up

Defina \(t = F_X(x_0)\)

  1. Gerar um número aleatório ao acaso, \(U\) e faça \(X=x_0\), \(F_X = t\) e \(s = F_X\).

  2. Se \(U \leq s\), então \(X\) é o valor desejado.

  3. Caso contrário, faça \(X = X + 1\), \(F_X = f(F_X)\) e \(s = s + F_X\) e retorne ao passo 2.

Algoritmo chop-down para distribuições discretas

O cálculo recursivo de probabilidades para acelerar a busca, é utilizado também no método chop-down, no qual a variável uniforme gerada é diminuída por um valor igual à função de distribuição acumulada de interesse. O algoritmo abaixo resumo esse tipo de procedimento:

Procedimento de busca chop-down

Defina \(t = F_X(x_0)\)

  1. Gerar um número aleatório ao acaso, \(U\) e faça \(X=x_0\), \(F_X = t\).

  2. Se \(U \leq F_X\), então \(X\) é o valor desejado.

  3. Caso contrário, faça \(U = U - F_X\), \(X = X + 1\), \(F_X = f(F_X)\) e retorne ao passo 2.

Gerador de Variável Aleatória de Poisson:

Seja a variável aleatória \(X\) com distribuição de Poisson de média \(\lambda\). Sua função de probabilidade é:

\[ p_i = \operatorname{P}\{X = i \} = \operatorname{e}^{-\lambda}\frac{\lambda^i}{i!} \text{, } i = 0, 1, \dots \]

Salienta-se que, para calcular as probabilidades de Poisson, pode ser usada a seguinte fórmula recursiva:

\[ \frac{p_{i+1}}{p_i} = \frac{\operatorname{e}^{-\lambda}\frac{\lambda^{i+1}}{(i+1)!}}{\operatorname{e}^{-\lambda}\frac{\lambda^i}{i!}}=\frac{\lambda}{i+1} \]

ou, equivalentemente, \[ p_{i+1} = \frac{\lambda}{i+1}p_i \text{, } i > 0 \text{, com } p_0 = \operatorname{e}^{-\lambda}. \]

Gerador de Variável Aleatória de Poisson - Procedimento A

A quantidade \(i\) refere-se ao valor em consideração da variável, \(p = p_i\) é a probabilidade de que \(X\) é igual a \(i\) e F.x é a probabilidade de que \(X\) é menor ou igual a \(i\).

  1. Gerar um número aleatório ao acaso, \(U\).

  2. Faça \(i = 0\), \(p = \operatorname{e}^{-\lambda}\), \(Fx = p\).

  3. Se \(U < F\), faça \(X = i\) e pare.

  4. Faça \(p = \frac{\lambda p}{i+1}\), \(F = F + p\), \(i = i+1\) (tomar o sucessor de \(i\)).

  5. Retorne ao passo 3.

De maneira geral esse algoritmo, que usa o princípio da transformação inversa pode ser usado para qualquer distribuição discreta. Assim, para gerar \(X \sim P_\theta\), em que \(P_\theta\) tem suporte nos números inteiros, pode-se calcular as probabilidades (assumindo que é possível armazená-las):

\[ p_0 = P_\theta\{X\leq 0 \}, \quad p_1 = P_\theta\{X\leq 1 \}, \quad p_2 = P_\theta\{X\leq 2 \}, \quad \dots, \]

em seguida gera-se \(U \sim \mathcal{U}(0, 1)\) e escolhe-se:

\[ X = k \text{, se } p_{k-1} < U < p_k. \]

Exemplo

Para gerar \(X \sim \mathcal{P}(7)\), tomam-se os valores:

\[ p_0 = 0,0009; p_1 = 0,0073; p_2 = 0,0296; \dots ;, \]

concluindo a sequência quando atinge-se 1, com uma dada quantidade de casas decimais. (Por exemplo \(p_{20} = 0,999985\).)

Fonte: Example 2.4, de ROBERT e CASELLA (2010)

Esse algoritmo (Procedimento A) verifica sucessivamante (a partir de \(X=0\)) qual o valo da função de distribuição acumulada de uma Poisson contempla a probabilidade sorteada (\(U\)). Assim, a quantidade necessária de comparações será uma unidade maior que o valor gerado da Poisson. Portanto, em média, serão necessária \(\lambda + 1\) buscas. Esse fato pode ser razoável quando \(\lambda\) é pequeno. Muitas vezes podemos melhorar o algoritmo acima por uma escolha criteriosa de quais probabilidades devem ser calculadas primeiro. Por exemplo, se quisermos gerar variáveis aleatórias de uma distribuição de Poisson com média \(\lambda = 100\), o algoritmo acima é muito ineficiente. Isso ocorre porque esperamos que a maioria das observações esteja no intervalo \(\lambda \pm 3 \sqrt{\lambda}\). Para \(\lambda = 100\) esse intervalo é \((70, 130)\). Assim, começar em \(0\) quase sempre produzirá 70 verificações inúteis de \(p_{k-1} < U < p_k\), porque serão quase certamente rejeitadas. Uma primeira solução seria descartar o que está fora de um intervalo altamente prováve. No exemplo, tudo que estiver fora do intervalo \((70, 130)\), ou seja:

\[ \operatorname{P}\{X <70 \}+ \operatorname{P}\{X >100 \} = 0,00268. \] Formalmente, devemos encontrar um limite inferior e um limite superior para tornar essa probabilidade pequena o suficiente, mas informalmente \(\pm 3 \signma\) funciona bem.

Exemplo

Fonte: Example 2.5, de ROBERT e CASELLA (2010)

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

Nsim <- 1E4
lambda <- 100

raio <- 3*sqrt(lambda)

t <- round(seq(max(0, lambda - raio), lambda + raio, 1))

prob <- ppois(t, lambda)

X <- rep(0, Nsim)

for (i in 1:Nsim){
    U <- runif(1)
    X[i] = t[1] + sum(prob < U) 
}

A última linha do código verifica a quais intervalos da Poisson pertence o valor do número aletório gerado \(U\) e atribui o valor correto da variável \(X\) de Poisson.

Uma maneira mais formal para a evitar a ineficiência de iniciar em \(p_0\) as probabilidades acumuladas é começar a partir da moda da distribuição discreta \(P_k\), explorando os valores vizinhos até que a probabilidade acumulada seja 1, considerado um erro de aproximação. Usa-se então o fato de que uma variável aleatória de Poisson pode ter uma ou duas modas próximos a \(\lambda\). Um algoritmo mais eficiente verificaria primeiro um desses valores, ao invés de ser iniciado em \(0\). Esse algoritmo, faz \(I= \lfloor i \rfloor\) e calcula \(p_I\) (tomando os logaritmos, elevando o resultado à potência de \(\operatorname{e}\). Ele então usa a fórmula recursiva para determinar \(F(I)\). Agora, pode-se gerar uma variável aleatória \(X\), com distribuição de Poisson com média \(\lambda\), gerando um número ao acaso \(U\), verificando se \(X \leq I\), se \(U\leq F(I)\). Procurando o valor de \(X\) em sentido decrescente, no caso em que \(X\leq I\) e, caso contrário, crescente, partindo de \(I+1\). Esse algorimo também evita a potencial dificuldade que o computador, quando \(\lambda\) é grande, possa tornar \(\operatorname{e}^{-\lambda}\) igual a zero, devido a erros de arredondamento.

A quantidade de buscas do algoritmo é uma unidade maior que a diferença absoluta entre a variável aleatória \(X\) e sua média \(\lambda\). Pelo Teorema central do Limite, para \(\lambda\) grande \(X\) é aproximadamente normal com média e variância iguais a \(\lambda\). Assim:

\[\begin{align} & \text{Qte. média de buscas} \simeq 1 + \operatorname{E}[|X - \lambda|] \text{, com } X \sim \mathcal{N}(\lambda, \lambda) \\ & = 1 + \sqrt{\lambda} \operatorname{E}\left[ \frac{|X - \lambda|}{\sqrt{\lambda}}\right] \\ & = 1 + \sqrt{\lambda} \operatorname{E}[|Z|] \text{, em que } Z \sim \mathcal{N}(0, 1) \\ $ = 1 + 0,798 \sqrt{\lambda} \end{align}\]

Usando esse aprimoramento do algoritmo, a quantidade média de buscas cresce com a raiz quadrada de \(\lambda\) ao invés de \(\lambda\), o que é preferível quando \(\lambda\) é grande,

Gerando Poisson por meio de variáveis aleatórias exponenciais;

Proposition [ROSS (2006), pg. 30/315]

O tempo entre chegadas de um processo de poisson, \(T_1, T_2, \dots\) são variáveis aleatórias independentes e exponencialmente distribuídas com taxa \(\lambda\).

Seja \(S_n = \sum_{i=1}^nT_i\) o instante da ocorrência do n-ésimo evento. Dessa maneira \(S_n\) será menor ou igual a \(t\) se, e somente se, tenha ocorrido pelo menos \(n\) eventos até o instante t, ou seja:

\[\begin{align} \operatorname{P}\{ S_n \leq t \} & = \operatorname{P}\{N(t) \geq n \} \\ & = \sum_{j=n}^\infty \operatorname{e}^{-\lambda t} \frac{(\lambda t)^j}{j!} \end{align}\]

Como a expressão acima é a função de distribuição acumulada de \(S_n\), sua derivada fornece a função de densidade de \(S_n\), ou seja:

\[\begin{align} f_{S_n}(t) & = \sum_{j=n}^\infty j\lambda \operatorname{e}^{-\lambda t} \frac{(\lambda t)^{j-1}}{j!} - \sum_{j = n}^\infty \lambda \operatorname{e}^{-\lambda t} \frac{(\lambda t)^{j}}{j!} \\ & = \lambda \operatorname{e}^{-\lambda t} \left[ \sum_{j=n}^\infty \frac{(\lambda t)^{j-1}}{(j-1)!} - \sum_{j = n}^\infty \frac{(\lambda t)^{j}}{j!} \right] \\ & = \lambda \operatorname{e}^{-\lambda t} \frac{(\lambda t)^{n -1}}{(n -1)!} \end{align}\]

Verifica-se então que \(f_{S_n}\) é a função de densidade de probabilidade de uma distribuição gama, com parâmetros \(n\) e \(lambda\). Por outro lado, essa é a distribuição da soma de \(n\) variáveis aleatórias independentes e exponencialmente distribuídas com média \(1/\lambda\).

Dessa maneira, pode-se construir outro algoritmo para gerar uma variável aleatória de Poisson. Em um processo de Poisson com taxa λ, tem-se que os tempos entre eventos sucessivos são variáveis aleatórias independentes, distribuídas exponencialmente com taxa λ. (Consulte a Seção 2.9 do Capítulo 2.) Para tal processo de Poisson, \(N(1)\) denota o número de eventos no tempo \(1\), é distribuído de Poisson com média λ. Seja a sequência de variáveis aleatória \(\{ T_i\}, \, i = 1, 2, \dots\), associada aos tempos sucessivos entre chegadas, independentes e exponencialmente distribuídos com taxa \(\lambda\). Dessa maneira, o n-ésimo evento ocorrerá no instante \(\sum_{i=1}^nT_i\) e a quantidade de eventos em uma unidade de tempo \(N(1)\) pode ser expressa como:

\[ N(1) = \max\left\{ n; \sum_{i=1}^n T_i \leq 1 \right\} \]

Isso é, o número de eventos no instante 1 é igual ao maior inteiro \(n\), relativo ao último evento que tenha ocorrido até o instante 1. Por exemple, se o quarto evento ocorreu no intervalo de tempo até o instante 1, mas o quinto evento ocorreu após esse instante, então ocorreram 4 eventos até o instante 1).

Dessa maneira, podemos utilizar o procedimento descrito a seguir para gerar uma variável aleatória de Poisson com média \(\lambda\), \(N(1)\):

  1. Gerar \(U_1, \dots, U_n, \dots\), uma sequência de números escolhidos ao acaso no intervalo \((0,1)\).

  2. Determinar o valor da variável aleatória \(N\), que conta eventos ocorridos ao acaso em uma unidade de tempo (variável aleatória de Poisson).

\[\begin{align} N & = \max\left\{ n; \, \sum_{i=1}^n - \frac{1}{\lambda} \log U_i \leq 1 \right\} \\ & = \max\left\{ n; \, \sum_{i=1}^n \log U_i \geq - \lambda \right\} \\ & = \max \left\{n; \, \log\left(\prod_{i=1}^n U_i \right) \geq - \lambda \right\} \\ & = \max\left\{n; \, \prod_{i=1}^n U_i \geq \operatorname{e}^{-\lambda} \right\} \end{align}\]

Portanto, os valores aleatórios de \(N\), variável com distribuição Poisson, com média \(\lambda\) podem ser determinados por meio da geração de números aleatórios escolhidos ao acaso no intervalo \((0,1)\) até que seu produto seja acima de \(\operatorname{e}^{-\lambda}\), tomando-se o valor de \(N\) como a quantidade de números aleatórios utilizados (nas gerações sucessivas) menos um, ou seja:

\[ N = \min\left\{ \prod_{i=1}^n U_i < \operatorname{e}^{-\lambda}\right\} - 1. \]

Esse método é simples, mas, na realidade, ele é prático apenas para pequenos valores de \(\lambda\). Em média, são necessárias \(\lambda\) variáveis aleatórias exponenciais e o custo dessa geração pode ser proibitivo para grandes valores de \(\lambda\). DEVROYE (1981) propôs um método cujo tempo computacional é uniformemente limitado.

Será apresentada posteriormente uma abordagem pelo Método da Aceitação/Rejeição para gerar variável aleatória de Poisson para grandes valores de \(\lambda\).

Gerador de Variável Aleatória Binomial

Suponha que \(X\) tenha distribuição binomial com parâmetros \(n\) e \(p\). Sua função de probabilidade é:

\[ p_i = \operatorname{P}\{X = i \}= \binom{n}{i} p^i(1-p)^{n-i}\text{, } i = 1, 2, \dots, n. \]

A fórmula recursiva abaixo, expressando \(p_{i+1}\) em termos de \(p_i\) é útil quando são necessários cálculos de probabilidades binomiais;

\[\begin{align} p_{i+1} & = \binom{n}{i+1} p^{i+1}(1-p)^{n-i-1} \\ & = \frac{n!}{(n-i-1)!(i+1)!}p^{i+1}(1-p)^{n-i-1} \\ & = \frac{n!(n-i)}{(n-i)! i! (i+1)}p^{i}(1-p)^{n-i} \frac{p}{1-p} \\ & = \frac{n-i}{i+1} \frac{p}{1-p} \frac{n!}{(n-i)! i!}p^{i}(1-p)^{n-i} \\ & = \frac{n-i}{i+1} \frac{p}{1-p} p_i \end{align}\]

Com auxílio dessa identidade recursiva, podemos empregar o método da transformação inversa.

Gerador de Variável Aleatória Binomial - Procedimento

A quantidade \(i\) refere-se ao valor em consideração da variável, \(p = p_i\) é a probabilidade de que \(X\) é igual a \(i\) e \(F = F(i)\) é a probabilidade de que \(X\) é menor ou igual a \(i\).

  1. Gerar um número aleatório ao acaso, \(U\).

  2. Faça \(c= p/(1-p)\), \(i = 0\), \(p = (1-p)^n\), \(F = p\).

  3. Se \(U < F\), faça \(X = i\) e pare.

  4. Faça \(p = \frac{c(n-1)}{i+1}p\), \(F = F + p\), \(i = i+1\) (tomar o sucessor de \(i\)).

  5. Retorne ao passo 3.

Esse procedimento verifica inicialmente se \(X = 0\), depois se \(X = 1\) e assim por diante. Portanto, o número de buscas que ele faz é uma unidade a mais que o valor de \(X\). Dessa maneira, serão necessárias em média \(1 + np\) buscas para gerar \(X\). A variável aleatória binomial \(X\) pode ser associada ao número de sucessos em \(n\) tentativas independentes, cada uma com mesma probabilidade de sucesso \(p\). Salienta-se que essa variável aleatória também pode ser gerada subtraindo de \(n\) o valor de uma variável aleatória binomial \((n, 1 − p)\) Esse fato, pode ser usado quando \(p > 1/2\), usando o procedimento acima para gerar uma variável aleatória binomial \((n, 1 − p)\) e subtraindo seu valor de \(n\) para obter a o número aleatório \(X\) desejado.

Observações:

  • Outra maneira de gerar uma variável aleatória binomial \(X\) é utilizando sua associação à quantidade de sucessos em \(n\) tentativas de Bernoulli independentes, com mesma probabilidade de sucesso \(p\). Assim, podemos simular \(X\), contando a quantidade de sucessos na geração de uma sequência de \(n\) ensaios de Bernoulli.

  • Como no caso de Poisson, quando a média \(np\) é grande, é mais eficiente determinar inicialmente se o valor \(U\) gerado é menor ou igual a \(I ≡ \lfloor np \rfloor\), ou se \(U > I\) . No primeiro caso a busca é iniciada com \(I\) , a seguir \(I − 1\) e assim por diante. No último caso inicia-se a busca com \(I + 1\), em sentido crescente.

Método da Aceitação/Rejeição:

Suponha que haja um método eficiente para simular uma variável aleatória \(Y\), com função de probabilidade \(q_j, j\geq 0\). Essa função de probabilidade pode ser usada como base para simular variável aleatória \(X\), com função de probabilidade \(p_j, j \geq 0\), gerando inicialmente um valor de \(Y\) e então aceitando ou rejeitando esse valor com uma probabilidade \(\frac{p_Y}{q_y}\).

Especificamente, seja \(M\) uma constante tal que:

\[ \frac{p_j}{q_j} \leq M \text{, } \forall \, j \text{ tal que } p_j >0. \]

Com essas considerações, pode-se simular a variável aleatória \(X\), com função de probabilidade \(p_j\) por meio de procedimento denominado Método de Aceitação/Rejeição.

Método de Aceitação/Rejeição (Discretas) - Procedimento

  1. Simular o valor de \(Y\), com função de probabilidade \(q_j\).

  2. Gerar um número aleatório ao acaso, \(U\).

  3. Se \(U < \frac{p_Y}{M q_Y}\), faça \(X = Y\) e pare. Caso contrário, retorne ao passo 2.

Teorema:

Sejam \(X\) e \(Y\) variáveis aleatórias discretas, com função de probabilidade \(p_j\) e \(q_j\), respectivamente, para as quais \(R_X \subset R_Y\), tal que \(\frac{p_j}{q_j} \geq M\). \(\forall j \in R_Y\), com \(M\) sendo uma constante. Então:

\[ \operatorname{P}\left\{ Y = k | U < \frac{p_Y}{M q_Y} \right\}=p_k \text{, em que } U\sim \mathcal{U}(0,1) \text{ e } M \geq 1. \]

\(U < \frac{p_Y}{M q_Y}\) é chamada condição de aceitação.

Prova:

A probabilidade de aceitação é:

\[\begin{align} % \operatorname{P}\{\text{aceitação} \} & = \operatorname{P} \left\{ U < \frac{p_Y}{M q_Y} \right\} \\ & = \sum_{k \in R_Y} \operatorname{P} \left\{ U < \frac{p_Y}{M q_Y} \biggr| Y = k \right\} \operatorname{P}\{ Y = k\} \\ & = \sum_{k \in R_Y} F_U \left( \frac{p_k}{M q_k} \right) q_k \\ & = \sum_{k \in R_Y} \frac{p_k}{M q_k} q_k \\ & = \frac{1}{M} \sum_{k \in R_Y} p_k = \frac{1}{M}, % \end{align}\]

ou seja,

\[ \operatorname{P} \left\{ U < \frac{p_Y}{M q_Y} \right\} = \frac{1}{M} \]

A probabilidade envolvida com o Método da Aceitação/Rejeição é então:

\[\begin{align} % \operatorname{P} \left\{ Y = k \biggr| U < \frac{p_Y}{M q_Y} \right\} & = \frac{\operatorname{P} \left\{ U < \frac{p_Y}{M q_Y} \biggr| Y = k \right\} \operatorname{P}\{ Y = k\}}{ \operatorname{P} \left\{ U < \frac{p_Y}{M q_y}\right\} } \\ & = \frac{\operatorname{P}\left\{ U < \frac{p_Y}{M q_Y} \biggr| Y = k \right\} \operatorname{P}\{ Y = k\}} {\frac{1}{M}} \\ & = M F_U\left(\frac{p_k}{M q_k} \right) q_k \\ & = M \frac{p_k}{M q_k} q_k = p_k. % \end{align}\]

Observações:

  • Para o algoritmo funcionar, deve ser fácil geral um número aleatório de Y.

  • O algoritmo é uma variável aleatória geométrica com sucesso associado à aceitação do valor candidato \(Y\), ou seja, \(\operatorname{P}\{\text{sucesso} \} = \frac{1}{M}\). Portanto, o número esperado de tentativas até a aceitação (sucesso) para gerar um número aleatório da variável \(X\) é \(M\). Quanto menor seu valor, mais rápido será o algoritmo para gerar a amostra desejada.

Exemplo

Fonte: Example 4e, de ROSS (2013)

Suponha que seja necessário simular o valor de uma variável aleatória \(X\), com contradomínio \(R_X = \{ 1, 2, \dots, 10 \}\), cujas probabilidades são, respectivamente, \(0,11; 0,12; 0,09; 0,08; 0,12; 0,10; 0,09; 0,09; 0,10; 0,10\). Uma alternativa mais adequada ao Método da Transformação Inversa é usar o Método da Rejeição/Aceitação, com a distribuição uniforme discreta como a distribuição candidata \(q\), ou seja, \(q_j = \frac{1}{10}\), para \(j=1, 2, \dots, 10\). Para essa escolha de \(q_j\), podemos escolher

\[ M = \max\left\{\frac{p_j}{q_j} \right\} = 1,2 \]

O algoritmo terá então os seguintes passos:

Método de Aceitação/Rejeição (Discretas) - Exemplo

  1. Gerar um número aleatório ao acaso, \(U_1\) e fazer $Y= 10U_1 .

  2. Gerar um segundo número aleatório ao acaso, \(U_2\).

  3. Se \(U_2 \leq \frac{p_Y}{0,12}\), faça \(X = Y\) e pare. Caso contrário, retorne ao passo 1.

No passo 3, a constante é \(0,12\) já que \(M q_Y = \frac{1,2}{10}= 0,12\). Em média, esse algoritmo executará apenas 1,2 iterações para gerar um valor de \(X\).

Exemplo

Fonte: Example 2.3.9, de ROBERT e CASELLA (1999)

Aqui é apresentada uma alternativa que usa a relação entre a distribuição de Poisson e a distribuição logística, cuja função de densidade de probabilidade é:

\[ f(x) = \frac{1}{\beta}\frac{ \exp \left\{-\frac{(x-\alpha)}{\beta} \right\} }{ \left[ 1 + \exp \left\{-\frac{(x-\alpha)}{\beta} \right\} \right]^2 }, \]

e função de distribuição acumulada

\[ F(x)= \frac{1}{ 1 + \exp \left\{-\frac{(x-\alpha)}{\beta} \right\}}, \]

a qual é, portanto, analiticamente invertível.

Para relacionar melhor as duas distribuições, uma contínua e a outra discreta, consideramos que \(N = \lfloor x + 0,5 \rfloor\), ou seja, \(N\) é a parte inteira de \(x + 0,5\). Embora o suporte da logística seja \((-\infty, infty)\), ela será restrita ao intervalo \((-1/2, \infty)\) para se ajustar à Poisson. Portanto, a variável aleatória \(N\) tem função de distribuição:

\[ F_N(n) = \operatorname{P}\{N = n \} = \frac{1}{1+\operatorname{e}^{-(n+0,5-\alpha)/\beta}} - \frac{1}{1+\operatorname{e}^{-(n- 0,5-\alpha)/\beta}} \]

Se \(x > \frac{1}{2}\) e

\[ F_N(n) = \left( \frac{1}{1+\operatorname{e}^{-(n+0,5-\alpha)/\beta}} - \frac{1}{1+\operatorname{e}^{-(n- 0,5-\alpha)/\beta}} \right) \frac{1+\operatorname{e}^{-(0,5+\alpha)/\beta}}{\operatorname{e}^{-(0,5+\alpha)/\beta}} \]

se \(-\frac{1}{2} < x \leq \frac{1}{2}\) e a razão das densidades for :: \[ \frac{\lambda^n}{F_N(n)\operatorname{e}^\lambda n!} \] Embora seja difícil calcular um limite apara a expressão acima e, portanto, para otimizá-la em \((\alpha, \beta)\), ATKINSO (1979) propôs a escolha de \(\alpha = \lambda\) e \(\beta = \frac{\pi}{\sqrt{3\lambda}}\). Essas escolhas identificam os dois primeiros momentos de \(X\) com os de \(F_N(n)\). Embora essas escolhas de \(\alpha\) e \(\beta\) não possibilitem a otimização analítica do limite da expressão acima, por meio de maximização numérica e interpolação, é possível estabelecer que aquele limite é \(c = 0.767 - \frac{3.36}{\lambda}\), resultando então no procedimento abaixo:

Variável aleatória de Poisson - Procedimento de Atkinson

Defina \(\beta = \frac{\pi}{\sqrt{3\lambda}}\), \(\alpha = \lambda \beta\), e \(k = \log c - \lambda - \log \beta\)

  1. Gerar um número aleatório ao acaso, \(U_1\) e calcular: \[ x = \frac{\alpha - \log\left\{ \frac{1-U_1}{U_1} \right\} }{\beta} \text{, até } X > -0.5. \]

  2. Defina \(N= \lfloor X + 0.5 \rfloor\)

  3. Gerar um número aleatório ao acaso, $U_2.

  4. Aceitar \(N\) se

\[ \alpha - \beta X + \log \left( \frac{U_2}{ \left[ 1 + \exp (\alpha - \beta x) \right]^2} \right) \leq k + N\log \lambda - \log N!. \]

Embora a simulação resultante desse algoritmo seja exata, ele é baseado em um número de aproximações, tanto pela escolha de \((\alpha, \beta)\) quanto pelo cálculo dos limites e das razões das densidades. Perceba entretanto que esse procedimento exige o cálculo de fatoriais que pode ser um fator que aumente consideravelmente o custo computacional. Embora esse procedimento seja razoavelmente eficiente, pode ser preferível algoritmos mais complexos como aqueles estabelecidos em DEVROYE (1985).

Gerador de Variável Aleatória Binomial Negativa

Seja \(X\) uma variável aleatória binomial negativa com parâmetros \(r\) e \(p\). Sua função de probabilidade é:

\[ p_n = \operatorname{P}\{X = n \} = \binom{n - 1}{r - 1} p^r (1 - p)^{n-r}\text{, para } n = r, r+1, \dots, \]

onde \(r > 0\) e \(0 < p < 1\), em que \(p\) é a probabilidade de sucesso em uma única tentativa de Bernoulli. A distribuição binomial negativa também é denominada distribuição de Pascal. A distribuição binomial negativa é uma extensão da distribuição geométrica. Ela pode ser interpretada como o número de tentativas independentes antes de serem obtidos \(r\) sucessos.

Se \(r \frac{p}{1 - p}\) for relativamente pequeno e \((1 − π)^r\) não for muito pequeno, o método da transformação inversa funciona bem. Caso contrário, pode ser gerada uma variável aleatória gama, com parâmetros \((r, \frac{p}{1 - p})\), usando o número gerado como parâmetro para gerar um variável aleatória de Poisson. Esse valor é então uma variável aleatória binomial negativa. Note que esse procedimento baseia-se no fato que, quando \(Y \sim \mathcal{Ga}(n, \frac{1-p}{p})\) e \(X|y \sim \mathcal{P}(y)\), então \(X \sim \mathcal{BN}(n, p)\).

Gerador de Variável Aleatória Hipergeométrica:

Seja \(X\) uma variável aleatória hipergeométrica com parâmetros \(r\) e \(p\). Sua função de probabilidade é:

\[ p_i = \operatorname{P}\{X = i \} = \frac{\binom{M}{i} \binom{L-M}{N-i}}{\binom{L}{N}} \text{, para } i = \max\{0, N - L + M \}, \dots, \min\{N, M \} \]

O método usual de desenvolver a distribuição hipergeométrica é com um modelo de amostragem finita: N itens devem ser selecionados sem reposição, independentemente e com igual probabilidade, de um lote com um total de \(L\) itens, dos quais \(M\) são especiais (por exemplo, defeituosos); a variável aleatória \(X\) conta a quantidade de itens especiais na amostra aleatória de tamanho \(N\).

Um bom método para gerar a variável aleatória hipergeométrica é o método da transformação inversa. A inversa da função de distribuição acumulada pode ser avaliada usando recursivamente a razão \(\frac{p_{i + 1}}{p_i}\), por meio de busca de build-up (Algoritmo AA) ou pela busca chop-down (Algoritmo BB)s. Em ambos os casos, a busca pode ser acelerada, iniciando-se o algoritmo na média, \(\frac{MN}{L}\), pode acelerar a busca.

Outro método simples e eficiente é o uso do modelo de amostragem finita que define a distribuição.

Kachitvichyanukul e Schmeiser (1985) fornecem um algoritmo baseado no Método da Aceitação/Rejeição de uma função de probabilidade decomposta como uma mistura, e Stadlober (1990) descreve um algoritmo baseado em um método de razão de uniformes. Ambos podem ser mais rápidos que o CDF inverso para valores maiores de N e M. Kachitvichyanukul e Schmeiser (1988b) fornecem um programa Fortran para amostragem da distribuição hipergeométrica. O programa usa a inversa da função de distribuição acumulada ou o Método da Aceitação/Rejeição, dependendo do modo,

\[ m = \left\lceil \frac{(N+1)(M+1)}{L+2} \right\rceil \]

Se m−max(0,N+M−L) < 10, então é usado o Método da Transformação Inversa; caso contrário, é usado o Método da Aceitação/Rejeição.

Mistura de Distribuições Discretas

Suponha que esteja disponível um método eficiente para simular o valor de uma variável aleatória com uma de duas funções de probabilidade \(\{p^1_j\}\) ou ${p^2_j} e que desejamos simular o valor da variável aleatória \(X\), cuja função de probabilidade é

\[ P_j = \operatorname{P}\{ X = j\} = \alpha p^1_j + (1-\alpha)p^2_j \text{, } j\geq 0, \]

em que \(0 < \alpha < 1\).

Sejam \(X_1\) e \(X_2\) variáveis aleatórias com funções de probabilidade \(\{p^1_j\}\) e ${p^2_j}, respectivamente. Então, a variável aleatória \(X\) pode ser definida como:

\[ X = \begin{cases} X_1, & \text{ com probabilidade } \alpha \\ X_2, & \text{ com probabilidade } 1 - \alpha \end{cases} \]

Dessa maneira, podemos gerar o valor dessa variável aleatória \(X\), gerando primeiro um número aleatório ao acaso \(U\) e então gerando um valor de \(X_1\), se \(U < \alpha\) ou de \(X_2\), se \(U > \alpha\)

Exemplo

Fonte: Example 4f, de ROSS (2013)

Deseja-se gerar um valor da variável aleatória \(X\) tal que:

\[ P_j = \operatorname{P}\{ X = j\} = \begin{cases} 0,05 & \text{, para } j = 1, 2, 3, 4, 5 \\ 0,15 & \text{, para } j = 6, 7, 8, 9, 10 \end{cases} \]

Perceba que \(p_j = 0,5 p_j^1 + 0,5 p_j^2\), onde \(p_j^1 = 0.1\), para \(j = 1, 2, \dots, 10\) e \(p_j^2 = \begin{cases} 0 & \text{, para } j = 1, 2, 3, 4, 5 \\ 0.2 & \text{, para } j = 6, 7, 8, 9, 10 \end{cases}\). Dessa maneira, incialmente geramos ao acaso um número aleatório \(U\) e então, se \(U < 0.5\), geramos um valor de uma variável aleatória uniforme discreta sobre o conjunto \(\{1, 2, \dots, 10\}\). Caso contrário, geramos um valor de uma variável aleatória uniforme discreta sobre o conjunto \(\{6, 7, 8, 9, 10 \}\). Ou seja, podemos simular \(X\) de acordo com o algoritmo abaixo:

  1. Gerar um número aleatório ao acaso \(U_1\).

  2. Gerar um número aleatório ao acaso \(U_2\).

  3. Se \(U_1 < 0.5\), faça \(X = \lceil 10 U_2 \rceil\). Caso contrário, faça \(X = \lceil 5 U_2 \rceil + 5\).

Método da Composição:

Seja uma mistura de distribuições, com função de distribuição acumulada \(F\), dada por:

\[ F(x)=\sum_{i=1}^n F_i(x), \]

em que, \(F_i\), \(1, 2, \dots, n\) são funções de distribuição acumulada e \(\alpha_i\), \(i= 1, 2, \dots, n\), probabilidades, tal que \(\sum_{i=1}^n \alpha_i = 1\).

Uma maneira de simular um valor a partir de \(F\) é primeiro simular uma variável aleatória \(I = i\), com probabilidade \(\alpha_i\), \(i= 1, 2, \dots, n\) e então simular a partir da distribuição \(F_I\). Por exemplo se \(I=j\), então simula-se o segundo valor de \(F_j\). Esse procedimento é denominado método da composição.

REFERÊNCIAS