El papel y el lápiz ejemplo de ajuste de la distribución normal de los datos con MCMC

Question

He estado tratando de entender la Cadena de Markov Monte Carlo métodos durante un tiempo y aunque yo un poco la idea, cuando me llega la aplicación de MCMC, no estoy seguro de lo que debo hacer. Muchas veces he recibido la respuesta "de usar un paquete" de parte de los profesores, pero no quiero usar un paquete, quiero hacerlo yo mismo, así que puedo entender! ;)

He estado tratando de buscar en internet para los buenos ejemplos, pero hasta ahora no he encontrado ninguno de esos ejemplos que explícitamente se muestra paso a paso lo que está sucediendo bajo el capó. Muchos de los ejemplos que he visto de tratar con un poco abstracto ejemplos que son difíciles de comprender por completo. En muchos de los ejemplos de lo que me queda es: "Bueno muy buen método, ahora quiero escribir un MCMC programa de mí mismo, que me encaja por ejemplo, la distribución normal a los datos". El problema es que no puedo programa de conceptos abstractos, necesito un ejemplo concreto. Algunos ejemplos que he encontrado perderse en los detalles del dominio del problema y se olvidan de explicar la MCMC el método en sí.

Así que mi pregunta es: ¿me Pueden brindar con un simple papel y lápiz ejemplo? Esto podría ser, por ejemplo, un problema en el que tiene tres puntos de datos, decir $(4,5), (3,3)$ $(4,2)$ en el plano xy y que necesita para adaptarse a estos datos a la distribución normal. Vamos a suponer también que, por alguna razón, necesitamos aplicar MCMC para este problema (de Metropolis-Hastings sería bueno). Me muestran los pasos detallados de cómo hacerlo, teniendo los derivados, igualando a cero, la elección de la propuesta de distribución, tomando una muestra aleatoria de los posteriores (si utiliza el método Bayesiano) et cetera. Ello, como una o dos iteraciones, sólo para mostrar a todos cómo los cálculos de proceder.

Creo que este ejemplo podría proporcionar información valiosa para los principiantes en métodos MCMC. Gracias!

"Oigo y olvido. Veo y recuerdo. Yo lo hago y entiendo."

Confucio

Answer 1

Aunque no es estrictamente de "lápiz y papel" me imagino que usted está planeando este código en algún lenguaje, así que voy a escribir un ejemplo en el código de python, pero esperemos que esta escrito de tal forma que he de sacrificar la eficiencia y la limpieza del código a cambio de algo que es fácil de seguir y traducir a un idioma de su elección - esto es no un buen código, pero espero que sea legible el código!

Así que primero básico de las importaciones, especificando la previa de las distribuciones y sus hyperparameters así como la carga de los datos y la especificación de la log-verosimilitud. En este bloque de código que nos

1. Especificar la función de densidad de probabilidad para la previa
2. Carga los datos de la muestra
3. Especificar el modelo loglikelihood

import numpy as np
from scipy.stats import multivariate_normal, wishart

#######
# Change this section to use your own priors
# with suitable hyperparameters

def mu_prior_logpdf(m):
    m0 = np.zeros(2)           
    C0 = np.diag(np.ones(2)) 
    return multivariate_normal.logpdf(m, mean=m0, cov=C0)

def S_prior_logpdf(S):
    df0 = 2
    V0 = np.diag(np.ones(2))
    return wishart.logpdf(S, df=df0, scale=V0)

#######

# Your sample data 
nData = 3
X = np.array([[4, 5], [3, 3], [4, 2] ])

###
# Your specified model
# 
#     X_i ~ N(m,S), with X_i i.i.d

def loglik(X, m, S):
    ll = 0.
    for n in range(nData):
        ll += multivariate_normal.logpdf( X[n, ], mean=m, cov=S )
    return ll

Paso 2: Especificar la propuesta de distribución de

Hay varias opciones en este punto con respecto a lo algoritmo MCMC para elegir, en particular, hay buenas razones para el uso de muestreo de Gibbs - sin embargo , porque usted quiere una manera bastante general método vamos a demostrar una Metropolis-Hastings algoritmo MCMC que no requiere ningún conocimiento especial de las distribuciones condicionales.

Recordar para MH-MCMC necesitamos especificar una propuesta de distribución de $q(\theta^* | \theta^{(n)} )$, lo que da la distribución condicional de una nueva muestra $\theta^*$ condicional en el valor actual de $\theta^{(n)}$.

Ahora algunos pseudo-código para la propuesta de las distribuciones, de nuevo debe ser claro cómo modificar todo esto, en este ejemplo en particular la propuesta de $q(\mu^* | \mu^{(n)})$ es administrado por un paseo aleatorio centra en el valor actual. Mientras que la propuesta de la covarianza es sólo una base independiente sampler $q(\Sigma^* | \Sigma^{(n)}) = q(\Sigma^*)$, esto es, probablemente, para mostrar la convergencia pobres y muestreo de las propiedades, pero de nuevo esto es sólo para demostrar que el conjunto básico de seguridad

# Proposal distributions: 
#   - return a random variable
#   - able to evaluate the pdf of q(z' | z)

# mean proposal is simple random walk proposal 
# centered at the current value so need to specify the 
# std. dev of the step size
scale = 0.1
mpropCov = np.array([[ scale*scale, 0.], [0., scale*scale]])

def mu_proposal_rvs(mcur, size=1):
    d = np.random.normal(size=2, scale=scale) # pair of ind. mean 0 steps
    return mcur + d

def mu_proposal_logpdf(mnew, mcur):
    return multivariate_normal.logpdf(mnew, mean=mcur, cov=mpropCov)

# proposal for the covariance, just going to put a simple independent 
# sampler here, so change this for something more suitable

def S_proposal_rvs(scur, size=1):
    return wishart.rvs(df=2, scale=np.diag(np.ones(2)))

def S_proposal_logpdf(Snew, Scur):
    return wishart.logpdf(Snew, df=2, scale=np.diag(np.ones(2)))

Paso 3: Solo MCMC paso

Y ahora, un solo MH paso para que sus datos pueden proceder de la siguiente manera; dibujar una variable aleatoria a partir de la propuesta y, a continuación, acepta con una probabilidad de $$Un(\theta^*, \theta^{(n)}) = \min \left(1, \frac{p(\mathbf{X}|\theta^*)\pi(\theta^{*}) p(\theta^{(n)}|\theta^*)}{p(\mathbf{X}|\theta^{(n)})\pi(\theta^{(n)})q(\theta^* | \theta^{(n)})} \right)$$

def MHstep(mcur, Scur):

    # propose a new value of the mean parameter
    mnew = mu_proposal_rvs(mcur)

    # evaluate the acceptance ratio
    # note that since we are not updating S at this point the prior 
    # for S and the proposal of S will cancel in the numerator and
    # denominator of A so we don't need to evaluate it here

    logAnum = loglik(X, mnew, Scur) + mu_prior_logpdf(mnew) + mu_proposal_logpdf(mcur, mnew)
    logAden = loglik(X, mcur, Scur) + mu_prior_logpdf(mcur) + mu_proposal_logpdf(mnew, mcur)
    A = np.exp( logAnum - logAden )

    # Now we accept the proposed value of the new mean
    # with probability A
    U = np.random.uniform()
    if U <= A:
        # Accept the proposed update
        mcur = mnew
    else:
        # do nothing
        pass

    # Now repeat the steps to update the covariance term with the 
    # obvious changes
    # 
    # ...
    # TO DO!
    # ... 

    return mcur, Scur

And repeat...

So to get a sample of size $$ N de la parte posterior de todo lo que necesitamos hacer ahora es el suministro de ciertas condiciones iniciales y de inicio de muestreo, también se puede modificar de la manera obvia para agregar una quemadura y mantener sólo un subconjunto de la parte posterior de muestras, etc.

# Initialise some parameters for the chain
mcur = np.array([0., 0.])
Scur = np.array([[1., 0.], [0., 1.]])

meanRVs = []
covRVs = []

nt = 0; N = 10
while nt < N:
    mcur, Scur = MHstep(mcur, Scur)
    meanRVs.append( mcur )
    covRVs.append( Scur )
    nt += 1

Usted puede ahora modificar el código anterior con algunos print ... declaraciones, o incluso agregar usuario de entrada/salida lenta cada paso hacia abajo y ver exactamente lo que está sucediendo.