Capítulo 2

Conceitos de Números Pseudoaleatórios

“The greatest value of a picture is when it forces us to notice what we never expected to see.”

– John Tukey, 1977

2.1. Introdução

Em aplicações computacionais, especialmente em simulações estocásticas, torna-se imprescindível dispor de uma fonte confiável de aleatoriedade. Entretanto, dado o caráter determinístico dos computadores, não é possível gerar números verdadeiramente aleatórios de forma autônoma. Para mitigar essa limitação, empregam-se os chamados números pseudoaleatórios — sequências numéricas produzidas por algoritmos determinísticos, projetados para reproduzir satisfatoriamente as propriedades estatísticas de sequências genuinamente aleatórias.

Definição 2.1 (Números Pseudoaleatórios - Knuth 2014). Seja \(\mathcal{S}\) o espaço das sementes possíveis e \(G: \mathcal{S} \to [0,1]^\mathbb{N}\) um algoritmo determinístico que, para uma semente \(s \in \mathcal{S}\), gera uma sequência \((u_n)_{n\ \in \ \mathbb{N}}\) tal que:

\[\begin{align}\\ u_n = f_n(s)\\\\ \end{align}\]

onde cada \(f_n\) é uma função computável. Diz-se que \((u_n)\) é uma sequência de números pseudoaleatórios se satisfaz, de modo aproximado, as propriedades estatísticas de uma sequência de variáveis aleatórias independentes e identicamente distribuídas segundo uma distribuição uniforme no intervalo \([0,1]\).

Os algoritmos \(G\) descritos na Definição 2.1 são conhecidos como geradores de números pseudoaleatórios (GNPAs ou Pseudo-Random Number Generators — PRNGs), e são frequentemente construídos a partir de relações recursivas do tipo:

\[\begin{align}\\ T(x_n) &= x_{n+1}\\\\ h(x_n) &= u_n \\\\ \end{align}\]

onde \(T\) é uma função de transição determinística e \(h\) é uma função de normalização que mapeia o domínio de \(x_n\) para o intervalo \([0,1]\). Os principais critérios para avaliação de um GNPA incluem: (1) uniformidade da distribuição dos \(u_n\); (2) independência aproximada entre os termos da sequência; e (3) período longo, isto é, um grande número de termos distintos antes da repetição.

Figura 2.1. Gerador de números aleatórios.

Fonte: Random Number Generator | Dilbert Cartoons. Acesso: Junho, 2025.

O propósito desses geradores, contudo, não se restringe à produção de números uniformemente distribuídos: eles fundamentam também a geração de variáveis aleatórias com distribuições arbitrárias, essencial para simulações estatísticas. Para isso, aplicam-se métodos de transformação que mapeiam números pseudoaleatórios uniformes em variáveis com a distribuição desejada. Uma das abordagens mais fundamentais nesse contexto é o método da transformação inversa, cuja ideia central consiste em aplicar a uma função \(F^{-1}\), chamada de função quantil, de uma distribuição contínua ao número pseudoaleatório uniforme \(U \sim \text{Uniforme}(0,1)\). Entretanto, muitas distribuições relevantes em estatística aplicada não dispõem de funções quantis analiticamente tratáveis, o que demanda o uso de métodos mais gerais. O método de rejeição (rejection sampling), por exemplo, constitui uma abordagem comumente adotada nessas situações. Esse procedimento gera candidatos a partir de uma distribuição auxiliar \(g(x)\), de fácil amostragem, aceitando-os com probabilidade proporcional à razão entre a densidade alvo \(f(x)\) e a densidade auxiliar \(g(x)\).

Esses métodos — transformação inversa e rejeição — formam os pilares da estatística computacional moderna, permitindo a geração controlada de variáveis aleatórias, e servindo como alicerces para algoritmos avançados, como amostragem de Gibbs e Metropolis-Hastings, amplamente utilizados em inferência Bayesiana. Este capítulo visa estabelecer os fundamentos matemáticos necessários para a implementação e avaliação destes métodos de geração de valores pseudoaleatórios, especialmente quando aplicados a modelagem estatśitica.

2.2. Método da Transformação Inversa

A geração de variáveis aleatórias com distribuições específicas, a partir de números pseudoaleatórios uniformemente distribuídos, constitui uma etapa essencial em procedimentos de simulação estocástica. Dentre os métodos disponíveis, destaca-se, tanto pela fundamentação teórica quanto pela simplicidade computacional, o método da transformação inversa. Essa abordagem baseia-se no princípio de que, sob condições adequadas de regularidade, é possível obter uma variável aleatória com a distribuição desejada por meio da aplicação da chamada função quantil — definida como a inversa generalizada da função de distribuição acumulada — a uma variável aleatória com distribuição uniforme no intervalo \((0,1)\). Embora frequentemente aplicado a variáveis aleatórias contínuas, o método da transformação inversa também pode ser estendido, com as devidas adaptações, ao caso de variáveis aleatórias discretas. Neste contexto, a função quantil assume uma forma discreta e não contínua, mas mantém sua interpretação como um mecanismo de correspondência entre probabilidades acumuladas e os valores do suporte da variável.

2.2.1. Função Quantil

---

Como discutido anteriormente, a função quantil é definida como a inversa generalizada da função de distribuição acumulada e é aplicável tanto a variáveis aleatórias contínuas quanto discretas. Em termos conceituais, essa função estabelece um mecanismo sistemático para associar níveis de probabilidade a valores específicos do suporte de uma distribuição de probabilidade. A Definição 2.2 a seguir formaliza matematicamente esse conceito no contexto de distribuições contínuas. O caso de variáveis aleatórias discretas, por sua vez, será tratado na Definição 2.3, com as devidas adaptações à descontinuidade inerente à função de distribuição acumulada nesse cenário.

Definição 2.2 (Função Quantil - Caso Contínuo - Devroye 1986). Seja \(F\) uma função de distribuição acumulada (FDA) contínua definida sobre \(\mathbb{R}\). A função quantil associada a \(F\), também denominada função inversa generalizada de \(F\) e denotada por \(Q\), é definida por:

\[\begin{align}\\ Q(u) = F^{-1}(u) = \inf\{x \in \mathbb{R}: F(x) \geqslant u, 0 < u < 1\}\\\\ \end{align} \]

Nesta definição, a função quantil \(Q\) associa a cada valor \(u\), no intervalo \((0,1)\), o menor valor \(x\) no suporte da distribuição, tal que a probabilidade acumulada até \(x\) seja, no mínimo, igual a \(u\).

Exemplo 2.1 (Distribuição Exponencial - Darmois 1935; Koopman 1936). Considere a variável aleatória contínua \(X\) com distribuição exponencial de parâmetro \(\lambda > 0\), denotada por \(X \sim \text{Exponencial}(\lambda)\). A FDA de \(X\) é dada por:

\[\begin{align}\\ F(x;\lambda) = 1 - e^{-\lambda x}, \quad x > 0\\\\ \end{align}\]

Suponha que nosso objetivo seja determinar a função quantil \(Q(u;\lambda)\), a qual, para cada \(u \in (0,1)\), fornece o valor \(x\) tal que \(F(x;\lambda) = u\). Para tal, resolve-se a equação:

\[\begin{align}\\ F(x;\lambda) = u \quad \Rightarrow \quad 1 - e^{-\lambda x} = u\\\\ \end{align}\]

Logo, isolando-se \(x\) nessa equação, obtém-se que a função quantil da distribuição exponencial é dada por:

\[\begin{align}\\ Q(u;\lambda) = F^{-1}(u;\lambda) = -\frac{1}{\lambda} \ln(1 - u)\\\\ \end{align}\]

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.2 (Distribuição Rayleigh - Hoffman e Karst 1975). Considere uma variável aleatória contínua \(X\) com distribuição Rayleigh de parâmetro \(\sigma > 0\), denotada por \(X \sim\) Rayleigh\((\sigma)\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(x; \sigma) = 1 - e^{-x^2/(2\sigma^2)}, \quad x > 0\\\\ \end{align}\]

Suponha que nosso objetivo seja determinar a função quantil \(Q(u;\sigma)\), a qual, para cada \(u \in (0,1)\), fornece o valor \(x\) tal que \(F(x;\sigma) = u\). Para tal, resolve-se a equação:

\[\begin{align}\\ F(x; \sigma) = u \quad \Rightarrow \quad 1 - e^{-x^2/(2\sigma^2)} = u\\\\ \end{align}\]

Logo, isolando-se \(x\) nessa equação, obtém-se que a função quantil da distribuição Rayleigh é dada por:

\[\begin{align}\\ Q(u; \sigma) = F^{-1}(u; \sigma) = \sigma \sqrt{-2 \ln(1 - u)}\\\\ \end{align}\]

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.3 (Distribuição de Pareto - Pareto 1964). Considere uma variável aleatória contínua \(X\) com distribuição de Pareto, caracterizada pelos parâmetros \(\alpha > 0\) (forma) e \(\beta > 0\) (escala), denotada por \(X \sim\) Pareto\((\alpha, \beta)\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(x; \alpha, \beta) = 1 - \left(\frac{\beta}{x}\right)^\alpha, \quad x \geqslant \beta\\\\ \end{align}\]

Suponha que nosso objetivo seja determinar a função quantil \(Q(u; \alpha, \beta)\), a qual, para cada \(u \in (0,1)\), fornece o valor \(x\) tal que \(F(x; \alpha, \beta) = u\). Neste caso, tem-se:

\[\begin{align}\\ F(x; \alpha, \beta) = u \quad \Rightarrow \quad 1 - \left(\frac{\beta}{x}\right)^\alpha = u\\\\ \end{align}\]

Logo, isolando-se \(x\) nessa equação, obtém-se que a função quantil da distribuição de Pareto é dada por:

\[\begin{align}\\ Q(u; \alpha, \beta) = F^{-1}(u; \alpha, \beta) = \frac{\beta}{(1 - u)^{1/\alpha}}\\\\ \end{align}\]

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.4 (Distribuição Lindley - Lindley 1958). Considere uma variável aleatória contínua \(X\) com distribuição Lindley de parâmetro \(\theta > 0\), denotada por \(X \sim\) Lindley\((\theta)\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(x; \theta) = 1 - \left(1 + \frac{\theta x}{\theta + 1}\right) e^{-\theta x}, \quad x > 0\\\\ \end{align}\]

Suponha que nosso objetivo seja determinar a função quantil \(Q(u; \theta)\), a qual, para cada \(u \in (0,1)\), fornece o valor \(x\) tal que \(F(x; \theta) = u\). Neste caso, tem-se que:

\[\begin{align}\\ F(x; \theta) = u \quad \Rightarrow \quad 1 - \left(1 + \frac{\theta x}{\theta + 1}\right) e^{-\theta x} = u\\\\ \end{align}\]

Para resolver essa equação em \(x\), defina \(y = \theta x\). Logo,

\[\begin{align}\\ \left(1 + \frac{y}{\theta + 1}\right) e^{-y} = 1 - u\\\\ \end{align}\]

Multiplicando ambos os lados por \(e^{y}\), a equação torna-se:

\[\begin{align}\\ 1 + \frac{y}{\theta + 1} = (1 - u) e^{y}\\\\ \end{align}\]

Rearranjando e isolando os termos, pode-se aplicar à função W de Lambert (Lambert 1758; Corless et al. 1996; Jodrá 2010) para resolver a equação em \(y\), uma vez que essa função é definida como a solução da equação \(z, e^{z} = c\), para \(z \in \mathbb{R}\). Então, após as manipulações algébricas necessárias, obtém-se, então, que a função quantil da distribuição Lindley é dada por:

\[\begin{align}\\ Q(u; \theta) = -1 -\frac{1}{\theta} - \frac{1}{\theta} W_{-1} \left\{- (1 + \theta)(1- u) e^{-(1+\theta)} \right\}\\\\ \end{align}\]

em que \(W_{-1}(\cdot)\) representa o ramo negativo da função W de Lambert.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Observação 2.1 (Aproximação Numérica da Função Quantil). Embora em diversos contextos clássicos — como os exemplos previamente apresentados das distribuições Exponencial, Rayleigh, de Pareto e Lindley — seja possível obter expressões analíticas fechadas para a função quantil, tal propriedade não é uma característica universal das distribuições contínuas. Em muitas situações práticas, a FDA apresenta uma estrutura funcional que inviabiliza a obtenção de sua inversa por métodos analíticos. Isso ocorre, por exemplo, quando a FDA envolve combinações não elementares de funções especiais, quando é definida implicitamente ou ainda quando sua forma depende de múltiplos parâmetros. Distribuições como, por exemplo, Nakagami (Nakagami 1960) e Exponencial-Poisson (Kuş 2007), bem como aquelas obtidas por procedimentos de ajuste empírico a dados observacionais, ilustram bem esse cenário. Para tais distribuições, a função quantil não pode ser expressa em termos analíticos. Nesses casos, recorre-se à aproximação numérica da função quantil, que consiste em resolver, para um dado número pseudoaleatório \(u \in (0,1)\), o problema de encontrar o valor de \(x\) que satisfaz:

\[\begin{align}\\ F(x; \theta) = u \quad \Longleftrightarrow \quad h(x) = F(x; \theta) - u = 0\\\\ \end{align}\]

em que \(\theta\) denota, de modo geral, o vetor de parâmetros que caracteriza a distribuição. Essa equação, em particular, configura um problema clássico de determinação de raízes de funções não-lineares. Para sua resolução, pode-se empregar algoritmos numéricos bem estabelecidos, como o método da bisseção, o método de Newton-Raphson (quando a derivada de \(F\) está disponível e é bem comportada) ou o método de Brent (Brent 1971), que combina estratégias de bisseção, interpolação e secante. A escolha do método apropriado, naturalmente, dependerá de considerações como, por exemplo, a suavidade da função, a disponibilidade de derivadas, e o controle de erros numéricos.

Exemplo 2.5 (Distribuição Exponencial - Darmois 1935; Koopman 1936). Considere uma variável aleatória contínua \(X\) com distribuição exponencial de parâmetro \(\lambda > 0\), denotada por \(X \sim\) Exponencial\((\lambda)\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(x;\lambda) = 1 - e^{-\lambda x}, \quad x > 0\\\\ \end{align}\]

Suponha que o objetivo aqui seja determinar, para um dado \(u \in (0,1)\), o valor do quantil \(Q(u;\lambda) = F^{-1}(u;\lambda)\), sem recorrer à expressão analítica explícita. Neste caso, devemos determinar a solução numérica da equação:

\[\begin{align}\\ F(x; \lambda) - u = 0\\\\ \end{align}\]

Para esse fim, será utilizado o ambiente R. Como primeiro passo, é necessário implementar a expressão analítica da função densidade de probabilidade da distribuição exponencial como um objeto do tipo function, conforme ilustrado no Código 2.1.

Código 2.1. Implementação em R da função densidade da distribuição exponencial.

# --------------------------------------------------------------
# Implementação da Função Densidade da Distribuição Exponencial
# --------------------------------------------------------------

# --- 0. Limpar todas as variáveis do console para evitar confrontos de nomes ---

rm(list = ls())    

# --- 1. Implmentação da densidade ---

my.dexp     <- function(x, lambda)
{
  dens      <- lambda * exp(-lambda * x)
  return(dens)
}

Em seguida, como segundo passo, é necessário implementar a FDA. Neste caso, uma vez que o interesse reside na obtenção numérica da função quantil, adota-se, também, a aproximação numérica da FDA por métodos de quadratura para integrais definidas, como a quadratura de Gauss-Kronrod. No ambiente R, tal procedimento pode ser realizado por meio da função integrate, a qual possui, entre seus argumentos:

• f: A função a ser integrada. Essa função deve ser unidimensional e contínua no intervalo de integração.
• lower: O limite inferior do intervalo de integração.
• upper: O limite superior do intervalo de integração.
• ...: Argumentos adicionais a serem passados para a função f.
• subdivisions: O número máximo de subdivisões permitidas para aproximar a função. O padrão é 100.
• rel.tol: A tolerância relativa desejada na solução.
• abs.tol: A tolerância absoluta desejada na solução.

O Código 2.2 apresenta a implementação desse procedimento para a distribuição exponencial.

Código 2.2. Implementação em R da função de distribuição acumulada da distribuição exponencial via aproximação numérica por métodos de quadratura.

# -----------------------------------------------------------------------------
# Implementação Numérica da Função de Distribuição da Distribuição Exponencial
# -----------------------------------------------------------------------------

# --- Aproximação numérica via quadratura de Gauss-Kronrod ---

my.pexp.numerical     <- function(q, lambda) 
{
  sapply(q, function(t) integrate(my.dexp, lower = 0, upper = t, lambda = lambda)$value)
}

Com a FDA definida numericamente, a obtenção da função quantil passa a ser possível pela determinação da raiz da equação \(F(x;\lambda) - u = 0\), para cada valor de \(u \in (0,1)\). Esse procedimento pode ser realizado utilizando a função uniroot que constitui uma procedure no ambiente R baseada em uma variação do método da bisseção. Seus principais argumentos incluem:

• f: A função para a qual se deseja encontrar a raiz. Essa função deve ser unidimensional e contínua em um intervalo específico.
• interval: Um vetor de dois elementos que define o intervalo inicial no qual a raiz deve ser procurada.
• ...: Argumentos adicionais a serem passados para a função f.
• tol: A tolerância relativa desejada na solução.
• maxiter: O número máximo de iterações permitidas.
• extendInt: Um valor lógico que indica se o intervalo de busca deve ser estendido automaticamente se a raiz não for encontrada no intervalo inicial. O padrão é TRUE.

O Código 2.3 apresenta a implementação desse procedimento para a distribuição exponencial.

Código 2.3. Implementação em R da função quantil da distribuição exponencial via aproximação numérica.

# ------------------------------------------------------------------------
# Implementação Numérica da Função de Quantil da Distribuição Exponencial
# ------------------------------------------------------------------------

# --- Uso do uniroot para a aproximação numérica da função quantil ---

my.qexp.numerical     <- function(p, lambda, lower = 0, upper = 100) 
{
  sapply(p, function(u) {root <- uniroot(function(x) my.pexp.numerical(x, lambda) - u, 
                                         interval = c(lower, upper))$root
                         return(root)})
}

Essa implementação nos leva a seguinte qustão: Será que a solução numérica obtida para o quantil constitui uma boa aproximação da solução analítica, considerando o mesmo valor de \(\lambda\)? Para responder a essa questão, será realizada uma comparação entre os valores do quantil obtidos numericamente e os valores exatos fornecidos pela função quantil analítica da distribuição exponencial. Neste experimento, adotam-se os parâmetros \(u = 0.9\) e \(\lambda = 2\). A implementação e os resultados dessa comparação estão apresentados no Código 2.4.

Código 2.4. Comparação entre as funções quantis, analítica e numérica, da distribuição exponencial.

# -----------------------------------------------------------------
# Quantil Analítico x Quantil Numérico da Distribuição Exponencial
# -----------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Cálculos comparativos ---

u_seq         <- seq(0.10, 0.90, by = 0.10)
qa_exp        <- qexp(u_seq, rate = 2)
qn_exp        <- my.qexp.numerical(u_seq, lambda = 2)
diff_q        <- round(qn_exp - qa_exp, 7)

# --- 2. Exibição dos resultados da comparação ---

tabela_res    <- data.frame("u" = u_seq,
                            "Quantil_Analitico" = qa_exp,
                            "Quantil_Numerico" = qn_exp,
                            "Diff" = diff_q)

kable(tabela_res,
      digits = 7,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("Valores de <br> u",
                    "Quantil Analítico <br> (λ = 2.0)",
                    "Quantil Numérico <br> (λ = 2.0)",
                    "Diferença Entre <br> Quantis")) 

Valores de u	Quantil Analítico (λ = 2.0)	Quantil Numérico (λ = 2.0)	Diferença Entre Quantis
0.1	0.0526803	0.0526802	0.00e+00
0.2	0.1115718	0.1115730	1.30e-06
0.3	0.1783375	0.1783515	1.40e-05
0.4	0.2554128	0.2554137	8.00e-07
0.5	0.3465736	0.3465773	3.70e-06
0.6	0.4581454	0.4581161	-2.93e-05
0.7	0.6019864	0.6019827	-3.70e-06
0.8	0.8047190	0.8047232	4.20e-06
0.9	1.1512925	1.1512995	6.90e-06

Os resultados obtidos evidenciam que a diferença entre os valores calculados da função quantil pelas abordagens analítica e numérica é da ordem mínima de \(10^{-5}\), o que corresponde a valor extremamente pequeno e numericamente desprezível. Tal resultado confirma a precisão e a robustez do procedimento de inversão numérica implementado, validando seu uso como alternativa viável em situações nas quais a inversa da FDA não é conhecida em forma fechada.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.6 (Distribuição Exponencial-Poisson - Kuş 2007). Considere uma variável aleatória contínua \(X\) com distribuição Exponencial-Poisson, parametrizada por \(\lambda > 0\) e \(\beta > 0\), cuja função densidade de probabilidade (FDP) é dada por:

\[\begin{align}\\ f(x;\lambda,\beta) = \frac{\lambda \beta}{1 - e^{-\lambda}} e^{-\lambda - \beta x + \lambda e^{-\beta x}}, \quad x > 0\\\\ \end{align}\]

com FDA descrita por:

\[\begin{align}\\ F(x;\lambda,\beta) = \left[e^{\lambda e^{-\beta x} - e^{\lambda}}\right]\left[{1 - e^{-\lambda}}\right]^{-1}, \quad x > 0\\\\ \end{align}\]

Note que a função quantil desta distribuição não admite expressão analítica fechada em termos elementares, o que inviabiliza sua inversão direta. Assim, para um dado valor \(u \in (0,1)\), a obtenção do quantil \(Q(u; \lambda, \beta)\) requer a resolução numérica da equação:

\[\begin{align}\\ F(x;\lambda,\beta) - u = 0\\\\ \end{align}\]

Para esse fim, será utilizado o ambiente R. Como primeiro passo, é necessário implementar a expressão analítica da FDP da distribuição exponencial-Poisson como um objeto do tipo function, conforme ilustrado no Código 2.5.

Código 2.5. Implementação em R da função densidade da distribuição exponencial-Poisson.

# ----------------------------------------------------------------------
# Implementação da Função Densidade da Distribuição Exponencial-Poisson
# ----------------------------------------------------------------------

# --- 0. Limpar todas as variáveis do console para evitar confrontos de nomes ---

rm(list = ls())    

# --- 1. Implmentação da densidade ---

my.dexp.poisson <- function(x, lambda, beta)
{
  d1    <- lambda * beta / (1 - exp(-lambda))
  d2    <- exp(-lambda - beta * x + lambda * exp(-beta * x))
  dens  <- d1 * d2
  return(dens)
}

Em seguida, como segundo passo, é necessário implementar a FDA. Neste caso, uma vez que o interesse reside na obtenção numérica da função quantil, adota-se, também, a aproximação numérica da FDA por métodos de quadratura para integrais definidas, como a quadratura de Gauss-Kronrod. No ambiente R, tal procedimento pode ser realizado por meio da função integrate. O Código 2.6 apresenta a implementação desse procedimento para a distribuição exponencial-Poisson.

Código 2.6. Implementação em R da função de distribuição acumulada da distribuição exponencial-Poisson via aproximação numérica por métodos de quadratura.

# -------------------------------------------------------------------------------------
# Implementação Numérica da Função de Distribuição da Distribuição Exponencial-Poisson
# -------------------------------------------------------------------------------------

# --- Aproximação numérica via quadratura de Gauss-Kronrod ---

my.pexp.poisson     <- function(q, lambda, beta) 
{
  sapply(q, function(t) integrate(my.dexp.poisson, 
                                  lower = 0, 
                                  upper = t, 
                                  lambda = lambda, 
                                  beta = beta)$value)
}

Com a FDA definida numericamente, a obtenção da função quantil pode ser aproximada numericamente resolvendo-se a equação \(F(x;\lambda,\beta) = u\) para cada valor \(u\). Para isso, utiliza-se o algoritmo de busca da raiz uniroot, conforme o Código 2.7.

Código 2.7. Implementação em R da função quantil da distribuição exponencial-Poisson via aproximação numérica.

# --------------------------------------------------------------------------------
# Implementação Numérica da Função de Quantil da Distribuição Exponencial-Poisson
# --------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Uso do uniroot para a aproximação numérica da função quantil ---

my.qexp.poisson     <- function(p, lambda, beta, lower = 0.001, upper = 100)
{
  sapply(p, function(u) {root <- uniroot(function(x) my.pexp.poisson(x, lambda, beta) - u, 
                                         interval = c(lower, upper))$root
                         return(root)})
}

# --- 2. Cálculo do quantil para diferentes valores de u, lambda e beta ---

u_seq         <- seq(0.10, 0.90, by = 0.10)

quantis_2_05  <- my.qexp.poisson(p = u_seq, lambda = 2, beta = 0.5, lower = 0, upper = 100)
quantis_05_2  <- my.qexp.poisson(p = u_seq, lambda = 0.5, beta = 2, lower = 0, upper = 100)
quantis_1_1   <- my.qexp.poisson(p = u_seq, lambda = 1, beta = 1, lower = 0, upper = 100)

tabela_res    <- data.frame("u" = u_seq,
                            "Quantil_2_05" = quantis_2_05,
                            "Quantil_05_2" = quantis_05_2,
                            "Quantil_1_1" = quantis_1_1)

kable(tabela_res,
      digits = 6,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("Valores de <br> u",
                    "Quantil Exponencial-Poisson <br> (λ = 2.0, β = 0.5)",
                    "Quantil Exponencial-Poisson <br> (λ = 0.5, β = 2.0)",
                    "Quantil Exponencial-Poisson <br> (λ = 1.0, β = 1.0)")) 

Valores de u	Quantil Exponencial-Poisson (λ = 2.0, β = 0.5)	Quantil Exponencial-Poisson (λ = 0.5, β = 2.0)	Quantil Exponencial-Poisson (λ = 1.0, β = 1.0)
0.1	0.092527	0.041845	0.067506
0.2	0.199498	0.089524	0.145213
0.3	0.325395	0.144648	0.236086
0.4	0.477075	0.209660	0.344586
0.5	0.665643	0.288252	0.477849
0.6	0.910577	0.386733	0.647888
0.7	1.250095	0.516913	0.877701
0.8	1.775380	0.705356	1.219450
0.9	2.796618	1.036895	1.841590

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Suponha que nosso interesse agora seja trabalhar com a função quantil no contexto de variáveis aleatórias discretas. Para esse tipo de variável, a definição da função quantil exige uma formulação específica, em virtude da natureza da FDA, sendo compreendida como o menor valor possível da variável aleatória cuja probabilidade acumulada seja maior ou igual a um dado nível \(u \in (0,1)\). Diferentemente do caso contínuo, em que a inversa da FDA constitui uma função contínua e estritamente crescente, no caso discreto essa inversa apresenta um formato em degraus, refletindo os saltos característicos da distribuição, o que impõe cuidados adicionais na análise e interpretação dos quantis nesse contexto.

Definição 2.3 (Função Quantil – Caso Discreto - Devroye 1986). Seja \(X\) uma variável aleatória discreta com suporte em \(\mathbb{Z}\) e FDA \(F: \mathbb{Z} \to [0,1]\), definida por:

\[\begin{align}\\ F(k) = P(X \leqslant k) = \sum_{i \ \leqslant \ k} P(X = i), \quad k \in \mathbb{Z}\\\\ \end{align}\]

Para um dado nível \(u \in (0,1)\), a função quantil \(Q(u)\) é definida como o menor valor inteiro \(k \in \mathbb{Z}\) que satisfaz:

\[\begin{align}\\ Q(u) = \min \left\{ x \in \mathbb{Z} : F(x) \geqslant u \right\}\\\\ \end{align}\]

Equivalentemente, \(Q(u)\) é o menor inteiro \(k\) tal que a soma acumulada das probabilidades até \(k\) ultrapassa ou atinge \(u\), enquanto a soma acumulada até \(k-1\) permanece estritamente menor que \(u\), isto é,

\[\begin{align}\\ \sum_{i \ < \ k} P(X = i) < u \leqslant \sum_{i \ \leqslant \ k} P(X = i)\\\\ \end{align}\]

Exemplo 2.7 (Distribuição de Poisson - Poisson 1837). Considere uma variável aleatória discreta \(X\) com distribuição de Poisson de parâmetro \(\lambda > 0\), denotada por \(X \sim \text{Poisson}(\lambda)\). A FDA de \(X\) é dada por:

\[\begin{align}\\ F(k; \lambda) = P(X \leqslant k) = \sum_{i=0}^{k} \frac{e^{-\lambda} \lambda^i}{i!}, \quad k \in \mathbb{N}\\\\ \end{align}\]

Nosso objetivo é determinar a função quantil \(Q(u; \lambda)\), que para um dado \(u \in (0,1)\) fornece o menor valor \(k \in \mathbb{N}\) tal que \(F(k; \lambda) \geqslant u\), conforme estabelecido na Definição 1.3. Assim, dado \(u \in (0,1)\), o procedimento para obter \(Q(u; \lambda)\) consiste em identificar o menor inteiro \(k\) tal que:

\[\begin{align}\\ \sum_{i=0}^{k-1} \frac{e^{-\lambda} \lambda^{i}}{i!} < u \leqslant \sum_{i=0}^{k} \frac{e^{-\lambda} \lambda^{i}}{i!}\\\\ \end{align}\]

Dessa forma, para cada \(u\), a função quantil \(Q(u; \lambda)\) é obtida encontrando o menor inteiro \(k\) para o qual a soma acumulada da função de probabilidade até \(k\) ultrapassa ou atinge \(u\). Por exemplo, suponha \(\lambda = 3\) e \(u = 0{.}65\). Calculando os valores acumulados:

\[\begin{align}\\ F(0;3) &= e^{-3} \approx 0{,}0498 \\\\ F(1;3) &= F(0) + \frac{e^{-3} \cdot 3^1}{1!} \approx 0{.}1991 \\\\ F(2;3) &= F(1) + \frac{e^{-3} \cdot 3^2}{2!} \approx 0{.}4232 \\\\ F(3;3) &= F(2) + \frac{e^{-3} \cdot 3^3}{3!} \approx 0{.}6472 \\\\ F(4;3) &= F(3) + \frac{e^{-3} \cdot 3^4}{4!} \approx 0{.}8153 \\\\ \end{align}\]

Como \(F(3;3) \approx 0{.}6472 < 0{.}65 < 0{.}8153 \approx F(4;3)\), temos que \(Q(0{,}65; 3) = 4\). Portanto, a função quantil de ordem \(u = 0{.}65\) da distribuição de Poisson com \(\lambda = 3\) é dado por \(Q(0{.}65; 3) = 4\). O Código 2.8 ilustra o cálculo da função quantil da distribuição de Poisson, assumindo o parâmetro \(\lambda = 3\), com base em duas abordagens complementares: (i) a aplicação direta da Definição 2.3, e (ii) o uso da função nativa qpois, disponível no ambiente R.

Código 2.8. Cálculo a função quantil da distribuição de Poisson, assumindo \(\lambda = 3\), por meio de duas abordagens: condições estabelecidas na Definição 2.3, e a utilização da função nativa qpois, disponível no ambiente R.

# -----------------------------------------------------------
# Implementação da Função Quantil da Distribuição de Poisson
# -----------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Sequência de valores para u e valor de lambda ---

u_seq           <- seq(0.10, 0.90, by = 0.10)
lambda          <- 3

# --- 2. Implementação da função quantil com base na Definição 2.3 ---

quantil_poisson <- function(u, lambda)
{
  k             <- 0:100
  cdf_values    <- ppois(k, lambda)
  quantil       <- k[min(which(cdf_values >= u))]
  return(quantil)
}

# --- 3. Cálculo do quantil para cada u usando função implementada e função nativa qpois ---

q_manual        <- sapply(u_seq, function(u) quantil_poisson(u, lambda))
q_nativa        <- qpois(u_seq, lambda)

tabela_res      <- data.frame("u" = u_seq,
                              "Quantil_Def23" = q_manual,
                              "Quantil_qpois" = q_nativa)

kable(tabela_res,
      digits = 6,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("Valores de <br> u",
                    "Quantil Poisson <br> (Definição 2.3)",
                    "Quantil Poisson <br> (Função qpois)"))

Valores de u	Quantil Poisson (Definição 2.3)	Quantil Poisson (Função qpois)
0.1	1	1
0.2	2	2
0.3	2	2
0.4	2	2
0.5	3	3
0.6	3	3
0.7	4	4
0.8	4	4
0.9	5	5

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.8 (Distribuição Geométrica - Moivre 1711). Considere uma variável aleatória discreta \(X\) com distribuição geométrica de parâmetro \(p \in (0,1)\), denotada por \(X \sim \text{Geom}(p)\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(k; p) = P(X \leqslant k) = \sum_{i=0}^{k} (1-p)^i p = 1 - (1-p)^{k+1}, \quad k \in \mathbb{N}\\\\ \end{align}\]

Nosso objetivo é determinar a função quantil \(Q(u; p)\), que para um dado \(u \in (0,1)\) fornece o menor valor \(k \in \mathbb{N}\) tal que \(F(k; p) \geq u\), conforme estabelecido na Definição 1.3. Isso equivale a identificar o menor inteiro \(k\) tal que:

\[\begin{align}\\ \sum_{i=0}^{k-1} (1-p)^i p < u \leqslant \sum_{i=0}^{k} (1-p)^i p\\\\ \end{align}\]

ou, de forma equivalente, utilizando a expressão fechada para a FDA:

\[\begin{align}\\ 1 - (1-p)^k < u \leqslant 1 - (1-p)^{k+1}\\\\ \end{align}\]

Assim, para cada \(u \in (0,1)\), a função quantil \(Q(u; p)\) é obtida determinando o menor inteiro \(k\) tal que:

\[\begin{align}\\ 1 - (1-p)^{k+1} \geqslant u \quad \Leftrightarrow \quad (1-p)^{k+1} \leqslant 1 - u\\\\ \end{align}\]

Tomando o logaritmo natural em ambos os lados:

\[\begin{align}\\ (k+1) \ln(1-p) \leqslant \ln(1 - u)\\\\ \end{align}\]

o que, como \(\ln(1-p) < 0\), implica:

\[\begin{align}\\ k \geqslant \left\lceil \frac{\ln(1 - u)}{\ln(1 - p)} - 1 \right\rceil\\\\ \end{align}\]

Por exemplo, suponha \(p = 0{.}3\) e \(u = 0{.}65\). Calculando os valores acumulados da distribuição geométrica:

\[\begin{align}\\ F(0; 0{.}3) &= 0{.}3 \\\\ F(1; 0{.}3) &= 0{.}3 + 0{.}3(0{.}7) = 0{.}51 \\\\ F(2; 0{.}3) &= F(1) + 0{.}3(0{.}7)^2 = 0{.}657 \\\\ F(3; 0{.}3) &= F(2) + 0{.}3(0{.}7)^3 = 0{.}7593\\\\ \end{align}\]

Como \(F(1) = 0{.}51 < 0{.}65 < 0{.}657 = F(2)\), segue que \(Q(0{.}65; 0{.}3) = 2\). Esse resultado também pode ser obtido diretamente pela equação de \(k\), isto é,

\[\begin{align}\\ k &\geqslant \left\lceil \frac{\ln(1 - 0{.}65)}{\ln(1 - 0{.}3)} - 1 \right\rceil = \left\lceil \frac{-1{.}0498}{-0{.}3567} - 1 \right\rceil \approx 2\\\\ \end{align}\]

Portanto, a função quantil de ordem \(u = 0{.}65\) da distribuição geométrica com \(p = 0{.}3\) é dado por \(Q(0{.}65; 0{.}3) = 2\). O Código 2.9 apresenta o cálculo da função quantil da distribuição geométrica, assumindo o parâmetro \(p = 0{.}3\), com base em duas abordagens complementares: (i) a aplicação direta da Definição 2.3, que conduz a uma ineequação para a determinação do valor de \(k\) correspondente, e (ii) o uso da função nativa qgeom, disponível no ambiente R.

Código 2.9. Cálculo da função quantil da distribuição geométrica, assumindo \(p = 0.3\), por meio de duas abordagens: a resolução da equação para \(k\), conforme estabelecido na Definição 2.3, e a utilização da função nativa qgeom, disponível no ambiente R.

# -----------------------------------------------------------
# Implementação da Função Quantil da Distribuição Geométrica
# -----------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Sequência de valores para u e valor de p ---

p             <- 0.3
u_seq         <- seq(0.10, 0.90, by = 0.10)

# --- 2. Implementação da função quantil via Definição 2.3 (busca discreta) ---

qdef_geom     <- function(u, p)
{
  k           <- 0:100
  cdf_values  <- pgeom(k, prob = p)
  quantil     <- k[min(which(cdf_values >= u))]
  return(quantil)
}

# --- 3. Implementação da função quantil via expressão analítica ---

qa_geom       <- function(u, p)
{
  k           <- ceiling(log(1 - u) / log(1 - p) - 1)
  return(k)
}

# --- 4. Cálculo do quantil para cada u usando as duas funções e a função nativa qgeom ---

q_manual      <- sapply(u_seq, function(u) qdef_geom(u, p))
q_analitico   <- sapply(u_seq, function(u) qa_geom(u, p))
q_nativa      <- qgeom(u_seq, prob = p)

tabela_res    <- data.frame("u" = u_seq,
                            "Quantil_Def23" = q_manual,
                            "Quantil_Analitico" = q_analitico,
                            "Quantil_qgeom)" = q_nativa)

kable(tabela_res,
      digits = 6,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("Valores de <br> u",
                    "Quantil Geométrica <br> (Definição 2.3)",
                    "Quantil Geométrica <br> (Versão Analítica)",
                    "Quantil Geométrica <br> (Função qgeom)"))

Valores de u	Quantil Geométrica (Definição 2.3)	Quantil Geométrica (Versão Analítica)	Quantil Geométrica (Função qgeom)
0.1	0	0	0
0.2	0	0	0
0.3	1	0	0
0.4	1	1	1
0.5	1	1	1
0.6	2	2	2
0.7	3	3	3
0.8	4	4	4
0.9	6	6	6

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Observação 2.2. Em aplicações práticas, quando a variável aleatória discreta \(X\) resulta de uma discretização de uma variável aleatória contínua \(Y\), ou seja, é construída por meio de métodos de discretização, torna-se viável empregar a aproximação \(Q(u) = \lfloor Q_c(u) \rfloor\), em que \(Q_c(u)\) representa a função quantil contínua associada a \(Y\) e \(\lfloor \cdot \rfloor\) representa a função parte inteira. Essa abordagem constitui uma alternativa conveniente para estimar quantis discretos a partir de funções contínuas relacionadas, assegurando maior coerência com a estrutura probabilística subjacente.

Exemplo 2.9 (Distribuição Lindley Discreta - Gómez-Déniz e Calderı́n-Ojeda 2011). Seja \(Y \sim\) Lindley\((\theta)\), com \(\theta > 0\), uma variável aleatória contínua cuja FDA é dada por:

\[\begin{align}\\ F_Y(y; \theta) = 1 - \left(1 + \frac{\theta y}{\theta + 1}\right)e^{-\theta y}, \quad y > 0\\\\ \end{align}\]

A partir dessa FDA, pode-se construir uma variável aleatória discreta \(X\) como versão discretizada de \(Y\), cuja função de probabilidade é definida por:

\[\begin{align}\\ P(X = x) &= F_Y(x + 1; \theta) - F_Y(x; \theta)\\\\ &= \dfrac{e^{-\theta x}}{1 + \theta}\left\{\theta\left(1 - 2e^{-\theta}\right) + \left(1 - e^{-\theta}\right)(1 + \theta x)\right\}, \quad x \in \mathbb{N}\\\\ \end{align}\]

Nesse caso, a FDA de \(X\) é definida por:

\[\begin{align}\\ F_X(k; \theta) = P(X \leqslant k) = P(Y < k+1) = F_Y(k+1; \theta), \quad k \in \mathbb{N}\\\\ \end{align}\]

Nosso objetivo é determinar a função quantil de \(X\), denotada por \(Q(u; \theta)\), para \(u \in (0,1)\). Conforme estabelecido na Observação 2.2, visto que \(X\) é obtida por discretização da distribuição Lindley contínua, uma aproximação natural para seu quantil é dada por:

\[\begin{align}\\ Q(u; \theta) \approx \lfloor Q_Y(u; \theta) \rfloor\\\\ \end{align}\]

em que \(Q_Y(u; \theta)\) representa a função quantil da distribuição Lindley contínua. Esta função pode ser expressa como:

\[\begin{align}\\ Q_Y(u; \theta) = -1 -\frac{1}{\theta} - \frac{1}{\theta} W_{-1} \left\{- (1 + \theta)(1- u) e^{-(1+\theta)} \right\}\\\\ \end{align}\]

em que \(W_{-1}(\cdot)\) denota o ramo negativo da função W de Lambert. Assim, a função quantil aproximada da distribuição Lindley discreta é descrita por:

\[\begin{align}\\ Q(u; \theta) \approx \lfloor Q_Y(u; \theta) \rfloor = \left\lfloor -1 -\frac{1}{\theta} - \frac{1}{\theta} W_{-1} \left\{- (1 + \theta)(1- u) e^{-(1+\theta)} \right\}\right\rfloor\\\\ \end{align}\]

O Código 2.10 apresenta o cálculo da função quantil da distribuição Lindley discreta, assumindo o parâmetro \(\theta = 0.7\), com base na expressão \(Q(u;\theta)\) e na função lambertWm1 do pacote lamW (Adler 2015) para o cálculo da função W de Lambert.

Código 2.10. Cálculo da função quantil da distribuição Lindley discreta, assumindo \(\theta = 0.7\), com base na expressão \(Q(u;\theta)\).

# -----------------------------------------------------------------
# Implementação da Função Quantil da Distribuição Lindley Discreta
# -----------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(lamW)
library(knitr)
library(kableExtra)

# --- 1. Sequência de valores para u e valor de p ---

theta         <- 0.7
u_seq         <- seq(0.10, 0.90, by = 0.10)

# --- 2. Implementação da função quantil contínuo ---

qlindley_cont <- function(u, theta)
{
  arg         <- - (1 + theta) * (1 - u) * exp(-(1 + theta))
  Wval        <- lambertWm1(arg)
  Q_continuo  <- -1 - (1/theta) - (1/theta) * Wval
  return(Q_continuo)
}

# --- 3. Implementação da função quantil discreto (aproximação via parte inteira) ---

qlindley_disc <- function(u, theta)
{
  Q_continuo  <- qlindley_cont(u, theta)
  Q_discreto  <- floor(Q_continuo)
  return(Q_discreto)
}

# --- 4. Cálculo do quantil para cada u ---

q_continuo    <- sapply(u_seq, function(u) qlindley_cont(u, theta))
q_discreto    <- sapply(u_seq, function(u) qlindley_disc(u, theta))

tabela_res    <- data.frame("u" = u_seq,
                            "Quantil_Lindley_Continua" = q_continuo,
                            "Quantil_Lindley_Discreto" = q_discreto)

kable(tabela_res,
      digits = 6,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("Valores de <br> u",
                    "Quantil Lindley <br> Contínua",
                    "Quantil Lindley <br> Discreto"))

Valores de u	Quantil Lindley Contínua	Quantil Lindley Discreto
0.1	0.335240	0
0.2	0.664115	0
0.3	1.003710	1
0.4	1.368048	1
0.5	1.773439	1
0.6	2.243788	2
0.7	2.821228	2
0.8	3.597463	3
0.9	4.859271	4

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.10 (Distribuição Gompertz Discreta - Hegazy et al. 2021). Seja \(Y\) uma variável aleatória contínua com distribuição Gompertz (Gompertz 1825) com parâmetros \(\eta > 0\) (escala) e \(\gamma > 0\) (forma), isto é, \(Y \sim \text{Gompertz}(\eta, \gamma)\). A FDA associada a \(Y\) é dada por:

\[\begin{align}\\ F_Y(y;\eta,\gamma) = 1 - \exp\left\{-\eta \left(e^{\gamma y} - 1\right)\right\} \quad y > 0\\\\ \end{align}\]

A partir dessa FDA, pode-se construir uma variável aleatória discreta \(X\) como versão discretizada de \(Y\), cuja função de probabilidade é definida por:

\[\begin{align}\\ P(X = x) &= F_Y(x + 1;\eta,\gamma) - F_Y(x;\eta,\gamma) \\\\ &= \exp\left\{-\eta \left(e^{\gamma x} - 1\right)\right\} - \exp\left\{-\eta \left(e^{\gamma (x+1)} - 1\right)\right\}, \quad x \in \mathbb{N}.\\\\ \end{align}\]

Nesse caso, a FDA de \(X\) é dada por:

\[\begin{align}\\ F_X(k;\eta,\gamma) = P(X \leqslant k) = P(Y < k+1) = F_Y(k+1;\eta,\gamma), \quad k \in \mathbb{N}\\\\ \end{align}\]

Nosso objetivo é determinar a função quantil de \(X\), denotada por \(Q(u;\eta,\gamma)\), para \(u \in (0,1)\). Conforme estabelecido na Observação 2.2, visto que \(X\) é obtida por discretização da distribuição Gompertz contínua, uma aproximação natural para seu quantil é dada por:

\[\begin{align}\\ Q(u;\eta,\gamma) \approx \lfloor Q_Y(u;\eta,\gamma) \rfloor\\\\ \end{align}\]

em que \(Q_Y(u;\eta,\gamma)\) representa a função quantil da distribuição Gompertz contínua. Esta função pode ser expressa por:

\[\begin{align}\\ Q_Y(u;\eta,\gamma) = \dfrac{1}{\gamma} \ln\left(1 - \dfrac{1}{\eta} \ln(1 - u)\right)\\\\ \end{align}\]

válida para \(u \in \left(0, 1 - e^{-\eta}\right)\), garantindo que o argumento do logaritmo permaneça positivo. Assim, a função quantil aproximada da distribuição Gompertz discreta é descrita por:

\[\begin{align}\\ Q(u;\eta,\gamma) \approx \lfloor Q_Y(u; \eta,\gamma) \rfloor = \left\lfloor \dfrac{1}{\gamma} \ln\left(1 - \dfrac{1}{\eta} \ln(1 - u)\right)\right\rfloor\\\\ \end{align}\]

O Código 2.11 apresenta o cálculo da função quantil da distribuição Gompertz discreta, assumindo os parâmetros \(\eta = 1.5\) e \(\gamma = 0.8\), com base na expressão aproximada de \(Q(u;\eta,\gamma)\).

Código 2.11. Cálculo da função quantil da distribuição Gompertz discreta, assumindo \(\eta = 1.5\) e \(\gamma = 0.8\), com base na expressão \(Q(u;\eta,\gamma)\).

# ------------------------------------------------------------------
# Implementação da Função Quantil da Distribuição Gompertz Discreta
# ------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Sequência de valores para u e parâmetros da Gompertz ---

eta            <- 1.5
gamma          <- 0.8
u_seq          <- seq(0.10, 0.90, by = 0.10)

# --- 2. Implementação da função quantil contínua Gompertz ---

qgompertz_cont <- function(u, eta, gamma)
{
  arg          <- 1 - (1/eta) * log(1 - u)
  
  # Se argumento ≤ 0, retorna NA (fora do domínio do log)
  
  Q_continuo   <- ifelse(arg > 0,
                         (1/gamma) * log(arg),
                         NA_real_)
  
  return(Q_continuo)
}

# --- 3. Implementação da função quantil discreta (aproximação via parte inteira) ---

qgompertz_disc <- function(u, eta, gamma)
{
  Q_continuo   <- qgompertz_cont(u, eta, gamma)
  Q_discreto   <- floor(Q_continuo)
  return(Q_discreto)
}

# --- 4. Cálculo do quantil para cada u ---

q_continuo     <- sapply(u_seq, function(u) qgompertz_cont(u, eta, gamma))
q_discreto     <- sapply(u_seq, function(u) qgompertz_disc(u, eta, gamma))

tabela_res     <- data.frame("u" = u_seq,
                             "Quantil_Gompertz_Contínua" = q_continuo,
                             "Quantil_Gompertz_Discreta" = q_discreto)

kable(tabela_res,
      digits     = 6,
      align      = 'c',
      format     = 'html',
      escape     = FALSE,
      col.names  = c("Valores de <br> u",
                     "Quantil Gompertz <br> Contínua",
                     "Quantil Gompertz <br> Discreta"))

Valores de u	Quantil Gompertz Contínua	Quantil Gompertz Discreta
0.1	0.084854	0
0.2	0.173356	0
0.3	0.266653	0
0.4	0.366350	0
0.5	0.474841	0
0.6	0.595961	0
0.7	0.736571	0
0.8	0.911221	0
0.9	1.162770	1

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

2.2.2. Geração de Números Pseudoaleatórios

---

Agora que se encontram devidamente estabelecidos os procedimentos para a obtenção da função quantil associada tanto a variáveis aleatórias contínuas — por meio de expressões analíticas explícitas ou técnicas numéricas — quanto a variáveis aleatórias discretas, é possível avançar para a aplicação do denominado método da transformação inversa. Esse método consiste na transformação de valores aleatórios gerados a partir de uma distribuição uniforme padrão por meio da função quantil da distribuição-alvo, assegurando, assim, a reprodução da estrutura probabilística desejada. O procedimento correspondente ao caso contínuo é formalizado no Algoritmo 2.1, enquanto a adaptação para o caso discreto será apresentada no Algoritmo 2.2.

Algoritmo 2.1 (Método da Transformação Inversa - Caso Contínuo).

• 1º Passo: Gerar um número pseudoaleatório \(u_i\) a partir de uma variável aleatória \(U \sim \text{Uniforme}(0,1)\), isto é, distribuída uniformemente no intervalo aberto \((0,1)\).

• 2º Passo: Determinar o valor correspondente \(x_i\) como solução da equação \(F(x) = u_i\), onde \(F\) representa a FDA da variável de interesse \(X\).

• 3º Passo: Definir \(x_i = Q(u_i)\), em que \(Q = F^{-1}\) denota a função quantil associada à distribuição de \(X\).

Exemplo 2.11 (Geração de Tempos de Falha - Distribuição Exponencial). Considere o problema de simular tempos de falha de um determinado componente eletrônico cuja vida útil, medida em horas, segue uma distribuição exponencial com taxa de falha \(\lambda = 0.5\). Nesse contexto, a variável aleatória \(X\), que representa o tempo até a falha, possui FDA dada por:

\[\begin{align}\\ F(x;\lambda) = 1 - e^{-\lambda x}, \quad x > 0\\\\ \end{align}\]

Suponha que se deseja gerar uma amostra de tamanho \(n = 1000\) de valores pseudoaleatórios da variável \(X\) utilizando o Algoritmo 2.1. Para isso, é necessário determinar a função quantil \(Q(u) = F^{-1}(u)\), a qual, neste caso, admite a forma analítica:

\[\begin{align}\\ Q(u;\lambda) = F^{-1}(u;\lambda) = -\frac{1}{\lambda} \ln(1 - u)\\\\ \end{align}\]

Dessa forma, cada valor simulado \(x_i\), correspondente ao tempo até a falha, pode ser obtido a partir da função quantil \(Q(u;\lambda)\), segundo a relação:

\[\begin{align}\\ x_i = -\frac{1}{\lambda} \ln(1 - u_i)\\\\ \end{align}\]

em que \(u_i\) representa um número pseudoaleatório gerado de uma distribuição uniforme padrão, \(\text{Uniforme}(0,1)\). Note que, devido à simetria da distribuição uniforme em \((0,1)\), a expressão acima pode, alternativamente, ser simplificada mediante o uso direto de \(u_i\) no lugar de \(1 - u_i\). O Código 2.12 apresenta a implementação computacional, no ambiente R, para a obtenção dos valores \(x_i\) a partir da forma analítica da função quantil, bem como, para fins comparativos, por meio de uma abordagem numérica alternativa, considerando o parâmetro \(\lambda = 0.5\). A Figura 2.2 exibe os histogramas das amostras geradas por ambas as metodologias, sobrepostos à curva teórica da densidade exponencial, permitindo avaliar a aderência dos dados simulados ao modelo teórico.

Código 2.12. Geração de uma amostra de tamanho \(n = 1000\) dos tempos de falha de um componente eletrônico, baseados em uma distribuição exponencial com \(\lambda = 0.5\), utilizando tanto a abordagem analítica quanto a abordagem numérica da função quantil.

# --------------------------------------------------------------------
# Geração de Amostras de Tempos de Falha via Distribuição Exponencial
# --------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)
library(patchwork)

# --- 1. Configurações iniciais ---

set.seed(123)           # Semente (para reprodutibilidade)
lambda          <- 0.5  # Parâmetro da distribuição exponencial
n               <- 1000 # Tamanho de amostra

# --- 2. Funções quantil da distribuição exponencial ---

qexp_analitica  <- function(u, lambda)
{
  -log(1 - u) / lambda
}

qexp_numerica   <- function(u, lambda, lower = 0.0001, upper = 100)
{
  sapply(u, function(p) {uniroot(function(x) pexp(q = x, rate = lambda) - p,
            interval = c(lower, upper))$root})
}

# --- 3. Geração de valores pseudoaleatórios para a distribuição exponencial ---

rexp_analitica  <- function(n, lambda)
{
  u             <- runif(n)
  qexp_analitica(u, lambda)
}

rexp_numerica   <- function(n, lambda)
{
  u             <- runif(n)
  qexp_numerica(u, lambda)
}

amostra_an      <- rexp_analitica(n, lambda)
amostra_nu      <- rexp_numerica(n, lambda)

# --- 4. Visualização gráfica dos valores gerados ---

p1              <- ggplot(data.frame(x = amostra_an), aes(x = x)) +
                     geom_histogram(aes(y = ..density..), bins = 30, fill = "grey80", color = "black") +
                     stat_function(fun = dexp, args = list(rate = lambda),
                                   aes(color = "Curva Teórica"), linewidth = 0.5) +
                     scale_color_manual(name = NULL, values = c("Curva Teórica" = "red")) +
                     theme_minimal() +
                     labs(x = "Tempo de Falha (x)",
                          y = "Densidade") +
                     theme(legend.position = "none", 
                           axis.title.x = element_text(margin = margin(t = 10)),
                           axis.title.y = element_text(margin = margin(r = 10)))

p2              <- ggplot(data.frame(x = amostra_nu), aes(x = x)) +
                     geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "black") +
                     stat_function(fun = dexp, args = list(rate = lambda),
                                   aes(color = "Curva Teórica"), linewidth = 0.5) +
                     scale_color_manual(name = NULL, values = c("Curva Teórica" = "magenta")) +
                     theme_minimal() +
                     labs(x = "Tempo de Falha (x)",
                          y = "Densidade") +
                     theme(legend.position = "none", 
                           axis.title.x = element_text(margin = margin(t = 10)),
                           axis.title.y = element_text(margin = margin(r = 10)))

p1 + p2

Figura 2.2. Histogramas das amostras geradas pelas funções quantil analítica (painel esquerdo) e numérica (painel direito) da distribuição exponencial com parâmetro \(\lambda = 0.5\), acompanhados da curva da densidade teórica.

Para finalizar, com o intuito de avaliar o desempenho computacional dos métodos de geração de valores pseudoaleatórios — analítico e numérico — será adotado um procedimento de benchmarking, consistente na execução sistemática e repetida de testes de desempenho, com o objetivo de mensurar o tempo necessário para a execução de uma função ou algoritmo. No ambiente R, essa análise pode ser realizada por meio do pacote microbenchmark, que fornece estimativas precisas dos tempos de execução, expressos em microssegundos ou nanossegundos, a partir de múltiplas repetições do código em análise. O Código 2.13 ilustra esse procedimento, no qual as funções rexp_analitica e rexp_numerica serão executadas \(B = 10\) vezes, cada uma gerando uma amostra de tamanho \(n = 100\) valores pseudoaleatórios.

Código 2.13. Benchmarking dos processos de geração de \(n = 100\) valores pseudoaleatórios de uma distribuição exponencial com \(\lambda = 0.5\).

# --------------------------------------------------------------------------------------------
# Benchmarking dos Métodos de Geração de Valores Pseudoaleatórios da Distribuição Exponencial
# --------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(microbenchmark)

# --- 1. Avaliação do método analítico ---

set.seed(1212)

microbenchmark(rexp_analitica(100, 2), times = 10, unit = "ms")

Unit: milliseconds
                   expr      min       lq      mean    median       uq      max
 rexp_analitica(100, 2) 0.002075 0.002097 0.1791936 0.0024765 0.002908 1.768289
 neval
    10

# --- 2. Avaliação do método numérico ---

set.seed(1212)

microbenchmark(rexp_numerica(100, 2), times = 10, unit = "ms")

Unit: milliseconds
                  expr      min       lq    mean  median       uq      max
 rexp_numerica(100, 2) 2.621617 2.720401 3.04089 2.79404 2.832428 5.515143
 neval
    10

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.12 (Distribuição Gompertz - Gompertz 1825). Seja \(Y\) uma variável aleatória contínua com distribuição Gompertz com parâmetros \(\eta > 0\) (escala) e \(\gamma > 0\) (forma), isto é, \(Y \sim \text{Gompertz}(\eta, \gamma)\). A FDA associada a \(Y\) é dada por:

\[\begin{align}\\ F_Y(y;\eta,\gamma) = 1 - \exp\left\{-\eta \left(e^{\gamma y} - 1\right)\right\} \quad y > 0\\\\ \end{align}\]

Suponha que se deseja gerar uma amostra de tamanho \(n = 1000\) de valores pseudoaleatórios da variável \(Y\) utilizando o Algoritmo 2.1. Para isso, é necessário determinar a função quantil \(Q(u) = F^{-1}(u)\), a qual, neste caso, admite a forma analítica:

\[\begin{align}\\ Q_Y(u;\eta,\gamma) = \dfrac{1}{\gamma} \ln\left(1 - \dfrac{1}{\eta} \ln(1 - u)\right)\\\\ \end{align}\]

válida para \(u \in \left(0, 1 - e^{-\eta}\right)\), garantindo que o argumento do logaritmo permaneça positivo. Dessa forma, cada valor simulado \(y_i\), correspondente ao tempo até a falha, pode ser obtido a partir da função quantil \(Q_Y(u;\eta,\gamma)\), segundo a relação:

\[\begin{align}\\ y_i = \dfrac{1}{\gamma} \ln\left(1 - \dfrac{1}{\eta} \ln(1 - u_i)\right)\\\\ \end{align}\]

em que \(u_i\) representa um número pseudoaleatório gerado de uma distribuição uniforme padrão, \(\text{Uniforme}(0,1)\). O Código 2.14 apresenta a implementação computacional, no ambiente R, para a obtenção dos valores \(y_i\) a partir da forma analítica da função quantil da distribuição Gompertz, considerando os parâmetros \(\eta = 1.0\) e \(\gamma = 1.5\). A Figura 2.3 exibe o histograma da amostra gerada, sobreposto à curva teórica da densidade Gompertz, permitindo avaliar a aderência dos dados simulados ao modelo teórico.

Código 2.14. Geração de uma amostra de tamanho \(n = 1000\) baseada em uma distribuição Gompertz com \(\eta = 1.0\) e \(\gamma = 1.5\).

# --------------------------------------------------------------
# Geração de Valores Pseudoaleatórios via Distribuição Gompertz
# --------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Configurações iniciais ---

set.seed(123)               # Semente (para reprodutibilidade)
eta           <- 1.0        # Parâmetro de escala da Gompertz
gamma         <- 1.5        # Parâmetro de forma da Gompertz
n             <- 1000       # Tamanho de amostra

# --- 2. Função quantil da distribuição Gompertz ---

qgompertz_an  <- function(u, eta, gamma)
{
  (1 / gamma) * log(1 - (1 / eta) * log(1 - u))
}

# --- 3. Geração de valores pseudoaleatórios para a distribuição Gompertz ---

rgompertz_an  <- function(n, eta, gamma)
{
  u           <- runif(n)
  qgompertz_an(u, eta, gamma)
}

amostra_an    <- rgompertz_an(n, eta, gamma)

# --- 4. Visualização gráfica dos valores gerados ---

# Função densidade da Gompertz

dgompertz <- function(y, eta, gamma)
{
  eta * gamma * exp(gamma * y) * exp(- eta * (exp(gamma * y) - 1))
}

# Gráfico 

ggplot(data.frame(y = amostra_an), aes(x = y)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "grey80", color = "black") +
  stat_function(fun = dgompertz, args = list(eta = eta, gamma = gamma),
                aes(color = "Curva Teórica"), linewidth = 0.5) +
  scale_color_manual(name = NULL, values = c("Curva Teórica" = "red")) +
  theme_minimal() +
  labs(x = expression(y[i]),
       y = "Densidade") +
  theme(legend.position = "none",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 2.3. Histograma da amostra gerada baseada em uma distribuição Gompertz com \(\eta = 1.0\) e \(\gamma = 1.5\), acompanhado da curva da densidade teórica.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Algoritmo 2.2 (Método da Transformação Inversa - Caso Discreto).

• 1º Passo: Gerar um número pseudoaleatório \(u_i\) a partir de uma variável aleatória \(U \sim \mathcal{U}(0,1)\), isto é, distribuída uniformemente no intervalo aberto \((0,1)\).

• 2º Passo: Inicializar as variáveis \(x_i = 0\) e \(S = p_0\), onde \(p_k = P(X = k)\) representa a função de probabilidade associada à variável discreta \(X\).

• 3º Passo: Enquanto \(u_i > S\), incrementar \(x_i \leftarrow x_i + 1\) e atualizar \(S \leftarrow S + p_{x_i}\).

• 4º Passo: Retornar \(x_i\) como a observação pseudoaleatória gerada da variável \(X\).

Exemplo 2.13 (Distribuição Lindley Discreta - Oliveira, Mazucheli, e Achcar 2017). Seja \(X\) uma variável aleatória com distribuição Lindley discreta, cuja formulação decorre do processo de discretização via séries infinitas (Good 1953) aplicado à distribuição Lindley contínua. A função de probabilidade de \(X\) é expressa por:

\[\begin{align}\\ P(X = x;\, \beta) = (1 + x)\, e^{-\beta(x + 2)}\, (e^{\beta} - 1)^2\\\\ \end{align}\]

em que \(x \in \mathbb{N}\) e \(\beta > 0\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(x;\, \beta) = P(X\leqslant x;\beta) =\sum_{k = 0}^{x} (1 + k)\, e^{-\beta(k + 2)}\, (e^{\beta} - 1)^2\\\\ \end{align}\]

Suponha que o objetivo seja gerar uma amostra de tamanho \(n = 10000\) de valores pseudoaleatórios provenientes da variável \(X\), utilizando o método da transformação inversa. No caso de distribuições discretas, tal método fundamenta-se no Algoritmo 2.2, específico para variáveis aleatórias de suporte discreto. O procedimento consiste, para cada observação, em gerar um número pseudoaleatório \(u_i\) proveniente da distribuição uniforme padrão, \(\text{Uniforme}(0,1)\), e determinar o menor inteiro \(x_i\) tal que \(F(x_i) \geqslant u_i\), onde \(F(x)\) corresponde à FDA da distribuição Lindley discreta. O Código 2.15 apresenta a implementação computacional deste procedimento no ambiente R, considerando o parâmetro \(\beta = 0.3\), mediante a criação da função rlindley, a qual retorna valores pseudoaleatórios de \(X\). Por fim, a Figura 2.4 exibe o gráfico da amostra gerada, sobreposto à função de probabilidade teórica, possibilitando a avaliação da adequação da amostra simulada em relação ao modelo teórico.

Código 2.15. Geração de uma amostra de tamanho \(n = 10000\) segundo a distribuição Lindley discreta com parâmetro \(\beta = 0{.}3\), por meio do método da transformação inversa para variáveis aleatórias discretas.

# ----------------------------------------------------------------------------------------------------
# Geração de Valores Pseudoaleatórios via Distribuição Lindley Discreta (Obtida via Séries Infinitas)
# ----------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Configurações iniciais ---

set.seed(1212)                      # Semente para reprodutibilidade
beta            <- 0.3              # Parâmetro da distribuição Lindley discreta
n               <- 10000            # Tamanho da amostra

# --- 2. Função de probabilidade da distribuição Lindley discreta ---

dlindley        <- function(x, beta)
{
  (1 + x) * exp(-beta * (x + 2)) * (exp(beta) - 1)^2
}

# --- 3. Geração de valores pseudoaleatórios para a distribuição Lindley discreta ---

rlindley        <- function(n, beta) 
{
  # Função auxiliar para o cálculo da FDA
  
  prob          <- function(q, beta)
  {
    sum(dlindley(0:q, beta))
  }
  
  U             <- runif(n)         # Amostra Uniforme(0,1)
  X             <- integer(n)       # Vetor para armazenar amostra gerada
  
  for (i in 1:n)
  {
    if(U[i] < prob(q = 0, beta)) {
      X[i]      <- 0
    } else {
      B         <- FALSE
      I         <- 0
      
      while (B == FALSE)
      {
        int     <- c(prob(q = I, beta), prob(q = I + 1, beta))
        
        if((U[i] > int[1]) & (U[i] < int[2])) {
          X[i]  <- I + 1
          B     <- TRUE
        } else {
          I     <- I + 1
        }
      }
    }
  }
  return(X)
}

amostra         <- rlindley(n, beta)

# --- 4. Visualização gráfica dos valores gerados ---

# Tabela de frequências relativas da amostra

probs           <- prop.table(table(amostra))
x               <- as.numeric(names(probs))
y               <- as.numeric(probs)

df              <- data.frame(x = x, 
                              y = y,
                              y_teorico = dlindley(x, beta = beta))

# Gráfico

ggplot(df, aes(x = x)) +
  geom_col(aes(y = y), width = 0.1, color = "black", fill = "grey70") +
  geom_point(aes(y = y_teorico, color = "Probabilidade Teórica"), size = 2) +
  scale_color_manual(name = NULL,
                     values = c("Probabilidade Teórica" = "red")) +
  scale_x_continuous(breaks = round(seq(min(x), max(x), length.out = 5), 0)) +
  scale_y_continuous(breaks = round(seq(min(y), max(y), length.out = 5), 2)) +
  labs(title = "",
       x = expression(x[i]),
       y = "Probabilidade") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 2.4. Gráfico da amostra gerada baseada em uma distribuição Lindley discreta com parâmetro \(\beta = 0{.}3\), acompanhado da função de probabilidade teórica (em vermelho).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.14 (Distribuição Burr Discreta - Krishna e Pundir 2009). Seja \(X\) uma variável aleatória com distribuição Burr discreta, cuja formulação se baseia na discretização da distribuição Burr contínua, conforme apresentado em Krishna e Pundir (2009). A função de probabilidade de \(X\) é expressa por:

\[\begin{align}\\ P(X = x; \theta,\alpha) = \theta^{\ln(1 + x^\alpha)} - \theta^{\ln(1 + (1 + x)^\alpha)}\\\\ \end{align}\]

em que \(x \in \mathbb{N}\), \(0 < \theta < 1\) e \(\alpha > 0\). A FDA associada a \(X\) é dada por:

\[\begin{align}\\ F(x;\theta,\alpha) = P(X\leqslant x; \theta,\alpha) = 1 - \theta^{\ln(1 + x^\alpha)}\\\\ \end{align}\]

Suponha que o objetivo, neste exemplo, seja gerar uma amostra de tamanho \(n = 1000\) de valores pseudoaleatórios provenientes da variável \(X\), utilizando o método da transformação inversa expresso no Algoritmo 2.2. O procedimento consiste, para cada observação, em gerar um número pseudoaleatório \(u_i\) proveniente da distribuição uniforme padrão, \(\text{Uniforme}(0,1)\), e determinar o menor inteiro \(x_i\) tal que \(F(x_i) \geqslant u_i\), onde \(F(x)\) corresponde à FDA da distribuição Burr discreta. O Código 2.16 apresenta a implementação computacional deste procedimento no ambiente R, considerando os parâmetros \(\theta = 0.2\) e \(\alpha = 1.5\), mediante a criação da função rburrdisc, a qual retorna valores pseudoaleatórios de \(X\). Por fim, a Figura 2.5 exibe o gráfico da amostra gerada, sobreposto à função de probabilidade teórica, possibilitando a avaliação da adequação da amostra simulada em relação ao modelo teórico.

Código 2.16. Geração de uma amostra de tamanho \(n = 1000\) segundo a distribuição Burr discreta com parâmetros \(\theta = 0{.}2\) e \(\alpha = 1{.}5\), por meio do método da transformação inversa para variáveis aleatórias discretas.

# -------------------------------------------------------------------
# Geração de Valores Pseudoaleatórios via Distribuição Burr Discreta
# -------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ggplot2)

# --- 1. Configurações iniciais ---

set.seed(123)                       # Semente para reprodutibilidade
theta          <- 0.2               # Parâmetro da distribuição Burr discreta (0 < theta < 1)
alpha          <- 1.5               # Parâmetro da distribuição Burr discreta (alpha > 0)
n              <- 1000              # Tamanho da amostra

# --- 2. Função de probabilidade da distribuição Burr discreta ---

dburrdisc      <- function(x, theta, alpha)
{
  theta^(log(1 + x^alpha)) - theta^(log(1 + (1 + x)^alpha))
}

# --- 3. Geração de valores pseudoaleatórios para a distribuição Burr discreta ---

rburrdisc      <- function(n, theta, alpha)
{
  # Função auxiliar para calcular a FDA acumulada até q

  prob         <- function(q, theta, alpha) 
  {
    sum(dburrdisc(0:q, theta, alpha))
  }
  
  U            <- runif(n)          # Amostra uniforme(0,1)
  X            <- integer(n)        # Vetor para armazenar amostra gerada
  
  for (i in 1:n)
  {
    if (U[i] < prob(0, theta, alpha)) {
      X[i]     <- 0
    } else {
      B        <- FALSE
      I        <- 0
      
      while (!B) 
      {
        int    <- c(prob(I, theta, alpha), prob(I + 1, theta, alpha))
        
        if ((U[i] > int[1]) & (U[i] <= int[2])) {
          X[i] <- I + 1
          B    <- TRUE
        } else {
          I    <- I + 1
        }
      }
    }
  }
  return(X)
}

amostra        <- rburrdisc(n, theta, alpha)

# --- 4. Visualização gráfica dos valores gerados ---

# Tabela de frequências relativas da amostra

probs          <- prop.table(table(amostra))
x              <- as.numeric(names(probs))
y              <- as.numeric(probs)

df             <- data.frame(x = x,
                             y = y,
                             y_teorico = dburrdisc(x, theta, alpha))

# Gráfico

ggplot(df, aes(x = x)) +
  geom_col(aes(y = y), width = 0.1, color = "black", fill = "grey10") +
  geom_point(aes(y = y_teorico, color = "Probabilidade Teórica"), size = 2) +
  scale_color_manual(name = NULL,
                     values = c("Probabilidade Teórica" = "red")) +
  scale_x_continuous(breaks = round(seq(min(x), max(x), length.out = 5), 0)) +
  scale_y_continuous(breaks = round(seq(min(y), max(y), length.out = 5), 4)) +
  labs(title = "",
       x = expression(x[i]),
       y = "Probabilidade") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 2.5. Gráfico da amostra gerada baseada em uma distribuição Burr discreta com parâmetros \(\theta = 0{.}2\) e \(\alpha = 1{.}5\), acompanhado da função de probabilidade teórica (em vermelho).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

2.3. Método da Rejeição

Quando a distribuição de interesse apresenta uma função densidade de probabilidade de difícil manipulação analítica — seja por ausência de uma forma fechada para a função quantil, seja por complexidade estrutural que inviabiliza a inversão da função de distribuição acumulada — torna-se necessário recorrer a estratégias alternativas para a geração de variáveis aleatórias. Dentre essas estratégias, destaca-se o método da rejeição (ou aceitação e rejeição), que consiste em introduzir uma distribuição proposta — da qual seja computacionalmente viável extrair amostras — e aplicar um critério probabilístico baseado em uma variável auxiliar uniformemente distribuída para decidir sobre a aceitação ou rejeição das amostras geradas. Especificamente, supõe-se a existência de uma constante \(c > 0\) tal que, para todo \(x\) no domínio de interesse, a densidade-alvo \(f(x)\) seja majorada por \(c \, g(x)\), onde \(g(x)\) é a densidade da distribuição proposta. O mecanismo de aceitação é então construído de forma a garantir que a distribuição marginal das amostras aceitas coincida com a densidade \(f\), o que é formalmente justificado por meio da propriedade fundamental de densidades, apresentada no Teorema 2.1 a seguir.

Teorema 2.1 (Propriedade Fundamental de Densidades - Devroye 1986; Givens e Hoeting 2012). Seja \(X\) um vetor aleatório com densidade \(f\) definida em \(\mathbb{R}^d\), e seja \(U\) uma variável aleatória com distribuição uniforme no intervalo \([0,1]\). Considere uma constante arbitrária \(k > 0\) tal que \(kf(x)\) seja finito para todo \(x \in \mathbb{R}^d\). Então, o vetor aleatório \((X, k U f(X))\) possui distribuição uniforme sobre o conjunto:

\[\begin{align}\\ A = \{(x, u) \in \mathbb{R}^d \times [0, \infty) : 0 \leqslant u \leqslant k f(x)\}\\\\ \end{align}\]

Ainda, se \((X, U) \in \mathbb{R}^{d+1}\) tem distribuição uniforme sobre \(A\), então a variável aleatória \(X\) tem densidade \(f\) em \(\mathbb{R}^d\).

Demonstração. Para a primeira parte, considere um conjunto de Borel \(\mathcal{B} \subseteq A\), e defina, para cada \(x\), a seção vertical:

\[\begin{align}\\ \mathcal{B}_x = \{u \in \mathbb{R} : (x, u) \in \mathcal{B}\}\\\\ \end{align}\]

Sabendo que \(X\) tem densidade \(f\) e \(U \sim \text{Uniform}(0,1)\), a probabilidade de \((X, k U f(X))\) pertencer a \(\mathcal{B}\) é definida como:

\[\begin{align}\\ P\left( (X, kUf(X)) \in \mathcal{B} \right) &= \int_{\mathbb{R}^d} \int_{\mathcal{B}_x} \frac{1}{k f(x)} \, du \cdot f(x) \, dx = \frac{1}{k} \int_{\mathcal{B}} du \, dx,\\\\ \end{align}\]

o que mostra que \((X, kUf(X))\) possui densidade uniforme em \(A\), cuja “área” (no sentido de volume \(d+1\) dimensional) é \(k\). Assim, a primeira parte está provada. Para a segunda parte, suponha que \((X, U)\) é uniformemente distribuído em \(A\). Então, dado um conjunto de Borel \(\mathcal{B} \subset \mathbb{R}^d\), defina:

\[\begin{align}\\ \mathcal{B}_1 = \{(x, u) \in A : x \in \mathcal{B} \} = \{(x,u) : x \in \mathcal{B}, 0 \leqslant u \leqslant kf(x)\}\\\\ \end{align}\]

Logo, a probabilidade de \(X \in \mathcal{B}\) é descrita por:

\[\begin{align}\\ P(X \in \mathcal{B}) &= P\big((X,U) \in \mathcal{B}_1\big) = \frac{\displaystyle\int_{\mathcal{B}_1} du \, dx}{\displaystyle\int_A du \, dx} = \frac{1}{k} \int_{\mathcal{B}} kf(x) \, dx = \int_{\mathcal{B}} f(x) \, dx,\\\\ \end{align}\]

o que implica que a variável aleatória \(X\) tem de fato densidade \(f\) em \(\mathbb{R}^d\), como queríamos demonstrar.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Algoritmo 2.3 (Método da Rejeição - Givens e Hoeting 2012). Seja \(f\) a densidade-alvo e \(g\) uma densidade auxiliar da qual seja possível amostrar diretamente. Suponha que exista uma constante \(c > 0\) tal que, para todo \(x\) no domínio de \(f\), tenha-se

\[\begin{align}\\ f(x) \leqslant e(x)\\\\ \end{align}\]

onde \(e(\cdot)\) é uma função, chamada de função envelope, definida por \(e(x) = cg(x)\). Então, o procedimento para gerar um valor pseudoaleatório de uma variável aleatória com densidade proporcional a \(f\) consiste nos seguintes passos:

• 1º Passo: Gerar um valor pseudoaleatório \(y\) de uma variável aleatória \(G\) com função densidade de probabilidade \(g(y)\).

• 2º Passo: Gerar um valor pseudoaleatório \(u\) de uma variável aleatória \(U \sim\) Uniform\((0,1)\), independente de \(G\).

• 3º Passo: Se \(u < \dfrac{f(y)}{c \, g(y)}\), então \(x = y\) é aceito como valor pseudoaleatório da distribuição alvo com densidade \(f\); caso contrário, rejeita-se o valor \(y\).

A constante \(c\) empregada no Algoritmo 2.3 é denominada constante de rejeição e decorre da escolha da função envelope, \(e(\cdot)\), que majoriza a densidade-alvo \(f\), ou seja, satisfaz \(f(x) \leqslant e(x)\) para todo \(x\) no domínio de interesse. Formalmente, define-se:

\[\begin{align}\\ c = \sup_x \frac{f(x)}{g(x)}\\\\ \end{align}\]

Na prática, é comum que a densidade auxiliar adotada no método da rejeição seja escolhida dentro de uma família paramétrica de densidades, denotada por \(\{ g_\theta : \theta \in \Theta \subseteq \mathbb{R}^k \}\), em que \(\theta\) representa um vetor de parâmetros de dimensão \(k\), pertencente ao espaço paramétrico \(\Theta\). Para essa família, define-se a constante de rejeição por:

\[\begin{align}\\ c_\theta = \sup_x \frac{f(x)}{g_\theta(x)}\\\\ \end{align}\]

A escolha ótima do parâmetro \(\theta\) é aquela que minimiza a constante de rejeição \(c_\theta\), isto é, \(\theta^\ast = \arg\min_{\theta \ \in \ \Theta}c_\theta\). Essa minimização é essecial, pois maximiza a eficiência do algoritmo, reduzindo a taxa de rejeição e, consequentemente, o custo computacional do procedimento. A Figura 2.6 ilustra graficamente essa construção: a densidade-alvo \(f(x)\) é coberta pela função envelope \(e(x) = c_\theta\, g_\theta(x)\), cuja forma depende do valor de \(\theta\). A área entre as curvas \(e(x)\) e \(f(x)\) representa o desperdício computacional decorrente das amostras rejeitadas. Assim, a escolha adequada de \(\theta\) — isto é, aquela que minimiza essa área — é essencial para obter um algoritmo eficiente.

Figura 2.6. Ilustração do método da rejeição para uma distribuição-alvo \(f\) utilizando uma função envelope \(e\).

Fonte: Rejection Method | Givens, G. H. & Hoeting, J. A. (2012). Computational Statistics. 2\(^{\text{nd}}\) Edition. John Wiley & Sons, Inc.

Exemplo 2.15 (Distribuição Normal Padrão). Considere a densidade da distribuição normal padrão dada por:

\[\begin{align}\\ f(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}\\\\ \end{align}\]

Nosso objetivo é gerar uma amostra pseudoaleatória de tamanho \(n = 10000\) dessa distribuição por meio do método da rejeição. Para tanto, é necessário selecionar uma densidade auxiliar \(g\) que satisfaça a condição de cobertura \(f(x) \leqslant cg(x)\), para algum valor finito de \(c > 0\). Uma escolha clássica na literatura é a distribuição de Cauchy com parâmetro de escala \(\theta > 0\), cuja densidade é dada por:

\[\begin{align}\\ g_\theta(x) = \frac{\theta}{\pi(\theta^2 + x^2)}, \quad x \in \mathbb{R}\\\\ \end{align}\]

A determinação da constante \(c\) que garante a cobertura de \(f\) por \(cg_\theta\) baseia-se na razão entre as funções \(f(x)\) e \(g_\theta(x)\), dada por:

\[\begin{align}\\ \frac{f(x)}{g_\theta(x)} = \frac{(1/\sqrt{2\pi}) e^{-x^2/2}}{(\theta/\pi)(\theta + x^2)^{-1}}\\\\ \end{align}\]

A análise gráfica dessa razão, para \(\theta = 1\), revela o comportamento da função de cobertura e orienta a escolha da constante de rejeição. O Código 2.17 apresenta a implementação, em ambiente R, da comparação entre as densidades da distribuição normal padrão e da distribuição Cauchy padrão, bem como a análise do comportamento da razão \(f/g_\theta\), considerando \(\theta = 1\). A Figura 2.7 ilustra graficamente essa comparação, evidenciando as diferenças entre as densidades e o perfil da razão entre elas.

Código 2.17. Implementação da comparação entre as densidades normal (\(f\)) e Cauchy (\(g_\theta\)), e, também, do comportamento da razão \(f/g_\theta\), com \(\theta = 1\).

# -------------------------------------------------------------
# Comparação: Distribuição Normal Padrão x Distribuição Cauchy
# -------------------------------------------------------------

# --- 0. Carregar os pacotes necessários ---

library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(tibble)

# --- 1. Gerar dados genéricos ---

x_vals    <- seq(-4, 4, length.out = 1000)

df        <- tibble(x       = x_vals,
                    Normal  = dnorm(x),
                    Cauchy  = dcauchy(x),
                    Ratio   = dnorm(x) / dcauchy(x))

# --- 2. Preparar dados para o gráfico de densidades ---

df_long   <- df %>%
              pivot_longer(cols       = c(Normal, Cauchy),
                           names_to   = "Distribuição",
                           values_to  = "Densidade")

# --- 3. Construir gráfico das densidades Normal e Cauchy ---

g1        <- ggplot(df_long, aes(x = x, y = Densidade, color = Distribuição, linetype = Distribuição)) +
              geom_line(linewidth = 0.5) +
              theme_minimal() +
              labs(title = '',
                   y     = "Densidade") +
              scale_color_manual(values = c("black", "gray50")) +
              theme(legend.position = "left",
                    axis.title.x = element_text(margin = margin(t = 10)),
                    axis.title.y = element_text(margin = margin(r = 10)))

# --- 4. Construir gráfico da razão f(x) / g_theta(x) ---

g2        <- ggplot(df, aes(x = x, y = Ratio)) +
              geom_line(color = "darkblue", linewidth = 0.5) +
              geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
              theme_minimal() +
              labs(title = '',
                   y     = expression(f(x) / g[theta](x)),
                   x     = "x") +
              theme(axis.title.x = element_text(margin = margin(t = 10)),
                    axis.title.y = element_text(margin = margin(r = 10)))

# --- 5. Visualização conjunta dos gráficos ---

g1 + g2

Figura 2.7. Ilustração da comparação entre as densidades normal e Cauchy (painel esquerdo), e, também, do comportamento da razão \(f/g_\theta\) (painel direito), com \(\theta = 1\).

Para obter a constante ótima de rejeição \(c_\theta\), note que \(f/g_\theta\) é máxima quando \(\ln(f/g_\theta)\) é máxima. Derivando essa função e igualando a zero, tem-se:

\[\begin{align}\\ \frac{d}{dx} \ln\left( \frac{f(x)}{g_\theta(x)} \right) = -x + \frac{2x}{\theta^2 + x^2} = 0\\\\ \end{align}\]

Resolvendo essa equação, obtém-se os pontos críticos \(x = 0\) e \(x = \pm \sqrt{2 - \theta^2}\), desde que \(\theta^2 \leqslant 2\). Para esses pontos, a razão \(f/g_\theta\) assume os seguintes valores:

\[\begin{align}\\ \frac{f(0)}{g_\theta(0)} = \theta \sqrt{\frac{\pi}{2}}, \quad \frac{f(\pm\sqrt{2 - \theta^2})}{g_\theta(\pm\sqrt{2 - \theta^2})} = \frac{\sqrt{2\pi}}{e\,\theta} e^{\theta^2/2}\\\\ \end{align}\]

Logo, a constante ótima de rejeição é dada por:

\[\begin{align}\\ c_\theta = \begin{cases} \dfrac{\sqrt{2\pi}}{e\,\theta} e^{\theta^2/2}, & \theta < \sqrt{2}, \\\\ \theta \sqrt{\pi/2}, & \theta \geqslant \sqrt{2} \end{cases}\\\\ \end{align}\]

Assim, o valor mínimo de \(c_\theta\) é obtido para \(\theta = 1\), resultando em \(c = \sqrt{2\pi/e}\). O Código 2.18 apresenta a implementação, em ambiente R, da função envelope \(e(x) = cg(x)\) cobrindo a densidade-alvo \(f(x)\), enquanto a Figura 2.8 ilustra graficamente essa cobertura da função envelope.

Código 2.18. Implementação da função envelope \(e(x) = cg(x)\) cobrindo a densidade-alvo \(f(x)\).

# --------------------------------------------------
# Função Envelope Para a Distribuição Normal Padrão
# --------------------------------------------------

# --- 0. Carregar os pacotes necessários ---

library(ggplot2)
library(dplyr)
library(tidyr)

# --- 1. Definir constante de normalização ---

c         <- sqrt(2 * pi / exp(1))

# --- 2. Calcular a função envelope e a densidade normal ---

df2       <- tibble(x         = x_vals,
                    Normal    = dnorm(x_vals),
                    Envelope  = c * dcauchy(x_vals))

# --- 3. Reestruturar base em formato longo ---

df2_long  <- df2 %>%
                   pivot_longer(cols      = c(Normal, Envelope), 
                                names_to  = "Curva", 
                                values_to = "Densidade")

# --- 4. Construir gráfico com curva normal e envelope ---

ggplot(df2_long, aes(x         = x, 
                     y         = Densidade, 
                     color     = Curva, 
                     linetype  = Curva)) +
       geom_line(linewidth     = 1) +
       theme_minimal() +
       labs(title              = '',
            y                  = "Densidade") +
       scale_color_manual(values = c("black", "gray50")) +
       theme(legend.position  = "right", 
             axis.title.x = element_text(margin = margin(t = 10)),
             axis.title.y = element_text(margin = margin(r = 10)))

Figura 2.8. Ilustração da função envelope \(e(x) = cg(x)\) cobrindo a densidade-alvo \(f(x)\).

Com a constante \(c\) adequadamente determinada, a geração da amostra pode ser realizada por meio do seguinte procedimento baseado no Algoritmo 2.3:

• 1º Passo: Gerar \(y_i \sim\) Cauchy\((0, 1)\), uma vez que \(\theta = 1\);
• 2º Passo: Gerar \(u_i \sim U(0, 1)\);
• 3º Passo: Calcular \(c = \sqrt{2\pi/e}\);
• 4º Passo: Aceitar \(y_i\) como \(x_i\) se \(u_i \leqslant f(y_i) / [c\,g(y_i)]\), caso contrário, repetir o processo.

O Código 2.19 apresenta, em ambiente R, a implementação dos procedimentos descritos anteriormente para a geração de uma amostra de tamanho \(n = 10000\) da distribuição normal padrão, empregando o método da rejeição. A Figura 2.9 exibe o histograma e a FDA da amostra gerada, sobrepostos, respectivamente, à função densidade de probabilidade teórica e à FDA teórica, permitindo avaliar a adequação da amostra simulada ao modelo teórico.

Código 2.19. Geração de uma amostra de tamanho \(n = 10000\) da distribuição normal padrão, utilizando o método da rejeição e a distribuição de Cauchy como proposta.

# --------------------------------------------------------------------------------------
# Método da Rejeição: Geração de Valores Pseudoaleatórios da Distribuição Normal Padrão
# --------------------------------------------------------------------------------------

# --- 0. Carregar os pacotes necessários ---

library(ggplot2)

# --- 1. Definir valores iniciais ---

set.seed(1212)
n         <- 10000
c         <- sqrt(2 * pi / exp(1))
x         <- numeric(n)
i         <- 1

# --- 2. Gerar amostra utilizando o método da rejeição ---

while (i < n)
{
  y       <- rcauchy(1)
  u       <- runif(1)
  
  if (u <= dnorm(y) / (c * dcauchy(y))) 
  {
    x[i]  <- y
    i     <- i + 1
  }
}

# --- 3. Construir histograma comparado à densidade normal padrão ---

p1        <- ggplot(data.frame(x = x), aes(x)) +
                    geom_histogram(aes(y = ..density.., fill = "Amostra Gerada"), 
                                   bins = 30,
                                   fill = "gray80", color = "black") +
                    stat_function(aes(color = "Curva Normal"), 
                                  fun = dnorm, 
                                  linewidth = 1) +
                    theme_minimal() +
                    labs(title = "",
                         y     = "Densidade") +
                    scale_color_manual(name = NULL, 
                                       values = c("Curva Normal" = "darkred")) +
                    theme(legend.position = "none",
                          axis.title.x    = element_text(margin = margin(t = 10)),
                          axis.title.y    = element_text(margin = margin(r = 10)))

# --- 4. Construir gráfico EFDA comparado à FDA normal padrão ---

p2        <- ggplot(data.frame(x = x), aes(x)) +
                    stat_ecdf(aes(color = "EFDA - Amostra Gerada"), 
                              geom = "step", 
                              linewidth = 0.8) +
                    stat_function(aes(color = "FDA - Normal Padrão"), 
                                  fun = pnorm, 
                                  linetype = "dashed", 
                                  linewidth = 0.8) +
                    theme_minimal() +
                    labs(title = "",
                         y     = "Probabilidade Acumulada") +
                    scale_color_manual(name = NULL,
                                       values = c("EFDA - Amostra Gerada" = "black",
                                                  "FDA - Normal Padrão" = "darkred")) +
                    theme(legend.position  = "right", 
                          axis.title.x     = element_text(margin = margin(t = 10)),
                          axis.title.y     = element_text(margin = margin(r = 10)))

# --- 5. Visualização conjunta dos gráficos ---

p1 + p2

Figura 2.9. Histograma (painel esquerdo) e FDA (painel direito) empírica da amostra gerada, sobrepostos à densidade e à FDA teóricas da distribuição normal padrão.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.16 (Distribuição Gama - Thom 1958). Considere a densidade da distribuição gama com parâmetros \(\alpha > 0\) e \(\beta > 0\), dada por:

\[\begin{align}\\ f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x}, \quad x > 0\\\\ \end{align}\]

Nosso objetivo é gerar uma amostra pseudoaleatória de tamanho \(n = 10000\) dessa distribuição por meio do método da rejeição. Para tanto, é necessário selecionar uma densidade auxiliar \(g\) que satisfaça a condição de cobertura \(f(x) \leqslant c g(x)\), para algum valor finito de \(c > 0\). Uma escolha clássica na literatura é a distribuição exponencial com parâmetro de taxa \(\lambda > 0\), cuja densidade é dada por:

\[\begin{align}\\ g_{\lambda}(x) = \lambda e^{-\lambda x}, \quad x > 0\\\\ \end{align}\]

A determinação da constante \(c\) que garante a cobertura de \(f\) por \(c g_{\lambda}\) baseia-se na razão entre as funções \(f(x)\) e \(g_{\lambda}(x)\), dada por:

\[\begin{align}\\ \frac{f(x)}{g_{\lambda}(x)} = \frac{\beta^\alpha}{\lambda \Gamma(\alpha)} x^{\alpha - 1} e^{-(\beta - \lambda) x}\\\\ \end{align}\]

A análise gráfica dessa razão, para \(\alpha = 3\), \(\beta = 2\) e \(\lambda = 0.5\), revela o comportamento da função de cobertura e orienta a escolha da constante de rejeição. O Código 2.20 apresenta a implementação, em ambiente R, da comparação entre as densidades gama e exponencial, bem como a análise do comportamento da razão \(f/g_\lambda\). A Figura 2.10 ilustra graficamente essa comparação, evidenciando as diferenças entre as densidades e o perfil da razão entre elas.

Código 2.20. Implementação da comparação entre as densidades gama e exponencial, e do comportamento da razão \(f/g_\lambda\), com \(\alpha = 3\), \(\beta = 2\) e \(\lambda = 0.5\).

# ---------------------------------------------------------
# Comparação: Distribuição Gama x Distribuição Exponencial
# ---------------------------------------------------------

# --- 0. Carregar os pacotes necessários ---

library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(tibble)

# --- 1. Gerar dados genéricos ---

alpha     <- 3
beta      <- 2
lambda    <- 0.5

x_vals    <- seq(0, 7, length.out = 1000)

df        <- tibble(x = x_vals,
                    Gama = dgamma(x, shape = alpha, rate = beta),
                    Exponencial = dexp(x, rate = lambda),
                    Ratio = dgamma(x, shape = alpha, rate = beta) / dexp(x, rate = lambda))

# --- 2. Preparar dados para o gráfico de densidades ---

df_long   <- df %>%
                  pivot_longer(cols = c(Gama, Exponencial), 
                               names_to = "Distribuicao", 
                               values_to = "Densidade")

# --- 3. Construir gráfico das densidades gama e exponencial ---

g1        <- ggplot(df_long, aes(x = x,
                                 y = Densidade, 
                                 color = Distribuicao, 
                                 linetype = Distribuicao)) +
                    geom_line(linewidth = 1) +
                    theme_minimal() +
                    labs(title = '',
                         y = "Densidade") +
                    scale_color_manual(values = c("black", "gray50")) +
                    theme(legend.position = "left",
                          axis.title.x = element_text(margin = margin(t = 10)),
                          axis.title.y = element_text(margin = margin(r = 10)))

# --- 4. Construir gráfico da razão f(x) / g_lambda(x) ---

g2        <- ggplot(df, aes(x = x, y = Ratio)) +
                    geom_line(color = "darkblue", linewidth = 1) +
                    theme_minimal() +
                    labs(title = '',
                         y = expression(f(x)/g[lambda](x)), x = "x") +
                    theme(axis.title.x = element_text(margin = margin(t = 10)),
                          axis.title.y = element_text(margin = margin(r = 10)))

# --- 5. Visualização conjunta dos gráficos ---

g1 + g2

Figura 2.10. Ilustração da comparação entre as densidades gama e exponencial (painel esquerdo), e, também, do comportamento da razão \(f/g_\lambda\) (painel direito).

Para obter a constante ótima de rejeição \(c\), note que \(f/g_\lambda\) é máxima quando \(\ln(f/g_\lambda)\) é máxima. Derivando essa função logarítmica e igualando a zero, tem-se:

\[\begin{align}\\ \frac{d}{dx} \ln\left( \frac{f(x)}{g_\lambda(x)} \right) = \frac{\alpha - 1}{x} - (\beta - \lambda) = 0\\\\ \end{align}\]

Resolvendo a equação acima, obtém-se o ponto crítico:

\[\begin{align}\\ x^* = \frac{\alpha - 1}{\beta - \lambda}, \quad \alpha > 1, \ \lambda < \beta\\\\ \end{align}\]

Substituindo esse valor de \(x^*\) na razão \(f(x) / g_\lambda(x)\), obtém-se a constante ótima de rejeição \(c\):

\[\begin{align}\\ c &=\frac{\beta^\alpha}{\lambda \Gamma(\alpha)} \left( \frac{\alpha - 1}{\beta - \lambda} \right)^{\alpha - 1} e^{-(\alpha - 1)}\\\\ \end{align}\]

Para minimizar \(c\), pode-se fixar \(\alpha\) e \(\beta\) (do alvo) e variar \(\lambda\) (da proposta), e, como \(x^* = (\alpha - 1)/(\beta - \lambda)\), a escolha ótima de \(\lambda\) está ligada ao valor que equilibra o decaimento exponencial da proposta e o pico da densidade-alvo. O Código 2.21 apresenta a implementação, em ambiente R, da função envelope \(e(x) = cg(x)\) cobrindo a densidade-alvo \(f(x)\) assumindo os valores ótimos \(\alpha = 3\), \(\beta = 2\) e \(\lambda = 0.5\), enquanto a Figura 2.11 ilustra graficamente essa cobertura da função envelope.

Código 2.21. Implementação da função envelope \(e(x) = cg(x)\) cobrindo a densidade-alvo \(f(x)\) para \(\alpha = 3\), \(\beta = 2\) e \(\lambda = 0.5\).

# -----------------------------------------
# Função Envelope Para a Distribuição Gama
# -----------------------------------------

# --- 0. Carregar os pacotes necessários ---

library(ggplot2)
library(dplyr)
library(tidyr)

# --- 1. Definir constante de normalização ---

alpha     <- 3
beta      <- 2
lambda    <- 0.5
x_vals    <- seq(0, 7, length.out = 1000)

c_1       <- (beta^alpha * (alpha - 1)^(alpha - 1) * exp(-(alpha - 1)))
c_2       <- (lambda * gamma(alpha) * (beta - lambda)^(alpha - 1))
c         <-  c_1 / c_2 

# --- 2. Calcular a função envelope e a densidade gama ---

df2       <- tibble(x = x_vals,
                    Gama = dgamma(x_vals, shape = alpha, rate = beta),
                    Envelope = c * dexp(x_vals, rate = lambda))

# --- 3. Reestruturar base em formato longo ---

df2_long  <- df2 %>%
                   pivot_longer(cols = c(Gama, Envelope),
                                names_to = "Curva",
                                values_to = "Densidade")

# --- 4. Construir gráfico com curva gama e envelope ---

ggplot(df2_long, aes(x = x, y = Densidade, color = Curva, linetype = Curva)) +
       geom_line(linewidth = 0.5) +
       theme_minimal() +
       labs(title = '',        
            y = "Densidade") +
       scale_color_manual(values = c("black", "gray50")) +
       theme(legend.position  = "right", 
             axis.title.x = element_text(margin = margin(t = 10)),
             axis.title.y = element_text(margin = margin(r = 10)))

Figura 2.11. Ilustração da função envelope \(e(x) = cg(x)\) cobrindo a densidade-alvo \(f(x)\).

Com a constante \(c\) adequadamente determinada, a geração da amostra pode ser realizada por meio do seguinte procedimento baseado no Algoritmo 2.3 adaptado para o caso da distribuição gama com proposta exponencial:

• 1º Passo: Gerar \(y_i \sim\) Exponencial\((\lambda)\), cuja densidade é \(g(y) = \lambda e^{-\lambda y}\);
• 2º Passo: Gerar \(u_i \sim U(0,1)\);
• 3º Passo: Calcular \(c = \dfrac{\beta^\alpha}{\lambda \Gamma(\alpha)} \left( \dfrac{\alpha - 1}{\beta - \lambda} \right)^{\alpha - 1} e^{-(\alpha - 1)}\);
• 4º Passo: Aceitar \(y_i\) como \(x_i\) se \(u_i \leqslant f(y_i) / [c\,g(y_i)]\), caso contrário, repetir o processo.

O Código 2.22 apresenta, em ambiente R, a implementação dos procedimentos descritos anteriormente para a geração de uma amostra de tamanho \(n = 10000\) da distribuição gama, empregando o método da rejeição. A Figura 2.12 exibe o histograma e a FDA da amostra gerada, sobrepostos, respectivamente, à função densidade de probabilidade teórica e à FDA teórica, permitindo avaliar a adequação da amostra simulada ao modelo teórico.

Código 2.22. Geração de uma amostra de tamanho \(n = 10000\) da distribuição gama, utilizando o método da rejeição e a distribuição de exponencial como proposta.

# -----------------------------------------------------------------------------
# Método da Rejeição: Geração de Valores Pseudoaleatórios da Distribuição Gama
# -----------------------------------------------------------------------------

# --- 0. Carregar os pacotes necessários ---

library(ggplot2)

# --- 1. Definir valores iniciais ---

set.seed(1212)
n         <- 10000
x         <- numeric(n)
i         <- 1
alpha     <- 3
beta      <- 2
lambda    <- 0.5

c_1       <- (beta^alpha * (alpha - 1)^(alpha - 1) * exp(-(alpha - 1)))
c_2       <- (lambda * gamma(alpha) * (beta - lambda)^(alpha - 1))
c         <-  c_1 / c_2

# --- 2. Gerar amostra utilizando o método da rejeição ---

while (i < n)
{
  y       <- rexp(1, rate = lambda)
  u       <- runif(1)
  
  if (u <= dgamma(y, shape = alpha, rate = beta) / (c * dexp(y, rate = lambda))) 
  {
    x[i]  <- y
    i     <- i + 1
  }
}

# --- 3. Construir histograma comparado à densidade gama ---

p1        <- ggplot(data.frame(x = x), aes(x)) +
                    geom_histogram(aes(y = ..density.., fill = "Amostra Gerada"), bins = 30,
                                   fill = "gray80", color = "black") +
                    stat_function(aes(color = "Curva Gama"), 
                                  fun = function(z) dgamma(z, shape = alpha, rate = beta), 
                                  linewidth = 0.5) +
                    theme_minimal() +
                    labs(title = "",
                         y = "Densidade") +
                    scale_color_manual(name = NULL, values = c("Curva Gama" = "darkred")) +
                    theme(legend.position = "none",
                          axis.title.x    = element_text(margin = margin(t = 10)),
                          axis.title.y    = element_text(margin = margin(r = 10)))

# --- 4. Construir gráfico EFDA comparado à FDA gama ---

p2        <- ggplot(data.frame(x = x), aes(x)) +
                    stat_ecdf(aes(color = "EFDA - Amostra Gerada"), 
                              geom = "step", linewidth = 0.8) +
                    stat_function(aes(color = "FDA - Gama"), 
                                  fun = function(z) pgamma(z, shape = alpha, rate = beta), 
                                  linetype = "dashed", 
                                  linewidth = 0.5) +
                    theme_minimal() +
                    labs(title = "",
                         y = "Probabilidade Acumulada") +
                    scale_color_manual(name = NULL,
                                       values = c("EFDA - Amostra Gerada" = "black",
                                                  "FDA - Gama" = "darkred")) +
                    theme(legend.position = "right",
                          axis.title.x    = element_text(margin = margin(t = 10)),
                          axis.title.y    = element_text(margin = margin(r = 10)))

# --- 5. Visualização conjunta dos gráficos ---

p1 + p2

Figura 2.12. Histograma (painel esquerdo) e FDA (painel direito) empírica da amostra gerada, sobrepostos à densidade e à FDA teóricas da distribuição gama.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Observação 2.3. Visto que o método da rejeição consiste em aceitar ou rejeitar amostras geradas a partir de uma distribuição auxiliar com base na razão entre a densidade-alvo e uma função envelope. Nesse contexto, duas medidas de desempenho são de particular interesse: a taxa de aceitação observada e a taxa de aceitação esperada. A primeira é obtida empiricamente, sendo definida como a razão entre o número de valores efetivamente aceitos pelo algoritmo e o número total de tentativas realizadas. Por outro lado, a taxa de aceitação esperada é uma quantidade teórica associada à constante de rejeição \(c\), definida na condição de cobertura \(f(x) \leqslant cg(x)\) para todo \(x\). Essa taxa esperada corresponde a \(1/c\) e representa o valor médio da proporção de amostras aceitas, assumindo que a razão \(f(x)/g(x)\) é corretamente limitada superiormente por \(c\).

Exemplo 2.17. Considere novamente a distribuição gama com parâmetros \(\alpha > 0\) e \(\beta > 0\), bem como a densidade auxiliar \(g\) definida pela distribuição exponencial com parâmetro de taxa \(\lambda > 0\), descritas no Exemplo 2.16. Suponha que o objetivo aqui seja determinar a taxa de aceitação observada e a taxa de aceitação esperada no contexto do método da rejeição. Para tal, será gerada uma amostra pseudoaleatória de tamanho \(n = 10000\), empregando-se a constante ótima de rejeição dada por:

\[\begin{align}\\ c &=\frac{\beta^\alpha}{\lambda \Gamma(\alpha)} \left( \frac{\alpha - 1}{\beta - \lambda} \right)^{\alpha - 1} e^{-(\alpha - 1)}\\\\ \end{align}\]

O Código 2.23 apresenta, em ambiente R, a implementação da geração da amostra por meio do método da rejeição, bem como o cálculo das taxas de aceitação observada e esperada para diferentes combinações de valores paramétricos.

Código 2.23. Cálculo das taxas de aceitação observada e esperada para uma amostra de tamanho \(n = 10000\) gerada da distribuição gama, utilizando o método da rejeição e adotando a distribuição exponencial como densidade proposta, para diferentes combinações de valores paramétricos.

# ----------------------------------------------------------------------------------------------------
# Taxas de Aceitação Observada e Esperada da Distribuição Gamma via Rejeição com Proposta Exponencial
# ----------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Função de rejeição para a distribuição Gamma ---

rejection.gamma <- function(n, alpha, beta, lambda)
{
  x             <- numeric(n)
  i             <- 1
  l             <- 1
  
  c_1           <- (beta^alpha * (alpha - 1)^(alpha - 1) * exp(-(alpha - 1)))
  c_2           <- (lambda * gamma(alpha) * (beta - lambda)^(alpha - 1))
  c             <- c_1 / c_2
  
  while (i < n)
  {
    y           <- rexp(1, rate = lambda)
    u           <- runif(1)
    
    if (u <= dgamma(y, shape = alpha, rate = beta) / (c * dexp(y, rate = lambda)))
    {
      x[i]      <- y
      i         <- i + 1
    }
    l           <- l + 1
  }
  
  tx_obs        <- round(i / l * 100, 2)
  tx_esp        <- round(1 / c * 100, 2)
  
  return(c(tx_obs, tx_esp))
}

# --- 2. Combinações viáveis: lambda < beta ---

parametros      <- expand.grid(alpha  = c(2, 3, 5),
                               beta   = c(1, 2),
                               lambda = c(0.5, 1))

parametros      <- subset(parametros, lambda < beta)

# --- 3. Construir das funções envelope para cada combinação válida ---

x_vals          <- seq(0.001, 10, length.out = 1000)

graficos        <- list()

for (k in 1:nrow(parametros)) 
{
  a             <- parametros$alpha[k]
  b             <- parametros$beta[k]
  l             <- parametros$lambda[k]
  
  c_1           <- (b^a * (a - 1)^(a - 1) * exp(-(a - 1)))
  c_2           <- (l * gamma(a) * (b - l)^(a - 1))
  c             <- c_1 / c_2
  
  df_env        <- tibble(x = x_vals,
                          Gama = dgamma(x_vals, shape = a, rate = b),
                          Envelope = c * dexp(x_vals, rate = l)) %>%
                   pivot_longer(cols = c(Gama, Envelope),
                                names_to = "Curva",
                                values_to = "Densidade")
  
  g             <- ggplot(df_env, aes(x = x, y = Densidade, color = Curva, linetype = Curva)) +
                   geom_line(linewidth = 0.5) +
                   theme_minimal() +
                   labs(title = bquote(paste(alpha == .(a), ", ", beta == .(b), ", ", lambda == .(l))),
                        y = "Densidade", x = "x") +
                   scale_color_manual(values = c("black", "gray50")) +
                   theme(legend.position = "none",
                         plot.title = element_text(size = 9),
                         axis.title.x = element_text(size = 8),
                         axis.title.y = element_text(size = 8))
  
  graficos[[k]] <- g
}

wrap_plots(graficos, ncol = 3)

# --- 3. Cálculo das taxas para cada configuração de parâmetros ---

set.seed(1212)
n               <- 100
resultados      <- data.frame()

for (k in 1:nrow(parametros)) 
{
  a             <- parametros$alpha[k]
  b             <- parametros$beta[k]
  l             <- parametros$lambda[k]
  
  taxas         <- rejection.gamma(n, a, b, l)
  
  resultados    <- rbind(resultados,
                         data.frame(alpha = a,
                                    beta = b,
                                    lambda = l,
                                    Taxa_Observada = taxas[1],
                                    Taxa_Esperada  = taxas[2]))
}

# --- 4. Exibir resultados ---

kable(resultados,
      digits = 2,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("&alpha;", "&beta;", "&lambda;",
                    "Taxa Observada (%)", 
                    "Taxa Esperada (%)"))

α	β	λ	Taxa Observada (%)	Taxa Esperada (%)
2	1	0.5	70.92	67.96
3	1	0.5	51.55	46.18
5	1	0.5	18.28	16.00
2	2	0.5	54.95	50.97
3	2	0.5	48.31	51.95
5	2	0.5	41.15	40.49
2	2	1.0	64.94	67.96
3	2	1.0	43.10	46.18
5	2	1.0	14.99	16.00

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

2.4. Testes de Verificação

Para concluir este capítulo, será abordado o estudo dos testes de verificação, com ênfase nos testes de bondade de ajuste para a distribuição alvo. Esses testes são essenciais para assegurar que os algoritmos empregados na geração de números pseudoaleatórios — especialmente os métodos da transformação inversa e da rejeição — produzem sequências de números pseudoaleatórios que simulam adequadamente a variável aleatória com a distribuição teórica desejada. Entre os testes mais utilizados destacam-se: o teste de Kolmogorov-Smirnov, o teste de Anderson-Darling e o teste de Cramér-von Mises, voltados para variáveis contínuas; e o teste qui-quadrado de aderência, apropriado para variáveis discretas.

Definição 2.4 (Teste de Kolmogorov-Smirnov - Kolmogorov 1933; Smirnov 1948). Seja \(X_1, X_2, \ldots, X_n\) uma amostra independente e identicamente distribuída (i.i.d.) de uma variável aleatória contínua com FDA definida por \(F(x)\). O teste de Kolmogorov-Smirnov (KS) tem como objetivo verificar se a amostra é derivada ou não da distribuição \(F(x)\). A estatística do teste é definida como a maior diferença absoluta entre a função distribuição empírica (FDE), \(F_n(x)\), e a FDA teórica, \(F(x)\), dada por:

\[\begin{align}\\ D_n = \sup_{x} \left| F_n(x) - F(x) \right|\\\\ \end{align}\]

em que \(F_n(x) = (1/n) \sum_{i=1}^n \mathbf{1}_{\{X_i \leqslant x\}}\) é a função distribuição empírica, e \(\mathbf{1}_{\{ \cdot \}}\) é a função indicadora. Neste teste, o valor observado de \(D_n\) é comparado com valores críticos derivados da distribuição de Kolmogorov para decidir sobre a rejeição ou não da hipótese nula, a um nível de significância pré-estabelecido \(\alpha\).

Exemplo 2.18 (Distribuição Weibull - Weibull 1951). Considere uma amostra aleatória \(X_1, X_2, \ldots, X_n\), com \(n=1000\), gerada a partir de uma distribuição Weibull com parâmetros forma \(\alpha = 2\) e escala \(\lambda = 1\). A FDA associada a distribuição é dada por:

\[\begin{align}\\ F(x; \lambda,\alpha) = 1 - \exp\left\{-\left(\frac{x}{\lambda}\right)^\alpha\right\}, \quad x > 0\\\\ \end{align}\]

em que \(\alpha > 0\) é o parâmetro de forma e \(\lambda > 0\) o parâmetro de escala. Para simular essa amostra, utiliza-se o método da transformação inversa, que consiste em aplicar a função quantil \(Q(u) = F^{-1}(u)\) a uma amostra \(u_1, u_2, \ldots, u_n\) de variáveis uniformemente distribuídas no intervalo \([0,1]\). A função quantil da distribuição Weibull é dada por:

\[\begin{align}\\ Q(u; \lambda,\alpha) = F^{-1}(u; \lambda,\alpha) = \lambda \left(-\log(1 - u)\right)^{1/\alpha}\\\\ \end{align}\]

Suponha que o objetivo deste exemplo seja aplicar o teste de Kolmogorov-Smirnov à amostra gerada para avaliar a hipótese de que os dados gerados seguem, de fato, a distribuição Weibull especificada pelos parâmetros \(\alpha = 2\) e \(\lambda = 1\) ao nível de significância de 5%. Para tal, será considerada a função ks.test do R, cujos argumentos são:

• x: corresponde ao vetor de dados amostrais a ser testado;
• name: especifica a FDA teórica correspondente à hipótese nula, especificada por uma string (por exemplo, "pweibull" para a distribuição Weibull);
• ...: corresponde à argumentos adicionais passados a função (por exemplo, os parâmetros da FDA - shape (\(\alpha\)) e scale (\(\lambda\)) no caso da Weibull).

O Código 2.24 apresenta, em ambiente R, a implementação da geração da amostra de tamanho \(n = 1000\) proveniente da distribuição Weibull, utilizando o método da transformação inversa, bem como a execução do teste de Kolmogorov-Smirnov para avaliar, ao nível de significância de 5%, a hipótese de que a amostra gerada segue a distribuição Weibull com parâmetros \(\alpha = 2\) e \(\lambda = 1\).

Código 2.24. Geração de amostra da distribuição Weibull com tamanho \(n = 1000\) e parâmetros \(\alpha = 2\) e \(\lambda = 1\), por meio do método da transformação inversa, seguida da aplicação do teste de bondade de ajuste de Kolmogorov-Smirnov.

# -----------------------------------------------------------------------------------------------------
# Teste de Kolmogorov-Smirnov para a Distribuição Weibull Com Amostra Gerada via Transformação Inversa
# -----------------------------------------------------------------------------------------------------

# --- 1. Definir parâmetros da distribuição Weibull e tamanho da amostra ---

alpha         <- 2     # Parâmetro de forma
lambda        <- 1     # Parâmetro de escala
n             <- 1000  # Tamanho da amostra

# --- 2. Gerar amostra da distribuição Weibull via transformação inversa ---

set.seed(123)
u             <- runif(n)
amostra       <- lambda * (-log(1 - u))^(1 / alpha)

# --- 3. Visualizar as primeiras observações geradas ---

head(amostra, n = 50)

 [1] 0.5823093 1.2460375 0.7251898 1.4648311 1.6796515 0.2159325 0.8666025
 [8] 1.4931549 0.8953780 0.7809844 1.7727625 0.7771213 1.0638944 0.9220158
[15] 0.3295686 1.5168508 0.5314878 0.2072912 0.6303800 1.7578746 1.4842829
[22] 1.0864011 1.0114643 2.2720037 1.0325981 1.1103242 0.8862321 0.9496062
[29] 0.5842153 0.3989097 1.8158999 1.5250718 1.0832639 1.2597730 0.1578660
[36] 0.8060378 1.1919389 0.4938286 0.6188627 0.5133015 0.3925354 0.7316886
[43] 0.7307292 0.6783838 0.4066931 0.3865690 0.5150854 0.7920159 0.5560656
[50] 1.3966802

# --- 4. Aplicar o teste de Kolmogorov-Smirnov ---

ks_result     <- ks.test(amostra, "pweibull", shape = alpha, scale = lambda)
print(ks_result)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  amostra
D = 0.014051, p-value = 0.9891
alternative hypothesis: two-sided

O resultado do teste de Kolmogorov-Smirnov aplicado à amostra revelou uma estatística observada \(D = 0.0140\), valor relativamente baixo que indica uma boa concordância entre a função distribuição empírica da amostra e a função distribuição acumulada teórica da distribuição Weibull com parâmetros \(\alpha = 2\) e \(\lambda = 1\). O valor-p associado, superior a 0.05, também mostra que não há evidência suficiente para rejeitar a hipótese nula \(H_0\), que postula que os dados amostrais seguem a distribuição Weibull especificada, considerando um nível de significância de 5%. Dessa forma, conclui-se que a amostra é proveniente com uma distribuição Weibull de parâmetros \(\alpha = 2\) e \(\lambda = 1\).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\ \end{align}\]

Exemplo 2.19 (Distribuição Gama). Considere uma amostra aleatória \(X_1, X_2, \ldots, X_n\), com \(n=1000\), gerada a partir de uma distribuição gama com parâmetros forma \(\alpha = 3\) e taxa \(\beta = 2\). A FDA associada a essa distribuição é dada por:

\[\begin{align}\\ F(x; \alpha, \beta) = \frac{1}{\Gamma(\alpha)} \gamma\left(\alpha, \beta x\right), \quad x > 0\\\\ \end{align}\]

em que \(\Gamma(\alpha)\) é a função gama e \(\gamma(\alpha, \cdot)\) é a função gama incompleta inferior. Para simular essa amostra, será utilizado o método da rejeição, com a distribuição exponencial com parâmetro \(\lambda = 0.5\) como densidade proposta. A constante ótima de rejeição, neste caso, é definida por:

\[\begin{align}\\ c &=\frac{\beta^\alpha}{\lambda \Gamma(\alpha)} \left( \frac{\alpha - 1}{\beta - \lambda} \right)^{\alpha - 1} e^{-(\alpha - 1)}\\\\ \end{align}\]

Suponha que o objetivo deste exemplo seja aplicar o teste de Kolmogorov-Smirnov para avaliar a hipótese de que os dados seguem a distribuição gama especificada pelos parâmetros \(\alpha = 3\) e \(\beta = 2\) ao nível de significância de 5%. Para tal, será considera a função ks.test. O Código 2.25 apresenta, em ambiente R, a implementação da geração de uma amostra de tamanho \(n = 1000\) proveniente da distribuição gama com parâmetros \(\alpha = 3\) e \(\beta = 2\), utilizando o método da rejeição com a distribuição exponencial de parâmetro \(\lambda = 0.5\) como densidade proposta. Além disso, o código executa o teste de Kolmogorov-Smirnov para avaliar, ao nível de significância de 5%, a hipótese de que a amostra gerada segue a distribuição gama especificada pelos parâmetros indicados, utilizando a função ks.test do ambiente R.

Código 2.25. Geração de amostra da distribuição gama com tamanho \(n = 1000\) e parâmetros \(\alpha = 3\) e \(\beta = 2\), por meio do método da rejeição com distribuição exponencial como densidade proposta, seguida da aplicação do teste de bondade de ajuste de Kolmogorov-Smirnov.

# -----------------------------------------------------------------------------------------------
# Teste de Kolmogorov-Smirnov para a Distribuição Gama com Amostra Gerada via Método da Rejeição
# -----------------------------------------------------------------------------------------------

# --- 1. Definir parâmetros iniciais e variáveis auxiliares ---

set.seed(1212)
n         <- 1000
x         <- numeric(n)
i         <- 1
alpha     <- 3
beta      <- 2
lambda    <- 0.5

# --- 2. Calcular a constante de normalização c ---

c_1       <- (beta^alpha * (alpha - 1)^(alpha - 1) * exp(-(alpha - 1)))
c_2       <- (lambda * gamma(alpha) * (beta - lambda)^(alpha - 1))
c         <- c_1 / c_2

# --- 3. Gerar amostra da distribuição gama via método da rejeição ---

while (i < n)
{
  y       <- rexp(1, rate = lambda)
  u       <- runif(1)
  
  if (u <= dgamma(y, shape = alpha, rate = beta) / (c * dexp(y, rate = lambda))) 
  {
    x[i]  <- y
    i     <- i + 1
  }
}

# --- 4. Visualizar as primeiras observações geradas ---

head(x, n = 50)

 [1] 1.5034659 1.2423327 1.4033567 2.3361018 1.3445968 1.8822605 1.6358332
 [8] 2.2705762 0.4197185 1.5154502 1.1434963 0.4764195 0.2937304 0.9712513
[15] 2.9406662 1.6990578 0.4721587 1.8244135 3.3208655 0.8754741 1.4324538
[22] 0.8002435 0.2901695 0.6815301 2.1593844 2.1871619 1.2581018 0.8633273
[29] 1.8582996 0.3019156 2.1228579 2.0397908 2.2844274 2.0018213 3.4929925
[36] 3.9847285 2.3428645 1.2473672 1.8908973 2.7921358 1.5843157 2.8963475
[43] 2.6562529 2.1419938 0.7978212 1.8898706 0.6100282 2.1555828 1.8219644
[50] 1.3151346

# --- 5. Aplicar o teste de Kolmogorov-Smirnov para a distribuição gama ---

ks_result <- ks.test(x, "pgamma", shape = alpha, rate = beta)
print(ks_result)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  x
D = 0.02643, p-value = 0.4872
alternative hypothesis: two-sided

O resultado do teste de Kolmogorov-Smirnov para a amostra apresenta uma estatística observada \(D = 0.0264\), valor pequeno que indica proximidade entre a função distribuição empírica da amostra e a função distribuição acumulada teórica da distribuição gama com parâmetros \(\alpha = 3\) (forma) e \(\beta = 2\) (taxa). O valor-p associado ao teste, igual a 0.4872, é significativamente maior ao nível de significância de 5%, indicando também que não há evidência estatística suficiente para rejeitar a hipótese nula \(H_0\), que postula que os dados amostrais seguem a distribuição Gama especificada. Dessa forma, conclui-se que a amostra gerada é proveniente de uma distribuição Gama de parâmetros \(\alpha = 3\) e \(\beta = 2\).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Definição 2.5 (Teste de Anderson-Darling - Anderson e Darling 1952, 1954). Seja \(X_1, X_2, \ldots, X_n\) uma amostra independente e identicamente distribuída (i.i.d.) de uma variável aleatória contínua com FDA definida por \(F(x)\). O teste de Anderson-Darling (AD) tem como objetivo avaliar a aderência dos dados à distribuição \(F(x)\), incorporando um peso maior às discrepâncias observadas nas caudas da distribuição. A estatística do teste é definida por:

\[\begin{align}\\ A^2 = -n - \frac{1}{n} \sum_{i=1}^n \left\{ (2i - 1) \left[ \log F(X_{(i)}) + \log \left\{1 - F(X_{(n+1-i)}) \right\} \right] \right\}\\\\ \end{align}\]

em que \(X_{(1)} \leqslant X_{(2)} \leqslant \cdots \leqslant X_{(n)}\) são as estatísticas de ordem da amostra, e \(F(x)\) é a FDA da distribuição sob a hipótese nula \(H_0\). A decisão sobre a rejeição da hipótese nula é feita pela comparação da estatística observada \(A^2\), a um nível de significância pré-estabelecido \(\alpha\), com os valores críticos apropriados, os quais dependem da distribuição assumida sob \(H_0\) e, em alguns casos, também da estimação dos parâmetros.

Exemplo 2.20 (Distribuição Lomax - Lomax 1954). Considere uma amostra aleatória \(X_1, X_2, \ldots, X_n\), com \(n = 1000\), gerada a partir de uma distribuição Lomax com parâmetros forma \(\alpha = 2\) e escala \(\lambda = 1\). A FDA associada a essa distribuição é dada por:

\[\begin{align}\\ F(x; \lambda,\alpha) = 1 - \left(1 + \frac{x}{\lambda}\right)^{-\alpha}, \quad x > 0\\\\ \end{align}\]

Para simular essa amostra, utiliza-se o método da transformação inversa, que consiste em aplicar a função quantil \(Q(u) = F^{-1}(u)\) a uma amostra \(u_1, u_2, \ldots, u_n\) de variáveis uniformemente distribuídas no intervalo \([0,1]\). A função quantil da distribuição Lomax é dada por:

\[\begin{align}\\ Q(u; \lambda,\alpha) = \lambda \left[(1 - u)^{-1/\alpha} - 1\right]\\\\ \end{align}\]

Suponha que o objetivo deste exemplo seja aplicar o teste de Anderson-Darling à amostra gerada para avaliar a hipótese de que os dados seguem a distribuição Lomax especificada pelos parâmetros \(\alpha = 2\) e \(\lambda = 1\), ao nível de significância de 5%. Para tal, será considerada a função ad.test do pacote ADGofTest, cujos argumentos são:

• x: corresponde ao vetor de dados amostrais a ser testado;
• distr.fun: especifica a FDA teórica correspondente à hipótese nula (por exemplo, plomax para a distribuição Lomax);
• ...: corresponde à argumentos adicionais passados a função (por exemplo, os parâmetros da FDA - shape (\(\alpha\)) e scale (\(\lambda\)) no caso da Lomax).

O Código 2.26 apresenta, em ambiente R, a implementação da geração de uma amostra de tamanho \(n = 1000\) proveniente da distribuição Lomax com parâmetros \(\alpha = 2\) e \(\lambda = 1\), utilizando o método da transformação inversa. Além disso, o código executa o teste de Anderson-Darling para avaliar, ao nível de significância de 5%, a hipótese de que a amostra gerada segue a distribuição Lomax especificada, por meio da função ad.test do pacote ADGofTest.

Código 2.26. Geração de amostra da distribuição Lomax com tamanho \(n = 1000\) e parâmetros \(\alpha = 2\) e \(\lambda = 1\), utilizando o método da transformação inversa, seguida da aplicação do teste de bondade de ajuste de Anderson-Darling.

# -------------------------------------------------------------------------------------------------
# Teste de Anderson-Darling para a Distribuição Lomax Com Amostra Gerada via Transformação Inversa
# -------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(ADGofTest)

# --- 1. Definir parâmetros da distribuição Lomax e tamanho da amostra ---

alpha     <- 2        # Parâmetro de forma
lambda    <- 1        # Parâmetro de escala
n         <- 1000     # Tamanho da amostra

# --- 2. Gerar amostra da distribuição Lomax via método da inversa da CDF ---

set.seed(123)
u         <- runif(n)
x         <- lambda * ((1 - u)^(-1 / alpha) - 1)

# --- 3. Visualizar as primeiras observações geradas ---

head(x, n = 50)

 [1]  0.18476221  1.17342590  0.30076182  1.92374420  3.09847370  0.02358730
 [7]  0.45571895  2.04882363  0.49309457  0.35658207  3.81311018  0.35250547
[13]  0.76109468  0.52967691  0.05580945  2.15951383  0.15170061  0.02171728
[19]  0.21980314  3.68826031  2.00881951  0.80422959  0.66784109 12.21034695
[25]  0.70425723  0.85226627  0.48097941  0.56968662  0.18608002  0.08281538
[31]  4.20045570  2.19926781  0.79809961  1.21115277  0.01253879  0.38382092
[37]  1.03472187  0.12967882  0.21105942  0.14081076  0.08008746  0.30693424
[43]  0.30601770  0.25872872  0.08621549  0.07757999  0.14185767  0.36840327
[49]  0.16719623  1.65211608

# --- 4. Definir a função de distribuição acumulada (FDA) da Lomax ---

plomax <- function(q, shape, scale)
{
  1 - (1 + q / scale)^(-shape)
}

# --- 5. Aplicar o teste de Anderson-Darling para a distribuição Lomax ---

ad_result <- ad.test(x, distr.fun = plomax, shape = alpha, scale = lambda)
print(ad_result)

    Anderson-Darling GoF Test

data:  x  and  plomax
AD = 0.16038, p-value = 0.9977
alternative hypothesis: NA

O resultado do teste de Anderson-Darling para a amostra apresenta uma estatística observada \(A^2 = 0{.}1604\), valor moderadamente pequeno que indica aderência razoável entre a função distribuição empírica da amostra e a função distribuição acumulada teórica da distribuição Lomax com parâmetros \(\alpha = 2\) (forma) e \(\lambda = 1\) (escala). O valor-p associado ao teste, igual a 0.9977, é maior que o nível de significância usual de 5%, sugerindo também que não há evidência estatística suficiente para rejeitar a hipótese nula \(H_0\), segundo a qual os dados seguem a distribuição Lomax especificada. Dessa forma, conclui-se que a amostra gerada é consistente com a hipótese de que foi extraída de uma distribuição Lomax com \(\alpha = 2\) e \(\lambda = 1\).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Definição 2.6 (Teste de Cramér-von Mises - Cramér 1928; Mises 1936). Seja \(X_1, X_2, \ldots, X_n\) uma amostra independente e identicamente distribuída (i.i.d.) de uma variável aleatória contínua com FDA definida por \(F(x)\). O teste de Cramér-von Mises tem como objetivo verificar a aderência da distribuição empírica dos dados à FDA \(F(x)\), considerando o erro quadrático médio entre ambas as distribuições ao longo de todo o domínio. A estatística do teste é dada por:

\[\begin{align}\\ W^2 = \int_{-\infty}^{\infty} \left[ F_n(x) - F(x) \right]^2 \, dF(x)\\\\ \end{align}\]

onde \(F_n(x) = (1/n) \sum_{i=1}^n \mathbf{1}_{\{X_i \leqslant x\}}\) é a função distribuição empírica da amostra. Na prática, a estatística \(W^2\) é computada, de forma equivalente, pela expressão:

\[\begin{align}\\ W^2 = \frac{1}{12n} + \sum_{i=1}^{n} \left[ F(X_{(i)}) - \frac{2i - 1}{2n} \right]^2\\\\ \end{align}\]

em que \(X_{(1)} \leqslant X_{(2)} \leqslant \cdots \leqslant X_{(n)}\) são as estatísticas de ordem da amostra. A decisão sobre a hipótese nula \(H_0\), que assume que os dados seguem a distribuição \(F(x)\), é feita por comparação da estatística \(W^2\) com valores críticos apropriados, a um nível de significância pré-estabelecido \(\alpha\).

Exemplo 2.21 (Distribuição Burr - Burr 1942). Considere uma amostra aleatória \(X_1, X_2, \ldots, X_n\), com \(n = 1000\), gerada a partir de uma distribuição Burr de parâmetros \(c = 3\) e \(k = 2\), cuja FDA é dada por:

\[\begin{align}\\ F(x; c, k) = 1 - \left(1 + x^c\right)^{-k}, \quad x > 0\\\\ \end{align}\]

onde \(c > 0\) é o parâmetro de forma e \(k > 0\). Para simular essa amostra, utiliza-se o método da transformação inversa, que consiste em aplicar a função quantil \(Q(u) = F^{-1}(u)\) a uma amostra \(u_1, u_2, \ldots, u_n\) de variáveis uniformemente distribuídas no intervalo \([0,1]\). A função quantil da distribuição Burr é dada por:

\[\begin{align}\\ Q(u; c, k) = \left[\left(1 - u\right)^{-1/k} - 1\right]^{1/c}\\\\ \end{align}\]

Suponha que o objetivo deste exemplo seja aplicar o teste de Cramér-von Mises à amostra gerada para avaliar a hipótese de que os dados seguem a distribuição Burr especificada com parâmetros \(c = 3\) e \(k = 2\), ao nível de significância de 5%. Para tal, será considerada a função cvm.test do pacote goftest, cujos argumentos são:

• x: corresponde ao vetor de dados amostrais a ser testado;
• null: especifica a FDA teórica correspondente à hipótese nula (por exemplo, pburr para a distribuição Lomax);
• ...: corresponde à argumentos adicionais passados a função (por exemplo, os parâmetros da FDA - c e k no caso da Burr).

O Código 2.27 apresenta, em ambiente R, a implementação da geração de uma amostra de tamanho \(n = 1000\) proveniente da distribuição Burr com parâmetros \(c = 3\) e \(k = 2\), utilizando o método da transformação inversa. Além disso, o código executa o teste de Cramér-von Mises para avaliar, ao nível de significância de 5%, a hipótese de que a amostra gerada segue a distribuição Burr especificada, por meio da função cvm.test do pacote goftest.

Código 2.27. Geração de amostra da distribuição Burr com tamanho \(n = 1000\) e parâmetros \(c = 3\) e \(k = 2\), utilizando o método da transformação inversa, seguida da aplicação do teste de bondade de ajuste de Cramér-von Mises.

# ------------------------------------------------------------------------------------------------
# Teste de Cramér-von Mises para a Distribuição Burr com Amostra Gerada via Transformação Inversa
# ------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(goftest)

# --- 1. Definir parâmetros da distribuição Burr e tamanho da amostra ---

c         <- 3       # Parâmetro de forma c
k         <- 2       # Parâmetro de forma k
n         <- 1000    # Tamanho da amostra

# --- 2. Gerar amostra da distribuição Burr via método da inversa da CDF ---

set.seed(123)
u         <- runif(n)
x         <- ((1 - u)^(-1 / k) - 1)^(1 / c)

# --- 3. Visualizar as primeiras observações geradas ---

head(x, n = 50)

 [1] 0.5695577 1.0547557 0.6699991 1.2437004 1.4578604 0.2867870 0.7695421
 [8] 1.2700911 0.7900297 0.7091202 1.5622833 0.7064075 0.9130185 0.8091028
[15] 0.3821518 1.2925638 0.5333297 0.2789985 0.6035010 1.5450428 1.2617703
[22] 0.9299509 0.8740931 2.3027281 0.8897004 0.9481094 0.7835057 0.8289825
[29] 0.5709086 0.4358834 1.6134870 1.3004471 0.9275821 1.0659405 0.2323193
[36] 0.7267352 1.0114425 0.5061622 0.5953901 0.5202498 0.4310439 0.6745515
[43] 0.6738794 0.6372085 0.4417689 0.4264976 0.5215360 0.7168712 0.5509035
[50] 1.1821707

# --- 4. Definir a função de distribuição acumulada (FDA) da Burr ---

pburr <- function(q, c, k)
{
  1 - (1 + q^c)^(-k)
}

# --- 5. Aplicar o teste de Cramér-von Mises para a distribuição Burr ---

cvm_result <- cvm.test(x, null = "pburr", c = c, k = k)
print(cvm_result)

    Cramer-von Mises test of goodness-of-fit
    Null hypothesis: distribution 'pburr'
    with parameters c = 3, k = 2
    Parameters assumed to be fixed

data:  x
omega2 = 0.023904, p-value = 0.9917

O resultado do teste de Cramér-von Mises para a amostra apresenta uma estatística observada \(W^2 = 0{.}0239\), valor pequeno que sugere boa concordância entre a função distribuição empírica da amostra e a função distribuição acumulada teórica da distribuição Burr com parâmetros \(c = 3\) e \(k = 2\). O valor-p obtido, igual a 0.9917, é superior ao nível de significância de 5%, indicando também que não há evidência estatística suficiente para rejeitar a hipótese nula \(H_0\). Conclui-se, portanto, que a amostra pode ser considerada proveniente de uma distribuição Burr com os parâmetros especificados.

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Definição 2.7 (Teste Qui-Quadrado de Aderência - Pearson 1900). Seja \(X_1, X_2, \ldots, X_n\) uma amostra de uma variável aleatória discreta que assume valores em um conjunto finito de categorias \(\{c_1, c_2, \ldots, c_k\}\). O teste qui-quadrado de aderência tem como objetivo avaliar se a distribuição empírica das frequências observadas nas \(k\) categorias é consistente com uma distribuição de probabilidades teórica especificada \(\{p_1, p_2, \ldots, p_k\}\), em que \(p_j = \mathbb{P}(X = c_j)\) sob a hipótese nula. A estatística do teste é definida por:

\[\begin{align}\\ X^2 = \sum_{j=1}^{k} \frac{(O_j - E_j)^2}{E_j}\\\\ \end{align}\]

em que \(O_j\) representa a frequência observada na categoria \(c_j\), \(E_j = n p_j\) corresponde à frequência esperada sob \(H_0\), e \(n\) é o tamanho da amostra. Sob a hipótese nula \(H_0\), que postula que os dados seguem a distribuição discreta com probabilidades \(\{p_1, \ldots, p_k\}\), a estatística \(X^2\) tem, assintoticamente, distribuição qui-quadrado com \(k - 1\) graus de liberdade. A rejeição de \(H_0\) ocorre quando o valor observado de \(X^2\) excede o valor crítico da distribuição qui-quadrado para um dado nível de significância \(\alpha\).

Exemplo 2.22 (Distribuição Geométrica). Considere uma amostra aleatória \(X_1, X_2, \ldots, X_n\), com \(n = 1000\), gerada a partir de uma distribuição geométrica com parâmetro \(p = 0{.}2\), cuja função de probabilidade é dada por:

\[\begin{align}\\ P(X = k) = p(1 - p)^k, \quad k = 0, 1, 2, \ldots\\\\ \end{align}\]

Para simular essa amostra, utiliza-se o método da transformação inversa (caso discreto) por meio da função rgeom(n, prob = p) em ambiente R. Suponha que o objetivo deste exemplo seja aplicar o teste de aderência do tipo qui-quadrado para avaliar a hipótese de que os dados seguem a distribuição Geométrica com parâmetro \(p = 0{.}2\), ao nível de significância de 5%. Para realizar o teste qui-quadrado, é necessário:

• 1º Passo: Agrupar os dados amostrais em categorias (valores inteiros observados);
• 2º Passo: Calcular as frequências observadas e as esperadas sob a hipótese nula;
• 3º Passo: Aplicar o teste qui-quadrado com base na estatística: \[\begin{align}\\ \chi^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i}\\\\ \end{align}\] onde \(O_i\) e \(E_i\) representam, respectivamente, as frequências observadas e esperadas para a \(i\)-ésima categoria. Para tal, será considerada a função chisq.test do R, cujos argumentos são:

• x: corresponde as frequências observadas dos dados amostrais a serem testados;
• p: corresponde ao vetor de probabilidades da distribuição teórica;
• rescale.p: corresponde à um valor lógico que se for igual TRUE, então p é reescalonado (se necessário) para que a soma seja 1.

O Código 2.28 apresenta, em ambiente R, a implementação da geração de uma amostra de tamanho \(n = 1000\) proveniente da distribuição geométrica com parâmetro \(p = 0.2\), utilizando o método da transformação inversa (caso discreto). Além disso, o código executa o teste qui-quadrado de aderência Mises para avaliar, ao nível de significância de 5%, a hipótese de que a amostra gerada segue a distribuição geométrica especificada, por meio da função chisq.test.

Código 2.28. Geração de amostra da distribuição geométrica com tamanho \(n = 1000\) e parâmetro \(p = 0.2\), utilizando o método da transformação inversa (caso discreto), seguida da aplicação do teste qui-quadrado de aderência.

# ----------------------------------------------------------------------------------------------------------
# Teste Qui-Quadrado para a Distribuição Geométrica Com Amostra Gerada via Transformação Inversa (Discreta)
# ----------------------------------------------------------------------------------------------------------

# --- 0. Pacotes necessários ---

library(knitr)
library(kableExtra)

# --- 1. Definir parâmetros da distribuição geométrica e tamanho da amostra ---

p         <- 0.2       # Probabilidade de sucesso
n         <- 1000      # Tamanho da amostra

# --- 2. Gerar amostra da distribuição geométrica ---

set.seed(123)
x         <- rgeom(n, prob = p)

# --- 3. Definir categorias e construir tabela de frequências observadas ---

max_cat   <- 10
obs_freq  <- table(factor(ifelse(x < max_cat, x, max_cat), 
                          levels = 0:max_cat))

# --- 4. Calcular probabilidades teóricas da distribuição geométrica agrupada ---

probs     <- c(dgeom(0:(max_cat - 1), prob = p), 
               1 - pgeom(max_cat - 1, prob = p))

# --- 5. Calcular frequências esperadas sob a hipótese geométrica ---

exp_freq  <- n * probs

# --- 6. Tabela com os valores observados e esperados, e as probabilidades teóricas ---

tab_res   <- data.frame("Probabilidades" = probs,
                        "Observado" = as.numeric(obs_freq),
                        "Esperado" = exp_freq)

kable(tab_res,
      digits = 0,
      align = 'c',
      format = 'html',
      escape = FALSE,
      col.names = c("Probabilidades <br> (Distribuição Geométrica)",
                    "Frequência <br> Observada",
                    "Frequência <br> Esperada"))

Probabilidades (Distribuição Geométrica)	Frequência Observada	Frequência Esperada
0	205	200
0	179	160
0	103	128
0	102	102
0	99	82
0	64	66
0	44	52
0	41	42
0	32	34
0	28	27
0	103	107

# --- 7. Aplicar o teste Qui-Quadrado para comparar frequências observadas e esperadas ---

test_result <- chisq.test(obs_freq, p = probs, rescale.p = TRUE)
print(test_result)

    Chi-squared test for given probabilities

data:  obs_freq
X-squared = 12.539, df = 10, p-value = 0.2506

O resultado do teste qui-quadrado para a amostra apresenta uma estatística observada \(\chi^2 = 12.5390\), com 10 graus de liberdade. Este valor sugere concordância entre a distribuição empírica da amostra e a distribuição teórica geométrica com parâmetro \(p = 0{.}2\). O valor-p correspondente, igual a 0.2506, é maior que o nível de significância usual de 5%, indicando que não há evidência estatística suficiente para rejeitar a hipótese nula \(H_0\). Dessa forma, conclui-se que a amostra pode ser considerada proveniente de uma distribuição Geométrica com parâmetro \(p = 0{.}2\).

\[\small \begin{align}\\ \tag*{$\blacksquare$}\\\\\\ \end{align}\]

Adler, Avraham. 2015. «lamW: Lambert-W Function (R package version 2.2.3)». https://CRAN.R-project.org/package=lamW.

Anderson, Theodore W, e Donald A Darling. 1952. «Asymptotic theory of certain" goodness of fit" criteria based on stochastic processes». The annals of mathematical statistics, 193–212.

———. 1954. «A test of goodness of fit». Journal of the American statistical association 49 (268): 765–69.

Brent, Richard P. 1971. «An algorithm with guaranteed convergence for finding a zero of a function». The computer journal 14 (4): 422–25.

Burr, Irving W. 1942. «Cumulative frequency functions». The Annals of mathematical statistics 13 (2): 215–32.

Cinlar, Erhan. 2013. Introduction to Stochastic Processes. Courier Corporation.

Corless, Robert M, Gaston H Gonnet, David EG Hare, David J Jeffrey, e Donald E Knuth. 1996. «On the Lambert W function». Advances in Computational mathematics 5: 329–59.

Cramér, Harald. 1928. «On the composition of elementary errors: First paper: Mathematical deductions». Scandinavian Actuarial Journal 1928 (1): 13–74.

Darmois, Georges. 1935. «Sur les lois de probabilitéa estimation exhaustive». Comptes rendus de l’Académie des Sciences 260 (1265): 85.

Devroye, Luc. 1986. Non-Uniform Random Variate Generation. New York: Springer-Verlag.

Dobrow, Robert P. 2016. Introduction to Stochastic Processes with R. John Wiley & Sons.

Givens, Geof H, e Jennifer A Hoeting. 2012. Computational Statistics. John Wiley & Sons.

Gómez-Déniz, Emilio, e Enrique Calderı́n-Ojeda. 2011. «The discrete Lindley distribution: properties and applications». Journal of statistical computation and simulation 81 (11): 1405–16.

Gompertz, Benjamin. 1825. «On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies». Philosophical Transactions of the Royal Society of London, n.º 115: 513–83.

Good, Irving J. 1953. «The population frequencies of species and the estimation of population parameters». Biometrika 40 (3-4): 237–64.

Hegazy, MA, REA El-Kader, AA El-Helbawy, e GR Al-Dayian. 2021. «Bayesian estimation and prediction of discrete gompertz distribution». Journal of Advances in Mathematics and Computer Science 36 (2): 1–21.

Hoel, Paul G, Sidney C Port, e Charles J Stone. 1986. Introduction to Stochastic Processes. Waveland Press.

Hoffman, Dan, e Otto J Karst. 1975. «The theory of the Rayleigh distribution and some of its applications». Journal of Ship Research 19 (03): 172–91.

Jodrá, Pedro. 2010. «Computer generation of random variables with Lindley or Poisson–Lindley distribution via the Lambert W function». Mathematics and Computers in Simulation 81 (4): 851–59.

Knuth, Donald E. 2014. The Art of Computer Programming: Seminumerical Algorithms. Addison-Wesley Professional.

Kolmogorov, Andrei Nikolaevich. 1933. «Sulla determinazione empirica di una legge di distribuzione». Giornale dell’Istituto Italiano degli Attuari 4: 83–91.

Koopman, Bernard Osgood. 1936. «On distributions admitting a sufficient statistic». Transactions of the American Mathematical society 39 (3): 399–409.

Krishna, Hare, e Pramendra Singh Pundir. 2009. «Discrete Burr and discrete Pareto distributions». Statistical methodology 6 (2): 177–88.

Kuş, Coşkun. 2007. «A new lifetime distribution». Computational Statistics & Data Analysis 51 (9): 4497–4509.

Lambert, Johann Heinrich. 1758. «Observationes variae in mathesin puram». Acta Helvetica 3 (1): 128–68.

Lawler, Gregory F. 2018. Introduction to Stochastic Processes. Chapman; Hall/CRC.

Lindley, Dennis V. 1958. «Fiducial distributions and Bayes’ theorem». Journal of the Royal Statistical Society. Series B (Methodological), 102–7.

Lomax, Kenneth S. 1954. «Business failures: Another example of the analysis of failure data». Journal of the American statistical association 49 (268): 847–52.

Markov, Andrei Andreevich. 1906. «Rasprostranenie zakona bol’shih chisel na velichiny, zavisyaschie drug ot druga». Izvestiya Fiziko-matematicheskogo obschestva pri Kazanskom universitete 15 (135-156): 18.

Meyn, Sean P, e Richard L Tweedie. 2012. Markov Chains and Stochastic Stability. Springer Science & Business Media.

Mises, Richard von. 1936. «Wahrscheinlichkeit, Statistik und Wahrheit: Einführung in die neue Wahrscheinlichkeitslehre und ihre Anwendung».

Moivre, Abraham de. 1711. «De mensura sortis, seu, de probabilitate eventuum in ludis a casu fortuito pendentibus». Philosophical transactions of the Royal Society of London 27 (329): 213–64.

Nakagami, Minoru. 1960. «The m-distribution—A general formula of intensity distribution of rapid fading». Em Statistical methods in radio wave propagation, 3–36. Elsevier.

Oliveira, Ricardo Puziol de, Josmar Mazucheli, e Jorge Alberto Achcar. 2017. «A comparative study between two discrete Lindley distributions». Ciência e Natura 39 (3): 539–52.

Pareto, Vilfredo. 1964. Cours d’économie Politique. Vol. 1. Librairie Droz.

Pearson, Karl. 1900. «On the Criterion that a Given System of Deviations from the Probable in the Case of a Correlated System of Variables is Such that it Can Be Reasonably Supposed to Have Arisen from Random Sampling». The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 50 (302): 157–75.

Poisson, Siméon Denis. 1837. «Probabilité des jugements en matière criminelle et en matière civile, précédées des règles générales du calcul des probabilitiés». Paris, France: Bachelier 1 (1837): 1837.

Smirnov, Nickolay. 1948. «Table for estimating the goodness of fit of empirical distributions». The annals of mathematical statistics 19 (2): 279–81.

Thom, Herbert CS. 1958. «A note on the gamma distribution». Monthly Weather Review 86 (4): 117–22.

Weibull, Waloddi. 1951. «A statistical distribution function of wide applicability». Journal of applied mechanics.