Introducción.

En esta publicación se muestra como se estiman los parámetros de un modelo de Regresión Lineal Normal Bayesiano mediante el algoritmo de Metropolis.

Modelo de Regresión Normal Bayesiano.

Considere el Modelo de Regresión Lineal Simple \[ y_i=a\,x_{i1}+b+\epsilon_i\] para \(i=1,\,2,\,3,\,\ldots\,n\). Donde \(\epsilon_i\sim N(0,\sigma_\epsilon^2)\). En este caso se asume que \(a=5\), \(b=10\), \(\sigma_\epsilon=5\). Para la generación de datos en este modelo en \(R\) se utiliza el código:

trueA <- 5
trueB <- 10
trueSd <- 5
n <- 31
# create independent x-values
x <- (-(n-1)/2):((n-1)/2)
# create dependent values according to ax+b+N(0,sd)
y <- trueA * x + trueB + rnorm(n=n, mean=0, sd=trueSd)

La función de verosimilitud está dada por:

\[ \begin{align} L(a,\,b,\,\sigma_\epsilon|Datos)&=\prod^{n}_{i=1}\frac{1}{\sigma_\epsilon\sqrt{2\pi}}\exp\left(-\frac{\left(y_i-a-b\,x_i\right)^2}{2\sigma^2_\epsilon}\right)\\ &=\frac{1}{{(2\pi\sigma^2_\epsilon)}^{n/2}}\exp\left(-\frac{\sum^n_{i=1}\left(y_i-a-b\,x_i\right)^2}{2\sigma^2_\epsilon}\right) \end{align} \] y su logaritmo es:

\[\log\left(L(a,\,b,\,\sigma_\epsilon|Datos)\right)=-\frac{n}{2}\log\left(2\pi\sigma^2_\epsilon\right)-\frac{\sum^n_{i=1}\left(y_i-a-b\,x_i\right)^2}{2\sigma^2_\epsilon} \]

#Log-Likelihood
ll <- function(param){
  a <- param[1]
  b <- param[2]
  sd <- param[3]
  pred <- a * x + b
  loglikelihoods <- dnorm(y, mean=pred, sd=sd, log=T)
  sumll <- sum(loglikelihoods)
  return(sumll)
}

Para la implementación Bayesiana del modelo se asumirán las siguientes distribuciones a priori para \(a\), \(b\) y \(\sigma^2_\epsilon\): \[ \begin{align} a\sim &\, N(\mu=5,\,\sigma^2=0.1^2)\\ b\sim &\, N(\mu=10,\,\sigma^2=0.1^2)\\ \sigma^2_\epsilon\sim &\, Gamma(\alpha=4,\,\beta=1.25) \end{align} \]

Denotando \(p(a)\), \(p(b)\) y \(p(\sigma_\epsilon^2)\) como las distribuciones a priori para \(a\), \(b\) y \(\sigma^2_\epsilon\) respectivamente. Se tiene que la distribución conjunta a priori está dada por:

\[p(a,\,b,\,\sigma^2_\epsilon)=p(a)p(b)p(\sigma_\epsilon^2) \] y su logaritmo es:

\[\log\left(p(a,\,b,\,\sigma^2_\epsilon)\right)=\log(p(a))+\log(p(b))+\log(p(\sigma_\epsilon^2)) \]

log_prior <- function(param){
  a <- param[1]
  b <- param[2]
  sd <- param[3]
  aprior <- dnorm(a,mean=a,sd=0.1,log=T)
  bprior <- dnorm(b,mean=b,sd=0.1, log=T)
  sdprior <- dgamma(sd, shape=2.5, scale=2, log=T)
  return(aprior + bprior + sdprior)
}

Por tanto, la función de distribución a posteriori está dada por: \[p(a,\,b,\,\sigma^{2}_\epsilon|Datos)\propto L(a,\,b,\,\sigma_\epsilon|Datos)p(a,\,b,\,\sigma^2_\epsilon)\] y su logaritmo:

\[\log\left(p(a,\,b,\,\sigma^{2}_\epsilon|Datos)\right)\propto \log(L(a,\,b,\,\sigma_\epsilon|Datos))+\log(p(a))+\log(p(b))+\log(p(\sigma_\epsilon^2))\]

post <- function(param){
  return (ll(param) + log_prior(param))
}

Estimación de los parámetros

Para obtener los estimadores de \(a\), \(b\) y \(\sigma_\epsilon^2\) se generaran muestras aleatorias a partir de la distribución a posterior y se calculan la Media y la mediana. Dado que estas muestras se generan mediante el algoritmo de Metropolis es necesario proponer una función de distribución auxiliar mediante la cual se generen los valores aleatorios. Asumiendo que la función de distribución propuesta está dada por: \[g(a,\,b,\,\sigma^2_{\epsilon}|p_1,\,p_2\,p_3)=g(a|p_1)g(b|p_2)g(\sigma^2_{\epsilon}|p_3) \] En este caso \(a\), \(b\) y \(\sigma^2_{\epsilon}\) son independientes y \(p_1\), \(p_2\) , \(p_3\) son parámetros de localización. Las distribuciones propuestas para estos parámetros están dadas por: \[ \begin{align} a\sim &\, \sqrt{\frac{1}{4}}t_{\nu=4}+p_1\\ b\sim &\, \sqrt{\frac{1}{4}}t_{\nu=4}+p_2\\ \sigma^2_\epsilon\sim &\, \sqrt{\frac{1}{4}}t_{\nu=4}+p_3 \end{align} \] donde \(a\), \(b\) y \(\sigma^2_\epsilon\) pueden ser generados en \(R\) mediante:

proposalfunction <- function(param){
  return(c(0.5^1/2*sqrt(2/4)*rt(3,df=4)+param))
}

El algoritmo de Metropolis en este caso está dado por: + Se parte de un vector arbitrario \((a_0, b_0, \sigma^2_{\epsilon0})\) + Se genera un vector aleatorio a partir de \(g(a,\,b,\,\sigma^2_{\epsilon}|a_0,\, b_0,\, \sigma^2_{\epsilon0})\). Este vector se denota como \(a^t,\, b^t,\, {\sigma^{2t}_{\epsilon}}\) + Se calcula \(k=p(a^t,\,b^t,\,\sigma^{2t}_\epsilon|Datos)/p(a_0,\,b_0,\,\sigma^{2}_{\epsilon0}|Datos)\) + El nuevo valor de la cadena se decide mediante:

\[(a^{t+1},\, b^{t+1},\, {\sigma^{2{(t+1)}}_{\epsilon}}) = \left\{ \begin{array}{ll} (a^t,\, b^t,\, {\sigma^{2t}_{\epsilon}}) & \mbox{con probabilidad $\min(k,1)$};\\ (a_0, b_0, \sigma^2_{\epsilon0}) & \mbox{con probabilidad $1-\min(k,1)$}.\end{array} \right. \]

se ejecuta mediante la siguiente linea de código:

run_metropolis_MCMC <- function(startvalue, iterations){
  chain <- array(dim=c(iterations + 1, 3))
  chain[1, ] <- startvalue
  for (i in 1:iterations){
    proposal <- proposalfunction(chain[i, ])
    probab <- exp(post(proposal) - post(chain[i, ]))
    if (runif(1) < probab) {
      chain[i+1, ] <- proposal
    }
    else chain[i+1, ] <- chain[i, ]
  }
  return(chain)
}


startvalue <- c(4.8, 9.8, 4.8)
chain <- run_metropolis_MCMC(startvalue, 10000)
burnIn <- 5000


mod <- lm(y ~ x)
res <- cbind(c(trueA, trueB, trueSd),
             colMeans(chain[-(1:burnIn),]),
             apply((chain[-(1:burnIn),]), 2, median),
             c(coef(mod)[2:1], summary(mod)$sigma))
colnames(res) <- c('True', 'Means', 'Medians', 'lm')
rownames(res) <- c('a', 'b', 'sd')

Y los estimadores bayesianos y clásicos están dados por:

res
##    True     Means   Medians        lm
## a     5  5.030606  5.028109  5.022448
## b    10 10.765338 10.787916 10.702239
## sd    5  6.823049  6.752740  6.710297

Se puede observar que la aproximación bayesiana es mejor que la clásica.