Capítulo 4

Método Monte Carlo

“The quiet statisticians have changed our world – not by discovering new facts or technical developments but by changing the ways we reason, experiment and form our opinions about it.”

– Ian Hacking, 1984

4.1. Introdução

O termo Monte Carlo possui uma conotação glamorosa por estar associado ao célebre Cassino Monte Carlo, situado no Principado de Mônaco, onde jogos de azar exploram fenômenos essencialmente aleatórios. Contudo, no âmbito científico, o método Monte Carlo está intrinsecamente ligado a Stanislaw Ulam (1909–1984) que, no final da década de 1940, atuava em projetos de armamentos nucleares no Laboratório Nacional de Los Alamos. Em 1946, físicos nucleares em Los Alamos investigavam a difusão de nêutrons no núcleo de uma arma nuclear. Embora dispusessem da maioria dos dados necessários — tais como a distância média percorrida por um nêutron em determinado material antes de colidir com um núcleo atômico, bem como a quantidade de energia que o nêutron tenderia a liberar após a colisão —, os físicos não conseguiam resolver o problema utilizando métodos matemáticos convencionais, de caráter determinístico. Diante dessa limitação, Ulam propôs recorrer a experimentos aleatórios como estratégia de solução. Sua principal contribuição consistiu em perceber que os fenômenos relacionados à difusão de nêutrons poderiam ser formulados como integrais e reinterpretados como valores esperados de variáveis aleatórias adequadamente definidas, possibilitando sua estimação por meio de procedimentos estatísticos baseados em simulações. Essa perspectiva abriu caminho para o desenvolvimento do método Monte Carlo (Metropolis e Ulam 1949; Metropolis 1987), uma técnica numérica para a resolução de problemas que, até então, eram considerados intratáveis pelas abordagens analíticas clássicas.

Figura 4.1.. Cassino Monte Carlo.

Fonte: Monte Carlo Cassino | iStock Photo. Acesso: Julho, 2025.

4.2. Definição do Método Monte Carlo

Considere uma integral da forma:

\[\begin{align}\\ \theta = \int_D g(x) \, dx\\\\ \end{align}\]

em que \(g: D \subseteq \mathbb{R}^p \to \mathbb{R}\) é uma função mensurável e integrável. Quando essa integral admite solução analítica, não há necessidade de recorrer a métodos de aproximação numérica. Entretanto, caso a solução analítica seja inviável — e sobretudo quando o domínio \(D\) possui uma ou duas dimensões (\(p = 1\) ou \(p = 2\)) —, é comum empregar métodos numéricos de quadratura, tais como a Regra do Trapézio ou Regra de Simpson. Contudo, à medida que a dimensão \(p\) aumenta, esses métodos tornam-se progressivamente ineficazes, em razão da chamada maldição da dimensionalidade, pela qual a complexidade computacional cresce de forma exponencial com \(p\).

Nessas circunstâncias, o uso do método Monte Carlo, conforme proposto por Ulam (Metropolis e Ulam 1949), revela-se vantajoso, pois sua taxa de convergência não depende da dimensionalidade do espaço de integração, uma vez que este método baseia-se em uma aproximação estocástica por meio de amostragem aleatória, substituindo a avaliação determinística da integral. Em essência, o método Monte Carlo estima o valor de integrais por meio da amostragem aleatória da função \(g\) sobre seu domínio de suporte \(p\)-dimensional, sem a necessidade de uma exploração sistemática de todo o domínio. Para isso, supõe-se que a função \(g\) possa ser decomposta como \(g(x) = h(x), f(x)\), em que \(f(x)\) é uma densidade de probabilidade definida sobre \(D\), isto é,

\[\begin{align}\\ \int_D f(x) = 1, \quad f(x) \geqslant 0\\\\ \end{align}\]

Então, a integral \(\theta\) pode ser reinterpretada como o valor esperado de uma variável aleatória \(Y = h(X)\), onde \(X \sim f\). Nesse contexto, \(Y\) é uma transformação da variável aleatória \(X\), e a integral original assume a forma:

\[\begin{align}\\ \theta = E[Y] = E[h(X)] = \int_D h(x) f(x) \, dx\\\\ \end{align}\]

Essa reformulação permite estimar \(\theta\) por meio da média amostral de observações independentes da variável aleatória \(Y\), isto é,

\[\begin{align}\\ \hat{\theta}_{MC} = \frac{1}{n} \sum_{i=1}^n Y_i = \frac{1}{n} \sum_{i=1}^n h(X_i)\\\\ \end{align}\]

em que \(X_1, \ldots, X_n\) são amostras independentes e identicamente distribuídas segundo a densidade \(f\). O valor \(\hat{\theta}_n\) é conhecido como estimador Monte Carlo para \(\theta\).

Definição 4.1 (Estimador Monte Carlo - Givens e Hoeting 2012). Seja \(X \sim f\) uma variável aleatória com densidade de probabilidade \(f: D \subseteq \mathbb{R}^p \to \mathbb{R}_+\), tal que \(\int_D f(x)\,dx = 1\), e seja \(h: D \to \mathbb{R}\) uma função mensurável para a qual \(\theta = E[h(X)]\) existe e é finito. Seja \(X_1, X_2, \ldots, X_n\) uma amostra aleatória extraída de forma independente e identicamente distribuída segundo \(f\). O estimador de Monte Carlo para \(\theta\) é definido por:

\[\begin{align}\\ \hat{\theta}_{MC} = \frac{1}{n} \sum_{i=1}^n h(X_i)\\\\ \end{align}\]

Note que \(\hat{\theta}_{MC}\) é um estimador não-viciado para \(\theta\). De fato, dado que os \(X_i\) são variáveis independentes e identicamente distribuídas segundo \(f\), com \(E[h(X_i)] = \theta\) para todo \(i\), tem-se que.

\[\begin{align}\\ E[\hat{\theta}_{MC}] = \frac{1}{n} \sum_{i=1}^n E[h(X_i)] = \frac{1}{n}\cdot n\theta = \theta\\\\ \end{align}\]

Exemplo 4.1 (Integral Unidimensional no Intervalo [0,1]). Seja \(X_1, \ldots, X_n\) uma amostra aleatória independente e identicamente distribuída segundo a distribuição uniforme \(\mathcal{U}(0,1)\). Considere o problema de estimar a integral:

\[\begin{align}\\ \theta = \int_{0}^{1} e^{-x} \, dx\\\\ \end{align}\]

Neste caso, como a densidade da distribuição uniforme no intervalo \([0,1]\) é \(f(x) = 1\), a função integranda \(g(x) = e^{-x}\) pode ser reescrita como \(g(x) = h(x) f(x)\), com \(h(x) = e^{-x}\). Assim, a integral \(\theta\) pode ser reinterpretada como o valor esperado:

\[\begin{align}\\ \theta = E[h(X)] = E[e^{-X}], \quad X \sim \mathcal{U}(0,1)\\\\ \end{align}\]

Dessa forma, utilizando o método Monte Carlo com \(n\) observações \(X_1, \ldots, X_n \sim \mathcal{U}(0,1)\), obtém-se o estimador:

\[\begin{align}\\ \hat{\theta}_{MC} = \frac{1}{n} \sum_{i=1}^n e^{-X_i}\\\\ \end{align}\]

O Código 4.1 apresenta, em ambiente R, uma rotina para o cálculo do estimador Monte Carlo de \(\theta\) com \(n = 10000\) observações, enquanto a Figura 4.2 ilustra a convergência do estimador ao longo das simulações.

Código 4.1. Implementação em R do estimador Monte Carlo para aproximação da integral \(\int_0^1 e^{-x} \, dx\), utilizando amostras da distribuição uniforme no intervalo [0,1].

# ----------------------------------
# Integração Monte Carlo via U(0,1)
# ----------------------------------

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

library(ggplot2)
library(dplyr)

# --- 1. Geração da amostra aleatória U(0,1) ---

set.seed(123)
n         <- 10000
x         <- runif(n)

# --- 2. Cálculo do estimador Monte Carlo ---

y         <- exp(-x)
theta_hat <- mean(y)
theta_hat

[1] 0.6333057

# --- 3. Ilustração da convergência do estimador ---

theta_acum     <- cumsum(y) / seq_along(y)
theta_exact    <- 1 - exp(-1)

df_convergence <- data.frame(n = 1:n, theta_acum = theta_acum)

ggplot(df_convergence, aes(x = n)) +
  geom_line(aes(y = theta_acum, color = "Estimador Monte Carlo"), size = 0.5) +
  geom_hline(aes(yintercept = theta_exact, color = "Valor Exato"), 
             linetype = "dashed", size = 0.5) +
  scale_color_manual(name = NULL,
                     values = c("Estimador Monte Carlo" = "blue", 
                                "Valor Exato" = "red")) +
  labs(title = "",
       x = "Tamanho da amostra (n)",
       y = expression(hat(theta)[MC])) +
  theme_minimal() +
  theme(axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)),
        legend.position = "right")

Figura 4.2.. Convergência do estimador Monte Carlo para o valor exato da integral \(\theta\). A linha contínua azul representa o valor estimado conforme o aumento do tamanho da amostra, enquanto a linha tracejada vermelha indica o valor exato da integral.

A execução do Código 4.1 retorna a estimativa \(\hat{\theta}_{MC} \approx 0{.}6333\), a qual se mostra bastante próxima do valor exato da integral, obtido analiticamente como:

\[\begin{align}\\ \theta = \int_0^1 e^{-x} \, dx = -e^{-x} \Big|_0^1 = 1 - e^{-1} \approx 0{.}6321\\\\ \end{align}\]

Esse resultado evidencia a acurácia do método Monte Carlo, mesmo com um número moderado de simulações.

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

Exemplo 4.2 (Integral Unidimensional no Intervalo [a,b]). Considere o caso em que a região de integração \(\mathcal{D}\) corresponde ao intervalo \([a, b]\), tal que \(X \sim \mathcal{U}(a,b)\). Neste contexto, a densidade de probabilidade associada é definida por \(f(x) = 1/(b-a)\), o que permite reescrever a integral de interesse na forma:

\[\begin{align}\\ \theta = \frac{1}{f(x)} \int_{a}^{b} h(x)f(x)\,dx = (b - a) \int_{a}^{b} h(x)\,\frac{1}{b - a}\,dx = (b - a)\, E[h(X)]\\\\ \end{align}\]

Deste modo, o problema de calcular a integral \(\theta\) é convertido, via técnica de Monte Carlo, em uma tarefa estatística de estimação do valor esperado da função \(h(X)\). Assim, considerando uma amostra aleatória de tamanho \(n\), a estimativa de \(\theta\) por meio do método do Monte Carlo é dada por:

\[\begin{align}\\ \hat{\theta}_{MC} = (b - a) \cdot \frac{1}{n} \sum_{i=1}^{n} h(X_i)\\\\ \end{align}\]

em que \(X_1, \ldots, X_n\) são variáveis aleatórias independentes e identicamente distribuídas segundo a distribuição uniforme \(\mathcal{U}(a,b)\). Para fins ilustrativos, considere o problema de estimar numericamente a seguinte integral definida:

\[\begin{align}\\ \theta = \int_{2}^{4} e^{-x} \, dx\\\\ \end{align}\]

Neste exemplo, a densidade da distribuição uniforme no intervalo \([2,4]\) é constante e igual a \(f(x) = 1/2\). A função integranda pode, então, ser interpretada como \(g(x) = h(x) f(x)\), com \(h(x) = 2 e^{-x}\), de modo que a integral \(\theta\) admite a seguinte representação:

\[\begin{align}\\ \theta = \int_2^4 e^{-x} \, dx = (b - a) \cdot E[e^{-X}], \quad X \sim \mathcal{U}(2,4)\\\ \end{align}\]

Consequentemente, a estimativa Monte Carlo para \(\theta\), baseada em \(n\) amostras independentes da distribuição \(\mathcal{U}(2,4)\), é expressa por:

\[\begin{align}\\ \hat{\theta}_{MC} = (b - a) \cdot \frac{1}{n} \sum_{i=1}^n e^{-X_i} = 2 \cdot \frac{1}{n} \sum_{i=1}^n e^{-X_i}\\\\ \end{align}\]

Embora a integral em questão admita solução analítica — neste caso,

\[\begin{align}\\ \theta = \int_{2}^{4} e^{-x}dx = - e^{-x} \Big|_{2}^{4} = e^{-2} - e^{-4} \approx 0{.}1170\\\\ \end{align}\]

— o propósito aqui é demonstrar a aplicação do método de integração simples Monte Carlo como ferramenta de aproximação numérica. Nesse sentido, o Código 4.2 apresenta, em ambiente R, uma rotina para o cálculo do estimador Monte Carlo de \(\theta\) com \(n = 10000\) observações e amostras geradas de uma distribuição uniforme \(\mathcal{U}(2,4)\), enquanto a Figura 4.3 ilustra a convergência do estimador ao longo das simulações.

Código 4.2. Implementação em R do estimador Monte Carlo para aproximação da integral \(\int_2^4 e^{-x} \, dx\), utilizando amostras da distribuição uniforme no intervalo [2,4].

# ----------------------------------
# Integração Monte Carlo via U(0,1)
# ----------------------------------

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

library(ggplot2)
library(dplyr)

# --- 1. Geração da amostra aleatória U(0,1) ---

set.seed(123)
n         <- 10000
x         <- runif(n, min = 2, max = 4)

# --- 2. Cálculo do estimador Monte Carlo ---

y         <- exp(-x)
theta_hat <- (4 - 2) * mean(y)
theta_hat

[1] 0.117334

# --- 3. Ilustração da convergência do estimador ---

theta_acum     <- 2 * cumsum(y) / seq_along(y)
theta_exact    <- exp(-2) - exp(-4)

df_convergence <- data.frame(n = 1:n, theta_acum = theta_acum)

ggplot(df_convergence, aes(x = n)) +
  geom_line(aes(y = theta_acum, color = "Estimador Monte Carlo"), size = 0.5) +
  geom_hline(aes(yintercept = theta_exact, color = "Valor Exato"), 
             linetype = "dashed", size = 0.5) +
  scale_color_manual(name = NULL,
                     values = c("Estimador Monte Carlo" = "blue", 
                                "Valor Exato" = "red")) +
  labs(title = "",
       x = "Tamanho da amostra (n)",
       y = expression(hat(theta)[MC])) +
  theme_minimal() +
  theme(axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)),
        legend.position = "right")

Figura 4.3.. Convergência do estimador Monte Carlo para o valor exato da integral \(\theta\). A linha contínua azul representa o valor estimado conforme o aumento do tamanho da amostra, enquanto a linha tracejada vermelha indica o valor exato da integral.

A execução do Código 4.2 retorna a estimativa \(\hat{\theta}_{MC} \approx 0{.}117334\), a qual se mostra bastante próxima do valor exato da integral, obtido analiticamente como \(\theta = 0{.}1170\), e evidenciando a acurácia do método Monte Carlo, mesmo com um número moderado de simulações.

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

Definição 4.2 (Variância do Estimador Monte Carlo - Givens e Hoeting 2012). Seja \(\nu(X) = [h(X) - \theta]^2\), onde \(\theta = E_f[h(X)]\), e suponha que \(E_f[h(X)^2] < \infty\). Sob tais condições, a variância do estimador Monte Carlo \(\hat{\theta}_{MC}\) pode ser expressa como:

\[\begin{align}\\ \text{Var}(\hat{\theta}_{MC}) = E_f\left\{ \frac{\nu(X)}{n} \right\} = \frac{\sigma^2}{n}\\\\ \end{align}\]

em que \(\sigma^2 = \text{Var}_f[h(X)]\), e o valor esperado é tomado com respeito à densidade de probabilidade \(f\). Na prática, o valor de \(\sigma^2\) é geralmente desconhecido, mas pode ser estimado a partir da própria amostra gerada. Para isso, recorre-se ao estimador:

\[\begin{align}\\ \hat{\sigma}^2 = \frac{1}{n - 1} \sum_{i=1}^n \left[ h(X_i) - \hat{\theta}_{MC} \right]^2\\\\ \end{align}\]

que corresponde a variância amostral. Esse estimador é não-viesado e converge quase certamente para \(\sigma^2\) quando \(n \to \infty\), desde que \(h(X)\) possua variância finita sob \(f\).

Exemplo 4.3. Considere o problema de estimar numericamente o valor da integral definida por:

\[\begin{align}\\ \theta = \int_{0}^{1} \sqrt{1 - x^2} \, dx\\\\ \end{align}\]

Este integral corresponde geometricamente à área de um quarto de círculo de raio 1. Seu valor exato é conhecido e pode ser calculado analiticamente como \(\theta = \pi/4 \approx 0,7854\). Para estimar \(\theta\) via o método Monte Carlo, considere uma variável aleatória \(X\) com distribuição uniforme contínua no intervalo \([0,1]\), ou seja, \(X \sim \mathcal{U}(0,1)\). Desta forma, a integral pode ser expressa como o valor esperado da função \(h(x) = \sqrt{1 - x^2},\), isto é,

\[\begin{align}\\ \theta = \int_0^1 h(x) \, dx = E[h(X)]\\\\ \end{align}\]

Com base em uma amostra aleatória independente e identicamente distribuída \(X_1, X_2, \ldots, X_n \sim \mathcal{U}(0,1)\), o estimador de Monte Carlo para \(\theta\) é dado por:

\[\begin{align}\\ \hat{\theta}_{MC} = \frac{1}{n} \sum_{i=1}^n \sqrt{1 - X_i^2}\\\\ \end{align}\]

Além da estimativa pontual, é importante quantificar a variabilidade deste estimador para avaliar a precisão da aproximação numérica. Assim, supondo que a variância da função \(h(X)\) seja finita, definimos \(\sigma^2 = \mathrm{Var}[h(X)] = \mathrm{Var}\left(\sqrt{1 - X^2}\right)\). Como \(\sigma^2\) é geralmente desconhecida, utiliza-se a variância amostral para sua estimação, isto é,

\[\begin{align}\\ \hat{\sigma}^2 = \frac{1}{n-1} \sum_{i=1}^n \left(h(X_i) - \hat{\theta}_{MC}\right)^2\\\\ \end{align}\]

de forma que uma estimativa da variância do estimador de Monte Carlo, \(\hat{\theta}_{MC}\), é dada por:

\[\begin{align}\\ \widehat{\mathrm{Var}}(\hat{\theta}_{MC}) = \frac{\hat{\sigma}^2}{n}\\\\ \end{align}\]

O Código 4.3 apresenta, em ambiente R, uma rotina para o cálculo do estimador Monte Carlo de \(\theta\) com \(n = 10000\) observações.

Código 4.3. Implementação em R do estimador Monte Carlo e cálculo da variância para a aproximação da integral \(\int_0^1 \sqrt{1 - x^2} \, dx\), usando amostras da distribuição uniforme no intervalo [0,1].

# -----------------------------------
# Variância do Estimador Monte Carlo
# -----------------------------------

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

library(knitr)

# --- 1. Geração da amostra aleatória U(0,1) ---

set.seed(2025)
n             <- 10000
X             <- runif(n, min = 0, max = 1)

# --- 2. Cálculo do estimador Monte Carlo e variância associada ---

hX            <- sqrt(1 - X^2)

theta_hat     <- mean(hX)          # Estimador Monte Carlo
sigma2_hat    <- var(hX)           # Variância amostral de h(X)
var_theta_hat <- sigma2_hat / n    # Variância do estimador Monte Carlo

# --- 3. Apresentação dos resultados ---

resultados    <- data.frame(row.names = c("Estimativa Monte Carlo de θ",
                                       "Estimativa da variância do estimador Monte Carlo de θ",
                                       "Valor exato da integral"),
                            Valor = c(round(theta_hat, 6),
                                      round(var_theta_hat, 6),
                                      round(pi / 4, 6)))
kable(resultados, align = 'c')

	Valor
Estimativa Monte Carlo de θ	0.785920
Estimativa da variância do estimador Monte Carlo de θ	0.000005
Valor exato da integral	0.785398

A estimativa de \(\theta\) obtida pela média amostral da função \(h(X) = \sqrt{1 - X^2}\) aproxima-se do valor exato \(\pi/4 \approx 0,7854\) à medida que o tamanho da amostra cresce, refletindo a consistência do estimador Monte Carlo. A variância estimada do estimador, dada por \(\hat{\sigma}^2 / n\), quantifica a incerteza associada à simulação, possibilitando a avaliação da precisão da aproximação numérica. Para \(n = 10\,000\) amostras, os resultados mostraram uma estimativa de \(\theta\) de aproximadamente 0,7859, com variância estimada em torno de \(5 \times 10^{-6}\), valores esses muito próximos do valor exato \(0,7854\), o que evidencia a eficácia e a precisão do método aplicado.

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

Observação 4.1. A partir dos conceitos apresentados nas Definições 4.1 e 4.2, é possível construir a distribuição amostral do estimador Monte Carlo \(\hat{\theta}_{MC}\). Para tal, considere:

\[\begin{align}\\ \hat{\theta}_{MC} = \frac{1}{n} \sum_{i=1}^n h(X_i)\\\\ \end{align}\]

em que \(X_i \sim f(x)\) independentemente e identicamente distribuídos, e \(h: \mathbb{R} \to \mathbb{R}\) é uma função tal que \(E[h(X)] = \theta\) e \(\mathrm{Var}[h(X)] = \sigma^2 < \infty\). Sob as condições do Teorema Central do Limite, e para \(n\) suficientemente grande, a distribuição assintótica de \(\hat{\theta}_{MC}\) é aproximadamente normal, com média \(\theta\) e variância \(\sigma^2/n\). Matematicamente, tem-se que:

\[\begin{align}\\ \frac{\hat{\theta}_{MC} - \theta}{\sqrt{\sigma^2 / n}} \xrightarrow{D} \mathcal{N}(0,1)\\\\ \end{align}\]

quando \(n \to \infty\). Assim, um intervalo de confiança de nível \(100(1 - \alpha)\%\) para o parâmetro \(\theta = E[h(X)]\) é dado por:

\[\begin{align}\\ IC_\theta = \left[ \hat{\theta}_{MC} - z_{\alpha/2} \sqrt{\frac{\hat{\sigma}^2}{n}}, \; \hat{\theta}_{MC} + z_{\alpha/2} \sqrt{\frac{\hat{\sigma}^2}{n}} \right]\\\\ \end{align}\]

em que \(z_{\alpha/2}\) representa o quantil superior de ordem \(1 - \alpha/2\) da distribuição normal padrão \(\mathcal{N}(0,1)\), e \(\hat{\sigma}^2\) é uma estimativa consistente de \(\sigma^2\). Esse intervalo baseia-se na normalidade assintótica de \(\hat{\theta}_{MC}\), permitindo quantificar a incerteza associada à estimativa computacional do parâmetro de interesse.

4.3. Simulação Monte Carlo

Quando se deseja estimar o valor de uma integral ou esperança matemática associada a uma função mensurável, uma estratégia natural é empregar o estimador Monte Carlo, definido como a média amostral de valores obtidos a partir de uma amostra aleatória extraída da distribuição de interesse. Esse estimador fornece uma aproximação pontual da quantidade desejada, com simplicidade computacional e independência da dimensionalidade do problema.

No entanto, para investigar as propriedades estatísticas desse estimador — isto é, aspectos teóricos que caracterizam seu desempenho — é necessário ir além de uma única aplicação. Nesse contexto, introduz-se a simulação Monte Carlo, que consiste em repetir o processo de estimação diversas vezes. Mais precisamente, geram-se \(R\) amostras independentes, cada uma com tamanho \(n\), provenientes da distribuição de probabilidade subjacente ao modelo. Em cada amostra, calcula-se o estimador de interesse, obtendo-se uma coleção de \(R\) valores simulados. Essa coleção constitui uma amostra empírica da distribuição amostral do estimador, cuja análise permite inferir propriedades teóricas tais como viés, variância, erro quadrático médio, distribuição assintótica e intervalos de confiança.

Definição 4.3 (Distribuição Amostral de \(\hat{\theta}\) via Simulação Monte Carlo - Givens e Hoeting 2012). Seja \(\hat{\theta}\) um estimador para um parâmetro populacional \(\theta\), baseado em uma amostra aleatória de tamanho \(n\), cujas observações são extraídas de uma distribuição \(f\). A simulação Monte Carlo da distribuição amostral de \(\hat{\theta}\) envolve os seguintes passos:

• 1º Passo: Repetir o experimento de amostragem \(R\) vezes, de forma independente, gerando amostras \(\{X_1^{(r)}, \ldots, X_n^{(r)}\}\), para \(r = 1, \ldots, R\), em que cada amostra é composta por \(n\) observações independentes e identicamente distribuídas segundo \(f\).

• 2º Passo: Calcular o estimador \(\hat{\theta}\) para cada uma das amostras geradas, obtendo assim \(R\) estimativas, denotadas por \(\hat{\theta}^{(r)}\).

A coleção de valores \(\{\hat{\theta}^{(1)}, \ldots, \hat{\theta}^{(R)}\}\) constitui, portanto, uma amostra empírica da distribuição amostral do estimador \(\hat{\theta}\). Sob condições de regularidade — tais como independência entre observações, identidade de distribuição, existência de momentos finitos de ordem superior e tamanho amostral \(n\) suficientemente grande —, o Teorema Central do Limite assegura que a distribuição amostral de \(\hat{\theta}\) converge, em distribuição, para uma distribuição normal padrão, isto é,

\[\begin{align}\\ \frac{\hat{\theta} - \theta}{\sqrt{\mathrm{Var}(\hat{\theta})}} \xrightarrow{D} \mathcal{N}(0,1)\\\\ \end{align}\]

quando \(n \to \infty\). Na prática, como o valor da variância \(\mathrm{Var}(\hat{\theta})\) é geralmente desconhecido, costuma-se estimá-la a partir das replicações obtidas na simulação. Nesse caso, o estimador da variância via simulação Monte Carlo é dado por:

\[\begin{align}\\ \widehat{\mathrm{Var}}_{MC}(\hat{\theta}) = \frac{1}{R - 1} \sum_{r=1}^R \left(\hat{\theta}^{(r)} - \bar{\theta}_{MC}\right)^2\\\\ \end{align}\]

em que \(\bar{\theta}_{MC}\) representa a média das estimativas simuladas, isto é, o estimador Monte Carlo:

\[\begin{align}\\ \bar{\theta}_{MC} = \frac{1}{R} \sum_{r=1}^R \hat{\theta}^{(r)} \approx E[\hat\theta]\\\\ \end{align}\]

Dessa forma, a simulação Monte Carlo fornece uma aproximação numérica da distribuição amostral de \(\hat{\theta}\), permitindo estimar, de maneira empírica, suas propriedades inferenciais, como viés, variância, erro quadrático médio e comportamento assintótico.

Exemplo 4.4 (Distribuição Amostral de \(\hat{\lambda}_{\text{MLE}}\) da Distribuição Exponencial via Simulação Monte Carlo). Seja \(X\) uma variável aleatória contínua tal que \(X \sim\) Exponencial\((\lambda)\), com parâmetro de taxa \(\lambda > 0\), cuja função densidade de probabilidade é dada por:

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

Considere uma amostra aleatória simples \(X_1, X_2, \ldots, X_n\), composta por observações independentes e identicamente distribuídas segundo a distribuição exponencial acima. O estimador de máxima verossimilhança (MLE) para o parâmetro \(\lambda\) é conhecido e assume a forma:

\[\begin{align}\\ \hat{\lambda}_{\text{MLE}} = \frac{1}{\overline{X}}\\\\ \end{align}\]

em que \(\overline{X} = (1/n) \sum_{i=1}^{n} X_i\). Nosso objetivo é estudar a distribuição amostral do estimador \(\hat{\lambda}_{\text{MLE}}\) via simulação Monte Carlo, conforme estabelecido na Definição 4.3. Para tal, procede-se da seguinte forma:

• 1º Passo (Especificação dos Parâmetros da Simulação). Define-se o número de replicações \(R \in \mathbb{N}\) e o tamanho amostral \(n\). Estipula-se, ainda, o valor verdadeiro do parâmetro populacional \(\lambda\), que será utilizado na geração das amostras.

• 2º Passo (Geração das Replicações Amostrais). Para cada replicação \(r = 1, \ldots, R\), simula-se uma amostra independente \(\left\{X_1^{(r)}, \ldots, X_n^{(r)}\right\}\), com \(X_i^{(r)} \sim\) Exponencial\((\lambda)\), independentes e identicamente distribuídas.

• 3º Passo (Cálculo do Estimador em Cada Replicação). Para cada uma das \(R\) amostras geradas, calcula-se o valor estimado do parâmetro \(\lambda\) via máxima verossimilhança: \[\begin{align}\\ \hat{\lambda}_{\mathrm{MLE}}^{(r)} = \frac{1}{\bar{X}^{(r)}}\\\\ \end{align}\] em que \(\bar{X}^{(r)} = (1/n) \sum_{i=1}^{n} X_i^{(r)}\).

• 4º Passo (Construção da Distribuição do Estimador). A coleção \(\left\{\hat{\lambda}_{\mathrm{MLE}}^{(1)}, \ldots, \hat{\lambda}_{\mathrm{MLE}}^{(R)}\right\}\) constitui uma amostra empírica da distribuição amostral do estimador \(\hat{\lambda}_{\mathrm{MLE}}\).

• 5º Passo (Estimação de Propriedades Inferenciais). A partir dessa amostra, podem-se estimar o valor esperado e a variância do estimador pelas equações: \[\begin{align}\\ \bar{\lambda}_{MC} &= \frac{1}{R} \sum_{r=1}^R \hat{\lambda}_{\mathrm{MLE}}^{(r)} \approx E[\hat{\lambda}_{\mathrm{MLE}}]\\\\ \widehat{\mathrm{Var}}_{MC}(\hat{\lambda}_{\mathrm{MLE}}) &= \frac{1}{R - 1} \sum_{r=1}^R \left( \hat{\lambda}_{\mathrm{MLE}}^{(r)} - \bar{\lambda}_{MC} \right)^2 \approx\mathrm{Var}(\hat{\lambda}_{\mathrm{MLE}})\\\\ \end{align}\]

Assim, sob condições de regularidade — como independência entre observações, existência de momentos finitos e tamanho amostral \(n\) suficientemente grande —, o Teorema Central do Limite assegura que a distribuição amostral de \(\hat{\lambda}_{\mathrm{MLE}}\) pode ser aproximada por uma distribuição normal, com média \(E[\hat{\lambda}_{\mathrm{MLE}}]\) e variância \(\mathrm{Var}(\hat{\lambda}_{\mathrm{MLE}})\). Na prática, essas quantidades podem ser estimadas empiricamente pelas estatísticas \(\bar{\lambda}_{MC}\) e \(\widehat{\mathrm{Var}}_{MC}(\hat{\lambda}_{\mathrm{MLE}})\), respectivamente.

O Código 4.4 apresenta, em ambiente R, a implementação da simulação Monte Carlo voltada ao estudo da distribuição amostral do estimador de máxima verossimilhança \(\hat{\lambda}_{\mathrm{MLE}}\), considerando dados gerados a partir de uma distribuição exponencial com parâmetro verdadeiro \(\lambda = 2\). A Figura 4.4, por sua vez, apresenta o histograma das estimativas obtidas nas replicações, ilustrando a distribuição empírica do estimador.

Código 4.4. Estudo da distribuição amostral de \(\hat{\lambda}_{\mathrm{MLE}}\) via simulação Monte Carlo, para dados gerados de uma distribuição exponencial com parâmetro verdadeiro \(\lambda = 2\).

# --------------------------------------------------------------
# Distribuição Amostral do MLE de λ da Distribuição Exponencial
# --------------------------------------------------------------

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

library(ggplot2)
library(knitr)

# --- 1. Definição dos parâmetros da simulação ---

lambda_true <- 2       # Valor verdadeiro do parâmetro λ
n           <- 1000    # Tamanho da amostra em cada replicação
R           <- 10000   # Número de replicações Monte Carlo

# --- 2. Inicialização do vetor para armazenar os estimadores ---

lambda_mle  <- numeric(R)

# --- 3. Simulação Monte Carlo ---

set.seed(2025)

for (r in 1:R)
{
  # Geração de uma amostra i.i.d. da Exp(λ)
  
  amostra       <- rexp(n, rate = lambda_true)
  
  # Cálculo do MLE de λ para a amostra gerada
  
  lambda_mle[r] <- 1 / mean(amostra)
}

# --- 4. Cálculo das estimativas via Monte Carlo ---

media_mc        <- mean(lambda_mle)        # Estimativa do valor esperado do MLE de λ
var_mc          <- var(lambda_mle)         # Estimativa da variância do MLE de λ

# --- 5. Apresentação dos resultados ---

resultados <- data.frame(row.names = c("Estimativa Monte Carlo do Valor Esperado do MLE de λ",
                                       "Estimativa Monte Carlo da Variância do MLE de λ"),
                         Valor = c(round(media_mc, 5), 
                                   round(var_mc, 5)))
kable(resultados, align = 'c')

	Valor
Estimativa Monte Carlo do Valor Esperado do MLE de λ	2.00150
Estimativa Monte Carlo da Variância do MLE de λ	0.00393

# --- 6. Visualização gráfica da distribuição amostral ---

df_dist     <- data.frame(lambda_mle = lambda_mle)

ggplot(df_dist, aes(x = lambda_mle)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "black") +
  geom_density(aes(color = "Densidade Empírica (Kernel)"), size = 0.5, key_glyph = "path") +
  stat_function(fun = dnorm,
                args = list(mean = media_mc, sd = sqrt(var_mc)),
                aes(color = "Aproximação Normal (Simulação MC)"),
                size = 0.5) +
  scale_color_manual(name = 'Curvas',
                     values = c("Densidade Empírica (Kernel)" = "blue",
                                "Aproximação Normal (Simulação MC)" = "red")) +
  labs(title = '',
       x = expression(hat(lambda)[MLE]),
       y = "Densidade") +
  theme_minimal() +
  theme(legend.position = "right", 
        axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)))

Figura 4.4.. Histograma da distribuição amostral do estimador de máxima verossimilhança para o parâmetro \(\lambda\) da distribuição exponencial, obtido por simulação Monte Carlo com \(R=10000\) replicações e tamanho amostral \(n=1000\). A curva azul representa a densidade empírica estimada por kernel, enquanto a curva vermelha corresponde à aproximação normal baseada na média e variância calculadas a partir das simulações.

Os resultados obtidos indicam que a média das estimativas simuladas, isto é, o estimador Monte Carlo de \(\hat{\lambda}_{\mathrm{MLE}}\) aproxima-se do valor verdadeiro do parâmetro, evidenciando a consistência do estimador no cenário analisado. Além disso, a variância amostral das estimativas apresenta-se muito próxima da variância assintótica teórica, dada por:

\[\begin{align}\\ \mathrm{Var}(\hat{\lambda}_{\text{MLE}}) \approx \dfrac{\lambda^2}{n} = \dfrac{4}{100} = 0{.}04\\\\ \end{align}\]

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

Definição 4.4 (Propriedades de \(\hat{\theta}\) via Simulação Monte Carlo - Givens e Hoeting 2012). Seja \(\hat{\theta}\) um estimador para o parâmetro populacional \(\theta\), construído a partir de uma amostra aleatória de tamanho \(n\), cujas observações são extraídas de uma distribuição \(f\). Suponha que são geradas \(R\) amostras independentes \(\{X_1^{(r)}, \ldots, X_n^{(r)}\}\), com \(r = 1, \ldots, R\), cada uma composta por \(n\) observações independentes e identicamente distribuídas segundo \(f\), e que o estimador \(\hat{\theta}\) é calculado em cada amostra, produzindo o conjunto de estimativas \(\{\hat{\theta}^{(1)}, \ldots, \hat{\theta}^{(R)}\}\). A partir desse conjunto, as seguintes propriedades inferenciais de \(\hat{\theta}\) podem ser estimadas numericamente via simulação Monte Carlo:

• Viés (B): \[\begin{align}\\ \widehat{B}(\hat{\theta}) = \bar{\theta}_{MC} - \theta = \frac{1}{R} \sum_{r=1}^R \hat{\theta}^{(r)} - \theta\\\\ \end{align}\] Quando \(\widehat{B}(\hat{\theta}) \approx 0\), diz-se que o estimador \(\hat{\theta}\) é aproximadamente não-viciado sob simulação Monte Carlo.

• Erro Quadrático Médio (EQM): \[\begin{align}\\ \widehat{\mathrm{EQM}}(\hat{\theta}) = \frac{1}{R} \sum_{r=1}^{R} \left( \hat{\theta}^{(r)} - \theta \right)^2 = \widehat{B}^2(\hat{\theta}) + \widehat{\mathrm{Var}}_{MC}(\hat{\theta})\\\\ \end{align}\] em que \(\widehat{\mathrm{Var}}_{MC}(\hat{\theta})\) é a variância estimada via simulação Monte Carlo, sendo definida por: \[\begin{align}\\ \widehat{\mathrm{Var}}_{MC}(\hat{\theta}) = \frac{1}{R - 1} \sum_{r=1}^{R} \left( \hat{\theta}^{(r)} - \bar{\theta}_{MC} \right)^2\\\\ \end{align}\]

• Probabilidade de Cobertura (CP): A proporção de intervalos de confiança de nível \(1 - \alpha\), construídos com base na aproximação normal, que contêm o valor verdadeiro de \(\theta\), é dada por: \[\begin{align}\\ \widehat{CP}(\theta) = \frac{1}{R} \sum_{r=1}^{R} \mathbb{I} \left\{ \hat{\theta}^{(r)} - z_{\alpha/2} \cdot \widehat{\mathrm{se}}^{(r)} < \theta < \hat{\theta}^{(r)} + z_{\alpha/2} \cdot \widehat{\mathrm{se}}^{(r)} \right\}\\\\ \end{align}\] em que \(\mathbb{I}\{\cdot\}\) é a função indicadora, \(z_{\alpha/2}\) é o quantil superior da distribuição normal padrão, e \(\widehat{\mathrm{se}}^{(r)}\) é o erro-padrão estimado do estimador \(\hat{\theta}^{(r)}\) calculado a partir da amostra \(r\).

• Comprimento Médio do Intervalo (CL): \[\begin{align}\\ \widehat{CL}(\theta) = \frac{1}{R} \sum_{r=1}^{R} 2 \cdot z_{\alpha/2} \cdot \widehat{\mathrm{se}}^{(r)}\\\\ \end{align}\] em que \(z_{\alpha/2}\) é o quantil superior da distribuição normal padrão, e \(\widehat{\mathrm{se}}^{(r)}\) é o erro-padrão estimado definido na probabilidade de cobertura (CP).

Exemplo 4.5 (Propriedades de \(\hat\alpha_{\text{MLE}}\) e \(\hat\beta_{\text{MLE}}\) da Distribuição Weibull via Simulação Monte Carlo). Seja \(X\) uma variável aleatória contínua tal que \(X\sim\) Weibull\((\alpha,\beta)\), cuja função densidade de probabilidade é dada por:

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

em que \(\alpha > 0\) é o parâmetro de forma (chamado de declividade de Weibull), e \(\beta > 0\) é o parâmetro de escala (medido na mesma unidade que \(X\); por exemplo, se \(X\) é uma variável de tempo, então \(\beta\) é interpretado como um tempo característico). Para uma amostra aleatória independente e identicamente distribuída, \(\mathbf{X} = (X_1, \ldots, X_n)\), a função de log-verossimilhança associada a essa distribuição é expressa por:

\[\begin{align}\\ \ell(\alpha, \beta; \mathbf{x}) = n \ln(\alpha) - n \alpha \ln(\beta) + (\alpha - 1) \sum_{i=1}^n \ln(x_i) - \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha\\\\ \end{align}\]

Para obter os estimadores de máxima verossimilhança, \(\hat\alpha_{\text{MLE}}\) e \(\hat\beta_{\text{MLE}}\), é necessário determinar, primeiramente, as derivadas de \(\ell(\alpha, \beta; \mathbf{x})\) com respeito a \(\alpha\) e \(\beta\), que compõe o vetor escore de Fisher, \(\mathbf{U}(\alpha, \beta; \mathbf{x})\), e resolver o sistema \(\mathbf{U}(\alpha, \beta; \mathbf{x}) = 0\). Para a distribuição Weibull, o vetor escore de Fisher é dado por:

\[\begin{align}\\ \mathbf{U}(\alpha, \beta; \mathbf{x}) = \begin{bmatrix} \displaystyle \frac{n}{\alpha} - n\ln (\beta) + \sum_{i=1}^n \ln(x_i) - \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha \ln\left( \frac{x_i}{\beta} \right) \\[10pt] \displaystyle -\frac{n\alpha}{\beta} + \frac{\alpha}{\beta} \sum_{i=1}^n \left( \frac{x_i}{\beta} \right)^\alpha \end{bmatrix}\\\\ \end{align}\]

Como não há forma analítica explícita para os estimadores de máxima verossimilhança \(\hat{\alpha}_{\text{MLE}}\) e \(\hat{\beta}_{\text{MLE}}\), sua obtenção deve ser realizada por meio de métodos numéricos, tais como o método de Newton-Raphson, simulated annealing ou algoritmos implementados nos otimizadores optim(), maxLik() e fitdist() do ambiente R. Assim, para estudar as propriedades desses estimadores, também se faz necessário o uso de abordagem numérica. Neste caso, pode-se trabalhar com simulação Monte Carlo, o que permite avaliar o comportamento dos estimadores sob diferentes condições amostrais. Para tanto, propõe-se o seguinte experimento:

• 1º Passo (Especificação dos Parâmetros da Simulação). Fixam-se o número de replicações \(R \in \mathbb{N}\) e o tamanho amostral \(n\). Estipulam-se ainda os valores verdadeiros dos parâmetros populacionais \(\alpha\) e \(\beta\), que serão utilizados na geração das amostras.

• 2º Passo (Geração das Replicações Amostrais). Para cada replicação \(r = 1, \ldots, R\), simula-se uma amostra independente e identicamente distribuída \(\left\{X_1^{(r)}, \ldots, X_n^{(r)}\right\}\), com \(X_i^{(r)} \sim\) Weibull\((\alpha,\beta)\).

• 3º Passo (Cálculo do Estimador em Cada Replicação). ara cada uma das \(R\) amostras geradas, obtêm-se numericamente os estimadores de máxima verossimilhança \(\hat{\alpha}_{\text{MLE}}^{(r)}\) e \(\hat{\beta}_{\text{MLE}}^{(r)}\) para os parâmetros \(\alpha\) e \(\beta\), respectivamente.

• 4º Passo (Avaliar as Propriedades dos MLEs). Com base nas \(R\) réplicas, estimam-se as propriedades inferenciais dos MLEs, tais como viés, erro quadrático médio, probabilidade de cobertura e largura média dos intervalos de confiança assintóticos.

A seguir, apresenta-se a implementação, em ambiente R, da estrutura computacional para os três cenários previamente delineados, voltada ao estudo das propriedades assintóticas dos estimadores de máxima verossimilhança (MLEs) \(\hat{\alpha}_{\mathrm{MLE}}\) e \(\hat{\beta}_{\mathrm{MLE}}\) da distribuição Weibull\((\alpha, \beta)\). A estimação será realizada por meio do otimizador fitdist(), do pacote fitdistrplus. O objetivo é analisar o desempenho desses estimadores considerando distintas combinações de valores paramétricos, tamanhos amostrais e quantidades de replicações da simulação, conforme especificado a seguir.

1º Cenário (Código 4.5):

• Número de Simulações \((R)\): 10, 100 e 1000.
• Tamanho Amostral \((n)\): 30.
• Valores Nominais (Verdadeiros) dos Parâmetros: \(\alpha = 1.5\) e \(\beta = 2.0\).
• Propriedades Avaliadas: B, EQM, CP, CL.

Código 4.5. Estudo das propriedades dos estimadores de máxima verossimilhança \(\hat{\alpha}_{\mathrm{MLE}}\) e \(\hat{\beta}_{\mathrm{MLE}}\) da distribuição Weibull\((\alpha, \beta)\) via simulação Monte Carlo, considerando o primeiro cenário experimental.

# -------------------------------------------------------------------------------------
# Propriedades dos MLEs da Distribuição Weibull via Simulação Monte Carlo - 1º Cenário
# -------------------------------------------------------------------------------------

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

library(fitdistrplus)
library(ggplot2)
library(tidyr)
library(dplyr)
library(knitr)

# --- 1. Criação da função para realizar as simulações e os cálculos das propriedades ---

avaliar_mle_weibull   <- function(R, n, par, sig.level = 0.05, seed = 123)
{
  # Semente (para reprodutibilidade)
  
  set.seed(seed)
  
  # Gerar dados: R amostras independentes de tamanho n da Weibull(α, β)
  
  x                   <- rweibull(n * R, shape = par[1], scale = par[2]) 
  dados               <- data.frame(matrix(sample(x), ncol = R)) 
  
  # Estimação dos parâmetros via fitdist()
  
  tempo_exec          <- system.time(mles <- lapply(dados, fitdist, distr = "weibull", 
                                                    start = list(shape = par[1], scale = par[2])))
  
  # Extração dos MLEs e dos erros-padrão (SE)
  
  mle                 <- do.call("rbind", lapply(mles, "[[", "estimate")) 
  mle.se              <- do.call("rbind", lapply(mles, "[[", "sd"))  
  
  # Cálculo das propriedades inferenciais para α
  
  vicio.alpha         <- mean(mle[,1] - par[1], na.rm = TRUE) 
  EQM.alpha           <- mean((mle[,1] - par[1])^2, na.rm = TRUE) 
  cp.alpha            <- mean((par[1] > (mle[,1] - qnorm(1 - sig.level/2) * mle.se[,1])) &
                                (par[1] < (mle[,1] + qnorm(1 - sig.level/2) * mle.se[,1])), na.rm = TRUE) 
  cl.alpha            <- mean(2 * qnorm(1 - sig.level/2) * mle.se[,1], na.rm = TRUE)
  
  # Cálculo das propriedades inferenciais para β
  
  vicio.beta          <- mean(mle[,2] - par[2], na.rm = TRUE) 
  EQM.beta            <- mean((mle[,2] - par[2])^2, na.rm = TRUE) 
  cp.beta             <- mean((par[2] > (mle[,2] - qnorm(1 - sig.level/2) * mle.se[,2])) &
                                (par[2] < (mle[,2] + qnorm(1 - sig.level/2) * mle.se[,2])), na.rm = TRUE) 
  cl.beta             <- mean(2 * qnorm(1 - sig.level/2) * mle.se[,2], na.rm = TRUE)
  
  # Organização dos resultados
  
  res.alpha           <- round(c(vicio.alpha, 
                                 EQM.alpha, 
                                 cp.alpha, 
                                 cl.alpha, 
                                 tempo_exec[3]), 4)
  
  res.beta            <- round(c(vicio.beta, 
                                 EQM.beta, 
                                 cp.beta, 
                                 cl.beta, 
                                 tempo_exec[3]), 4)
  
  res.final           <- data.frame(alpha = res.alpha, beta = res.beta)
  rownames(res.final) <- c("B", "EQM", "CP", "CL", "Tempo de Execução (s)")
  return(res.final)
}

# --- 2. Execução das simulações Monte Carlo ---

param                 <- c(1.5, 2.0) # Parâmetros da distribuição Weibull
n_amostra             <- 30          # Tamanho amostral

resultado_10          <- avaliar_mle_weibull(R = 10  , n = n_amostra, par = param, seed = 1212)
resultado_100         <- avaliar_mle_weibull(R = 100 , n = n_amostra, par = param, seed = 1212)
resultado_1000        <- avaliar_mle_weibull(R = 1000, n = n_amostra, par = param, seed = 1212)

# --- 3. Organização dos resultados ---

res_f                 <- data.frame(R_10_alpha = resultado_10[,1],
                                    R_100_alpha = resultado_100[,1],
                                    R_1000_alpha = resultado_1000[,1],
                                    R_10_beta = resultado_10[,2],
                                    R_100_beta = resultado_100[,2],
                                    R_1000_beta = resultado_1000[,2],
                                    row.names = c('Viés', 
                                                  'Erro Quadrático Médio', 
                                                  'Probabilidade de Cobertura', 
                                                  'Comprimento Médio do Intervalo', 
                                                  'Tempo de Execução (s)'))

kable(res_f, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('R = 10 (α)', 'R = 100 (α)', 'R = 1000 (α)', 
                    'R = 10 (β)', 'R = 100 (β)', 'R = 1000 (β)'))

	R = 10 (α)	R = 100 (α)	R = 1000 (α)	R = 10 (β)	R = 100 (β)	R = 1000 (β)
Viés	0.1785	0.1169	0.0750	0.1758	0.0520	0.0144
Erro Quadrático Médio	0.1315	0.0779	0.0588	0.0724	0.0761	0.0638
Probabilidade de Cobertura	0.9000	0.9300	0.9530	0.9000	0.9000	0.9430
Comprimento Médio do Intervalo	0.9404	0.9091	0.8859	0.9989	0.9724	0.9794
Tempo de Execução (s)	0.0170	0.1680	1.5070	0.0170	0.1680	1.5070

# --- 4. Geração dos gráficos ---

# Organização dos dados

lista_resultados      <- list(`R = 10`   = resultado_10,
                              `R = 100`  = resultado_100,
                              `R = 1000` = resultado_1000)

df                    <- bind_rows(lapply(names(lista_resultados),
                                          function(nome){
                                              res <- lista_resultados[[nome]]
                                              res$Metrica <- rownames(res)
                                              res$R <- nome
                                              res}), .id = NULL)

df_long               <- df %>%
                          pivot_longer(cols = c("alpha", "beta"),
                                       names_to = "Parametro",
                                       values_to = "Valor")

df_long$R             <- factor(df_long$R, levels = c("R = 10", "R = 100", "R = 1000"))

# Definir siglas para as propriedades

siglas                <- c("B" = "B",
                           "EQM" = "EQM",
                           "CP" = "CP",
                           "CL" = "CL",
                           "Tempo de Execução (s)" = "Tempo")

df_long$Metrica         <- factor(df_long$Metrica, levels = names(siglas))
levels(df_long$Metrica) <- siglas[levels(df_long$Metrica)]

# Filtrar para gráficos (excluir tempo de execução)

df_plot                 <- df_long %>% filter(Metrica != "Tempo")

# Função auxiliar: gera o gráfico das propriedades por parâmetro

plot_parametro          <- function(param, titulo_grafico) 
{
  df_sub                <- df_plot %>% filter(Parametro == param)
  linhas_ref            <- data.frame(Metrica = factor(c("B", "EQM", "CP", "CL"), 
                                                       levels = levels(df_sub$Metrica)),
                                      yintercept = c(0, 0, 0.95, 0))
  cores_metricas        <- c(B = "#1b9e77", EQM = "#d95f02", CP = "#7570b3", CL = "#e7298a")
  
  ggplot(df_sub, aes(x = R, y = Valor, group = 1, color = Metrica)) +
    geom_point(size = 3) +
    geom_line(linewidth = 0.5) +
    facet_wrap(~Metrica, scales = "free_y", strip.position = "top", ncol = 2, nrow = 2,
               labeller = as_labeller(c(B = "Viés (B)",
                                        EQM = "Erro Quadrático Médio (EQM)",
                                        CP = "Probabilidade de Cobertura (CP)",
                                        CL = "Comprimento Médio do Intervalo (CL)"))) +
    geom_hline(data = linhas_ref, aes(yintercept = yintercept), color = "red", 
               linetype = "dashed", linewidth = 0.7) +
    scale_color_manual(values = cores_metricas) +
    labs(title = titulo_grafico,
         x = "Número de simulações (R)",
         y = "Valores") +
    theme_minimal() +
    theme(strip.background = element_blank(),
          strip.text = element_text(face = "bold", size = 11, margin = margin(b = 8)),
          axis.text.x = element_text(angle = 45, hjust = 1),
          panel.spacing = unit(1.2, "lines"),
          axis.title.y = element_text(margin = margin(r = 15)),
          axis.title.x = element_text(margin = margin(t = 15)),
          legend.position = "none")
}

# --- 5. Ilustração gráfica para o parâmetro α ---

plot_parametro("alpha", '')

Figura 4.5.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\alpha\) da distribuição Weibull, obtidas via simulação Monte Carlo para diferentes números de replicações \((R)\).

# --- 6. Ilustração gráfica para o parâmetro β ---

plot_parametro("beta", '')

Figura 4.6.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\beta\) da distribuição Weibull, obtidas via simulação Monte Carlo para diferentes números de replicações \((R)\).

2º Cenário (Código 4.6):

• Número de Simulações \((R)\): 100.
• Tamanhos Amostrais \((n)\): 10, 20, 30, 40, 50, 60, 70, 80, 90.
• Valores Nominais (Verdadeiros) dos Parâmetros: \(\alpha = 1.5\) e \(\beta = 2.0\).
• Propriedades Avaliadas: B, EQM, CP, CL.

Código 4.6. Estudo das propriedades dos estimadores de máxima verossimilhança \(\hat{\alpha}_{\mathrm{MLE}}\) e \(\hat{\beta}_{\mathrm{MLE}}\) da distribuição Weibull\((\alpha, \beta)\) via simulação Monte Carlo, considerando o segundo cenário experimental.

# -------------------------------------------------------------------------------------
# Propriedades dos MLEs da Distribuição Weibull via Simulação Monte Carlo - 2º Cenário
# -------------------------------------------------------------------------------------

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

library(fitdistrplus)
library(ggplot2)
library(tidyr)
library(dplyr)
library(knitr)

# --- 1. Criação da função para realizar as simulações e os cálculos das propriedades ---

avaliar_mle_weibull   <- function(R, n, par, sig.level = 0.05, seed = 123)
{
  # Semente (para reprodutibilidade)
  
  set.seed(seed)
  
  # Gerar dados: N amostras independentes de tamanho n da Weibull(α, β)
  
  x                   <- rweibull(n * R, shape = par[1], scale = par[2]) 
  dados               <- data.frame(matrix(sample(x), ncol = R)) 
  
  # Estimação dos parâmetros via fitdist()
  
  tempo_exec          <- system.time(mles <- lapply(dados, fitdist, distr = "weibull", 
                                                    start = list(shape = par[1], scale = par[2])))
  
  # Extração dos MLEs e dos erros-padrão (SE)
  
  mle                 <- do.call("rbind", lapply(mles, "[[", "estimate")) 
  mle.se              <- do.call("rbind", lapply(mles, "[[", "sd"))  
  
  # Cálculo das propriedades inferenciais para α
  
  vicio.alpha         <- mean(mle[,1] - par[1], na.rm = TRUE) 
  EQM.alpha           <- mean((mle[,1] - par[1])^2, na.rm = TRUE) 
  cp.alpha            <- mean((par[1] > (mle[,1] - qnorm(1 - sig.level/2) * mle.se[,1])) &
                                (par[1] < (mle[,1] + qnorm(1 - sig.level/2) * mle.se[,1])), na.rm = TRUE) 
  cl.alpha            <- mean(2 * qnorm(1 - sig.level/2) * mle.se[,1], na.rm = TRUE)
  
  # Cálculo das propriedades inferenciais para β
  
  vicio.beta          <- mean(mle[,2] - par[2], na.rm = TRUE) 
  EQM.beta            <- mean((mle[,2] - par[2])^2, na.rm = TRUE) 
  cp.beta             <- mean((par[2] > (mle[,2] - qnorm(1 - sig.level/2) * mle.se[,2])) &
                                (par[2] < (mle[,2] + qnorm(1 - sig.level/2) * mle.se[,2])), na.rm = TRUE) 
  cl.beta             <- mean(2 * qnorm(1 - sig.level/2) * mle.se[,2], na.rm = TRUE)
  
  # Organização dos resultados
  
  res.alpha           <- round(c(vicio.alpha, 
                                 EQM.alpha, 
                                 cp.alpha, 
                                 cl.alpha, 
                                 tempo_exec[3]), 4)
  
  res.beta            <- round(c(vicio.beta, 
                                 EQM.beta, 
                                 cp.beta, 
                                 cl.beta, 
                                 tempo_exec[3]), 4)
  
  res.final           <- data.frame(alpha = res.alpha, beta = res.beta)
  rownames(res.final) <- c("B", "EQM", "CP", "CL", "Tempo de Execução (s)")
  return(res.final)
}

# --- 2. Execução das simulações Monte Carlo ---

param                 <- c(1.5, 2)              # Parâmetros da distribuição Weibull
R_simu                <- 100                    # Número de simulações (fixo)
amostras_n            <- seq(10, 90, by = 10)   # Tamanhos amostrais

resultados            <- vector("list", length(amostras_n))
names(resultados)     <- paste0("n = ", amostras_n)

for (i in seq_along(amostras_n))
{
  resultados[[i]]     <- avaliar_mle_weibull(R = R_simu, 
                                             n = amostras_n[i], 
                                             par = param, 
                                             seed = 1212)
}

# --- 3. Organização dos resultados ---

# Parâmetro α

res_alpha             <- data.frame(n_10_alpha = resultados[[1]][,1],
                                    n_20_alpha = resultados[[2]][,1],
                                    n_30_alpha = resultados[[3]][,1],
                                    n_40_alpha = resultados[[4]][,1],
                                    n_50_alpha = resultados[[5]][,1],
                                    n_60_alpha = resultados[[6]][,1],
                                    n_70_alpha = resultados[[7]][,1],
                                    n_80_alpha = resultados[[8]][,1],
                                    n_90_alpha = resultados[[9]][,1],
                                    row.names = c('Viés', 
                                                  'Erro Quadrático Médio', 
                                                  'Probabilidade de Cobertura', 
                                                  'Comprimento Médio do Intervalo', 
                                                  'Tempo de Execução (s)'))
kable(res_alpha, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('n = 10', 'n = 20', 'n = 30', 
                    'n = 40', 'n = 50', 'n = 60',
                    'n = 70', 'n = 80', 'n = 90'))

	n = 10	n = 20	n = 30	n = 40	n = 50	n = 60	n = 70	n = 80	n = 90
Viés	0.3241	0.1964	0.1169	0.0797	0.0603	0.0325	0.0327	0.0225	0.0152
Erro Quadrático Médio	0.3526	0.1700	0.0779	0.0503	0.0270	0.0223	0.0200	0.0179	0.0147
Probabilidade de Cobertura	0.9800	0.8900	0.9300	0.9400	0.9900	0.9700	0.9800	0.9500	0.9600
Comprimento Médio do Intervalo	1.8008	1.1807	0.9091	0.7693	0.6775	0.6073	0.5623	0.5223	0.4889
Tempo de Execução (s)	0.1580	0.1510	0.1580	0.1560	0.1670	0.1640	0.1600	0.1630	0.1750

# Parâmetro β

res_beta              <- data.frame(n_10_beta = resultados[[1]][,2],
                                    n_20_beta = resultados[[2]][,2],
                                    n_30_beta = resultados[[3]][,2],
                                    n_40_beta = resultados[[4]][,2],
                                    n_50_beta = resultados[[5]][,2],
                                    n_60_beta = resultados[[6]][,2],
                                    n_70_beta = resultados[[7]][,2],
                                    n_80_beta = resultados[[8]][,2],
                                    n_90_beta = resultados[[9]][,2],
                                    row.names = c('Viés', 
                                                  'Erro Quadrático Médio', 
                                                  'Probabilidade de Cobertura', 
                                                  'Comprimento Médio do Intervalo', 
                                                  'Tempo de Execução (s)'))

kable(res_beta, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('n = 10', 'n = 20', 'n = 30', 
                    'n = 40', 'n = 50', 'n = 60',
                    'n = 70', 'n = 80', 'n = 90'))

	n = 10	n = 20	n = 30	n = 40	n = 50	n = 60	n = 70	n = 80	n = 90
Viés	0.1028	0.0493	0.0520	0.0336	0.0151	0.0107	0.0149	0.0110	0.0039
Erro Quadrático Médio	0.2054	0.0994	0.0761	0.0554	0.0396	0.0263	0.0261	0.0208	0.0202
Probabilidade de Cobertura	0.8800	0.9300	0.9000	0.9400	0.9400	0.9700	0.9400	0.9500	0.9300
Comprimento Médio do Intervalo	1.5854	1.1504	0.9724	0.8504	0.7585	0.7034	0.6529	0.6128	0.5780
Tempo de Execução (s)	0.1580	0.1510	0.1580	0.1560	0.1670	0.1640	0.1600	0.1630	0.1750

# --- 4. Geração dos gráficos ---

# Organização dos dados

df                    <- bind_rows(lapply(names(resultados),
                                          function(nome){
                                              res <- resultados[[nome]]
                                              res$Metrica <- rownames(res)
                                              res$Tamanho <- nome
                                              res}), .id = NULL)

df_long               <- df %>%
                          pivot_longer(cols = c("alpha", "beta"),
                                       names_to = "Parametro",
                                       values_to = "Valor")

df_long$Tamanho       <- factor(df_long$Tamanho, levels = paste0("n = ", amostras_n))

# Definir siglas para as propriedades

siglas                <- c("B" = "B",
                           "EQM" = "EQM",
                           "CP" = "CP",
                           "CL" = "CL",
                           "Tempo de Execução (s)" = "Tempo")

df_long$Metrica         <- factor(df_long$Metrica, levels = names(siglas))
levels(df_long$Metrica) <- siglas[levels(df_long$Metrica)]

# Filtrar para gráficos (excluir tempo de execução)

df_plot                 <- df_long %>% filter(Metrica != "Tempo")

# Função auxiliar: gera o gráfico das propriedades por parâmetro

plot_parametro          <- function(param, titulo_grafico) 
{
  df_sub                <- df_plot %>% filter(Parametro == param)
  linhas_ref            <- data.frame(Metrica = factor(c("B", "EQM", "CP", "CL"), 
                                                       levels = levels(df_sub$Metrica)),
                                      yintercept = c(0, 0, 0.95, 0))
  cores_metricas        <- c(B = "#1b9e77", EQM = "#d95f02", CP = "#7570b3", CL = "#e7298a")
  
  ggplot(df_sub, aes(x = Tamanho, y = Valor, group = 1, color = Metrica)) +
    geom_point(size = 3) +
    geom_line(linewidth = 0.5) +
    facet_wrap(~Metrica, scales = "free_y", strip.position = "top", ncol = 2, nrow = 2,
               labeller = as_labeller(c(B = "Viés (B)",
                                        EQM = "Erro Quadrático Médio (EQM)",
                                        CP = "Probabilidade de Cobertura (CP)",
                                        CL = "Comprimento Médio do Intervalo (CL)"))) +
    geom_hline(data = linhas_ref, aes(yintercept = yintercept), color = "red", 
               linetype = "dashed", linewidth = 0.7) +
    scale_color_manual(values = cores_metricas) +
    labs(title = titulo_grafico,
         x = "Tamanho Amostral (n)",
         y = "valores") +
    theme_minimal() +
    theme(strip.background = element_blank(),
          strip.text = element_text(face = "bold", size = 11, margin = margin(b = 8)),
          axis.text.x = element_text(angle = 45, hjust = 1),
          panel.spacing = unit(1.2, "lines"),
          axis.title.y = element_text(margin = margin(r = 15)),
          axis.title.x = element_text(margin = margin(t = 15)),
          legend.position = "none")
}

# --- 5. Ilustração gráfica para o parâmetro α ---

plot_parametro("alpha", '')

Figura 4.7.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\alpha\) da distribuição Weibull, obtidas via simulação Monte Carlo para diferentes tamanhos amostrais com \(R = 100\) replicações.

# --- 6. Ilustração gráfica para o parâmetro β ---

plot_parametro("beta", '')

Figura 4.8.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\beta\) da distribuição Weibull, obtidas via simulação Monte Carlo para diferentes tamanhos amostrais com \(R = 100\) replicações.

3º Cenário (Código 4.7):

• Número de Simulações \((R)\): 1000.
• Tamanhos Amostrais \((n)\): 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.
• Valores Nominais (Verdadeiros) dos Parâmetros: \(\alpha \in \{0{.}5,\ 1{.}0,\ 1{.}5\}\) e \(\lambda \in \{0{.}5,\ 1{.}0,\ 1{.}5\}\).
• Propriedades Avaliadas: B, EQM, CP, CL.

Código 4.7. Estudo das propriedades dos estimadores de máxima verossimilhança \(\hat{\alpha}_{\mathrm{MLE}}\) e \(\hat{\beta}_{\mathrm{MLE}}\) da distribuição Weibull\((\alpha, \beta)\) via simulação Monte Carlo, considerando o terceiro cenário experimental.

# -------------------------------------------------------------------------------------
# Propriedades dos MLEs da Distribuição Weibull via Simulação Monte Carlo - 3º Cenário
# -------------------------------------------------------------------------------------

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

library(fitdistrplus)
library(ggplot2)
library(tidyr)
library(dplyr)
library(knitr)

# --- 1. Função para o cálculo do MLE da distribuição Weibull ---

mle.weibull                   <- function(x, par)
{
  p                           <- as.numeric(par)
  fit                         <- fitdist(x, 
                                         'weibull', 
                                         start = list(shape = p[1], scale = p[2]))
  mles                        <- fit$estimate
  ses                         <- fit$sd
  if(!is.numeric(mles)) mles  <- NA
  if(!is.numeric(ses)) ses    <- NA
  return(rbind('MLEs' = mles, 'SEs' = ses))
}

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

set.seed(1212) 

shapes                        <- seq(0.5, 1.5, 0.5) # shape = 𝛼
scales                        <- seq(0.5, 1.5, 0.5) # scale = β
param                         <- expand.grid(alpha = shapes, 
                                             beta = scales) # Combinações dos parâmetros
n_param                       <- nrow(param)

R                             <- 1000                 # Número de simulações
nmin                          <- 10                   # Menor tamanho amostral
nmax                          <- 100                  # Maior tamanho amostral
enes                          <- seq(nmin, nmax, 10)  # Tamanhos amostrais
n_enes                        <- length(enes)

# --- 3. Criação das matrizes para armazenar os resultados ---

v.shape                       <- matrix(nrow = n_enes, ncol = n_param)
eqm.shape                     <- matrix(nrow = n_enes, ncol = n_param)
cp.shape                      <- matrix(nrow = n_enes, ncol = n_param)
cl.shape                      <- matrix(nrow = n_enes, ncol = n_param)

v.scale                       <- matrix(nrow = n_enes, ncol = n_param)
eqm.scale                     <- matrix(nrow = n_enes, ncol = n_param)
cp.scale                      <- matrix(nrow = n_enes, ncol = n_param)
cl.scale                      <- matrix(nrow = n_enes, ncol = n_param)

time_mat                      <- matrix(nrow = n_enes, ncol = n_param)

# --- 4. Execução das simulações Monte Carlo ---

for(i in 1:n_param)
{
  X                           <- rweibull(nmax * R, 
                                          shape = param$alpha[i], 
                                          scale = param$beta[i])
  X                           <- matrix(X, nrow = nmax, ncol = R)
  
  for(k in 1:n_enes)
  {
    n                         <- enes[k]
    x                         <- data.frame(X[1:n,])
    
    tempo                     <- system.time({fit <- sapply(x, 
                                                            mle.weibull, 
                                                            par = c(param$alpha[i], 
                                                                    param$beta[i]))})["elapsed"]
    
    # Parâmetro 𝛼 (shape)
    
    v.shape[k, i]             <- mean(fit[1,] - param$alpha[i], na.rm = TRUE)
    eqm.shape[k, i]           <- mean((fit[1,] - param$alpha[i])^2, na.rm = TRUE)
    cp.shape[k, i]            <- mean(((fit[1,] - qnorm(1 - 0.05/2) * fit[2,]) < param$alpha[i]) &
                                        ((fit[1,] + qnorm(1 - 0.05/2) * fit[2,]) > param$alpha[i]), na.rm = TRUE)
    cl.shape[k, i]            <- mean(2 * qnorm(1 - 0.05/2) * fit[2,], na.rm = TRUE)
    
    # Parâmetro β (scale)
    
    v.scale[k, i]             <- mean(fit[3,] - param$beta[i], na.rm = TRUE)
    eqm.scale[k, i]           <- mean((fit[3,] - param$beta[i])^2, na.rm = TRUE)
    cp.scale[k, i]            <- mean(((fit[3,] - qnorm(1 - 0.05/2) * fit[4,]) < param$beta[i]) &
                                        ((fit[3] + qnorm(1 - 0.05/2) * fit[4,]) > param$beta[i]), na.rm = TRUE)
    cl.scale[k, i]            <- mean(2 * qnorm(1 - 0.05/2) * fit[4,], na.rm = TRUE)
    
    # Tempo computacional
    
    time_mat[k, i]            <- tempo
  }
}

# --- 5. Prepração dos dados para os gráficos ---

mat_to_df                     <- function(mat, medida_name, enes, n_param)
{
  df                          <- as.data.frame(mat)
  df$n                        <- enes
  df_long                     <- pivot_longer(df, cols = -n, names_to = "parametro", 
                                              values_to = "valor")
  df_long$medida              <- medida_name
  return(df_long)
}

# --- 6. Função de geração dos gráficos ---

df_plot                       <- function(v_mat, eqm_mat, cp_mat, cl_mat, enes, param, titulo)
{
  n_param                     <- nrow(param)
  
  df_all                      <- bind_rows(mat_to_df(v_mat, "B", enes, n_param),
                                           mat_to_df(eqm_mat, "EQM", enes, n_param),
                                           mat_to_df(cp_mat, "CP", enes, n_param),
                                           mat_to_df(cl_mat, "CL", enes, n_param))
  
  df_all$parametro            <- factor(df_all$parametro, levels = paste0("V", 1:n_param))
  levels(df_all$parametro)    <- paste0("param_", 1:n_param)
  
  df_all                      <- df_all %>%
                                          left_join(data.frame( parametro = paste0("param_", 1:n_param),
                                                                alpha = param$alpha,
                                                                beta = param$beta,
                                                                param_num = 1:n_param),
                                                    by = "parametro") %>%
                                                                        mutate(param_label = paste0("α = ", alpha, ", β = ", beta))
  
  ref_lines                   <- data.frame(medida = c("B", "EQM", "CP", "CL"),
                                            yintercept = c(0, 0, 0.95, 0))
  
  x_labels                    <- data.frame(n = enes, label = paste0("n = ", enes))
  
  p                           <- ggplot(df_all, aes(x = n, y = valor, group = param_num, color = param_label)) +
                                        geom_line() +
                                        geom_point() +
                                        geom_hline(data = ref_lines, aes(yintercept = yintercept),
                                                   linetype = "dashed", color = "red", size = 0.5) +
                                        geom_text(data = x_labels, aes(x = n, y = -Inf, label = label),
                                                  inherit.aes = FALSE, angle = 0, hjust = 0.5, vjust = 1.2, size = 3, 
                                                  color = "black") +
                                        facet_wrap(~ medida, nrow = 2, ncol = 2, scales = "free_y",
                                                   labeller = as_labeller(c(B = "Viés (B)",
                                                                            EQM = "Erro Quadrático Médio (EQM)",
                                                                            CP = "Probabilidade de Cobertura (CP)",
                                                                            CL = "Comprimento Médio do Intervalo (CL)"))) +
                                        labs(title = titulo,
                                             x = "Tamanho Amostral (n)",
                                             y = 'Valores',
                                             color = "Cenários") +
                                        scale_x_continuous(limits = c(min(enes), max(enes)), 
                                                           expand = expansion(mult = c(0.05, 0.05))) +
                                        theme_minimal() +
                                        theme(legend.position = "bottom",
                                              legend.title = element_text(size = 11, face = "bold"),
                                              legend.text = element_text(size = 11),
                                              strip.text = element_text(face = "bold", size = 11),
                                              plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
                                              axis.text.x = element_text(angle = 0, hjust = 0.5),
                                              axis.title.y = element_text(margin = margin(r = 15)),
                                              axis.title.x = element_text(margin = margin(t = 15))) +
                                        guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
                                        scale_color_discrete()
    return(p)
}

# --- 7. Ilustração gráfica para o parâmetro α ---

df_plot(v.shape, eqm.shape, cp.shape, cl.shape, enes, param, '')

Figura 4.9.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\alpha\) da distribuição Weibull, obtidas via simulação Monte Carlo para diferentes tamanhos amostrais e valores paramétricos com \(R = 1000\) replicações.

# --- 8. Ilustração gráfica para o parâmetro β ---

df_plot(v.scale, eqm.scale, cp.scale, cl.scale, enes, param, '')

Figura 4.10.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\beta\) da distribuição Weibull, obtidas via simulação Monte Carlo para diferentes tamanhos amostrais e valores paramétricos com \(R = 1000\) replicações.

# --- 9. Ilustração gráfica do tempo de execução das simulações Monte Carlo em cada cenário ---

time_df                       <- as.data.frame(time_mat)
time_df$n                     <- enes
time_long                     <- pivot_longer(time_df, cols = -n, names_to = "parametro", 
                                              values_to = "tempo")
time_long$parametro           <- factor(time_long$parametro, levels = paste0("V", 1:n_param))
levels(time_long$parametro)   <- paste0("param_", 1:n_param)

time_long                     <- time_long %>%
                                             left_join(data.frame(parametro = paste0("param_", 1:n_param),
                                                                  alpha = param$alpha,
                                                                  beta = param$beta,
                                                                  param_num = 1:n_param),
                                                       by = "parametro") %>%
                                                                     mutate(param_label = paste0("α = ", alpha, ", β = ", beta))


ggplot(time_long, aes(x = n, y = tempo, group = param_num, color = param_label)) +
  geom_line() +
  geom_point() +
  geom_text(data = data.frame(n = enes, label = paste0("n = ", enes)),
            aes(x = n, y = -Inf, label = label),
            inherit.aes = FALSE, angle = 0, hjust = 0.5, vjust = 1.2, size = 3, color = "black") +
  labs(title = "",
       x = "Tamanho Amostral (n)",
       y = "Tempo (segundos)",
       color = "Cenários") +
  scale_x_continuous(limits = c(min(enes), max(enes)), expand = expansion(mult = c(0.05, 0.05))) +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 11, face = "bold"),
        legend.text = element_text(size = 11),
        plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
        axis.text.x = element_text(angle = 0, hjust = 0.5),
        axis.title.y = element_text(margin = margin(r = 15)),
        axis.title.x = element_text(margin = margin(t = 15))) +
  guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
  scale_color_discrete()

Figura 4.11.. Gráfico do tempo de execução (em segundos) em função do tamanho amostral (n), estratificado por cenários definidos por diferentes valores dos parâmetros da distribuição Weibull com \(R = 1000\) replicações.

Exemplo 4.6 (Propriedades de \(\hat\beta_{\text{MLE}}\) da Distribuição Lindley via Simulação Monte Carlo) A distribuição de Lindley foi originalmente introduzida por Lindley (1958) no contexto de distribuições fiduciais. Durante várias décadas, sua utilização restringiu-se, sobretudo, a processos de composição em conjunto com a distribuição de Poisson. Contudo, nas últimas décadas, essa distribuição tem despertado crescente interesse, notadamente nas áreas de análise de sobrevivência e confiabilidade, onde se apresenta como uma alternativa competitiva à distribuição Weibull. De modo geral, diz-se que uma variável aleatória contínua \(X\) segue uma distribuição de Lindley com parâmetro \(\beta > 0\) se sua função densidade de probabilidade é expressa por:

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

em que \(\beta > 0\) representa o parâmetro de forma. No ambiente R, diferentemente do que ocorre com a distribuição Weibull, a distribuição de Lindley não possui funções básicas pré-implementadas. Assim, para trabalhar com essa distribuição, torna-se necessário implementar manualmente tais funções. Então, para iniciar este processo, apresenta-se a função densidade de probabilidade, a qual será implementada de acordo com o Código 4.8. Por definição e padronização no ambiente R, tal função deve necessariamente conter três argumentos:

• x: vetor numérico correspondente aos possíveis valores assumidos pela variável aleatória \(X\).
• beta: parâmetro de forma da distribuição Lindley.
• log: argumento lógico; caso definido como TRUE, a função retorna os valores das probabilidades no formato logarítmico.

Código 4.8. Implementação da função densidade da distribuição Lindley com parâmetro de forma \(\beta\).

# ----------------------------------------------------------
# Função Densidade de Probabilidade da Distribuição Lindley
# ----------------------------------------------------------

dlindley <- function(x, beta, log = FALSE) 
{
    stopifnot(beta > 0)
    if(log) 
    {
        t1 <- log(beta)
        t4 <- log1p(beta)
        t6 <- log1p(x)
        pdf <- -beta * x + 2 * t1 - t4 + t6
        return(pdf)
    }
    else 
    {
        t1 <- beta^2
        t7 <- exp(-beta * x)
        pdf <- t1/(1 + beta) * (1 + x) * t7
        return(pdf)
    }
}

Por outro lado, tem-se que a função de distribuição acumulada da distribuição Lindley é descrita por:

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

O Código 4.9 ilustra a implementação dessa função no ambiente R, com base nos seguintes argumentos:

• q: vetor numérico correspondente aos valores dos quantis nos quais se deseja avaliar a função de distribuição acumulada.
• beta: parâmetro de forma da distribuição Lindley.
• log.p: valor lógico; se TRUE, retorna o logaritmo da função de distribuição acumulada.
• lower.tail: valor lógico; se TRUE (padrão), retorna \(F(x)\), isto é, \(\mathrm{P}(X \leqslant x)\); caso contrário, retorna \(\mathrm{P}(X > x)\)

Código 4.9. Implementação da função de distribuição acumulada da distribuição de Lindley com parâmetro de forma \(\beta\).

# ---------------------------------------------------------
# Função de Distribuição Acumulada da Distribuição Lindley
# ---------------------------------------------------------

plindley <- function(q, beta, lower.tail = TRUE, log.p = FALSE) 
{
    stopifnot(beta > 0)
    if (lower.tail)
    {
        t1 <- beta * q
        t6 <- exp(-t1)
        cdf <- 1 - (1 + t1/(1 + beta)) * t6
    }
    else
    {
        t1 <- beta * q
        t6 <- exp(-t1)
        cdf <- (1 + t1/(1 + beta)) * t6
    }
    if (log.p) 
    {
        return(log(cdf))
    }
    else
    {
        return(cdf)
    } 
}

No que tange à função quantil, é necessário adotar maior cautela, dado que nem todas as distribuições admitem expressão analítica fechada para essa função. Contudo, no caso da distribuição de Lindley, a função quantil possui forma analítica, a qual é definida por:

\[\begin{align}\\ Q(p;\beta) = -1 - \dfrac{1}{\beta} - \dfrac{1}{\beta}W_{-1}\left[(1 + \beta)(p - 1) e^{-(1 + \beta)}\right]\\\\ \end{align}\]

em que \(W_{-1}(\cdot)\) representa o ramo negativo da função W de Lambert (Lambert 1758; Jodrá 2010) , a qual é uma função multivalorada, porém pode ser definida de forma unívoca como função real no intervalo \((-1/e, \infty)\). Assim, para implementar a função quantil da distribuição de Lindley no ambiente R, faz-se necessário utilizar a função lambertWm1() do pacote lamW, que retorna o valor da função W de Lambert correspondente ao ramo negativo. O Código 4.10 ilustra a implementação dessa função, com base nos seguintes argumentos:

• p: vetor numérico de probabilidades, com valores no intervalo \([0,1]\).
• beta: parâmetro de forma da distribuição Lindley.
• log.p: valor lógico; se TRUE, retorna o logaritmo da função de distribuição acumulada.
• lower.tail: valor lógico; se TRUE (padrão), as probabilidades referem-se a \(\mathrm{P}(X \leqslant x)\); caso contrário, a \(\mathrm{P}(X > x)\)

Código 4.10. Implementação da função quantil da distribuição de Lindley com parâmetro de forma \(\beta\).

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

qlindley <- function(p, beta, lower.tail = TRUE, log.p = FALSE) 
{
    stopifnot(beta > 0)
    if (lower.tail)
    {
        t1 <- 1 + beta
        t4 <- exp(-t1)
        t6 <- lambertWm1(t1 * (p - 1) * t4)
        qtf <- -(t6 + 1 + beta)/beta
    }
    else 
    {
        t1 <- 1 + beta
        t3 <- exp(-t1)
        t5 <- lambertWm1(-p * t1 * t3)
        qtf <- -(t5 + 1 + beta)/beta
    }
    if (log.p)
    {
        return(log(qtf))
    }
    else
    {   
        return(qtf)
    }
}

Por fim, resta implementar a função geradora de valores pseudoaleatórios segundo a distribuição de Lindley. Conforme Ghitany, Atieh, e Nadarajah (2008), existem duas formas principais para gerar amostras a partir dessa distribuição: o algoritmo de mistura e o método da transformação inversa. Ambos os algoritmos são descritos a seguir:

• Algoritmo I (Forma de Mistura da Distribuição Lindley):
  • Gerar \(U_i \sim Uniforme(0,1), i = 1, \ldots, n;\)
  • Gerar \(V_i \sim Exponencial(\beta), i = 1, \ldots, n;\)
  • Gerar \(W_i \sim Gama(2,\beta), i = 1, \ldots, n;\)
  • Se \(U_i \leqslant \dfrac{\beta}{\beta + 1}\), então defina \(X_i = V_i\), caso contrário, defina \(X_i = W_i, i = 1, \ldots, n.\)

• Algoritmo II (Transformação Inversa):
  • Gerar \(U_i \sim Uniforme(0,1), i = 1, \ldots, n;\)
  • Definir: \[\begin{align}\\ X_i = -1 - \dfrac{1}{\beta} - \dfrac{1}{\beta}W_{-1}\left[(1 + \beta)(p - 1) e^{-(1 + \beta)}\right]\\\\ \end{align}\] em que \(W_{-1}(\cdot)\) é o ramo negativo da função W de Lambert.

Com base nestes algoritmos, o Código 4.11 ilutra a implementação da função geradora de valores pseudoaleatórios da distribuição de Lindley, que recebe três argumentos básicos:

• n: número de valores pseudoaleatórios a gerem gerados.
• beta: parâmetro de forma da distribuição Lindley.
• mixture: argumento lógico; se TRUE, utiliza-se o Algoritmo I; caso contrário, emprega-se o Algoritmo II.

Código 4.11. Implementação da função geradora de valores pseudoaleatórios da distribuição de Lindley com parâmetro de forma \(\beta\).

# --------------------------------------------------------------------
# Função Geradora de Valores Pseudoaleatórios da Distribuição Lindley
# --------------------------------------------------------------------

rlindley <- function(n, beta, mixture = TRUE) 
{
    stopifnot(beta > 0)
    if (mixture) 
    {
        x  <- rbinom(n, size = 1, prob = beta/(1 + beta))
        rv <- x * rgamma(n, shape = 1, rate = beta) + (1 - x) * rgamma(n, shape = 2, rate = beta)
        return(rv)
    }
    else 
    {
        rv <- qlindley(p = runif(n), beta, lower.tail = TRUE, log.p = FALSE)
        rv
    }
}

Portanto, no ambiente R, a distribuição de Lindley passa agora a ser acessível por meio das seguintes funções implementadas:

• dlindley: função densidade de probabilidade;
• plindley: função de distribuição acumulada;
• qlindley: função quantil;
• rlindley: função geradora de pseudovalores aleatórios.

Em relação ao processo de estimação da distribuição de Lindley, considerando sua função densidade de probabilidade, a função de verossimilhança para uma amostra \(\mathbf{X} = (X_1, \ldots, X_n)\), composta por variáveis independentes e identicamente distribuídas, é dada por:

\[\begin{align}\\ L(\beta;\mathbf{x}) = \prod_{i=i}^{n} \frac{\beta^2}{1 + \beta} (1 + x_i) e^{-\beta x_i}\\\\ \end{align}\]

O estimador de máxima verossimilhança (MLE) do parâmetro \(\beta\) consiste em maximizar a função \(L(\beta)\). Na prática, é preferível trabalhar com o logaritmo da função de verossimilhança, uma vez que a maximização da função original e de seu logaritmo conduzem ao mesmo estimador, mas a segunda abordagem é computacionalmente mais estável e conveniente. Dessa forma, o logaritmo da função de verossimilhança, denotado por \(\ell(\beta;;\mathbf{x})\), é expresso como:

\[\begin{align}\\ \ell(\beta;\mathbf{x}) = 2n \ln \beta - n \ln (1 - \beta) + \sum_{i=1}^{n} \ln(1 + x_{i}) - n\beta \bar{x}\\\\ \end{align}\]

Para encontrar o estimador \(\hat{\beta}\), deriva-se \(\ell(\beta)\) em relação a \(\beta\) e resolve-se a equação de primeira ordem \(\partial \ell(\beta; \mathbf{x})/\partial \beta = 0\). No caso da distribuição de Lindley, o estimador de máxima verossimilhança do parâmetro \(\beta\) é dado pela expressão explícita:

\[\begin{align}\\ \hat\beta = \dfrac{(1 - \overline{X}) + \sqrt{(\overline{X} - 1)^2 + 8\overline{X}}}{2\overline{X}}\\\\ \end{align}\]

em que \(\overline{X} = \frac{1}{n} \sum_{i=1}^n X_i\) representa a média amostral. Considerando que esse estimador não possui uma forma trivial, propõe-se a realização de um estudo de simulação Monte Carlo com o objetivo de avaliar as propriedades estatísticas do estimador \(\hat{\beta}\), verificando se apresenta características desejáveis para ser considerado um bom estimador do parâmetro \(\beta\). Para tanto, o estudo será conduzido com um número fixo de simulações, e variação dos tamanhos amostrais e dos valores nominais do parâmetro da distribuição de Lindley. Especificamente, serão realizadas simulações considerando a seguinte configuração:

• Número de Simulações \((R)\): 1000.
• Tamanhos Amostrais \((n)\): 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.
• Valores Nominais (Verdadeiros) do Parâmetro \(\beta\): \(0{.}5\), \(1{.}0\), \(1{.}5\), \(2{.}0\), \(2{.}5\) e \(3{.}0\).
• Propriedades Avaliadas: B, EQM.

Como o estimador de máxima verossimilhança da distribuição de Lindley possui expressão fechada, não se faz necessário o emprego da função fitdist do pacote fitdistrplus. Além disso, a rotina será otimizada com o uso da função sapply, integrante da família apply, permitindo aplicar a função de estimação de forma vetorizada. Ainda assim, permanecerão necessários dois laços de repetição: um correspondente à variação dos tamanhos amostrais e outro referente aos diferentes valores nominais do parâmetro \(\beta\) considerados no experimento. Por fim, os resultados obtidos serão apresentados também sob a forma gráfica, facilitando a análise visual do desempenho do estimador \(\hat{\beta}\) em diferentes cenários simulados. Todo o experimento será detalhado e implementado no Código 4.12).

Código 4.12. Estudo das propriedades do estimador de máxima verossimilhança \(\hat{\beta}_{\mathrm{MLE}}\) da distribuição de Lindley via simulação Monte Carlo, considerando diferentes tamanhos amostrais e valores nominais do parâmetro \(\beta\).

# -----------------------------------------------------------------------
# Propriedades dos MLE da Distribuição Lindley via Simulação Monte Carlo 
# -----------------------------------------------------------------------

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

library(fitdistrplus)
library(ggplot2)
library(tidyr)
library(dplyr)
library(knitr)
library(lamW)

# --- 1. Função para o cálculo do MLE da distribuição Lindley ---

mle.lindley <- function(x, par)
{
  m_x       <- mean(x)
  p1        <- sqrt((m_x - 1)^2 + 8 * m_x)
  p2        <- (m_x - 1)
  beta_hat  <- (p1 - p2) / (2 * m_x)
  
  return(beta_hat)
}

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

set.seed(1212)

betas       <- seq(0.5, 3.0, 0.5)   # Valores do parâmetro beta
param       <- data.frame(beta = betas)
n_param     <- nrow(param)

R           <- 1000                  # Número de simulações
nmin        <- 10                    # Menor tamanho amostral
nmax        <- 100                   # Maior tamanho amostral
enes        <- seq(nmin, nmax, 10)   # Tamanhos amostrais
n_enes      <- length(enes)

# --- 3. Criação das matrizes para armazenar os resultados ---

v.beta      <- matrix(nrow = n_enes, ncol = n_param)
eqm.beta    <- matrix(nrow = n_enes, ncol = n_param)

time_mat    <- matrix(nrow = n_enes, ncol = n_param)

# --- 4. Execução das simulações Monte Carlo ---

for(i in 1:n_param)
{
  
  X         <- rlindley(nmax * R,
                        beta     = param$beta[i],
                        mixture  = FALSE)
  X         <- matrix(X, nrow = nmax, ncol = R)
  
  for(k in 1:n_enes)
  {
    n       <- enes[k]
    x       <- data.frame(X[1:n, ])
    
    tempo   <- system.time({fit <- sapply(x,
                                          mle.lindley,
                                          par = param$beta[i])})["elapsed"]
    
    # Cálculo do viés e do EQM do estimador beta
    
    v.beta[k, i]    <- mean(fit - param$beta[i], na.rm = TRUE)
    eqm.beta[k, i]  <- mean((fit - param$beta[i])^2, na.rm = TRUE)
    
    # Armazenamento do tempo computacional
    
    time_mat[k, i]  <- tempo
  }
}

# --- 5. Prepração dos dados para os gráficos ---

mat_to_df         <- function(mat, medida_name, enes, n_param)
{
  df              <- as.data.frame(mat)
  df$n            <- enes
  df_long         <- pivot_longer(df, cols = -n, names_to = "parametro", values_to = "valor")
  df_long$medida  <- medida_name
  return(df_long)
}

# --- 6. Função de geração dos gráficos ---

df_plot            <- function(v_mat, eqm_mat, enes, param, titulo)
{
  n_param          <- nrow(param)
  
  df_all           <- bind_rows(mat_to_df(v_mat, "B", enes, n_param),
                                mat_to_df(eqm_mat, "EQM", enes, n_param))
  
  df_all$parametro <- factor(df_all$parametro,
                             levels = paste0("V", 1:n_param),
                             labels = paste0("param_", 1:n_param))
  
  df_all           <- df_all %>%
                        left_join(data.frame(parametro = paste0("param_", 1:n_param),
                                             beta      = param$beta,
                                             param_num = 1:n_param),
                                  by = "parametro") %>%
                          mutate( param_label = paste0("β = ", beta))
  
  ref_lines        <- data.frame(medida = c("B", "EQM"), yintercept = c(0, 0))
  
  x_labels         <- data.frame(n = enes, label = paste0("n = ", enes))
  
  
  p                <- ggplot(df_all, aes(x = n, y = valor, group = param_num, color = param_label)) +
                        geom_line() +
                        geom_point() +
                        geom_hline(data = ref_lines, 
                                   aes(yintercept = yintercept), 
                                   linetype = "dashed",
                                   color = "red",
                                   size = 0.5) +
                        geom_text(data = x_labels,
                                  aes(x = n, y = -Inf, label = label),
                                  inherit.aes = FALSE,
                                  angle = 0,
                                  hjust = 0.5,
                                  vjust = 1.2,
                                  size = 3,
                                  color = "black") +
                        facet_wrap(~ medida, nrow = 1, ncol = 2, scales = "free_y",
                                   labeller = as_labeller(c(B = "Viés (B)", 
                                                            EQM = "Erro Quadrático Médio (EQM)"))) +
                        labs(title = titulo, 
                             x = "Tamanho Amostral (n)",
                             y = "Valores",
                             color = "Cenários") +
                        scale_x_continuous(limits = c(min(enes), max(enes)),
                                           expand = expansion(mult = c(0.05, 0.05))) +
                        theme_minimal() +
                        theme(legend.position = "bottom",
                              legend.title = element_text(size = 11, face = "bold"),
                              legend.text = element_text(size = 11),
                              strip.text = element_text(face = "bold", size = 11),
                              plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
                              axis.text.x = element_text(angle = 0, hjust = 0.5),
                              axis.title.y = element_text(margin = margin(r = 15)),
                              axis.title.x = element_text(margin = margin(t = 15))) +
                        guides(color = guide_legend(nrow = 1, byrow = TRUE)) +
                        scale_color_discrete()
  return(p)
}

# --- 7. Ilustração gráfica para o parâmetro β ---

df_plot(v.beta, eqm.beta, enes, param, '')

Figura 4.12.. Gráfico das propriedades inferenciais do estimador de máxima verossimilhança do parâmetro \(\beta\) da distribuição Lindley, obtidas via simulação Monte Carlo para diferentes tamanhos amostrais e valores paramétricos com \(R = 1000\) replicações.

# --- 8. Ilustração gráfica do tempo de execução das simulações Monte Carlo em cada cenário ---

time_df                       <- as.data.frame(time_mat)
time_df$n                     <- enes
time_long                     <- pivot_longer(time_df, cols = -n, names_to = "parametro", 
                                              values_to = "tempo")
time_long$parametro           <- factor(time_long$parametro, levels = paste0("V", 1:n_param))
levels(time_long$parametro)   <- paste0("param_", 1:n_param)

time_long                     <- time_long %>%
                                    left_join(data.frame(parametro = paste0("param_", 1:n_param),
                                                         beta      = param$beta,
                                                         param_num = 1:n_param),
                                              by = "parametro") %>%
                                                                mutate(param_label = paste0("β = ", beta))


ggplot(time_long, aes(x = n, y = tempo, group = param_num, color = param_label)) +
  geom_line() +
  geom_point() +
  geom_text(data = data.frame(n = enes, label = paste0("n = ", enes)),
            aes(x = n, y = -Inf, label = label),
            inherit.aes = FALSE, angle = 0, hjust = 0.5, vjust = 1.2, size = 3, color = "black") +
  labs(title = "",
       x = "Tamanho Amostral (n)",
       y = "Tempo (segundos)",
       color = "Cenários") +
  scale_x_continuous(limits = c(min(enes), max(enes)), expand = expansion(mult = c(0.05, 0.05))) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 11, face = "bold"),
        legend.text = element_text(size = 11),
        plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
        axis.text.x = element_text(angle = 0, hjust = 0.5),
        axis.title.y = element_text(margin = margin(r = 15)),
        axis.title.x = element_text(margin = margin(t = 15))) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE)) +
  scale_color_discrete()

Figura 4.13.. Gráfico do tempo de execução (em segundos) em função do tamanho amostral (n), estratificado por cenários definidos por diferentes valores dos parâmetros da distribuição Lindley com \(R = 1000\) replicações.

4.4. Técnicas de Redução de Variância

Embora o método Monte Carlo se destaque como uma abordagem computacional versátil para a solução de problemas complexos, sua aplicação enfrenta um desafio fundamental: a alta variância intrínseca às estimativas produzidas. Essa limitação torna-se particularmente acentuada em dois cenários: (i) problemas de alta dimensionalidade, onde o fenômeno da “maldição da dimensionalidade” se manifesta, e (ii) situações envolvendo distribuições de probabilidade com características não-triviais, como multimodalidade ou caudas pesadas. Para enfrentar esse desafio, desenvolveram-se técnicas sofisticadas de redução de variância, que permitem melhorar substancialmente a eficiência computacional do método. Dentre estas técnicas, destacam-se:

• Amostragem por Importância (Importance Sampling): Reduz a variância ao concentrar a amostragem em regiões mais relevantes do espaço de parâmetros, utilizando pesos de importância para corrigir o viés;

• Amostragem Antitética (Antithetic Sampling): Gera variáveis correlacionadas negativamente para reduzir a variância, mantendo o estimador não viciado;

• Variáveis de Controle (Control Variables): Aproveita a correlação entre uma variável auxiliar (cuja média é conhecida) e a variável de interesse para ajustar a estimativa, reduzindo o erro.

4.4.1. Amostragem por Importância

---

A amostragem por importância constitui uma técnica de redução de variância, fundamentada na reexpressão do valor esperado de uma função de interesse \(h(\mathbf{X})\), cujo objetivo é a estimação do parâmetro \(\theta = E_f[h(\mathbf{X})]\), sendo \(\mathbf{X} \in \mathbb{R}^d\) uma variável aleatória com densidade de probabilidade \(f(\mathbf{x})\). O método baseia-se na introdução de uma densidade auxiliar \(g(\mathbf{x})\), denominada densidade de importância (ou função de amostragem de importância; ou, ainda, envelope), a qual satisfaz a condição \(g(\mathbf{x}) > 0\) sempre que \(f(\mathbf{x}) h(\mathbf{x}) \neq 0\). Nessas condições, o valor esperado pode ser reescrito na forma:

\[\begin{align}\\ \theta = E_f[h(\mathbf{X})] &= \int_{\mathbb{R}^d} h(\mathbf{x}) f(\mathbf{x}) \, d\mathbf{x} \\\\ &= \int_{\mathbb{R}^d} h(\mathbf{x}) \frac{f(\mathbf{x})}{g(\mathbf{x})} g(\mathbf{x}) \, d\mathbf{x} \\\\ &= E_g\left[ h(\mathbf{X}) \frac{f(\mathbf{X})}{g(\mathbf{X})} \right]\\\\ \end{align}\]

A partir dessa formulação, pode-se construir um estimador de Monte Carlo baseado em amostras extraídas de \(g\). Neste caso, seja \(X_1, \dots, X_n \overset{\text{i.i.d.}}{\sim} g\), então o estimador Monte Carlo via amostragem por importância é definido como:

\[\begin{align}\\ \hat{\theta}_{MC_{AI}}^* = \frac{1}{n} \sum_{i=1}^n h(X_i) w^*(X_i)\\\\ \end{align}\]

em que \(w^*(X_i) = f(X_i)/g(X_i)\) definem os pesos de importância não-padronizados (também chamados de razões de importância). Esses pesos, em particular, corrigem o viés introduzido pela mudança da densidade de referência, permitindo que a média ponderada preserve as propriedades do valor esperado sob a distribuição original. Por outro lado, definindo a variável auxiliar \(T(X) = h(X) w^*(X)\), tem-se que o estimador \(\hat{\theta}_{MC_{AI}}^*\) pode ser reescrito como uma média amostral de \(T(X)\) sob \(g\), isto é,

\[\begin{align}\\ \hat{\theta}_{MC_{AI}}^* = \frac{1}{n} \sum_{i=1}^{n} T(X_i) \\\\ \end{align}\]

Sob condições de regularidade, \(\hat{\theta}_{MC_{AI}}^*\) é não-viciado, uma vez que:

\[\begin{align}\\ E_g\left[ \hat{\theta}_{MC_{AI}}^* \right] = \frac{1}{n} \sum_{i=1}^{n}E_g\left[ T(X_i) \right] = E_g\left[ T(X) \right] = \theta\\\\ \end{align}\]

E sua variância é dada por:

\[\begin{align}\\ \text{Var}_g\left( \hat{\theta}_{MC_{AI}}^* \right) = \frac{1}{n^2} \sum_{i=1}^{n} \text{Var}_g\left( T(X_i) \right) = \frac{1}{n} \text{Var}_g\left( T(X) \right)\\\\ \end{align}\]

Essa expressão, em particular, destaca a importância da escolha adequada de \(g\), de modo a minimizar a variância dos pesos e, consequentemente, a variância do estimador. Por fim, o erro-padrão de \(\hat{\theta}_{MC_{AI}}^*\) pode ser estimado por:

\[\begin{align}\\ \widehat{SE}\left( \hat{\theta}_{MC_{AI}}^* \right) = \frac{1}{\sqrt{n}} \, S_T\\\\ \end{align}\]

em que \(S_T\) representa o desvio-padrão amostral dos valores \(T(X_1), \dots, T(X_n)\). Suponha agora que os pesos de importância sejam padronizados de acordo com

\[\begin{align}\\ w(X_i) = \frac{w^*(X_i)}{\sum_{j=1}^{n} w^*(X_j)}\\\\ \end{align}\]

em que \(w^*(X_i) = \dfrac{f(X_i)}{g(X_i)}\) são os pesos de importância não-padronizados previamente definidos. Neste caso, o estimador Monte Carlo correspondente para \(\theta = E_f[h(\mathbf{X})]\) é então dado por:

\[\begin{align}\\ \hat{\theta}_{MC_{AI}} = \frac{1}{n} \sum_{i=1}^n h(X_i) w(X_i)\\\\ \end{align}\]

em que \(w(X_i)\) são os pesos de importância padronizados, satisfazendo \(\sum_{i=1}^n w(X_i) = 1\). Este estimador, em particular, pode ser reescrito como a razão entre médias, isto é,

\[\begin{align}\\ \hat{\theta}_{MC_{AI}} = \frac{\bar{t}}{\bar{w}^*}\\\\ \end{align}\]

em que \(\bar{t} = \frac{1}{n} \sum_{i=1}^n t(X_i)\), com \(t(X_i) = h(X_i) w^*(X_i)\), e \(\bar{w}^* = \frac{1}{n} \sum_{i=1}^n w^*(X_i)\). Neste caso, para avaliar o viés introduzido pela padronização, aplica-se uma expansão de Taylor de ordem superior em torno de \(\bar{w}^* = 1\), o que nos leva:

\[\begin{align}\\ E\left[ \hat{\theta}_{MC_{AI}} \right] &= E\left[ \bar{t} \left( 1 - (\bar{w}^* - 1) + (\bar{w}^* - 1)^2 + \dots \right) \right] \\\\ &= E\left[ \bar{t} - (\bar{t} - \theta)(\bar{w}^* - 1) - \theta(\bar{w}^* - 1) + \bar{t}(\bar{w}^* - 1)^2 + \dots \right] \\\\ &= \theta - \frac{1}{n} \text{Cov}\{t(X), w^*(X)\} + \frac{\theta}{n} \text{Var}\{w^*(X)\} + \mathcal{O}\left( \frac{1}{n^2} \right)\\\\ \end{align}\]

A variância do estimador com pesos padronizados, \(\hat{\theta}_{MC_{AI}}\), é obtida de forma análoga à do estimador com pesos não-padronizados, resultando em:

\[\begin{align}\\ \text{Var}\{\hat{\theta}_{MC_{AI}}\} &= \frac{1}{n} \left[ \text{Var}\{t(X)\} + \theta^2 \text{Var}\{w^*(X)\} - 2\theta \text{Cov}\{t(X), w^*(X)\} \right] + \mathcal{O}\left( \frac{1}{n^2} \right)\\\\ \end{align}\]

Esta variância, em particular, pode ser estimada empiricamente a partir das amostras simuladas. Os resultados acima implicam que: (i) a padronização dos pesos de importância introduz um viés de ordem \(\mathcal{O}(1/n)\) no estimador, e (ii) a comparação entre os erros quadráticos médios (EQM) dos estimadores com e sem padronização fornece:

\[\begin{align}\\ \text{EQM}\{\hat{\theta}_{MC_{AI}}\} - \text{EQM}\{\hat{\theta}^*_{MC_{AI}}\} &= \frac{1}{n} \left[ \theta^2 \text{Var}\{w^*(X)\} - 2\theta \text{Cov}\{t(X), w^*(X)\} \right] + \mathcal{O}\left( \frac{1}{n^2} \right)\\\\ \end{align}\]

Supondo, sem perda de generalidade, que \(\theta > 0\), os termos dominantes indicam que essa diferença será negativa — isto é, a padronização resulta em menor EQM — sempre que:

\[\begin{align}\\ \text{Corr}\{t(X), w^*(X)\} > \frac{\text{CV}\{w^*(X)\}^2}{\text{CV}\{t(X)\}}\\\\ \end{align}\]

em que \(\text{CV}\{\cdot\}\) denota o coeficiente de variação da variável considerada. Essa condição expressa uma relação entre a correlação e a dispersão relativa das variáveis ponderadas, indicando quando a padronização dos pesos é vantajosa em termos de desempenho médio-quadrático.

Exemplo 4.7 (Probabilidade de Falha de Redes - Givens e Hoeting 2012). Muitos sistemas podem ser representados por “gráficos conectados” (também conhecidos como grafos, ou gráficos de rede). Este tipo de gráfico é formado por nós (círculos) e rotas (segmentos de reta ou arestas), tais que um sinal de A para B deve seguir o caminho por qualquer rota disponível. A Figura 4.14 mostra um exemplo deste tipo gráfico.

Figura 4.14.. Ilustração de um gráfico de rede de um ponto A até um ponto B.

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

Uma falha de rede, por outro lado, acontece quando não há rotas ou há rotas “quebradas” na transmissão entre dois nós, o que nos implica que o sinal falha ao tentar ser transmitido corretamente entre qualquer par de nós conectados. Um exemplo deste tipo de situação é expresso na Figura 4.15, onde há poucas rotas disponíveis para a transmissão de sinal.

Figura 4.15.. Ilustração de um gráfico de rede com rotas “quebradas”.

Fonte: Broken Network Graph - Importance Sampling | Givens, G. H. & Hoeting, J. A. (2012). Computational Statistics. 2\(^{\text{nd}}\) Edition. John Wiley & Sons, Inc.

Grafos de rede, em geral, são utilizados modelar a transmissão de diversos tipos de sinais, como transmissão de voz analógica, sinais digitais eletromagnéticos e transmissão óptica de dados digitais. O modelo também pode ser mais conceitual, com cada aresta representando diferentes máquinas ou pessoas cuja participação pode ser necessária para alcançar algum resultado. Usualmente, uma quantidade de interesse é a probabilidade de falha da rede, dadas probabilidades específicas para a falha de cada aresta. Para ilustrar este conceito, seja \(\mathbf{X}\) uma rede qualquer tal que, para cada rota, tenha-se como resultado: rota intacta ou rota quebrada. Considere, para este exemplo, que a rede tenha 20 rotas, isto é, \(\mathbf{X} = (X_1, \dots, X_{20})\), e defina \(b(\mathbf{X})\) como o número de rotas quebradas em \(\mathbf{X}\), tal que \(\mu = h(\mathbf{X})\) seja o um parâmetro que indica a falha da rede de modo que:

• \(h(\mathbf{X}) = 1\) se \(A\) não está conectado com \(B\).

• \(h(\mathbf{X}) = 0\) se \(A\) está conectado com \(B\).

Logo, a probabilidade de falha da rede será descrita por \(\mu = E[h(\mathbf{X})]\). Naturalmente, estimar \(\mu\) para uma rede de qualquer tamanho pode ser um problema combinatório complexo devido as questões de convergência, custo, etc dos algoritmos computacionais. Porém, este problema pode ser resolvido utilizando o estimador de Monte Carlo. Para tal, considere uma amostra aleatória \(X_1, \dots, X_n\) tal que a rota falha independentemente com probabilidade \(\rho\). Assim, o estimador Monte Carlo é dado por:

\[\begin{align}\\ \hat{\mu}_{MC} = \frac{1}{n} \sum_{i=1}^n h(X_i)\\\\ \end{align}\]

Neste caso, note que este estimador tem variância dada por \(\text{Var}(\hat{\mu}_{MC}) = [\mu(1 - \mu)]/n\). No entanto, o problema do estimador \(\hat{\mu}_{MC}\) é que \(h(\mathbf{X})\) dificilmente é igual a 1, a menos que \(\rho\) seja extremamente grande. Isto nos diz que um grande número de redes pode necessitar de um estimador com mais precisão para identificar o número de falhas na rede, que pode ser obtido usando a amostragem por importância. Neste caso, iremos considerar os pesos de importância não-padronizados (apenas por questões de praticidade neste problema, em específico), com foco na simulação de \(\mathbf{X}\) de modo que \(h(\mathbf{X}) = 1\). Logo, a precisão do estimador fica compensada pela atribuição de pesos de importância. Isto é, será simulada uma amostra \(\mathbf{X}_1^*, \dots, \mathbf{X}_n^*\) por meio da geração de configurações de redes formadas pelas rotas quebradas, assumindo falhas de rota independentes com probabilidade \(\rho^* > \rho\). Logo, os pesos de importância para \(\mathbf{X}_i^*\) podem ser escritos como:

\[\begin{align}\\ w^*(\mathbf{X}_i^*) = \left[ \frac{1 - \rho}{1 - \rho^*} \right]^{20} \left[ \frac{\rho (1 - \rho^*)}{\rho^* (1 - \rho)} \right]^{b(\mathbf{X}_i^*)}\\\\ \end{align}\]

em que \(\rho\) é um parâmetro que indica a falha da rede. Deste modo, tem-se que o estimador Monte Carlo para \(\mu\), via amostragem por importância, é definido por:

\[\begin{align}\\ \hat{\mu}_{MC_{AI}}^* = \frac{1}{n} \sum_{i=1}^n h(\mathbf{X}_i^*) w^*(\mathbf{X}_i^*)\\\\ \end{align}\]

Seja, agora, \(\mathcal{C}\) o conjunto de todas as possíveis configurações de rede, e seja \(\mathcal{F}\) o subconjunto de todas as configurações de rede em que A e B não estão conectados. Então, a variância do estimador de Monte Carlo, via amostragem de importância, é descrita por:

\[\begin{align}\\ \text{Var}(\hat{\mu}_{MC_{AI}}^*) &= \frac{1}{n} \, \text{Var}\big{\{}h(\mathbf{X}_i^*) w^*(\mathbf{X}_i^*)\big{\}}\\\\ & = \frac{1}{n}\left(E\left\{[h(X_{i}^*) w^*(x^*_{i})]^2\right\} - \left[E\left\{h(X_{i}^*) w^*(x^*_{i})\right\}\right]^2\right)\\\\ & = \frac{1}{n} \left(\sum_{\mathbf{x} \in \mathcal{F}} w^*(\mathbf{x}) \, \rho^{b(\mathbf{x})} (1 - \rho)^{20 - b(\mathbf{x})} - \mu^2\right)\\\\ \end{align}\]

Considere, por exemplo, que as falhas ocorrem apenas quando \(b(\mathbf{X}) \geqslant 4\). Então, obtém-se que:

\[\begin{align}\\ w^*(\mathbf{X}_i^*) \leqslant \left[ \frac{1 - \rho}{1 - \rho^*} \right]^{20} \left[ \frac{\rho (1 - \rho^*)}{\rho^* (1 - \rho)} \right]^{b(\mathbf{X}_i^*)}\\\\ \end{align}\]

Note que, quando \(\rho^* = 0.25\) e \(\rho = 0.05\), obtém-se que \(w^*(\mathbf{X}) \leqslant 0.07\). Neste caso, tem-se que a variância deste estimador é dada por \(\text{Var}(\hat{\mu}_{MC_{AI}}^*) = (0.07\mu - \mu^2)/n\). Isto nos diz que \(\text{Var}(\hat{\mu}_{MC_{AI}}^*)\) é bem menor do que \(\text{Var}(\hat{\mu}_{MC})\). Em particular, se for considerado a aproximação \(c\mu - \mu^2 \approx c\mu\) para valores pequenos de \(\mu\) e valores grandes de \(c\), obtém-se que:

\[\begin{align}\\ \frac{\text{Var}(\hat{\mu}_{MC})}{\text{Var}(\hat{\mu}_{MC_{AI}}^*)} \approx 14\\\\ \end{align}\]

isto é, \(\text{Var}(\hat{\mu}_{MC})\) é muito maior que \(\text{Var}(\hat{\mu}_{MC_{AI}}^*)\). Para clarificar as ideias deste exemplo, o Código 4.13 uma ilustração em ambiente R utilizando uma rede de comunicação genérica com numedges = 20, onde numedges representa o número total de arestas na rede, considerando-se os valores \(\rho = 0{.}05\) e \(\rho^* = 0{.}2\), com base em \(n = 100000\) iterações aplicadas ao estimador Monte Carlo.

Código 4.13. Implementação de rede de comunicação genérica com numedges = 20, e considerando-se os valores \(\rho = 0{.}05\) e \(\rho^{*} = 0{.}2\), com base em \(n = 100000\) iterações aplicadas ao estimador Monte Carlo.

#| message: false
#| warning: false
#| eval: true
#| echo: true

# ---------------------------------------------------------
# Amostragem por Importância: Falha em Rede de Comunicação
# ---------------------------------------------------------

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

library(knitr)

# --- 1. Construção da rede ---

# As letras indicam os vínculos. Os caracteres são atribuídos em ordem alfabética,
# da esquerda para a direita e, em seguida, de cima para baixo. Assim, a, b, c, d 
# representam os quatro vínculos de A para um dos nós diretamente à sua direita. 
# O nó superior na segunda coluna possui os vínculos e, f, g, h para os quatro nós
# na terceira coluna. Em seguida, os rótulos são atribuídos ao segundo nó na segunda 
# coluna. E assim por diante.

connecteds <- c("ahnr", "ahnos", "ahnopt", "ahnolmt", "ahnokgmt", "ahnjfks", 
                "ahnjfkpt", "ahnjfklmt", "ahnjfgls", "ahnjfglpt", "ahnjfgmt", 
                "ahnjfgmps", "ahq", "ainq", "air", "aios", "aiopt", "aiolmt", 
                "aiokfmt", "aijfgls", "aijfglpt", "aijfgmt", "aijfgmps", "bjr",
                "aijfks", "aijfkpt", "aijfklmt", "aejnq", "aejr", "aejos",
                "aejopt", "aejolmt", "aejokgmt", "aefkonq", "aefkor", "aefks",
                "aefkpt", "aefklmt", "aefglonq", "aefgls", "aefglpt", "bfgmt",
                "aefgmt", "aefgmps", "aefgmpor", "aefgmponq", "bfks", "bjopt",
                "behq", "behnr", "behnos", "behnopt", "behnokgmt", "behnolmt", 
                "beinq", "beir", "beios", "beiopt", "beiolmt", "beiokgmt", 
                "bjnq", "bjos", "bjolmt", "bjokgmt", "bfgmponq", "bfgmpoihq",
                "bjihq", "bfkonq", "bfkor", "bfkoihq", "bfkpt", "bfgmpor",
                "bfklmt", "bfgloihq", "bfglonq", "bfglor", "bfgls", "bfglpt",
                "bfgmps", "cfehq", "cfehnr", "cfehnos", "cfehnolmt", "cfehnopt", 
                "cfeinq", "cfeir", "cfeios", "cfeiopt", "cfeiolmt", "cfjnq",
                "cfjr", "cfjos", "cfjopt", "cfjolmt", "cfjihq", "ckonq", "ckor",
                "ckojehq", "ckoihq", "cks", "ckpt", "cklmt", "cglonq", "cglor",
                "cglojehq", "cgloihq", "cgls", "cglpt", "cgmponq", "cgmpor",
                "cgmpojehq", "cgmpoihq", "cgmps", "cgmt", "dgfehq", "dgfehnr",  
                "dgfehnos", "dgfehnopt", "dgfeinq", "dgfeir", "dgfeios",
                "dgfeiopt", "dgfjnq", "dgfjr", "dgfjos", "dgfjopt", "dgfjihq",
                "dgkonq", "dgkor", "dgkojehq", "dgks", "dgkpt", "dlonq", "dlor",
                "dlojehq", "dloihq", "dls", "dlpt", "dlkfehq", "dlkfehnr", 
                "dlkfeir", "dlkfeinq", "dlkfjr", "dlkfjnq", "dlkfjihq", "dmponq",
                "dmpor", "dmpojehq", "dmpoihq", "dmps", "dmpkfehq", "dmpkfehnr",
                "dmpkfeinq", "dmpkfeir", "dmpkfjihq", "dmpkfjnq", "dmpkfjr", "dmt") 

head(connecteds, n = 30)

 [1] "ahnr"      "ahnos"     "ahnopt"    "ahnolmt"   "ahnokgmt"  "ahnjfks"  
 [7] "ahnjfkpt"  "ahnjfklmt" "ahnjfgls"  "ahnjfglpt" "ahnjfgmt"  "ahnjfgmps"
[13] "ahq"       "ainq"      "air"       "aios"      "aiopt"     "aiolmt"   
[19] "aiokfmt"   "aijfgls"   "aijfglpt"  "aijfgmt"   "aijfgmps"  "bjr"      
[25] "aijfks"    "aijfkpt"   "aijfklmt"  "aejnq"     "aejr"      "aejos"    

# --- 2. Valores iniciais ---

# n = nº de iteraçõse Monte Carlo
# p = probabilidade de falha  
# pstar = probabilidade de falha no método de amostragem por importância  
# numedges = número de arestas na rede completa
  
p             <- 0.05
pstar         <- 0.2
n             <- 100000
numedges      <- 20
failed.orig   <- rep(FALSE, n)
failed.is     <- rep(FALSE, n)

# --- 3. Matriz de rota e conexões ---

# Aqui é gerada a estrutura de matriz lógica sinalizando os vínculos necessários
# em cada possível configuração conectada.

temp        <- matrix(0, nrow = 158, ncol = 9)

for(i in 1:9)
{
  temp[, i] <- substring(connecteds, i, i)
}

mylet       <- c("", letters[1:numedges])
linkmat     <- matrix(0, 158, 9)

for(i in 1:9)
{
  for(j in 1:158)
  {
    linkmat[j, i] <- (0:numedges)[mylet == temp[j, i]]
  }
}

link.logic.mat    <- matrix(F, nrow = 158, ncol = numedges)
for(i in 1:158) 
{
    for(j in 1:9) 
    {
        if(linkmat[i,j] > 0) {link.logic.mat[i,linkmat[i,j]] = TRUE} 
    } 
}

list('link.mat' = head(linkmat), 
     'link.logic.mat' = head(link.logic.mat))

$link.mat
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,]    1    8   14   18    0    0    0    0    0
[2,]    1    8   14   15   19    0    0    0    0
[3,]    1    8   14   15   16   20    0    0    0
[4,]    1    8   14   15   12   13   20    0    0
[5,]    1    8   14   15   11    7   13   20    0
[6,]    1    8   14   10    6   11   19    0    0

$link.logic.mat
     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7] [,8]  [,9] [,10] [,11] [,12]
[1,] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
[2,] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
[3,] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
[4,] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE  TRUE
[5,] TRUE FALSE FALSE FALSE FALSE FALSE  TRUE TRUE FALSE FALSE  TRUE FALSE
[6,] TRUE FALSE FALSE FALSE FALSE  TRUE FALSE TRUE FALSE  TRUE  TRUE FALSE
     [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,] FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE
[2,] FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE
[3,] FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE
[4,]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE
[5,]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE
[6,] FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE

# --- 4. Estimador de Monte Carlo (Clássico) ---

# Aqui é calculado a probabilidade de falha da rede e o número de falhas na rede via 
# método Monte Carlo clássico.

set.seed(126)
breaks        <- sample(c(T,F), numedges*n, prob = c(p, 1 - p), replace = T)
brokenlinks   <- matrix(breaks, nrow = n, ncol = numedges)

for(j in 1:n)
{
    config      <- brokenlinks[j,]
        if(sum(config) > 0)
        {
          tested = t(t(link.logic.mat) & (!config))
        failed.orig[j] =! any(!apply(xor(tested, link.logic.mat), 1, any))
      } 
}

mu_mc         <- mean(failed.orig)
var_mu_mc     <- var(failed.orig)/n
falhas_mu_mc  <- sum(failed.orig)

res_mu_mc     <- data.frame(row.names = c('Estimativa Monte Carlo (Clássico)',
                                          'Variância Monte Carlo (Clássico)', 
                                          'Falhas Identificadas (Clássico)'),
                            valores = c(mu_mc, var_mu_mc, falhas_mu_mc))

kable(res_mu_mc, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      digits = 15,
      col.names = c('Valores')) 

	Valores
Estimativa Monte Carlo (Clássico)	1e-05
Variância Monte Carlo (Clássico)	1e-10
Falhas Identificadas (Clássico)	1e+00

# --- 5. Estimador de Monte Carlo (Amostragem por importância) ---

# Aqui é calculado a probabilidade de falha da rede e o número de falhas na rede via 
# método Monte Carlo por amostragem de importância.

set.seed(127)
breaks        <- sample(c(T,F), numedges*n, prob = c(pstar, 1 - pstar), replace = T)
brokenlinks   <- matrix(breaks, nrow = n, ncol = numedges)

for(j in 1:n)
{
    config      <- brokenlinks[j,]
        if(sum(config) > 0)
        {
          tested = t(t(link.logic.mat) & (!config))
        failed.is[j] =! any(!apply(xor(tested, link.logic.mat), 1, any))
      } 
}

failures.in.config  <- apply(brokenlinks,1,sum)
weights             <- failures.in.config * (log(p) - log(pstar)) +
                          (numedges - failures.in.config) * (log(1 - p) - log(1 - pstar))
weights             <- exp(weights)

mu_ai         <- mean(failed.is*weights)
var_mu_ai     <- var(failed.is*weights)/n
falhas_mu_ai  <- sum(failed.is)

res_mu_ai     <- data.frame(row.names = c('Estimativa Monte Carlo (Amostragem por Importância)',
                                          'Variância Monte Carlo (Amostragem por Importância)', 
                                          'Falhas Identificadas (Amostragem por Importância)'),
                            valores = c(mu_ai, var_mu_ai, falhas_mu_ai))

kable(res_mu_ai, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      digits = 15,
      col.names = c('Valores')) 

	Valores
Estimativa Monte Carlo (Amostragem por Importância)	1.022175e-05
Variância Monte Carlo (Amostragem por Importância)	2.458000e-12
Falhas Identificadas (Amostragem por Importância)	4.910000e+02

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

4.4.2. Amostragem Antitética

---

A amostragem antitética constitui um segundo método de redução de variância, sendo fundamentada na construção de dois estimadores não-viciados, identicamente distribuídos, denotados por \(\hat{\mu}_1\) e \(\hat{\mu}_2\), e negativamente correlacionados, os quais são combinados para formar um novo estimador com menor variância. Nessas condições, define-se o estimador Monte Carlo via amostragem antitética como:

\[\begin{align}\\ \hat{\mu}_{AA} = \frac{\hat{\mu}_1 + \hat{\mu}_2}{2}\\\\ \end{align}\]

isto é, o estimador antitético é a média aritmética de \(\hat{\mu}_1\) e \(\hat{\mu}_2\). A variância desse novo estimador é descrita por:

\[\begin{align}\\ \text{Var}(\hat{\mu}_{AA}) = \frac{1}{4} \left[ \text{Var}(\hat{\mu}_1) + \text{Var}(\hat{\mu}_2) \right] + \frac{1}{2} \text{Cov}(\hat{\mu}_1, \hat{\mu}_2) = \frac{(1 + \rho)\sigma^2}{2n}\\\\ \end{align}\]

em que \(\rho\) representa o coeficiente de correlação entre \(\hat{\mu}_1\) e \(\hat{\mu}_2\), e \(\sigma^2/n\) é a variância de cada estimador individual, assumindo que ambos são construídos a partir de amostras de tamanho \(n\). Essa expressão, em particular, nos mostra que, sob \(\rho < 0\), \(\text{Var}(\hat{\mu}_{AA})\) será inferior àquela dos estimadores originais. No entanto, a dificuldade prática dessa técnica reside na seguinte questão: Dado um estimador inicial \(\hat{\mu}_1\), como se pode construir um segundo estimador, igualmente distribuído, \(\hat{\mu}_2\), que seja negativamente correlacionado com \(\hat{\mu}_1\)? Para responder a essa pergunta, considere um conjunto de variáveis aleatórias independentes e identicamente distribuídas (iid), denotado por \(X = \{X_1, \dots, X_n\}\), e defina:

\[\begin{align}\\ \hat{\mu}_1(X) = \frac{1}{n} \sum_{i=1}^n h_1(X_i)\\\\ \end{align}\]

em que \(X_i = (X_{i1}, \dots, X_{im})\), e \(h_1: \mathbb{R}^m \to \mathbb{R}\) é uma função definida sob \(X_i\), tal que \(E[h_1(X_i)] = \mu\). Isso nos diz que \(\hat{\mu}_1\) é um estimador não-viesado de \(\mu\). Em paralelo, defina o segundo estimador como:

\[\begin{align}\\ \hat{\mu}_2(X) = \frac{1}{n} \sum_{i=1}^n h_2(X_i)\\\\ \end{align}\]

onde \(h_2\) é uma função de mesma natureza que \(h_1\). Naturalmente, a correlação entre \(\hat{\mu}_1\) e \(\hat{\mu}_2\) dependerá da covariância entre as funções \(h_1(X_i)\) e \(h_2(X_i)\). Se, por exemplo, ambas funções forem estritamente crescentes (ou estritamente decrescentes) em cada argumento, então:

\[\begin{align}\\ \text{Cov}\{h_1(X_i), h_2(X_i)\} > 0\\\\ \end{align}\]

isto é, a covariância é positiva. Portanto, para se obter covariância negativa, é necessário empregar argumentos complementares. Neste caso, suponha que as funções \(h_1\) e \(h_2\) sejam aplicadas sobre vetores \(U_i = (U_{i1}, \dots, U_{im})\), com \(U_{ij} \sim \mathcal{U}(0, 1)\) independentes. Se essas funções forem monótonas em cada argumento, então:

\[\begin{align}\\ \text{Cov}\{h_1(U_{i1}, \dots, U_{im}), h_2(1 - U_{i1}, \dots, 1 - U_{im})\} \leqslant 0 \\\\ \end{align}\]

Este resultado estabelece a base para a construção de \(\hat{\mu}_2\) a partir de uma transformação antitética dos argumentos de \(\hat{\mu}_1\). Para tornar explícita essa construção, seja \(F_j\) a função de distribuição acumulada da componente \(X_{ij}\). A transformação inversa de \(F_j\) nos permite expressar \(\hat{\mu}_1\) como:

\[\begin{align}\\ \hat{\mu}_1(X) = \frac{1}{n} \sum_{i=1}^{n} h_1 \left( F_1^{-1}(U_{i1}), \dots, F_m^{-1}(U_{im}) \right)\\\\ \end{align}\]

Como as funções \(F_j^{-1}\) são não decrescentes e \(h_1\) é monótona, a composição preserva a monotonicidade. A partir da simetria da distribuição uniforme, define-se então:

\[\begin{align}\\ h_2(1 - U_{i1}, \dots, 1 - U_{im}) = h_1 \left( F_1^{-1}(1 - U_{i1}), \dots, F_m^{-1}(1 - U_{im}) \right) \\\\ \end{align}\]

o que mantém a propriedade de monotonicidade por argumento, além de assegurar que \(h_2\) possui a mesma distribuição que \(h_1\), garantindo a identidade em distribuição entre os estimadores resultantes. Dessa forma, obtém-se um segundo estimador, \(\hat{\mu}_2\), como:

\[\begin{align}\\ \hat{\mu}_2(X) = \frac{1}{n} \sum_{i=1}^{n} h_1 \left( F_1^{-1}(1 - U_{i1}), \dots, F_m^{-1}(1 - U_{im}) \right)\\\\ \end{align}\]

O par \((\hat{\mu}_1, \hat{\mu}_2)\) satisfaz:

• \(\hat{\mu}_1 \overset{d}{=} \hat{\mu}_2\),
• \(E[\hat{\mu}_1] = E[\hat{\mu}_2] = \mu\),
• \(\text{Cov}(\hat{\mu}_1, \hat{\mu}_2) \leqslant 0\).

Logo, o estimador Monte Carlo via amostragem antitética:

\[\begin{align}\\ \hat{\mu}_{AA} = \frac{\hat{\mu}_1 + \hat{\mu}_2}{2} \\\\ \end{align}\]

possui variância inferior àquela de um estimador baseado em \(n\) observações independentes, desde que a correlação entre \(\hat{\mu}_1\) e \(\hat{\mu}_2\) seja negativa. A eficácia do método reside, portanto, na exploração de simetrias e propriedades de monotonicidade para induzir correlação negativa entre réplicas estocásticas de um mesmo estimador, promovendo ganho de eficiência sem aumento do custo computacional.

Exemplo 4.8 (Média da Distribuição Normal). Seja \(X\) uma variável aleatória contínua com distribuição normal padrão, isto é, \(X \sim \mathcal{N}(0,1)\), e suponha que o objetivo seja estimar o valor esperado \(\mu = E[h(X)]\) em que a função de interesse é dada por:

\[\begin{align}\\ h(x) = \frac{x}{2^x - 1} \\\\ \end{align}\]

Observe que o cálculo analítico da integral correspondente à \(\mu\) não admite solução analítica, o que motiva o uso de técnicas de simulação, em particular o método Monte Carlo, para obtenção de uma estimativa numérica. Neste caso, pela abordagem tradicional, o estimador Monte Carlo é construído como a média amostral de \(N = 100000\) observações independentes de \(h(X_i)\), com \(X_i \sim \mathcal{N}(0,1)\), geradas de forma independente. Isto é,

\[\begin{align}\\ \hat{\mu}_{MC} = \frac{1}{100000} \sum_{i=1}^{100000} h(X_i) \\\\ \end{align}\]

Entretanto, o foco deste exemplo é a construção de um estimador mais eficiente por meio da técnica de amostragem antitética, cujo propósito é reduzir a variância da estimativa introduzindo correlação negativa entre pares de observações. Essa abordagem é particularmente eficaz em distribuições simétricas, como é o caso da normal padrão. A ideia do método consiste em utilizar, para cada amostra \(X_i\), seu correspondente antitético \(-X_i\), aproveitando a simetria em torno da origem. Dessa forma, ao invés de gerar \(100000\) observações independentes, gera-se apenas \(n = 50000\) observações \(X_1, \dots, X_n \sim \mathcal{N}(0,1)\), e constrói-se o seguinte estimador antitético:

\[\begin{align}\\ \hat{\mu}_{MC_{AA}} = \frac{1}{100000} \sum_{i=1}^{50000} \left[ h(X_i) + h(-X_i) \right]\\\\ \end{align}\]

Este estimador preserva o mesmo custo computacional do estimador tradicional, no sentido de que trabalha com o mesmo número total de avaliações da função \(h(x)\), porém ele explora a dependência negativa induzida pela simetria da distribuição. Neste caso, foi observada uma correlação empírica negativa de \(\widehat{\text{Cor}}(h(X_i), h(-X_i)) = -0{.}95\). Esse valor elevado de correlação negativa indica que a estratégia antitética é adequada neste caso, e deve produzir redução significativa na variância do estimador. Esse fato pode ser verificado nos resultados obtidos, em que:

• Estimativa com Amostragem Antitética: \(\hat{\mu}_{MC_{AA}} = 1.4992\), com variância estimada igual a \(1 \times 10^{-7}\);
• Estimativa tradicional de Monte Carlo: \(\hat{\mu}_{MC} = 1.4993\), com variância estimada igual a \(2.6 \times 10^{-6}\).

Embora as estimativas pontuais sejam praticamente idênticas, a variância do estimador antitético é cerca de 26 vezes menor, demonstrando a eficácia da técnica na prática. Essa melhoria justifica seu uso sempre que a estrutura da função \(h\) e a simetria da distribuição permitirem a construção de pares antitéticos com forte correlação negativa. O Código 4.14 apresenta a rotina em ambiente R utilizada para implementação dos procedimentos descritos neste exemplo.

Código 4.14. Implementação da técnica de amostragem antitética para estimar o valor esperado \(\mu = E[h(X)]\), com \(X \sim \mathcal{N}(0,1)\) e \(h(x) = x / (2^x - 1)\).

# ----------------------------------------------------
# Amostragem Antitética: Média da Distribuição Normal
# ----------------------------------------------------

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

library(knitr)

# --- 1. Definição dos valores iniciais ---

set.seed(1212)
n         <- 50000
N         <- 2 * n

# --- 2. Geração dos valores de X ~ N(0, 1) ---

X         <- rnorm(n)
X_ant     <- -X  # Antitéticos

data.frame('Normais' = X[1:20], 'Antitéticos' = X_ant[1:20])

       Normais Antitéticos
1  -0.62908580  0.62908580
2   1.83742241 -1.83742241
3   0.34368151 -0.34368151
4  -1.36009553  1.36009553
5   0.53581450 -0.53581450
6   0.03719874 -0.03719874
7  -1.92311898  1.92311898
8   0.41087971 -0.41087971
9  -2.24738050  2.24738050
10  2.06471825 -2.06471825
11  0.40559869 -0.40559869
12 -0.34958717  0.34958717
13  1.34170670 -1.34170670
14  1.67346194 -1.67346194
15  0.25818520 -0.25818520
16 -0.96876431  0.96876431
17 -0.36444202  0.36444202
18 -1.37055387  1.37055387
19  0.03709340 -0.03709340
20  0.96219532 -0.96219532

# --- 3. Definição da função h(x) ---

h         <- function(x)
{
  values  <- x / (2^x - 1)
  return(values)
}

# --- 4. Cálculo da correlação ---

h_X       <- h(X)
h_X_ant   <- h(X_ant)

rho         <- cor(h_X, h_X_ant)
rho

[1] -0.9523687

# --- 5. Estimador de Monte Carlo (Clássico) ---

mu.mc     <- mean(c(h_X, h_X_ant))
var.mc    <- var(c(h_X, h_X_ant)) / N

res.mu.mc <- data.frame(row.names = c('Estimativa Monte Carlo (Clássico)',
                                      'Variância Monte Carlo (Clássico)'),
                        valores = c(mu.mc, var.mc))
kable(res.mu.mc, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('Valores')) 

	Valores
Estimativa Monte Carlo (Clássico)	1.4992782
Variância Monte Carlo (Clássico)	0.0000026

# --- 6. Estimador de Monte Carlo (Amostragem antitética) ---

mu.aa       <- mean((h_X + h_X_ant) / 2)
var.aa    <- ((1 + rho) * var(h_X)) / N

res.mu.aa <- data.frame(row.names = c('Estimativa Monte Carlo (Amostragem Antitética)',
                                      'Variância Monte Carlo (Amostragem Antitética)'),
                        valores = c(mu.aa, var.aa))

kable(res.mu.aa, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('Valores')) 

	Valores
Estimativa Monte Carlo (Amostragem Antitética)	1.4992782
Variância Monte Carlo (Amostragem Antitética)	0.0000001

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

4.4.3. Variáveis de Controle

---

A técnica de variável de controle constitui um terceiro método de redução de variância cujo objetivo é aprimorar a precisão da estimação de uma quantidade de interesse, por meio do uso de informações auxiliares disponíveis a priori. A ideia central consiste em utilizar uma variável auxiliar correlacionada com a função de interesse, cujo valor esperado seja conhecido analiticamente. Com base nessa informação adicional, é possível ajustar a estimativa original, explorando a correlação existente de modo a reduzir a variabilidade do estimador final. Para ilustrar a ideia, seja \(\mu = E[h(X)]\) a quantidade de interesse a ser estimada, onde \(h\) é uma função mensurável e \(X\) é uma variável aleatória com distribuição conhecida ou simulável. Suponha que exista uma função auxiliar \(c\) aplicada a uma variável \(Y\), de tal forma que \(\theta = E[c(Y)]\) tenha expressão analítica, e que \(h(X)\) e \(c(Y)\) sejam correlacionadas. Para fins de estimação, considere que se dispõe de uma amostra de tamanho \(n\) de pares independentes \((X_i, Y_i)\), com \(i = 1, \ldots, n\), gerados via simulação de forma tal que:

\[\begin{align}\\ \text{Cov}\{X_i, X_j\} = \text{Cov}\{Y_i, Y_j\} = \text{Cov}\{X_i, Y_j\} = 0 \\\\ \end{align}\]

para \(i \neq j\), isto é, os pares \((X_i, Y_i)\) são mutuamente independentes entre si. Nesse contexto, os estimadores de Monte Carlo para \(\mu\) e \(\theta\) são dados por:

\[\begin{align}\\ \hat{\mu}_{MC} &= \frac{1}{n} \sum_{i=1}^n h(X_i) \\\\ \hat{\theta}_{MC} &= \frac{1}{n} \sum_{i=1}^n c(Y_i)\\\\ \end{align}\]

A utilidade da técnica de variáveis de controle surge da possibilidade de que \(\hat{\mu}_{MC}\) e \(\hat{\theta}_{MC}\) estejam correlacionados, ou seja, \(\text{Corr}(h(X_i), c(Y_i)) \neq 0\). Por exemplo, quando essa correlação é positiva, valores elevados de \(\hat{\theta}_{MC}\) tendem a estar associados a valores igualmente elevados de \(\hat{\mu}{MC}\), de modo que um ajuste “para baixo” (negativo) em \(\hat{\mu}_{MC}\) pode ser apropriado, caso \(\hat{\theta}_{MC}\) supere o valor verdadeiro \(\theta\). Por outro lado, uma correlação negativa entre as funções implica ajuste “para cima” (positivo). Essa lógica fundamenta a construção do estimador Monte Carlo por variável de controle, definido por:

\[\begin{align}\\ \hat{\mu}_{CV} = \hat{\mu}_{MC} + \lambda (\hat{\theta}_{MC} - \theta)\\\\ \end{align}\]

em que \(\lambda\) é um parâmetro escalar a ser escolhido com o intuito de minimizar a variância do estimador ajustado. A variância de \(\hat{\mu}_{CV}\) é dada por:

\[\begin{align}\\ \text{Var}\{\hat{\mu}_{CV}\} &= \text{Var}\{\hat{\mu}_{MC}\} + \lambda^2\, \text{Var}\{\hat{\theta}_{MC}\} + 2\lambda\, \text{Cov}\{\hat{\mu}_{MC}, \hat{\theta}_{MC}\}\\\\ \end{align}\]

A escolha ótima de \(\lambda\) corresponde àquele valor que minimiza a expressão acima. Assim, derivando com respeito a \(\lambda\) e igualando a zero, obtém-se:

\[\begin{align}\\ \lambda^* &= -\frac{\text{Cov}(\hat{\mu}_{MC}, \hat{\theta}_{MC})}{\text{Var}(\hat{\theta}_{MC})}\\\\ \end{align}\]

Substituindo \(\lambda = \lambda^*\) na equação da variância de \(\hat{\mu}_{CV}\), obtém-se a variância mínima associada ao estimador com variável de controle é dada por:

\[\begin{align}\\ \min_\lambda\, \text{Var}(\hat{\mu}_{CV}) &= \text{Var}(\hat{\mu}_{MC}) - \frac{\text{Cov}^2(\hat{\mu}_{MC}, \hat{\theta}_{MC})}{\text{Var}(\hat{\theta}_{MC})}\\\\ \end{align}\]

Essa expressão, em particular, pode ser reescrita em termos do coeficiente de correlação de Pearson \(\rho = \text{Corr}(\hat{\mu}_{MC}, \hat{\theta}_{MC})\), de modo que:

\[\begin{align}\\ \rho^2 = \frac{\text{Cov}^2(\hat{\mu}_{MC}, \hat{\theta}_{MC})}{\text{Var}(\hat{\mu}_{MC}) \cdot \text{Var}(\hat{\theta}_{MC})}\\\\ \end{align}\]

e, portanto,

\[\begin{align}\\ \text{Var}(\hat{\mu}_{CV}) = \text{Var}(\hat{\mu}_{MC}) \cdot (1 - \rho^2)\\\\ \end{align}\]

Essa equação evidencia o ganho em eficiência proporcionado pelo uso da variável de controle, que é tanto maior quanto mais próxima de \(\pm \, 1\) for a correlação entre os estimadores. No caso ideal em que \(\rho = \pm 1\), obtém-se \(\text{Var}(\hat{\mu}_{CV}) = 0\), o que significa que o estimador torna-se determinístico — cenário que representa a máxima eficiência possível. Na prática, a aplicação da técnica exige o conhecimento da covariância e da variância envolvidas, o que geralmente não está disponível a priori. No entanto, essas quantidades podem ser estimadas empiricamente a partir da amostra por meio das expressões:

\[\begin{align}\\ \widehat{\text{Var}}\{\hat{\theta}_{MC}\} &= \frac{1}{n(n - 1)} \sum_{i=1}^n [c(Y_i) - \bar{c}]^2\\\\ \widehat{\text{Cov}}\{\hat{\mu}_{MC}, \hat{\theta}_{MC}\} &= \frac{1}{n(n - 1)} \sum_{i=1}^n [h(X_i) - \bar{h}] \cdot [c(Y_i) - \bar{c}]\\\\ \end{align}\]

em que \(\bar{c} = \frac{1}{n} \sum_{i=1}^n c(Y_i)\) e \(\bar{h} = \frac{1}{n} \sum_{i=1}^n h(X_i)\) denotam as médias amostrais correspondentes. A partir dessas estimativas, obtém-se um valor empírico \(\hat{\lambda}\) que pode ser utilizado na construção prática do estimador \(\hat{\mu}_{CV}\). Uma situação particularmente comum ocorre quando \(X_i = Y_i\) e, portanto, as funções \(h\) e \(c\) são avaliadas sobre as mesmas variáveis simuladas. Nesses casos, é possível generalizar o procedimento para incorporar múltiplas variáveis de controle. Suponha que estejam disponíveis \(m\) funções auxiliares \(c_1, \ldots, c_m\), com valores esperados conhecidos \(\theta_1, \ldots, \theta_m\). O estimador generalizado por variáveis de controle assume a forma:

\[\begin{align}\\ \hat{\mu}_{CV} = \hat{\mu}_{MC} + \sum_{j=1}^m \lambda_j (\hat{\theta}_{MC,j} - \theta_j)\\\\ \end{align}\]

onde \(\hat{\theta}_{MC,j} = \frac{1}{n} \sum{i=1}^n c_j(X_i)\) é o estimador de Monte Carlo da \(j\)-ésima quantidade auxiliar. Essa formulação possui uma estrutura análoga à de um modelo de regressão linear múltipla, e os coeficientes \(\lambda_j\) podem ser estimados por mínimos quadrados ordinários. Esse enfoque regressivo não apenas simplifica a implementação computacional, como também permite incorporar eficientemente diversas informações auxiliares, ampliando a redução da variância do estimador de interesse.

Exemplo 4.9 (Option Pricing - Givens e Hoeting 2012). Uma opção de compra (call option) é um instrumento financeiro derivativo que concede ao titular o direito — mas não a obrigação — de adquirir uma quantidade previamente estabelecida de um ativo subjacente, tal como ações, commodities ou moedas, em uma data de vencimento específica ou até ela, por um valor predeterminado, denominado preço de exercício (strike price). Seja, por exemplo, \(S(t)\) o preço do ativo subjacente (por exemplo, uma ação) no instante de tempo \(t\). Denote-se por \(K\) o preço de exercício e por \(T\) a data de vencimento da opção. No vencimento, o titular da opção não terá incentivo econômico para exercê-la se \(K > S(T)\), pois o ativo poderia ser adquirido por um valor inferior no mercado. No entanto, se \(K < S(T)\), o titular poderá comprar o ativo por \(K\) e revendê-lo imediatamente ao valor de mercado \(S(T)\), obtendo um lucro instantâneo.

Diante desse cenário, surge uma questão fundamental: qual seria o valor justo a ser pago no instante \(t = 0\) por uma opção de compra com preço de exercício \(K\) e vencimento em \(T\)? O modelo clássico introduzido por Black, Scholes e Merton em 1973 propõe uma estrutura matemática para responder essa pergunta, baseando-se em uma equação diferencial estocástica para modelar a dinâmica do preço do ativo. Nesse modelo, o preço justo da opção é interpretado como o valor presente do retorno esperado, ajustado ao risco de mercado. Inicialmente, consideremos o caso mais simples: uma opção de compra europeia sobre uma ação que não distribui dividendos. Embora o modelo de Black-Scholes permita o cálculo analítico do preço justo dessa opção, a estimativa via simulação de Monte Carlo oferece uma alternativa instrutiva para a compreensão do problema. Segundo o modelo de Black-Scholes, o preço do ativo no vencimento \(T\) é modelado por:

\[\begin{align}\\ S(T) = S(0) \exp\left\{\left(r - \frac{\sigma^2}{2}\right)\frac{T}{365} + \sigma Z \sqrt{\frac{T}{365}}\right\}\\\\ \end{align}\]

em que \(r\) representa a taxa de retorno livre de risco (comumente associada à taxa de juros de títulos públicos com vencimento próximo a \(T\)), \(\sigma\) é a volatilidade anualizada do ativo (estimada como o desvio padrão dos log-retornos diários, isto é, \(\log(S(t+1)/S(t))\)), e \(Z \sim \mathcal{N}(0,1)\) é uma variável aleatória com distribuição normal padrão, ou seja, \(Z \sim \mathcal{N}(0,1)\). Assim, se o preço da ação no instante \(T\) fosse conhecido e igual a \(S(T)\), o valor justo da opção seria:

\[\begin{align}\\ C = \exp\left(-\frac{rT}{365}\right) \max(0, S(T) - K)\\\\ \end{align}\]

representando o retorno descontado ao valor presente. No entanto, como \(S(T)\) é desconhecido em \(t = 0\), o preço justo a ser pago corresponde ao valor esperado desse retorno descontado, ou seja, \(E[C]\). Assim, uma estimativa Monte Carlo para o preço justo da opção pode ser dada por:

\[\begin{align}\\ \bar{C} = \frac{1}{n} \sum_{i=1}^n C_i\\\\ \end{align}\]

onde os valores \(C_i\) são obtidos por simulações independentes do modelo, com uma amostra iid de variáveis \(Z_i \sim \mathcal{N}(0,1)\), para \(i = 1, \dots, n\). Como neste caso \(E[C]\) pode ser obtido analiticamente, o uso de Monte Carlo não é imprescindível. Todavia, em contextos mais complexos, como no caso das opções asiáticas (Asian Options), a simulação torna-se essencial.

As opções asiáticas têm como característica distintiva o fato de que seu retorno depende do preço médio do ativo durante o período de vigência, e não apenas do preço no vencimento. Além disso, devido à menor variância associada à média, tais opções tendem a ter preços inferiores aos das opções europeias correspondentes. Nessa configuração, técnicas de redução de variância, como o uso de variáveis de controle, desempenham papel importante na precificação via simulação de Monte Carlo. Neste caso, para precificar uma opção asiática, simula-se a trajetória do ativo ao longo dos \(T\) períodos (tipicamente dias úteis), por meio da forma discreta do modelo de Black-Scholes:

\[\begin{align}\\ S(t+1) = S(t) \exp\left\{\left(r - \frac{\sigma^2}{2}\right)\frac{1}{365} + \sigma Z(t) \sqrt{\frac{1}{365}}\right\}\\\\ \end{align}\]

em que \(Z(t) \sim \mathcal{N}(0,1)\) são variáveis normais independentes para \(t = 0, \dots, T-1\). Então, o retorno descontado da opção asiática no vencimento é dado por:

\[\begin{align}\\ A = \exp\left(-\frac{rT}{365}\right) \max\{0, \bar{S} - K\}\\\\ \end{align}\]

com \(\bar{S} = \frac{1}{T} \sum_{t=1}^T S(t)\) sendo a média aritmética dos preços simulados. Neste caso, o preço justo inicial é \(\mu = E[A]\), o qual não possui solução analítica conhecida. Portanto, a estimativa de \(\mu\) deve ser obtida via simulação de Monte Carlo, por exemplo, por meio do estimador Monte Carlo (clássico):

\[\begin{align}\\ \hat \mu_{MC} = \bar{A} = \frac{1}{n} \sum_{i=1}^n A_i\\\\ \end{align}\]

em que os valores \(A_i\) são obtidos por simulações independentes. Entretanto, ao considerar a média geométrica dos preços simulados no lugar da média aritmética, é possível derivar uma expressão analítica para \(\mu = E[A]\). O preço justo \(\theta\) da opção asiática baseada na média geométrica é então dado por:

\[\begin{align}\\ \theta &= S(0) \Phi(c_1) \exp\left\{-T\left(r - \frac{c_3\sigma^2}{6}\right)\frac{1 - 1/N}{730}\right\} - K \Phi(c_1 - c_2) \exp\left(-\frac{rT}{365}\right)\\\\ \end{align}\]

em que \(\Phi(\cdot)\) representa a função de distribuição acumulada da normal padrão e \(N\) é o número de preços incluídos na média. Os coeficientes \(c_1\), \(c_2\) e \(c_3\) são definidos como:

\[\begin{align}\\ c_1 &= \frac{1}{c_2}\left[\log\left(\frac{S(0)}{K}\right) + \left(\frac{c_3 T}{730}\right)\left(r - \frac{\sigma^2}{2}\right) + \frac{c_3\sigma^2T}{1095} + \left(1 + \frac{1}{2N}\right)\right] \\\\ c_2 &= \sigma \sqrt{\frac{c_3 T}{1095} \left(1 + \frac{1}{2N}\right)} \\\\ c_3 &= 1 + \frac{1}{N}\\\\ \end{align}\]

É importante destacar que o preço justo \(\theta\) da opção com média geométrica também pode ser estimado por simulação de Monte Carlo, análoga àquela utilizada para a média aritmética. Seja \(\hat{\theta}_{MC}\) essa estimativa. Neste caso, como os preços das opções com médias aritmética e geométrica tendem a ser altamente correlacionados, \(\hat{\theta}_{MC}\) pode ser empregada como uma variável de controle na estimação de \(\mu = E[A]\). Define-se, então, o estimador com variável de controle por:

\[\begin{align}\\ \hat{\mu}_{CV} = \hat{\mu}_{MC} + \lambda(\hat{\theta}_{MC} - \theta)\\\\ \end{align}\]

em que um valor inicial razoável para o parâmetro de controle é \(\lambda = -1\). Para ilustrar essa abordagem, considere uma opção de compra asiática com payoff baseado na média aritmética dos preços do ativo subjacente ao longo do período de detenção. Suponha que o preço atual do ativo seja \(S(0) = 100\), o preço de exercício seja \(K = 102\), a taxa de juros livre de risco seja \(r = 0{.}05\), a volatilidade anualizada seja \(\sigma = 0{.}3\) e que restam \(N = 50\) dias úteis até o vencimento, o que implica a necessidade de 50 simulações diárias para a trajetória do preço do ativo. Com base na expressão analítica para a média geométrica, obtém-se o preço justo estimado da opção como \(\theta = 1{.}83\) (Código 4.15).

Código 4.15. Cálculo do preço justo de uma opção de compra asiática com base na média geométrica dos preços do ativo subjacente ao longo de 50 dias úteis até o vencimento.

# --------------------------------------------------
# Solução Analítica: Preço Justo (Média Geométrica)
# --------------------------------------------------

# --- 1. Valores iniciais ---

S0      <- 100     # Preço em t = 0
K       <- 102     # Preço de Exercício
sigma   <- 0.3     # Volatilidade
t       <- 50      # Tempo até o Vencimento
r       <- 0.05    # Taxa de Retorno Livre de Risco

# --- 2. Cálculo do preço justo ---

N       <- t
c3      <- 1 + 1/N
c2      <- sigma * ((c3 * t/1095) * (1 + 1/(2 * N)))^0.5
c1      <- (1/c2) * (log(S0/K) + (c3 * t/730) * (r - (sigma^2)/2) +
           (c3 * (sigma^2) * t/1095) * (1 + 1/(2 * N)))

theta   <- S0 * pnorm(c1) * exp(-t * (r + c3 * (sigma^2)/6) * (1- 1/N)/730) -
           K * pnorm(c1 - c2) * exp(-r * t/365)
theta

[1] 1.827964

Utilizando agora a abordagem de Monte Carlo, é possível estimar o preço justo da opção asiática (\(\mu = E[A]\)) por meio dos estimadores \(\hat{\mu}_{MC}\) (clássico) e \(\hat{\mu}_{CV}\) (com variável de controle). Neste caso, considerando \(n = 100000\) simulações, obtêm-se as seguintes estimativas aproximadas:

• \(\hat{\mu}_{MC} = 1{.}79\) (Estimador Clássico)
• \(\hat{\mu}_{CV} = 1{.}87\) (Estimador com Variável de Controle)

Vamos agora analisar o erro-padrão — e, consequentemente, a variância — associado a cada estimador, pois essas quantidades fornecem uma medida da variabilidade das estimativas em torno do valor esperado, isto é, uma medida de precisão. Para tal, mantendo \(n = 1000\) simulações em cada execução, o experimento de Monte Carlo será replicado 100 vezes, considerando ambos os estimadores: o clássico, \(\hat{\mu}_{MC}\), e o com variável de controle, \(\hat{\mu}_{CV}\), conforme descrito no Código 4.16.

Código 4.16. Estudo comparativo entre os estimadores de Monte Carlo clássico (\(\hat{\mu}_{MC}\)) e com variável de controle (\(\hat{\mu}_{CV}\)) para o cálculo do preço justo de uma opção de compra asiática.

# ----------------------------------------------------
# Simulação Monte Carlo: Preço Justo (Opção Asiática)
# ----------------------------------------------------

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

library(knitr)

# --- 1. Valores iniciais ---

S0      <- 100     # Preço em t = 0
K       <- 102     # Preço de exercício
sigma   <- 0.3     # Volatilidade
t       <- 50      # Tempo até o vencimento
r       <- 0.05    # Taxa de retorno livre de risco
n       <- 1000    # Número de iterações
m       <- 100     # Número de repetições do experimento Monte Carlo

# --- 2.Simulação Monte Carlo ---

mu.mc       <- NULL
theta.mc    <- NULL

for(j in 1:m)
{
    A       <- NULL
    theta0  <- NULL
    
    for(i in 1:n)
    {
        ST      <- NULL
        ST[1]   <- S0
        
        for(k in 2:t)
        {
            ST[k] <- ST[k-1]*exp(((r-(sigma^2)/2)/365) + sigma*rnorm(1)/sqrt(365))
        }
        
        A[i]      <- exp(-r*t/365)*max(c(0,mean(ST) - K))
        theta0[i] <- exp(-r*t/365)*max(c(0,exp(mean(log(ST))) - K))
    }
    
    mu.mc[j]      <- mean(A)
    theta.mc[j] <- mean(theta0)
}

# --- 3. Estimadores Monte Carlo (Clássico) ---

mu_mc     <- mean(mu.mc)
sd_mu_mc  <- sd(mu.mc)

res_mu_mc <- data.frame(row.names = c('Estimativa Monte Carlo (Clássico)',
                                      'Erro-Padrão Monte Carlo (Clássico)'),
                        valores = c(mu_mc, sd_mu_mc))

kable(res_mu_mc, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('Valores')) 

	Valores
Estimativa Monte Carlo (Clássico)	1.802671
Erro-Padrão Monte Carlo (Clássico)	0.111593

# --- 4. Estimador Monte Carlo (Variável de Controle) ---

mu.cv     <- mu.mc - 1 * (theta.mc - theta)
mu_cv     <- mean(mu.cv)
sd_mu_cv  <- sd(mu.cv)

res_mu_cv <- data.frame(row.names = c('Estimativa Monte Carlo (Variáveis de Controle)',
                                      'Erro-Padrão Monte Carlo (Variáveis de Controle)'),
                        valores = c(mu_cv, sd_mu_cv))
kable(res_mu_cv, 
      align = 'c', 
      escape = FALSE, 
      format = 'html',
      col.names = c('Valores')) 

	Valores
Estimativa Monte Carlo (Variáveis de Controle)	1.8757249
Erro-Padrão Monte Carlo (Variáveis de Controle)	0.0032923

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

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

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

Ghitany, Mohamed E, Barbra Atieh, e Saralees Nadarajah. 2008. «Lindley distribution and its application». Mathematics and computers in simulation 78 (4): 493–506.

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

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

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.

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.

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.

Metropolis, N. 1987. «The Beginning of the Monte Carlo Method». Los Alamos Science, n.º 15: 125–30.

Metropolis, N., e Stanislaw Ulam. 1949. «The Monte Carlo method». Journal of the American Statistical Association 44 (247): 335–41.

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