1ª Escola de Verão - GET/IME/UFF
Parte 1

1 Modelagem Estatística

  • Pode ser usada em muitas áreas diferentes, incluindo finanças, medicina, marketing, ciências sociais e muitas outras.

  • É o processo de usar técnicas estatísticas para descrever, explicar, identificar padrões e tendências e realizar previsões de fenômenos para ajudar na tomada de decisões.

  • Cria-se algum modelo matemático para explicar e analisar dados observados.

  • Para propor o(s) modelo(s), realiza-se análises exploratórias de dados, conhecimento prévio do pesquisador e estudos já realizados na literatura.

  • Em geral, propõe-se mais um modelo e então há a necessidade de selecionar o melhor modelo de acordo com o ajuste e/ou com a previsão.

  • Exemplo de um modelo estatístico para uma variável de interesse:

\[Y_i \sim N(\beta, V), \quad i=1, \ldots, N, \quad V>0.\]

Considere que \(\boldsymbol{\theta} = (\beta, V)\) seja o vetor paramétrico desconhecido. Para estimar este vetor paramétrico, utiliza-se alguma amostra desse fenômeno.

Com base no modelo acima, acredita-se que o histograma dos dados observados e que a densidade da variável de interesse possuem as seguintes formas:

2 Inferência Bayesiana

  • É um método estatístico que permite incorporar crença ao parâmetro desconhecido do modelo.

  • É uma abordagem alternativa à inferência clássica (frequentista).

  • Na inferência bayesiana, a estimação sobre os parâmetros é realizada usando o teorema de Bayes:

\[\begin{equation} p(\boldsymbol{\theta}|\boldsymbol{y}) = \frac{p(\boldsymbol{y}|\boldsymbol{\theta}) p(\boldsymbol{\theta})}{p(\boldsymbol{y})}, \label{eq:teobayes} \end{equation}\] sendo \(\boldsymbol{y} = (y_1, \ldots, y_n)\) uma amostra aleatória, \(\boldsymbol{\theta}\) o vetor paramétrico a ser estimado e

\(p(\boldsymbol{y}|\boldsymbol{\theta})\) a função de verossimilhança,

\(p(\boldsymbol{y})\) a distribuição marginal dos dados e

\(p(\boldsymbol{\theta})\) a distribuição a priori dos parâmetros desconhecidos.

  • No exemplo dado anteriormente, tem-se que atribuir uma distribuição a priori para \(\boldsymbol{\theta} = (\beta, V)\) e obter-se então a distribuição a posteriori deste vetor.

Considere que \(\beta\) e \(V\) sejam independentes, \(\beta \sim N(m_\beta, V_\beta)\) e \(V \sim GI(a_V, b_V)\). Tem-se então que a distribuição a posteriori é dada por

\[\begin{align} p(\beta, V|\boldsymbol{y}) &= \frac{p(\boldsymbol{y}|\beta, V) p(\beta, V)}{p(\boldsymbol{y})}\\ &\propto \prod_{i=1}^n\left\{\frac{1}{\sqrt{2\pi V}} \exp\left[-\frac{1}{2V}(y_i - \beta)^2\right]\right\} \times\\ &\frac{1}{\sqrt{2\pi V_\beta}} \exp\left[-\frac{1}{2V_\beta}(\beta - m_\beta)^2\right] \frac{b_V}{\Gamma(a_V)} V^{-(a+1)} \exp\left\{-b/V\right\} \\ &\propto V^{-(\frac{n}{2}+a_V+1)} \exp\left[-\frac{1}{V}\left( b_V + \frac{\sum_{i=1}^n {(y_i - \beta)^2}}{2}\right)\right] \exp\left[-\frac{1}{2V_\beta}(\beta - m_\beta)^2\right] \label{eq:post} \end{align}\]

A distribuição a posteriori de \((\beta, V)\) é desconhecida e neste caso pode-se recorrer aos Métodos de Monte Carlo via cadeias de Markov (MCMC) (Gamerman e Lopes (2006)).

Usando o Amostrador de Gibbs, proposto por Geman e Geman (1984) e introduzido a comunidade estatística por Gelfand e Smith (1990), realiza-se o seguinte algoritmo:

  1. Inicialize o contador em \(j = 1\) e determine valores arbitrários para \[\begin{eqnarray} \boldsymbol{\theta} ^{(1)} = (\beta^{(1)}, V^{(1)}). \end{eqnarray}\]

  2. Modifique o contador de \(j\) para \(j + 1\);

  3. Obtenha um novo valor para \(\boldsymbol{\theta}^{(j)}\) a partir de \(\boldsymbol{\theta}^{(j-1)}\) sequencialmente da seguinte forma \[\begin{eqnarray} \beta^{(j)} &\sim & p(\beta \mid V^{(j-1)},\boldsymbol{y})\nonumber \\ V^{(j)} & \sim & p(V \mid \beta^{(j)},\boldsymbol{y})\nonumber \end{eqnarray}\]

  4. Repita os passos (2) e (3) até que a cadeia convirja.

As distribuições \(p(\beta \mid V,\boldsymbol{y})\) e \(p(V \mid \beta,\boldsymbol{y})\) são chamadas de distribuições condicionais completas a posteriori.

As distribuições condicionais completas a posteriori são: \[\begin{align} \beta |V, \boldsymbol{y} &\sim N\left(m^{*}, V^{*}\right),\\ V |\beta, \boldsymbol{y} &\sim GI\left(\frac{n}{2}+a_V+1, b_V + \frac{\sum_{i=1}^n {(y_i - \beta)^2}}{2}\right) \label{eq:dccp} \end{align}\] sendo \(V^{*} = \left(\frac{n}{V}+\frac{1}{V_\beta}\right)^{-1}\) e \(m^{*} = V^{*}\left(\frac{m_\beta}{V_\beta} + \frac{\sum_{i=1}^n{y_i}}{V}\right)\).

Caso alguma distribuição condicional completa a posteriori fosse desconhecida, poderia-se recorrer ao Metropolis-Hastings (Metropolis et al. (1953)).

Para mostrar a inferência deste modelo no R, antes criaremos um dado simulado usando os seguintes comandos:

#Simulando um conjunto de dados
n     = 1000  #tamanho da amostra
beta  = 10    #média populacional
V     = 4     #variância populacional
y     = rnorm(n, beta, sqrt(V)) #gerando a amostra

#Visualizando os valores amostrados e sua respectiva densidade
hist(y, freq=F,bty="n", main="", ylab="densidade",xlab="y",cex.lab=1.4,cex.axis=1.4)
curve(dnorm(x, beta, sqrt(V)), lwd=2, col="red", add=T)

Suponha agora que os valores de \(\beta\) e \(V\) são desconhecidos. Precisamos atribuir valores para os hiperparâmetros \(m_\beta, V_\beta, a_V\) e \(b_V\). Considere os seguintes valores:

\[m_\beta = 0, \quad V_\beta = 1000, \quad a_V=2, \quad b_V=1.\]

Dessa forma, para obter uma amostra deste vetor paramétrico usando o algoritmo acima, pode-se usar os seguintes comandos no R:

#Definindo as variaveis necessarias para o MCMC
n.aquec   = 100     #numero de iteracoes para o aquecimento
n.espac   = 1       #numero do espacamento
n.post    = 1000    #tamanho da amostra
n.ite     = n.aquec + n.post*n.espac #total de iteracoes

#Definindo os vetores que armazenarão a amostra da distribuição a posteriori
beta.post = NULL
V.post    = NULL

#Inicializando a cadeia
beta.post[1]  = 0
V.post[1]     = 1

#Definindo os hiperparametros (= parâmetros da distribuição a priori)
m_beta  = 0
V_beta  = 1000
a_V     = 2
b_V     = 1

#Executando o MCMC
for(i in 2:n.ite){
  #p(beta|V^(i-1), y1, ..., yn)
  Vast          = 1/(n/V.post[i-1] + 1/V_beta)
  mast          = Vast * (sum(y)/V.post[i-1] + m_beta/V_beta)
  beta.post[i]  = rnorm(1, mast, sqrt(Vast))

  #p(V|beta^(i), y1, ..., yn)
  V.post[i]     = 1/rgamma(1, n/2 + a_V + 1, b_V + 0.5*sum((y - beta.post[i])^2))
}


#Rodando uma segunda cadeia
#Definindo os vetores que armazenarão a amostra da distribuição a posteriori
beta.post.2 = NULL
V.post.2    = NULL

#Inicializando a cadeia
beta.post.2[1]  = 100
V.post.2[1]     = 100

#Executando o MCMC
for(i in 2:n.ite){
  #p(beta|V^(i-1), y1, ..., yn)
  Vast          = 1/(n/V.post.2[i-1] + 1/V_beta)
  mast          = Vast * (sum(y)/V.post.2[i-1] + m_beta/V_beta)
  beta.post.2[i]  = rnorm(1, mast, sqrt(Vast))

  #p(V|beta^(i), y1, ..., yn)
  V.post.2[i]     = 1/rgamma(1, n/2 + a_V + 1, b_V + 0.5*sum((y - beta.post.2[i])^2))
}

Analisando a convergência

Para avaliar a convergência, analisamos o traço da cadeia e o gráfico de autocorrelação.

O traço da cadeia mostra como os valores da cadeia variam ao longo das iterações. Ele permite que você verifique se a cadeia está se movendo em torno da região de maior densidade de probabilidade e se está estabilizada.

O gráfico de autocorrelação mostra como as amostras da cadeia estão autocorrelacionadas entre si. Pode ser usado como um indicador da independência das amostras. Se as amostras estiverem altamente correlacionadas, pode ser necessário aumentar o tamanho da amostra ou reduzir o comprimento do espaçamento.

#Para decidir sobre o espaçamento e sobre o aquecimento
#Plotando o traço da cadeia
par(mfrow=c(2,2))
plot(beta.post, ylab=expression(beta), cex.lab=1.4, cex.axis=1.4, lwd=2,
     type="l", xlab="iterações")
lines(beta.post.2, col="red",lwd=2)

plot(V.post, ylab="V", cex.lab=1.4, cex.axis=1.4, lwd=2, type="l", 
     xlab="iterações")
lines(V.post.2, col="red",lwd=2)

#Gráfico de autocorrelação
par(mfrow=c(1,2))

acf(beta.post, cex.lab=1.4, cex.axis=1.4, main="",
     xlab=expression(beta), ylab="ACF")

acf(V.post, cex.lab=1.4, cex.axis=1.4, main="",
     xlab="V", ylab="ACF")

#Após definirmos o período de aquecimento e espaçamento, avaliamos 
#novamente o traço da cadeia e a autocorrelação. Em seguida,
#fazemos a inferência com a amostra resultante.

# Criando uma sequência de índices para remover as iterações
# referentes ao período de aquecimento e as do espaçamento
seq.post = seq(n.aquec,n.ite,by=n.espac)
#seq.post

#autocorrelação
par(mfrow=c(1,2))
acf(beta.post[seq.post], cex.lab=1.4, cex.axis=1.4, main="",
     xlab=expression(beta), ylab="ACF")

acf(V.post[seq.post], cex.lab=1.4, cex.axis=1.4, main="",
     xlab="V", ylab="ACF")

#Plotando o traço da cadeia
par(mfrow=c(2,2))
plot(beta.post[seq.post], ylab=expression(beta), cex.lab=1.4, cex.axis=1.4, lwd=2,
     type="l", xlab="iterações")
lines(beta.post.2[seq.post], col="red",lwd=2)

plot(V.post[seq.post], ylab="V", cex.lab=1.4, cex.axis=1.4, lwd=2, type="l", 
     xlab="iterações")
lines(V.post.2[seq.post], col="red",lwd=2)

#Gráfico de densidade
#Mostra a distribuição da amostra. Ele permite que você verifique 
#se a amostra está se aproximando da distribuição de equilíbrio 
#e se há vários picos na distribuição.

hist(beta.post[seq.post], freq=F, cex.lab=1.4, cex.axis=1.4, main="",
     xlab=expression(beta), ylab="densidade")
abline(v=beta,lwd=2, lty=3,col="red") #plotando o valor verdadeiro
abline(v=quantile(beta.post[seq.post],c(0.025,0.975)),lwd=2, lty=3,col="blue") #plotando o IC de 95%
abline(v=mean(beta.post[seq.post]),lwd=2, lty=3,col="purple") #plotando a média a posteriori

hist(V.post[seq.post], freq=F, cex.lab=1.4, cex.axis=1.4, main="",
     xlab="V", ylab="densidade")
abline(v=V,lwd=2, lty=3,col="red") #plotando o valor verdadeiro
abline(v=quantile(V.post[seq.post],c(0.025,0.975)),lwd=2, lty=3,col="blue") #plotando o IC de 95%
abline(v=mean(V.post[seq.post]),lwd=2, lty=3,col="purple") #plotando a média a posteriori

Diagnóstico de Gelman-Rubin

Também conhecido como diagnóstico de convergência de cadeias múltiplas, é uma técnica para avaliar se as cadeias de Markov geradas por um algoritmo MCMC convergiram para a mesma distribuição estacionária. O diagnóstico envolve o cálculo de uma estimativa do fator de escala da variância entre cadeias e dentro de cadeias. Se as cadeias convergiram para a mesma distribuição estacionária, então a variância entre cadeias deve ser aproximadamente igual à variância dentro das cadeias, e o fator de escala da variância deve ser próximo a 1. Caso contrário, se as cadeias não convergiram, então o fator de escala da variância será maior do que 1.

#Analisando a convergência usando o pacote CODA
if(!require(coda)){ install.packages("coda"); require(coda)}
## Carregando pacotes exigidos: coda
mcmc_obj = mcmc.list(as.mcmc(beta.post[seq.post]), as.mcmc(beta.post.2[seq.post]))

gelman.diag(mcmc_obj)
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1       1.01

Podemos recorrer a algum pacote para ajustar o modelo acima. Usando o pacote rstanarm, fazemos isso da seguinte forma:

if(!require(rstanarm)){ install.packages("rstanarm"); require(rstanarm)}

x       = rep(1,n)
dados   = data.frame(y=y, x=x)
modelo  = stan_glm(y ~ x -1, family = gaussian, data=dados)

modelo$coefficients
summary(modelo)
if(!require(ggplot2)){ install.packages("ggplot2"); require(ggplot2)}
if(!require(bayesplot)){ install.packages("bayesplot"); require(bayesplot)}
mcmc_trace(modelo)
mcmc_acf(modelo)

3 Exercício

  1. Simular um conjunto de dados do modelo de regressão linear simples dado pela seguinte expressão \[Y_i \sim N(\beta_0 + x_i \times \beta_1, V), i=1, \ldots, n.\] Para isso suponha que \(x_i \sim N(0,1), \forall i\), \(\beta_0 = 10\), \(\beta_1 = 2\) e \(V=4\).

  2. Estimar \(\boldsymbol{\theta} = (\beta_0, \beta_1, V)\).

As DCCPs são \[\begin{align} \boldsymbol{\beta} |V, \boldsymbol{y} &\sim N\left(\boldsymbol{m}^{*}, \boldsymbol{V}^{*}\right),\\ V |\boldsymbol{\beta}, \boldsymbol{y} &\sim GI\left(\frac{n}{2}+a_V+1, b_V + \frac{(\boldsymbol{y} - \boldsymbol{x}\boldsymbol{\beta})^{T}(\boldsymbol{y} - \boldsymbol{x}\boldsymbol{\beta})}{2}\right) \label{eq:dccp2} \end{align}\] sendo \(\boldsymbol{V}^{*} = \left(\frac{\boldsymbol{x}^T\boldsymbol{x}}{V}+\boldsymbol{V_\beta}^{-1}\right)^{-1}\) e \(m^{*} = \boldsymbol{V}^{*}\left(\boldsymbol{V_\beta}^{-1} \; m_\beta + \frac{\boldsymbol{x}^{T}\boldsymbol{y}}{V}\right)\).

Referências

Gamerman, Dani, e Hedibert F Lopes. 2006. Markov chain Monte Carlo: stochastic simulation for Bayesian inference. CRC Press.
Gelfand, A. E., e A. F. M. Smith. 1990. "Samping-Based Approaches to Calculating Marginal Densities". Journal of the American Statistical Association 85 (410): 398–409.
Geman, Stuart, e Donald Geman. 1984. "Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images". IEEE Transactions on pattern analysis and machine intelligence, nº 6: 721–41.
Metropolis, Nicholas, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, e Edward Teller. 1953. "Equation of state calculations by fast computing machines". The journal of chemical physics 21 (6): 1087–92.