9  Simulações de Monte Carlo

using PythonCall, RCall

Tem como objetivo replicar artificialmente os dados de um modelo estatístico para estudar o comportamento de estatísticas e estimadores (ou qualquer procedimento estatístico)

9.1 Método

  1. Defina o modelo estatístico: “Seja \((X_{1},\dots,X_{n})\)” a.a. de \(X\sim f_{\theta,}\theta \in \Theta\).
  2. Escolha \(\theta_{0} \in \Theta\) e considere-o fixado daqui em diante.
  3. Para \(n\) fixado, gere (\(x_{1}, x_{2},\dots,x_{n}\)) a amostra observada de \(X\sim f_{\theta_{0}}\)
  4. Armazene a amostra observada
  5. Repita 3. e 4. \(M=10000\) vezes (ou quantas vezes desejar, a depender do caso)

9.2 Exemplo

Vamos estudar o comportamento do estimador para a variância. Como estudamos no últio capítulo sobre víes, vimos no exemplo que o estimador para a variância é viesado e que, por isso, usamos o estimador corrigido.

Vamos testar esses estimadores por simulação pelo método de Monte Carlo para entendermos a necessidade de correção.

Julia é uma linguagem de programação voltada para computação científica (CTEM). É caracterizada por sua velocidade (comparável ao C) e sintaxe simples. É a linguagem de programação favorecida por este livrete. Contudo, por ser relativamente nova, carece de alguns recursos importantes para estatística aplicada, por isso, fornecemos implementações dos algoritmos em Python e R.

using Distributions, StatsBase, Random, Printf
Random.seed!(123)

"""
  simula_montecarlo(theta0, M, n)

  Seja Xn a.a. de X sim Poisson(theta), theta > 0,
  gera M amostras aleatórias de tamanho n e armazena suas médias e variâncias
  exibe no final a média dessas estimativas. Espera-se que se aproximem de theta0
"""
function simula_montecarlo(theta0, M, n, corrigido = false)
    # Inicializa os vetores que receberão as estimativas a cada iteração
    # Ao inicializarmos o vetor previamente, evitamos alocações desnecessárias.
    estimativas_mean = zeros(M)
    estimativas_var  = zeros(M)

    # Inicializa o vetor das amostras
    amostra = zeros(n)
    dist = Poisson(theta0)

    for i in 1:M
        rand!(dist, amostra) # Funções termindas em i são mutadoras
        estimativas_mean[i] = mean(amostra)
        estimativas_var[i]  = var(amostra, corrected = corrigido) # Corrected = true divide por n-1.
    end

    m_esperada = mean(estimativas_mean)
    v_esperada = mean(estimativas_var)
    vies_observado = v_esperada - theta0
    vies_teorico = corrigido ? 0.0 : -theta0/n # O ? é o operador de tersor, como um if.
    println("\n" * "="^30)
    @printf("Config: M=%d, n=%d, Corrigido=%s\n", M, n, corrigido)
    @printf("Theta real (θ): %.2f\n", theta0)
    @printf("Média das Médias:  %.4f (Erro: %.4f)\n", m_esperada, m_esperada - theta0)
    @printf("Média das Vars:    %.4f\n", v_esperada)
    @printf("Viés Observado:    %.4f\n", vies_observado)
    @printf("Viés Teórico:      %.4f\n", vies_teorico)
end

simula_montecarlo(2, 10_000, 50)
simula_montecarlo(2, 10_000, 1000)
simula_montecarlo(2, 10_000, 50, true)

==============================
Config: M=10000, n=50, Corrigido=false
Theta real (θ): 2.00
Média das Médias:  2.0006 (Erro: 0.0006)
Média das Vars:    1.9629
Viés Observado:    -0.0371
Viés Teórico:      -0.0400

==============================
Config: M=10000, n=1000, Corrigido=false
Theta real (θ): 2.00
Média das Médias:  2.0000 (Erro: -0.0000)
Média das Vars:    1.9993
Viés Observado:    -0.0007
Viés Teórico:      -0.0020

==============================
Config: M=10000, n=50, Corrigido=true
Theta real (θ): 2.00
Média das Médias:  1.9984 (Erro: -0.0016)
Média das Vars:    1.9922
Viés Observado:    -0.0078
Viés Teórico:      0.0000

Python é uma linguagem de programação de uso geral. É extremamente popular por sua sintaxe simples e curva de aprendizado amigável. Dessa popularidade, destaca-se pelo seu ecossistema rico, contando com pacotes para aplicações diversas, para além da área de estatítisca e ciência de dados.

# O NumPy é uma biblioteca indispensável para computação matemática no Python.
import numpy as np
np.random.seed(123)

def simula_montecarlo(theta0, M, n, corrigido=False):
    # Inicializa os vetores que receberão as estimativas
    estimativas_mean = np.zeros(M)
    estimativas_var  = np.zeros(M)

    correcao = 1 if corrigido else 0

    for i in range(M):
        # Gera a amostra de tamanho n de uma poisson de parâmetro theta0
        amostra = np.random.poisson(theta0, n)
        
        estimativas_mean[i] = np.mean(amostra)
        # No NumPy, o parâmetro 'ddof' (Delta Degrees of Freedom) 
        # controla o denominador para var() e sd().
        # ddof=1 divide por n-1, ddof=0 divide por n.
        estimativas_var[i]  = np.var(amostra, ddof = correcao)

    m_esperada = np.mean(estimativas_mean)
    v_esperada = np.mean(estimativas_var)
    vies_observado = v_esperada - theta0
    
    vies_teorico = 0.0 if corrigido else -theta0 / n

    print("\n" + "=" * 30)
    print(f"Config: M={M}, n={n}, Corrigido={corrigido}")
    print(f"Theta real (θ): {theta0:.2f}")
    print(f"Média das Médias:  {m_esperada:.4f} (Erro: {m_esperada - theta0:.4f})")
    print(f"Média das Vars:    {v_esperada:.4f}")
    print(f"Viés Observado:    {vies_observado:.4f}")
    print(f"Viés Teórico:      {vies_teorico:.4f}")

simula_montecarlo(2, 10_000, 50)
simula_montecarlo(2, 10_000, 1000)
simula_montecarlo(2, 10_000, 50, corrigido=True)
# O NumPy é uma biblioteca indispensável para computação matemática no Python.
import numpy as np
np.random.seed(123)

def simula_montecarlo(theta0, M, n, corrigido=False):
    # Inicializa os vetores que receberão as estimativas
    estimativas_mean = np.zeros(M)
    estimativas_var  = np.zeros(M)

    correcao = 1 if corrigido else 0

    for i in range(M):
        # Gera a amostra de tamanho n de uma poisson de parâmetro theta0
        amostra = np.random.poisson(theta0, n)
        
        estimativas_mean[i] = np.mean(amostra)
        # No NumPy, o parâmetro 'ddof' (Delta Degrees of Freedom) 
        # controla o denominador para var() e sd().
        # ddof=1 divide por n-1, ddof=0 divide por n.
        estimativas_var[i]  = np.var(amostra, ddof = correcao)

    m_esperada = np.mean(estimativas_mean)
    v_esperada = np.mean(estimativas_var)
    vies_observado = v_esperada - theta0
    
    vies_teorico = 0.0 if corrigido else -theta0 / n

    print("\n" + "=" * 30)
    print(f"Config: M={M}, n={n}, Corrigido={corrigido}")
    print(f"Theta real (θ): {theta0:.2f}")
    print(f"Média das Médias:  {m_esperada:.4f} (Erro: {m_esperada - theta0:.4f})")
    print(f"Média das Vars:    {v_esperada:.4f}")
    print(f"Viés Observado:    {vies_observado:.4f}")
    print(f"Viés Teórico:      {vies_teorico:.4f}")

simula_montecarlo(2, 10_000, 50)
simula_montecarlo(2, 10_000, 1000)
simula_montecarlo(2, 10_000, 50, corrigido=True)

==============================
Config: M=10000, n=50, Corrigido=False
Theta real (θ): 2.00
Média das Médias:  2.0004 (Erro: 0.0004)
Média das Vars:    1.9619
Viés Observado:    -0.0381
Viés Teórico:      -0.0400

==============================
Config: M=10000, n=1000, Corrigido=False
Theta real (θ): 2.00
Média das Médias:  1.9998 (Erro: -0.0002)
Média das Vars:    1.9979
Viés Observado:    -0.0021
Viés Teórico:      -0.0020

==============================
Config: M=10000, n=50, Corrigido=True
Theta real (θ): 2.00
Média das Médias:  1.9986 (Erro: -0.0014)
Média das Vars:    1.9977
Viés Observado:    -0.0023
Viés Teórico:      0.0000

R é uma linguagem de programação voltada para a estatística e ciência de dados. É a sucessora da antes dominante S, por isso, e por ser fortemente influenciada pela filosofia GNU, é a linguagem de programação dominante para estatística na academia. Dessa forma, possui um ecossistema ainda mais rico que o Python para a estatística, com implementação de diversos modelos de ponta.

set.seed(123)

simula_montecarlo <- function(theta0, M, n, corrigido = FALSE) {
  
  # Pré-alocação de vetores: exatamente como o X = [zeros(M), zeros(M)]
  medias <- numeric(M)
  vars   <- numeric(M)
  
  for (i in 1:M) {
    # rpois gera a amostra de uma distribuição poisson de parâmetro lambda
    amostra <- rpois(n, lambda = theta0)
    
    medias[i] <- mean(amostra)
    
    # var() no R é a corrigida (n-1). 
    # Se corrigido = FALSE, ajustamos para a populacional (n).
    v <- var(amostra)
    if (!corrigido) {
      v <- v * (n - 1) / n
    }
    vars[i] <- v
  }
  
  m_esperada     <- mean(medias)
  v_esperada     <- mean(vars)
  vies_observado <- v_esperada - theta0
  vies_teorico   <- if (corrigido) 0 else -theta0 / n

  cat("\n", strrep("=", 30), "\n", sep = "")
  cat(sprintf("Config: M=%d, n=%d, Corrigido=%s\n", M, n, as.character(corrigido)))
  cat(sprintf("Theta real (theta): %.2f\n", theta0))
  cat(sprintf("Média das Médias:  %.4f (Erro: %.4f)\n", m_esperada, m_esperada - theta0))
  cat(sprintf("Média das Vars:    %.4f\n", v_esperada))
  cat(sprintf("Viés Observado:    %.4f\n", vies_observado))
  cat(sprintf("Viés Teórico:      %.4f\n", vies_teorico))
}

# Chamadas originais
simula_montecarlo(2, 10000, 50)
simula_montecarlo(2, 10000, 1000)
simula_montecarlo(2, 10000, 50, corrigido = TRUE)
set.seed(123)

simula_montecarlo <- function(theta0, M, n, corrigido = FALSE) {
  
  # Pré-alocação de vetores: exatamente como o X = [zeros(M), zeros(M)]
  medias <- numeric(M)
  vars   <- numeric(M)
  
  for (i in 1:M) {
    # rpois gera a amostra de uma distribuição poisson de parâmetro lambda
    amostra <- rpois(n, lambda = theta0)
    
    medias[i] <- mean(amostra)
    
    # var() no R é a corrigida (n-1). 
    # Se corrigido = FALSE, ajustamos para a populacional (n).
    v <- var(amostra)
    if (!corrigido) {
      v <- v * (n - 1) / n
    }
    vars[i] <- v
  }
  
  m_esperada     <- mean(medias)
  v_esperada     <- mean(vars)
  vies_observado <- v_esperada - theta0
  vies_teorico   <- if (corrigido) 0 else -theta0 / n

  cat("\n", strrep("=", 30), "\n", sep = "")
  cat(sprintf("Config: M=%d, n=%d, Corrigido=%s\n", M, n, as.character(corrigido)))
  cat(sprintf("Theta real (theta): %.2f\n", theta0))
  cat(sprintf("Média das Médias:  %.4f (Erro: %.4f)\n", m_esperada, m_esperada - theta0))
  cat(sprintf("Média das Vars:    %.4f\n", v_esperada))
  cat(sprintf("Viés Observado:    %.4f\n", vies_observado))
  cat(sprintf("Viés Teórico:      %.4f\n", vies_teorico))
}

# Chamadas originais
simula_montecarlo(2, 10000, 50)
simula_montecarlo(2, 10000, 1000)
simula_montecarlo(2, 10000, 50, corrigido = TRUE)

==============================
Config: M=10000, n=50, Corrigido=FALSE
Theta real (theta): 2.00
Média das Médias:  2.0008 (Erro: 0.0008)
Média das Vars:    1.9665
Viés Observado:    -0.0335
Viés Teórico:      -0.0400

==============================
Config: M=10000, n=1000, Corrigido=FALSE
Theta real (theta): 2.00
Média das Médias:  1.9996 (Erro: -0.0004)
Média das Vars:    1.9969
Viés Observado:    -0.0031
Viés Teórico:      -0.0020

==============================
Config: M=10000, n=50, Corrigido=TRUE
Theta real (theta): 2.00
Média das Médias:  2.0005 (Erro: 0.0005)
Média das Vars:    2.0024
Viés Observado:    0.0024
Viés Teórico:      0.0000