Reunião PPAG

1 Simulação do número de stirling

Aproximação de stirling:

\(n! \approx \sqrt{2\pi n}(\cfrac{n}{e})^n\)

stirling <- function (n)((2*pi*n)^.5)*(n/exp(1))^n

amostra = c(seq(1,3,0.001))
plot(stirling(amostra),type="l",pch=19,ylab="Valor aproximado por Stirling",xlab="Número",col="blue")
lines(factorial(amostra),type="l",pch=19,col="red")
legend(1, 5.5, legend=c("Aproximação de Stirling", "Fatorial"),
       col=c("blue", "red"), lty=1, cex=0.8)

A aproximação confirma-se visualmente, embora não tenhamos ideia numericamente do quanto bem está sendo aproximado. Para tal, utilizou-se o cálculo de um erro quadrático médio para dois tamanhos de amostra diferente.

#Criando uma amostra
amostra = seq(1,3,0.01)

#Aplicando as funcoes
st=stirling(amostra)
fc=factorial(amostra)
#Erro quadratico medio
eqm1=mean((st-fc)^2)

#Supondo N grande
amostra = seq(1,100,0.01)
st=log(stirling(amostra))
fc=log(factorial(amostra))
#Erro quadratico medio
eqm2=mean((st-fc)^2)

eqm1;eqm2
## [1] 0.009414479
## [1] 6.841241e-05

Observa-se um erro que diminui conforme o tamanho da amostra aumenta, o que é o comportamento esperado. Aqui comparamos uma amostra de 200 contra uma de 9 mil.

2 Uma possível aplicação teórica

Como visto na reunião anterior, a aproximação de laplace é utilizada para o cálculo da seguinte família de integrais:

\(\int_{a}^{b} e^{Mf(x)}\)

Exceto a distribuição normal e a exponencial (onde esta aproximação não é aconselhável), as demais distribuições de probabilidade conhecidas não possuem um núcleo unicamente dependente de um termo exponencializado.

Para além disso, uma transformação faz-se necessária. Em geral, uma transformação redundante que coloca todos os termos da expressão na base exponencial é a já conhecida \(exp(ln(f(x)))\)

Essa manipulação é muito utilizada em um assunto específico da inferência e que envolve um conjunto de modelos sob uma mesma escrita, a família exponencial.

Tese: Toda distribuição pertencente a família exponencial pode ser bem aproximada por Laplace ?

Para tal, é preciso reescrever a forma geral da família exponencial para a forma em que a aproximação de Laplace é bem vinda.

Utiliza-se como sendo pertencente a família exponencial, para uma distribuição contínua univariada qualquer, aquela que pode ser reescrita como:

\(e^{a(\theta)b(x)+c(\theta)+d(x)}\)

Tomando a integral do termo acima, e retirando as constantes, pois estamos integrando em relação da x, temos o seguinte termo:

\(\propto \int_{a}^{b}e^{a(\theta)b(x)+d(x)}dx\)

Fazendo a seguinte transformação \(x=(b-a)y \implies dx = (b-a)dy\)

Tem-se que:

\(\propto (b-a)\int_{a}^{b}e^{a(\theta)b[(b-a)y]+d[(b-a)y]}dy\)

Se e somente se \(b(y)\) e \(d(y)\) forem funções lineares ou combinações lineares, pode-se reescrever como:

\(\propto (b-a)\int_{a}^{b}e^{a(\theta)(b-a)b[y]+(b-a)d[y]}dy\)

\(\propto (b-a)\int_{a}^{b}e^{a(\theta)(b-a)(b[y]+d[y])}dy\)

Que é uma forma condizente com a aplicação da aproximação de Laplace.

Onde \(M = a(\theta)(b-a)\) e \(f(y) = -[ b(y) + d(y) ]\)

Problema ?

\(b(x)\) e \(d(x)\) dificilmente serão linares. O usual é esperar que tenha alguma relação logarítmica.

3 Laplace para distribuições a posteriori

Have you noticed that posterior distributions often are heap shaped and symmetric? What other thing, that statisticians love, is heap shaped and symmetric? The normal distribution of course! Turns out that, under most conditions, the posterior is asymptotically normally distributed as the number of data points goes to \(\infty\).

De uma forma intuitiva, a aproximação de Laplace assume uma normal com a média sendo a moda da distribuição investigada e a variância de acordo com a “curvatura” da posteriori, curvatura esta baseada na matrix Hessiana.

3.1 Modelo

\(Yi|\mu,\sigma \sim N(\mu,\sigma)\)

\(\mu \sim N(0,100)\)

\(\sigma \sim LogN(0,4)\)

Nota: A variância a log normal é \((e^{\sigma^2}-1)\times e^{2\mu-\sigma^2}\)

Outra nota: O parâmetro de dispersão da LogNormal (\(\sigma^2\)) é diferente da sua variância.

#Setando uma semente referente a data 
set.seed(10032019)

#Amostrando para y (valores reais)
y <- rnorm(n = 20, mean = 10, sd = 5)

#Média e Desvio Padrão (queremos chegar aqui)
c(mean = mean(y), sd = sd(y))
##     mean       sd 
## 8.994611 5.137709
#Criando um modelo "na mão"
model <- function(p, y) {
    log_lik <- sum(dnorm(y, p["mu"], p["sigma"], log = T))  #log likelihood
    log_post <- log_lik + dnorm(p["mu"], 0, 100, log = T) + dlnorm(p["sigma"],
        0, 4, log = T) #formato analítico da posteriori (proporcional), a princípio não reconhecível
    log_post
}

#Valores iniciais
inits <- c(mu = 0, sigma = 1)

#Ajuste baseado na otimização 
#control = list(fnscale = -1) para maximizar ao invés de minimizar
#hessian = T para retornar a matriz hessiana
fit <- optim(inits, model, control = list(fnscale = -1), hessian = TRUE, y = y)

#Valores que precisaremos para obter
#Fit$par é o máximo da distribuição de interesse e logo a média da distribuição normal que será utilizada para aproximar
param_mean <- fit$par
#A inversa da negativa da hessiana retorna a matriz de covariância (tem a ver com função score)
param_cov_mat <- solve(-fit$hessian)

round(param_mean, 2)
##    mu sigma 
##  8.99  4.87
round(param_cov_mat, 3)
##          mu sigma
## mu    1.187 0.000
## sigma 0.000 0.561
#Gerando amostras da distribuição a posteriori
samples <- rmvnorm(10000, param_mean, param_cov_mat)
samples <- cbind(samples, pred = rnorm(n = nrow(samples), samples[, "mu"], samples[,
    "sigma"]))
samples <- mcmc(samples)
summary(samples)
## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean     SD Naive SE Time-series SE
## mu    8.985 1.0868 0.010868       0.010868
## sigma 4.877 0.7396 0.007396       0.007255
## pred  8.816 5.0483 0.050483       0.050483
## 
## 2. Quantiles for each variable:
## 
##         2.5%   25%   50%    75% 97.5%
## mu     6.879 8.244 8.982  9.719 11.12
## sigma  3.392 4.382 4.880  5.374  6.33
## pred  -1.235 5.528 8.822 12.163 18.86
mu    = samples[,1]
sigma = samples[,2]

plot(mu)

plot(sigma)

3.2 Uma função útil

laplace_approx <- function(model, inits, no_samples, ...) {
    fit <- optim(inits, model, control = list(fnscale = -1), hessian = TRUE,
        ...)
    param_mean <- fit$par
    param_cov_mat <- solve(-fit$hessian)
    mcmc(rmvnorm(no_samples, param_mean, param_cov_mat))
}

y <- rnorm(n = 20, mean = 10, sd = 5)

model <- function(p, y) {
    log_lik <- sum(dnorm(y, p["mu"], p["sigma"], log = T))
    log_post <- log_lik + dnorm(p["mu"], 0, 100, log = T) + dlnorm(p["sigma"],
        0, 4, log = T)
    log_post
}


inits <- c(mu = 0, sigma = 1)
samples <- laplace_approx(model, inits, 10000, y = y)
densplot(samples)

4 Funciona mesmo?

Vamos testar a funcionalidade deste método

4.1 Exemplo prático

Geração dos dados

sigma2_verdadeiro=5
beta0_verdadeiro=-2
beta1_verdadeiro=2
n=100
x=rnorm(n,0,1)
e=rnorm(n,0,sqrt(sigma2_verdadeiro))
y=beta0_verdadeiro+beta1_verdadeiro*x+e
plot(x,y)

Aplicação

model <- function(p,y,x) {
    log_lik <- sum(dnorm(y,p["beta0"]+p["beta1"]*x,sqrt(p["sigma2"]), log = T))
    log_post <- log_lik + dnorm(p["beta0"], 0, 10, log = T) + dnorm(p["beta1"], 0, 10, log = T)       + dgamma(p["sigma2"],1,0.01, log = T) 
    log_post
}

#Valores iniciais
inits <- c(beta0 = 1, beta1 = 1, sigma2 = 1)

#Ajuste baseado na otimização 
#control = list(fnscale = -1) para maximizar ao invés de minimizar
#hessian = T para retornar a matriz hessiana
fit <- optim(inits, model, control = list(fnscale = -1), hessian = TRUE, y = y,x = x)

#Valores que precisaremos para obter
#Fit$par é o máximo da distribuição de interesse e logo a média da distribuição normal que será utilizada para aproximar
param_mean <- fit$par
#A inversa da negativa da hessiana retorna a matriz de covariância (tem a ver com função score)
param_cov_mat <- solve(-fit$hessian)

round(param_mean, 2)
##  beta0  beta1 sigma2 
##  -1.86   1.86   4.71
round(param_cov_mat, 3)
##         beta0  beta1 sigma2
## beta0   0.047 -0.002  0.000
## beta1  -0.002  0.048  0.000
## sigma2  0.000  0.000  0.442
#Gerando amostras da distribuição a posteriori
samples <- rmvnorm(10000, param_mean, param_cov_mat)
samples <- cbind(samples, pred = rnorm(n = nrow(samples), samples[, "beta0"], samples[,
    "beta1"]), samples[, "sigma2"])
samples1 <- mcmc(samples)
summary(samples1)
## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean     SD Naive SE Time-series SE
## beta0  -1.855 0.2179 0.002179       0.002179
## beta1   1.865 0.2170 0.002170       0.002170
## sigma2  4.712 0.6672 0.006672       0.006672
## pred   -1.835 1.8949 0.018949       0.018949
##         4.712 0.6672 0.006672       0.006672
## 
## 2. Quantiles for each variable:
## 
##          2.5%    25%    50%     75%  97.5%
## beta0  -2.280 -2.003 -1.857 -1.7068 -1.432
## beta1   1.444  1.716  1.863  2.0103  2.288
## sigma2  3.415  4.257  4.716  5.1605  6.018
## pred   -5.565 -3.073 -1.844 -0.5819  1.856
##         3.415  4.257  4.716  5.1605  6.018
beta0    = samples[,1]
beta1    = samples[,2]
sigma2   = samples[,3]


plot(beta0,type="l")
abline(-2,0,col="red")

plot(beta1,type="l")
abline(2,0,col="red")

plot(sigma2,type="l")
abline(5,0,col="red")

plot(density(beta0))

plot(density(beta1))

plot(density(sigma2))

4.2 Simulação

set.seed(2019)
beta0=-2
beta1=2
beta2=10
sigma2=5
tau=1/sigma2

yt=list()

for (i in 1:100){
  n=100
  x=rnorm(n,0,1)
  e=rnorm(n,0, sd = sqrt(sigma2))  
  yt[[i]]=list(beta0+beta1*x+e,x)
}


fit=list()
param_mean = list()
param_cov_mat = list()
samples = list()
amostra = list()
final = list()
b0  = list()
b1  = list()
tau = list()
media_b0 = vector()
media_b1 = vector()
media_tau = vector()
mediana_b0 = vector()
mediana_b1 = vector()
mediana_tau = vector()
var_b0 = vector()
var_b1 = vector()
var_tau = vector()
eqm_b0  = vector()
eqm_b1  = vector()
eqm_tau = vector()
hdi = list()


tic()
for (i in 1:100){
y =  yt[[i]][[1]]
x =  yt[[i]][[2]]

model <- function(p,y,x) {
  log_lik <- sum(dnorm(y,p["beta0"]+p["beta1"]*x,sqrt(p["sigma2"]), log = T))
  log_post <- log_lik + dnorm(p["beta0"], 0, 10, log = T) + dnorm(p["beta1"], 0, 10, log = T)       + dgamma(p["sigma2"],1,0.01, log = T) 
  log_post
}


inits <- c(beta0 = 0, beta1 = 0, sigma2 = 2)

fit[[i]] <- optim(inits, model, control = list(fnscale = -1), hessian = TRUE, y = y,x = x)

param_mean[[i]] <- fit[[i]]$par
param_cov_mat[[i]] <- solve(-fit[[i]]$hessian)
samples[[i]] <- rmvnorm(10000, param_mean[[i]], param_cov_mat[[i]])



b0[[i]]  = samples[[i]][,1]
b1[[i]]  = samples[[i]][,2]
tau[[i]] = 1/samples[[i]][,3]

media_b0[i]  = mean(b0[[i]]) 
media_b1[i]  = mean(b1[[i]]) 
media_tau[i] = mean(tau[[i]]) 

mediana_b0[i]  = median(b0[[i]]) 
mediana_b1[i]  = median(b1[[i]]) 
mediana_tau[i] = median(tau[[i]])  

var_b0[i]  = var(b0[[i]])
var_b1[i]  = var(b1[[i]])
var_tau[i] = var(tau[[i]])

eqm_b0[i]  = mean((b0[[i]]+2)^2)
eqm_b1[i]  = mean((b1[[i]]-2)^2)
eqm_tau[i] = mean((tau[[i]]-.2)^2)

hdi[[i]] = cbind(hdi(b0[[i]]),hdi(b1[[i]]),hdi(tau[[i]]))
}
toc()


media=c(mean(media_b0),mean(media_b1),mean(media_tau))
mediana=c(mean(mediana_b0),mean(mediana_b1),mean(mediana_tau))
var=c(mean(var_b0),mean(var_b1),mean(var_tau))
eqm=c(mean(eqm_b0),mean(eqm_b1),mean(eqm_tau))
medidas=cbind(media,mediana,var,eqm)
rownames(medidas) = c("b0","b1","tau");medidas

int = as.data.frame(cbind(t(as.data.frame((hdi))),rep(c(-2,2,0.2),100)))

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

  1. 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.

  1. 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

  1. \(\widetilde{\pi}(x_i|\mathbf{y}) = \int \widetilde{\pi}(x_i|\mathbf{\theta},\mathbf{y})\widetilde{\pi}(\mathbf{\theta}|\mathbf{y})d\theta\)

  2. \(\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})\)

e

\(\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.

Rafael Cabral Fernandez

2019-03-26