Interpretación de Algoritmo Metropolis-Hasting
Antero Flores/Alfredo Zea
07/01/2023
El Algoritmo de Metropolis-Hasting es procedimiento ampliamente utilizado para extraer muestras aproximadas de la distribución de probabilidad \(\pi\) (univariado o multivariado) encontrado un proceso de Markov \(M(x,y)\) con \(\pi\) como distribución estacionaria y utilizando el hecho conocido de que después de ejecutar el proceso de Markov durante un largo periodo de tiempo se distribuye aproximadamente como \(\pi\).
El Algoritmo de Metropolis-Hasting para obtener una muestra de \(\pi\) se resumen en los siguientes pasos:
- Proponer un estado arbitrario \(x\) como estado inicial para iniciar el proceso \(M(x,y)\)
- Proponer un nuevo estado \(y\) utilizando una matriz estocática de transición \(K(x,y)\) arbitraria
- Calcular una tasa aceptación \(R(x,y)\) tomado en cuenta el estado inicial \(x\) y el estado propuesto \(y\)
- Aceptar o rechazar el estado propuesto \(y\) comparando el valor \(R(x,y)\) con una probabilidad \(F(x,y)\) arbitraria
- Volver al paso (2) hasta culminar el número de iteraciones
Al aplicar el algoritmo se verifica la condición o ecuación de Balance Detallado, concluyendose que \(\pi\) es la distribución estacionaria el proceso de Markov \(M(x,y)\).
\[\pi(x)M(x,y) = \pi(y)M(y,x)\]
1. Núcleo o función de probabilidad \(\pi(x)\)
Considerando un caso univariado, la función de distribución de probabilidad \(\pi(x)\) cumple las siguientes condiciones:
- \(\pi(x) > 0\)
- \(\sum{\pi(x)} = 1\)
Una caraterística importante del algoritmo de Metropolis-Hasting es que se puede aplicar conocimiento sólo del núcleo de la distribución \(\pi(x)\) en situaciones en la que la constantes de normalización no es posible calcularlo.
A continuación de propone el nucleo de una distribución \(\pi(x;\theta)\) que corresponde a la suma de los núcleos de dos distribuciones normales.
# Eliminar los objetos de Global Environment
rm(list = ls())
set.seed(123)
# ------------------------------------------------------------------------------
# I. Núcleo o función de distribución π(x)
# ------------------------------------------------------------------------------
# Función o un Nucleo de distribución Normal(mu/2,sigma/2) + Normal(mu,sigma)
# Constante de Normalización K = 1/(2*sqrt(2*pi))
fKnormal2 = function(x,mu,sigma)
{
# Normal(mu,sigma) con x en R
fx = exp( - ((x - mu/2)^2/(2*(sigma/2)^2)))/(sigma/2) +
exp( - ((x - mu)^2/(2*(sigma)^2)))/(sigma)
return(fx)
}Dependiendo de los valores que asignados a los parámetros podemos obtener una distribución concreta \(\pi(x)\)
# Funciones de densidad π(x) con parametro θ
x=seq(-16,16,0.01)
par(mfrow=c(2,2))
plot(x,fKnormal2(x,12,2), type="l",main="π(x) para θ=(6,2)",xlab="x",ylab="f(x)")
plot(x,fKnormal2(x,-8,2),type="l",main="π(x) para θ=(-4,2)",xlab="x",ylab="f(x)")
plot(x,fKnormal2(x,0,2), type="l",main="π(x) para θ=(0,2)",xlab="x",ylab="f(x)")
plot(x,fKnormal2(x,6,1), type="l",main="π(x) para θ=(3,1)",xlab="x",ylab="f(x)")
par(mfrow=c(1,1))Para efectos de aplicar el algoritmo obtener la muestra de la distribución \(\pi(x)\) consideraremos \(\theta(16,4)\)
# Funcion de distribución especifica π(x)
gmu = 16
gsigma= 4
fPI = function(x){
return(fKnormal2(x,gmu,gsigma))
}
x=seq(-10 ,42,0.01)
plot(x,fPI(x), type="l",main="π(x) α N(8,2) + N(16,4)",xlab="x",ylab="π(x)")
2. Matriz estocástica de transición \(K(x,y)\)
Para aplicar el algorítmo de Metropolis-Hasting es necesario contar con una matriz estocástica \(K(x,y)\) arbitraría que permitirá mediante una transformación encontrar la matriz \(M(x,y)\) cuya distribución estacionaria es \(\pi(x)\).
Las condiciones que debe cumplir la matriz \(K(x,y)\) son las siguientes:
- El espacio de estados de \(K(x,y)\) de ser igual al dominio de la distribución \(\pi(x)\)
- \(K(x,y) > 0\) pues representa \(P(y|x)\)
- \(\sum_{i}{K(x,y_{i})} = 1\) La suma de probabilidades de pasar de x a cualquier y es 1
- \(\sum_{i}{K(x_{i},y)} = 1\) Proceso estocástico reversible
Para efectos de aplicar el algoritmo vamos a considerar una matriz \(K(x,y)\) generada por una distribución de Cauchy
\[f(x;x_{o},\gamma) = \frac{1}{\pi\gamma [ \frac{\gamma}{1 + (\frac{x - x_{0}}{\gamma})^2} ] }\]
- \(x \in \mathbb{R}\) representa a la variable aleatoria
- \(x_{0} \in \mathbb{R}\) se denomina parámetro de localización del pico de la distribución
- \(\gamma \in \mathbb{R}^{+}\) se denomina parámetro de escala y especifica el ancho medio
Creando una funcion que respresenta \(K(x,y)\) y verificando que cumpla las condiciones (a) y (b)
gscale = 5
fpK = function(x,y){
pK = dcauchy(y,location = x, scale = gscale)
return(pK)
}
fpK (10,5)## [1] 0.03183099fpK (3,-8)## [1] 0.01090102fpK (-5,-7) ## [1] 0.05488101Verificando que \(K(x,y)\) cumpla las condiciones (c) y (d)
# ------------------------------------------------------------------------------
par(mfrow=c(2,2))
y=seq(-30,30,0.01)
x = 5
plot(y,fpK(x,y), type="l",main="Densidad de K(5,y)",xlab="y",ylab="K(5,y)")
x = -4
plot(y,fpK(x,y), type="l",main="Densidad de K(-4,y)",xlab="y",ylab="K(-4,y)")
# ------------------------------------------------------------------------------
x=seq(-30,30,0.01)
y = 2
plot(x,fpK(x,y), type="l",main="Densidad de K(x,2)",xlab="x",ylab="K(x,2)")
y = -7
plot(x,fpK(x,y), type="l",main="Densidad de K(x,-7)",xlab="x",ylab="K(x,-7)")
par(mfrow=c(1,1))
# ------------------------------------------------------------------------------Asimismo se requiere una función que proponga un valor \(y\) de manera aleatoria dado un un valor \(x\) incluído en \(K(x,y)\)
# Función que retorna un valor aleatorio
frK = function(x){
rK = rcauchy(1,location = x, scale = gscale)
return(rK)
}
frK(5)## [1] 11.345653. Algoritmo de Metropolis Hasting
Encontrar una matriz estocástica \(M(x,y)\) cuya distribución de estacionaria sea \(\pi(x)\), satisfaciendose la Condición de Detalle Balanceado
\[\pi(x)M(x,y) = \pi(y)M(y,x)\]
Donde:
\[M(x,y) = \left\lbrace\begin{array}{l} K(x,y)min(1, R(x,y))~si~ x\neq y \\ K(x,y) + \sum_{z}K(x,z)\times(1-min(1,R(x,z))) ~otro~caso \end{array}\right. \] y \[R(x,y) = \frac{\pi(y)K(y,x)}{\pi(x)K(x,y)}\] Para aceptar o rechazar un valor propuesto \(y\) se utiliza una función de probabilidad \(F(x,y)\) arbitraria pudiendo ser \(Uniforme(0,1)\):
- Si \(R(x,y)\) es mayor a 1, aceptamos el valor propuesto \(y\)
- Si \(F(x,y)\) es menor que \(R(x,y)\) aceptamos el valor propuesto \(y\)
- Si \(F(x,y)\) es igual o mayor que \(R(x,y)\) rechazamos el valor propuesto \(y\)
# Numero de iteraciones
N = 500000
# periodo de quemado
L = 0.1*N
MCMC = matrix(data = 0, nrow = N, ncol = 12)
colnames(MCMC) =
c("x","y","PIx","PIy","Kxy","Kyx","Rxy","Ryx","Mxy","Myx","Fxy","Salto")
# 1. Inicial con un valor arbitrario de x del dominio de π(x)
x = runif(1,-50,50)
for (i in 1:N){
# 2. Generamos la propuesta con una distribucion arbitraria K(x,y)
# Distribución de Cauchy(location = x,scala = 10) x en R
y = frK(x)
#3. Tasa de Aceptación
PIx = fPI(x)
PIy = fPI(y)
Kxy = fpK(x,y)
Kyx = fpK(y,x)
Rxy = (PIy*Kyx) / (PIx*Kxy)
Ryx = (PIx*Kxy) / (PIy*Kyx)
# Matriz estocática cuya distribución estacinaria en π(x)
if (x!=y){
Mxy = Kxy*min(1,Rxy)
Myx = Kyx*min(1,Ryx)
}
else
{ Mxy = -1
Myx = -1
}
#4. Criterio de Aceptacion o Rechazo
#Probabilidad de aceptación,runif(1)
Fxy = runif(1,0,1)
MCMC[i,] = c(x,y,PIx,PIy,Kxy,Kyx,Rxy,Ryx,Mxy,Myx,Fxy,0)
if (Fxy < Rxy){
x = y
lsalto = 1
}
else
{ lsalto = 0
}
MCMC[i,12] = lsalto
}
mcmc = MCMC[(L+1):N,"x"]
head(mcmc,50)## [1] 7.094521 3.933804 8.484872 8.484872 8.484872 8.296230 8.432403
## [8] 8.432403 8.432403 8.432403 9.018744 9.018744 9.018744 9.018744
## [15] 11.505748 11.505748 12.877303 12.877303 12.877303 13.597334 14.745099
## [22] 14.662498 14.662498 14.003880 12.165367 12.165367 12.165367 10.025636
## [29] 10.025636 8.559446 8.559446 8.559446 8.559446 6.175836 6.175836
## [36] 6.175836 6.175836 6.175836 6.175836 10.163149 10.163149 10.163149
## [43] 10.163149 10.163149 10.163149 10.163149 9.469715 7.098403 8.537378
## [50] 9.7360854. Analisis de la muestra obtenida
par(mfrow=c(2,2))
x=seq(0 ,30,0.01)
plot(x,fPI(x), type="l",main="π(x) α N(8,2) + N(16,4)",xlab="x",ylab="π(x)")
abline(v=gmu, col='blue', lwd=3)
abline(v=gmu/2, col='red', lwd=3)
hist(mcmc,
prob = TRUE,
main = "Distribución de muestra MCMC",
xlab = "x",
ylab = "π(x)",
breaks = 200)
abline(v=gmu, col='blue', lwd=3)
abline(v=gmu/2, col='red', lwd=3)
plot( mcmc,type="l",xlab = "x", ylab = "y", main = "Traceplot de muestra MCMC")
abline(h=gmu, col='blue', lwd=3)
abline(h=gmu/2, col='red', lwd=3)
acf(mcmc,main = "Autocorrelación de muestra MCMC")
par(mfrow=c(1,1))5. Condición de Detalle Balanceado
Mostramos los valores intermedios obtenidos por el algoritmo
MCMC[(L+1):(L+12),]## x y PIx PIy Kxy Kyx
## [1,] 7.094521 3.933804 0.47226527 6.594104e-02 0.0454856662 0.0454856662
## [2,] 3.933804 8.484872 0.06594104 5.283205e-01 0.0348167163 0.0348167163
## [3,] 8.484872 12.349790 0.52832051 2.118288e-01 0.0398509083 0.0398509083
## [4,] 8.484872 36.507309 0.52832051 4.901900e-07 0.0019642527 0.0019642527
## [5,] 8.484872 8.296230 0.52832051 5.336733e-01 0.0635714877 0.0635714877
## [6,] 8.296230 8.432403 0.53367327 5.302047e-01 0.0636147927 0.0636147927
## [7,] 8.432403 16.581489 0.53020474 2.474225e-01 0.0174115664 0.0174115664
## [8,] 8.432403 74.755167 0.53020474 3.516578e-48 0.0003597772 0.0003597772
## [9,] 8.432403 4.926636 0.53020474 1.589498e-01 0.0426798655 0.0426798655
## [10,] 8.432403 9.018744 0.53020474 4.936772e-01 0.0627983843 0.0627983843
## [11,] 9.018744 39.346535 0.49367719 1.001240e-08 0.0016845806 0.0016845806
## [12,] 9.018744 19.286324 0.49367719 1.783885e-01 0.0122029634 0.0122029634
## Rxy Ryx Mxy Myx Fxy Salto
## [1,] 1.396271e-01 7.161933e+00 6.351032e-03 0.0454856662 0.03656945 1
## [2,] 8.012014e+00 1.248126e-01 3.481672e-02 0.0043455635 0.87053886 1
## [3,] 4.009474e-01 2.494093e+00 1.597812e-02 0.0398509083 0.57890195 0
## [4,] 9.278269e-07 1.077787e+06 1.822487e-09 0.0019642527 0.36720988 0
## [5,] 1.010132e+00 9.899700e-01 6.357149e-02 0.0629338639 0.79060375 1
## [6,] 9.935006e-01 1.006542e+00 6.320134e-02 0.0636147927 0.58272910 1
## [7,] 4.666547e-01 2.142912e+00 8.125189e-03 0.0174115664 0.75209307 0
## [8,] 6.632490e-48 1.507729e+47 2.386218e-51 0.0003597772 0.14798677 0
## [9,] 2.997895e-01 3.335674e+00 1.279498e-02 0.0426798655 0.51007350 0
## [10,] 9.311067e-01 1.073991e+00 5.847200e-02 0.0627983843 0.92779020 1
## [11,] 2.028128e-08 4.930656e+07 3.416545e-11 0.0016845806 0.91590227 0
## [12,] 3.613464e-01 2.767428e+00 4.409497e-03 0.0122029634 0.91531269 0Seleccionamos las columnas necesarias para verificar la condición de Detalle Balanceado
dfMCMC = as.data.frame(MCMC[(L+1):N,c("PIx","PIy","Mxy","Myx","Salto")])
head(dfMCMC)## PIx PIy Mxy Myx Salto
## 1 0.47226527 6.594104e-02 6.351032e-03 0.045485666 1
## 2 0.06594104 5.283205e-01 3.481672e-02 0.004345564 1
## 3 0.52832051 2.118288e-01 1.597812e-02 0.039850908 0
## 4 0.52832051 4.901900e-07 1.822487e-09 0.001964253 0
## 5 0.52832051 5.336733e-01 6.357149e-02 0.062933864 1
## 6 0.53367327 5.302047e-01 6.320134e-02 0.063614793 1Condición de Detalle Balanceado equivalente: \(\pi(x)M(x,y) - \pi(y)M(y,x) = 0\)
dfMCMC$Balance = dfMCMC$PIx*dfMCMC$Mxy - dfMCMC$PIy*dfMCMC$Myx
par(mfrow=c(2,1))
hist(dfMCMC$Balance,
prob = FALSE,
main = "Distribución π(x)M(x,y) - π(y)M(y,x)",
xlab = "x",
ylab = "π(x)",
breaks = 300)
plot( dfMCMC$Balance,type="l",xlab = "x", ylab = "y",
main = "π(x)M(x,y) - π(y)M(y,x)")
abline(h=0.00, col='red', lwd=3)
par(mfrow=c(1,1))Si obtenemos la media y varianzada de los valores obtenidos de \(\pi(x)M(x,y) - \pi(y)M(y,x) = 0\) podemos observar que es 0.00 con 15 decimales de precisión, verificándose la condición del Balance Detallado.
cat("π(x)M(x,y) - π(y)M(y,x)\n",
"media :" ,round(mean(dfMCMC$Balance), 15),"\n",
"varianza:" ,round(var(dfMCMC$Balance) ,15))## π(x)M(x,y) - π(y)M(y,x)
## media : 0
## varianza: 0cat("Tasa de aceptación \n",
"NumeroSaltos/TotalIteraciones :" , mean(dfMCMC$Salto) ,"\n") ## Tasa de aceptación
## NumeroSaltos/TotalIteraciones : 0.5074711