Regresión lineal simple con LaplacesDemon
Introducción
LaplacesDemon is an open-source statistical package that is intended to provide a complete environment for Bayesian inference. LaplacesDemon has been used in numerous fields. The user writes their own model specification function and selects a numerical approximation algorithm to update their Bayesian model. Some numerical approximation families of algorithms include Laplace’s method (Laplace approximation), numerical integration (iterative quadrature), Markov chain Monte Carlo (MCMC), and Variational Bayes.
For more information please consult the webpage.
Problema
El objetivo es estimar el intercepto, la pendiente y la desviación en un modelo de regresión lineal simple.
Los datos
Vamos a simular unos datos con una estructura apropiada.
b0 <- -5
b1 <- 3
desvi <- 5
n <- 31
x <- 1:n
set.seed(666) # Semilla del diablo
y <- rnorm(n=n, mean=b0 + b1 * x, sd=desvi)
library(ggplot2)
ggplot(data=NULL, aes(x=x, y=y)) + geom_point()
Distribuciones a priori
Vamos a usar las siguientes dos a priori para \(\beta_0\), \(\beta_1\) y \(\sigma\) respectivamente.
# Hyperparameters for b0 and b1
mu0 <- 0
sd0 <- 10
curve(dnorm(x, mean=mu0, sd=sd0), from=-200, to=200,
las=1, lwd=3, ylim=c(0, 0.4),
main='Prior to b0 and b1', ylab='Density', xla='beta')
# Hyperparameters for sigma
scale0 <- 25
library(LaplacesDemon)
curve(dhalfcauchy(x, scale=scale0), from=0.01, to=150, las=1, lwd=3,
xlim=c(-50, 150),
main='Prior to sigma', ylab='Density', xlab='sigma')
Definiendo el modelo
# Model function
Model <- function(parm, Data) {
b0 <- parm[1]
b1 <- parm[2]
desvi <- exp(parm[3])
loglik <- sum(dnorm(Data$y, mean=b0+b1*Data$x, sd=desvi, log=T))
postlik <- loglik + dnorm(b0, mean=mu0, sd=sd0, log = T)
postlik <- loglik + dnorm(b1, mean=mu0, sd=sd0, log = T)
postlik <- postlik + dhalfcauchy(desvi, scale=scale0, log=T)
Modelout<-list(LP=postlik,
Dev= -2 * loglik,
Monitor=c(postlik, loglik),
yhat=rnorm(n=length(Data$x), mean=b0 + b1 * Data$x, sd=desvi),
parm=parm)
return(Modelout)
}
# Data object
MyData <- list(N = length(y),
mon.names = c("postlik", "loglik"),
parm.names = c("intercept", "slope", "standard_deviation"),
y=y, x=x)
# Initial values for the chain
Initial.Values <- c(0, 0, 1)Corriendo el modelo
library(LaplacesDemon)
N <- 10000
fit <- LaplacesDemon(Model=Model, Data = MyData,
Iterations = N, Status=1000,
Algorithm = "HARM", Thinning = 1,
Initial.Values=Initial.Values)##
## Laplace's Demon was called on Thu Mar 05 11:36:18 2020
##
## Performing initial checks...
## Algorithm: Hit-And-Run Metropolis
##
## Laplace's Demon is beginning to update...
## Iteration: 1000Iteration: 1000, Proposal: Multivariate, LP: -106.1
## Iteration: 2000Iteration: 2000, Proposal: Multivariate, LP: -106.5
## Iteration: 3000Iteration: 3000, Proposal: Multivariate, LP: -106.3
## Iteration: 4000Iteration: 4000, Proposal: Multivariate, LP: -107
## Iteration: 5000Iteration: 5000, Proposal: Multivariate, LP: -106.3
## Iteration: 6000Iteration: 6000, Proposal: Multivariate, LP: -106.8
## Iteration: 7000Iteration: 7000, Proposal: Multivariate, LP: -105.8
## Iteration: 8000Iteration: 8000, Proposal: Multivariate, LP: -107.9
## Iteration: 9000Iteration: 9000, Proposal: Multivariate, LP: -106.9
## Iteration: 10000Iteration: 10000, Proposal: Multivariate, LP: -106.3
##
## Assessing Stationarity
## Assessing Thinning and ESS
## Creating Summaries
## Estimating Log of the Marginal Likelihood
## Creating Output
##
## Laplace's Demon has finished.Para ver un resumen del modelo usamos la siguiente instrucción. En esa salida aparecerá información valiosa.
##
## #############################################################
## # Consort with Laplace's Demon #
## #############################################################
## Call:
## LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
## Iterations = N, Status = 1000, Thinning = 1, Algorithm = "HARM")
##
## Acceptance Rate: 0.2006
## Algorithm: Hit-And-Run Metropolis
## Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
## intercept slope standard_deviation
## 4.48632244 0.02587917 0.01909145
##
## Covariance (Diagonal) History: (NOT SHOWN HERE)
## Deviance Information Criterion (DIC):
## All Stationary
## Dbar 204.466 200.035
## pD 8125.768 2.850
## DIC 8330.234 202.884
## Initial Values:
## [1] 0 0 1
##
## Iterations: 10000
## Log(Marginal Likelihood): -98.73895
## Minutes of run-time: 0.01
## Model: (NOT SHOWN HERE)
## Monitor: (NOT SHOWN HERE)
## Parameters (Number of): 3
## Posterior1: (NOT SHOWN HERE)
## Posterior2: (NOT SHOWN HERE)
## Recommended Burn-In of Thinned Samples: 8000
## Recommended Burn-In of Un-thinned Samples: 8000
## Recommended Thinning: 40
## Specs: (NOT SHOWN HERE)
## Status is displayed every 1000 iterations
## Summary1: (SHOWN BELOW)
## Summary2: (SHOWN BELOW)
## Thinned Samples: 10000
## Thinning: 1
##
##
## Summary of All Samples
## Mean SD MCSE ESS LB
## intercept -3.399344 2.1179272 0.263473800 7.049428 -7.370342
## slope 2.823722 0.1583806 0.015873785 99.232765 2.609493
## standard_deviation 1.802927 0.1379452 0.008564484 457.775407 1.543450
## Deviance 204.465809 127.4815100 5.565760894 953.893960 197.774945
## postlik -109.224322 63.7391307 2.783602226 953.348937 -110.119251
## loglik -102.232904 63.7407550 2.782880447 953.893960 -103.131191
## Median UB
## intercept -3.446412 0.2279303
## slope 2.830322 3.0623950
## standard_deviation 1.801555 2.0717088
## Deviance 199.943625 206.2623821
## postlik -106.956178 -105.8765113
## loglik -99.971813 -98.8874724
##
##
## Summary of Stationary Samples
## Mean SD MCSE ESS LB
## intercept -3.489536 1.07651112 0.262467049 5.97658 -5.020017
## slope 2.835080 0.07230744 0.009768449 116.49350 2.684371
## standard_deviation 1.772194 0.14586770 0.019425666 97.25725 1.500590
## Deviance 200.034746 2.38734456 0.335057786 94.45446 197.673131
## postlik -107.005683 1.19665426 0.167751425 94.47573 -109.987894
## loglik -100.017373 1.19367228 0.167528893 94.45446 -103.021217
## Median UB
## intercept -3.633930 -0.8286269
## slope 2.843312 2.9549208
## standard_deviation 1.770342 2.0825562
## Deviance 199.088786 206.0424334
## postlik -106.526237 -105.8186965
## loglik -99.544393 -98.8365654
##
## Demonic Suggestion
##
## Due to the combination of the following conditions,
##
## 1. Hit-And-Run Metropolis
## 2. The acceptance rate (0.2006) is within the interval [0.15,0.5].
## 3. At least one target MCSE is >= 6.27% of its marginal posterior
## standard deviation.
## 4. At least one target distribution has an effective sample size
## (ESS) less than 100. The worst mixing chain is: intercept (ESS=7.049428).
## 5. Each target distribution became stationary by
## 8001 iterations.
##
## Laplace's Demon has not been appeased, and suggests
## copy/pasting the following R code into the R console,
## and running it.
##
## Initial.Values <- as.initial.values(fit)
## fit <- LaplacesDemon(Model, Data=MyData, Initial.Values,
## Covar=NULL, Iterations=4e+05, Status=632, Thinning=40,
## Algorithm="HARM", Specs=list(alpha.star=NA, B=NULL)
##
## Laplace's Demon is finished consorting.Como se obtuvo el mensaje Laplace’s Demon has not been appeased, vamos a usar la __Demonic Suggestion_ anterior y vamos a correr el modelo nuevamente con las condiciones sugeridas. El valor de Status se cambió a 100000 para no llenar la pantalla de tantos resultados. En este nuevo ajuste se usó Thinning=40 para hacer un adelgazamiento de la cadena, como Iterations=4e+05, la longitud de la cadena resultante será 4e+05/40=10000.
Initial.Values <- as.initial.values(fit)
fit <- LaplacesDemon(Model, Data=MyData, Initial.Values,
Covar=NULL, Iterations=4e+05, Status=100000, Thinning=40,
Algorithm="HARM")##
## Laplace's Demon was called on Thu Mar 05 11:36:18 2020
##
## Performing initial checks...
## Algorithm: Hit-And-Run Metropolis
##
## Laplace's Demon is beginning to update...
## Iteration: 100000Iteration: 100000, Proposal: Multivariate, LP: -110.5
## Iteration: 200000Iteration: 200000, Proposal: Multivariate, LP: -107.3
## Iteration: 300000Iteration: 300000, Proposal: Multivariate, LP: -106.9
## Iteration: 400000Iteration: 400000, Proposal: Multivariate, LP: -105.8
##
## Assessing Stationarity
## Assessing Thinning and ESS
## Creating Summaries
## Estimating Log of the Marginal Likelihood
## Creating Output
##
## Laplace's Demon has finished.Para ver un resumen del modelo usamos la siguiente instrucción. En esa salida aparecerá información valiosa.
##
## #############################################################
## # Consort with Laplace's Demon #
## #############################################################
## Call:
## LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
## Covar = NULL, Iterations = 4e+05, Status = 1e+05, Thinning = 40,
## Algorithm = "HARM")
##
## Acceptance Rate: 0.19724
## Algorithm: Hit-And-Run Metropolis
## Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
## intercept slope standard_deviation
## 5.23518541 0.01554570 0.01736806
##
## Covariance (Diagonal) History: (NOT SHOWN HERE)
## Deviance Information Criterion (DIC):
## All Stationary
## Dbar 200.621 200.621
## pD 3.271 3.271
## DIC 203.892 203.892
## Initial Values:
## [1] -5.547942 2.922841 1.846178
##
## Iterations: 4e+05
## Log(Marginal Likelihood): -100.095
## Minutes of run-time: 0.26
## Model: (NOT SHOWN HERE)
## Monitor: (NOT SHOWN HERE)
## Parameters (Number of): 3
## Posterior1: (NOT SHOWN HERE)
## Posterior2: (NOT SHOWN HERE)
## Recommended Burn-In of Thinned Samples: 0
## Recommended Burn-In of Un-thinned Samples: 0
## Recommended Thinning: 40
## Specs: (NOT SHOWN HERE)
## Status is displayed every 1e+05 iterations
## Summary1: (SHOWN BELOW)
## Summary2: (SHOWN BELOW)
## Thinned Samples: 10000
## Thinning: 40
##
##
## Summary of All Samples
## Mean SD MCSE ESS LB
## intercept -3.505119 2.2880760 0.203657535 207.0179 -7.918293
## slope 2.835234 0.1246856 0.009720149 246.2066 2.595172
## standard_deviation 1.812345 0.1317941 0.001631110 7290.6138 1.570414
## Deviance 200.621073 2.5575417 0.101663531 1061.6139 197.723811
## postlik -107.302994 1.2849666 0.051042261 1062.4692 -110.672966
## loglik -100.310537 1.2787708 0.050831766 1061.6139 -103.663346
## Median UB
## intercept -3.512904 0.9819901
## slope 2.834937 3.0862555
## standard_deviation 1.808374 2.0850710
## Deviance 199.975808 207.3266921
## postlik -106.979379 -105.8480490
## loglik -99.987904 -98.8619053
##
##
## Summary of Stationary Samples
## Mean SD MCSE ESS LB
## intercept -3.505119 2.2880760 0.203657535 207.0179 -7.918293
## slope 2.835234 0.1246856 0.009720149 246.2066 2.595172
## standard_deviation 1.812345 0.1317941 0.001631110 7290.6138 1.570414
## Deviance 200.621073 2.5575417 0.101663531 1061.6139 197.723811
## postlik -107.302994 1.2849666 0.051042261 1062.4692 -110.672966
## loglik -100.310537 1.2787708 0.050831766 1061.6139 -103.663346
## Median UB
## intercept -3.512904 0.9819901
## slope 2.834937 3.0862555
## standard_deviation 1.808374 2.0850710
## Deviance 199.975808 207.3266921
## postlik -106.979379 -105.8480490
## loglik -99.987904 -98.8619053
##
## Demonic Suggestion
##
## Due to the combination of the following conditions,
##
## 1. Hit-And-Run Metropolis
## 2. The acceptance rate (0.197245) is within the interval [0.15,0.5].
## 3. At least one target MCSE is >= 6.27% of its marginal posterior
## standard deviation.
## 4. Each target distribution has an effective sample size (ESS)
## of at least 100.
## 5. Each target distribution became stationary by
## 1 iteration.
##
## Laplace's Demon has not been appeased, and suggests
## copy/pasting the following R code into the R console,
## and running it.
##
## Initial.Values <- as.initial.values(fit)
## fit <- LaplacesDemon(Model, Data=MyData, Initial.Values,
## Covar=NULL, Iterations=4e+05, Status=632, Thinning=40,
## Algorithm="HARM", Specs=list(alpha.star=NA, B=NULL)
##
## Laplace's Demon is finished consorting.Medias de las distribuciones a posteriori
res <- fit$Summary1
post_mean_b0 <- res[1, 1]
post_mean_b1 <- res[2, 1]
post_mean_desvi <- exp(res[3, 1])
w <- cbind(c(b0, b1, desvi),
c(post_mean_b0, post_mean_b1, post_mean_desvi))
colnames(w) <- c('True', 'Posterior mean')
rownames(w) <- c("intercept", "slope", "standard_deviation")
w## True Posterior mean
## intercept -5 -3.505119
## slope 3 2.835234
## standard_deviation 5 6.124795
