Reunião PPAG

1 Bugs
2 Jags
3 Stan
4 Nimble
5 INLA
- 5.1 Aproximação analítica

1 Bugs

set.seed(09042019)

x = rnorm(100,3,sqrt(4))
y = ifelse(x<0,-5,x)

sink("mod1.txt")        
cat("
    MODEL LR1 {
    #Verossimilhança 

    for(i in 1:N) {
    y[i] ~ dnorm(mu,tau)I(0,)
    }

    #Proris
     
    mu ~ dnorm(0,0.1)
    tau ~ dgamma(1,0.1)

    } ", fill = TRUE)

## 
##     MODEL LR1 {
##     #Verossimilhança 
## 
##     for(i in 1:N) {
##     y[i] ~ dnorm(mu,tau)I(0,)
##     }
## 
##     #Proris
##      
##     mu ~ dnorm(0,0.1)
##     tau ~ dgamma(1,0.1)
## 
##     }

sink()
#Obeservação: A dist normal no bugs está definido para a precisão. 
#Tamanho amostral
N = length(y)
#Conjunto de dados observados
data = list("N","y")
#Parâmetros estudados
params = c("mu", "tau")  
#Valores iniciais
inits <- function () {list(mu = 0,
                           tau = 1)}

#result <- bugs(data = data, inits = inits, parameters.to.save = params,
#              model.file = "mod1.txt", n.chains = 3, n.iter = 10000,
#              n.burnin = 1000, n.thin=9, bugs.directory = "C:\\WinBUGS14",
#              debug = TRUE, save.history = TRUE, DIC = FALSE)
#result$summary

Parâmetros	Real	Média a Posteriori	2.5%	97.5%
mu	3	2,81	2,36	3,28
tau	0,25	0,18	0,13	0,23

2 Jags

model_code = 
'model{
for(i in 1:N){
x[i] ~ dinterval(y[i],0)
y[i] ~ dnorm(mu,tau)}
mu ~ dnorm(0,0.1)
tau ~ dgamma(1,0.1)}'
  
N = length(y)

model_data = list( y = y , N = N)

model_parameters =  c("mu", "tau")

model_run = jags(data = model_data,
                 parameters.to.save = model_parameters,
                 model.file=textConnection(model_code),
                 n.chains=3, # Number of different starting positions
                 n.iter=10000, # Number of iterations
                 n.burnin=1000, # Number of iterations to remove at start
                 n.thin=9,# Amount of thinning
                 DIC=FALSE)

## module glm loaded

## module dic loaded

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 2
##    Total graph size: 206
## 
## Initializing model

model_run$BUGSoutput

## Inference for Bugs model at "5", fit using jags,
##  3 chains, each with 10000 iterations (first 1000 discarded), n.thin = 9
##  n.sims = 3000 iterations saved
##     mean  sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## mu   2.8 0.2  2.4 2.7 2.8 3.0   3.3    1  2600
## tau  0.2 0.0  0.1 0.2 0.2 0.2   0.2    1  3000
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

hdi(model_run$BUGSoutput$sims.matrix[,2])

##     lower     upper 
## 0.1388523 0.2405011 
## attr(,"credMass")
## [1] 0.95

3 Stan

model = "data {
  int<lower=0> N_obs;
  int<lower=0> N_cens;
  real y_obs[N_obs];
  real<lower=max(y_obs)> U;
}
parameters {
  real<lower=U> y_cens[N_cens];
  real mu;
  real<lower=0> sigma;
}
model {
  y_obs ~ normal(mu, sigma);
  y_cens ~ normal(mu, sigma);
  mu ~ normal(0,10)
  sigma ~ gamma(0,0.1)
}"


data <- list(y = y, x = x, N = N)
fit <- stan(model_code = model,
            data = data, iter = 1000, warmup = 100, chains = 2)

4 Nimble

simpleCode1 <- nimbleCode({
for(i in 1:N){
x[i] ~ dinterval(y[i],0)
y[i] ~ dnorm(mu,tau)}
  
mu ~ dnorm(0,10)
tau ~ dgamma(1,0.1)})

simpleModel1 <- nimbleModel(simpleCode1,
                            data = list(y = y, x = x),
                            constants = list(N = N),
                            inits = list(mu = 0, tau = 1))

mcmc.out <- nimbleMCMC(code = simpleModel1, 
                       data = list(y = y, x = x), 
                       inits = list(mu = 0, tau = 1),
                       monitors=c("mu","tau"),
                       nchains = 2, 
                       niter = 1000,
                       summary = TRUE, 
                       progressBar = TRUE)

5 INLA

O MCMC é a metodologia mais usual na modelagem Bayesiana. Auferidas as condições necessárias, ganha-se propriedades de otimalidade que nenhum outro método, nesta seara, consegue oferecer. No entanto, tais propriedades carregam consigo um custo computacional que, em situações de ordem prática, são inalcançáveis. Com o surgimento de modelos estruturalmente mais complexos e com a incorporação de bases mais robustas, tem-se a necessidade de avaliar metodologias computacionalmente menos custosas.

Uma das propostas metodológicas é o Integrated Nested Laplace Approximation, doravante INLA. Parte-se de alguns pressupostos

A classe de modelos estudada tem estrutura aditiva, isto é:

\(\eta_i = \alpha + \sum_{j=1}^{n_f}f^{(j)}(u_{ij}) + \sum_{k=1}^{n_\beta}\beta_kZ_{ki} + \epsilon_i\)

Em particular, quando atribui-se priori aos parâmetros acima, tem-se a classe de modelos chamada Gaussiano Latente.

Para os casos particulares da Estrutura Aditiva, os modelos mais bem aplicáveis são: Regressão Linear, Modelos Dinâmicos, Modelos Espaciais, Modelos Espaço-Temporais.

Embora exista a limitação de modelos, isto é, o MCMC por exemplo é mais genérico, a classe supracitada é a mais amplamente utilizada em modelagem, de uma forma geral.

Gaussian Markov random field (GMRF)

O caso particular do Campo Latente Gaussiano que admite propriedade de independência condicional e possui uma matriz de precisão esparsa. É necessário assumir estas condições para a boa aproximação e menos custo computacional.

5.1 Aproximação analítica

A chave para esta metolodogia esta nas segunintes expressões

\(\widetilde{\pi}(x_i|\mathbf{y}) = \int \widetilde{\pi}(x_i|\mathbf{\theta},\mathbf{y})\widetilde{\pi}(\mathbf{\theta}|\mathbf{y})d\theta\)
\(\widetilde{\pi}(\theta_j|\mathbf{y}) = \int\widetilde{\pi}(\theta|\mathbf{y})d\theta_{-j}\)

Onde \(x\) é o Campo Gaussiano e \(\theta\) é o vetor hiperparamétrico, em geral contendo variâncias, precisões, correlações, etc.

O Campo Gaussiano é tão somente o preditor linear como ele é definido, por exemplo:

Seja \(\eta_i = \beta_0 + \beta_ic_t + u_t + v_t\)

Onde \(\beta_i\) para qualquer i é uma componente fixa, \(c_t\) é uma variável que varia no tempo, \(u_t\) é um processo autoregressivo(1) (\(u_t = \phi u_{t-1}+e_t\)) e \(v_t\) é um termo não estruturado, um ruído. \(u_t\) e \(v_t\) são chamados de componentes aleatórias.

Logo,

\(\mathbf{x} = (\beta_0,\beta_i,\mathbf{u},\mathbf{v},\mathbf{\eta})\)

\(\mathbf{\theta} = (\phi,\sigma^2,\sigma_{v}^{2})\)

Observe que em geral a dimensão do Campo Gaussiano é grande, e em contrapartida, a dimensão do vetor hiperparamétrico costuma ser pequena.

O algoritmo INLA funciona basicamnete atribuindo distribuições as conjuntas maiores e avaliando as marginais por meio de integração numérica.

A primeira epata é calcular as marginais hiperparamétricas, mas para tal, precisa-se definir uma distribuição para \(\mathbf{\theta|y}\), para que depois seja possível integral cada marginal separadamente.

A abordagem mais simples é assumir uma aproximação de Laplace, embora esta distribuição tende a divergir da gaussianidade, sendo assim não suficientemente acurado. A abordagem mais correta, embora computacionalmente mais custosa, baseia-se na seguinte aproximação.

\(\widetilde{\pi}(\theta|\mathbf{y}) \propto \cfrac{\pi(x,\theta,y)}{\widetilde{\pi_g}(x|\theta,y)}\)

Sendo o denominador apresentado uma aproximação gaussiana, auferida pela aproximação de Laplace, avaliado na moda dos conhecidos.

A próxima etapa é calcular as marginais para os elementos do Campo Gaussiano.

Tem-se que

\(\widetilde{\pi}(x_i|\mathbf{y},\theta) = N(\mu_i(\theta),\sigma^2_i(\theta))\)

Obtidas através de uma distribuição Normal Multivariada com média sendo as modas de cada parâmetro e a matriz de variância e covariância obtida através da Matriz Hessiana.

Por fim, obtem-se:

\(\widetilde{\pi}(x_i|\mathbf{y}) = \sum_k \widetilde{\pi}(x_i|\mathbf{y},\theta_k)\widetilde{\pi}(\theta_k|\mathbf{y})\Delta_k\)

A ideia do INLA é, como foi descrito anterior, aproximando distribuições maiores por Laplace, e ruduzindo as marginais através de integração numérica.