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
- Especificar la función de densidad de probabilidad para la previa
- Carga los datos de la muestra
- 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.