¿Cómo establecer las priores para la estimación bayesiana de la distribución normal multivariante cuando la matriz de correlación tiene valores pequeños?

Question

Estoy estimando los parámetros de la media y la covarianza en la distribución normal multivariante (MVN). Siguiendo Wikipedia Utilicé el MVN para la media y el Inverse-Wishart para la covarianza y probé el muestreo de Gibbs.

$$ \begin{align} x_i &\sim \mathcal{N}(\mu,\Sigma), \ i=1,\ldots,N \end{align} $$ Muestreo $\mu$ : $$ \begin{align} p(\mu | x, \Sigma, \mu_0, \Sigma_0, \nu, \Psi) &\propto p(x | \mu, \Sigma) p(\mu | \mu_0, \Sigma_0) \ \cdots\ (1) \end{align} $$ Muestreo $\Sigma$ : $$ \begin{align} p(\Sigma | x, \Sigma, \mu_0, \Sigma_0, \nu, \Psi) &\propto p(x | \mu, \Sigma) p(\Sigma | \nu, \Psi) \ \cdots\ (2) \end{align} $$

El lado derecho de (1) es MVN $\times$ MVN, por lo que dibujo un nuevo $\mu$ en el $t$ iteración siguiente $$ \begin{align} \mu^{(t)} &\sim \mathcal{N} (m, s)\\ s &= (\Sigma_{0}^{-1} + N \Sigma^{-1}_{(t)} )^{-1} \\ m &= s (\Sigma_{0}^{-1}\mu_0 + N \Sigma^{-1}_{(t)} \overline{x}) \\ \overline{x} &= \frac{1}{N} \sum_{i=1}^N x_i \end{align} $$

El lado derecho de (2) es $$ \begin{align} \Sigma^{(t)} &\sim InvW (N+\nu, \eta)\\ \eta &= \Psi + \sum_{i=1}^N (x_i - \mu^{(t)}) (x_i - \mu^{(t)})^T \end{align} $$

Implementé lo anterior en Python, pero no pude recuperar los valores verdaderos después de un número suficiente de iteraciones. ¿Esto se debe a los antecedentes?

JFYI, aquí está el código Python:

import numpy as np
import scipy as sp
import scipy.stats as sps
import numpy.random as npr
import seaborn as sns
import matplotlib.pyplot as plt
npr.seed(225)

# Create Data from Multivariate Normal
N = 1000 # number of data
D = 2 # dimensions
max_mean = 0.8
max_cov = 0.15

mean_vec = npr.normal(max_mean/2, 1, D)
cov_mat = npr.uniform(max_cov/2, max_cov, (D, D))
data = npr.multivariate_normal(mean_vec, cov_mat, N)

iter_num = 10000
show_num = 9500

# Initialization
# Prepare priors
## mean
mu_0 = np.repeat(0, D)
cov_0 = np.diag(np.repeat(0.5, D))

## cov
nu = D + 1
psi = np.identity(D)

mean_itr = npr.uniform(0, max_mean*2, D)
cov_itr = npr.uniform(0.01, max_cov*2, (D, D))

# Iteration
mean_chain = []
cov_chain = []
for i in range(iter_num):
    # Update mean
    cov0_inv = np.linalg.inv(cov_0)
    cov_inv = np.linalg.inv(cov_itr)
    cov_tmp = np.linalg.inv( cov0_inv  + N * cov_inv )
    mean_tmp = cov_tmp.dot( cov0_inv.dot(mu_0) + N * np.dot(cov_inv,  data.mean(axis=0)) )
    mean_itr = npr.multivariate_normal(mean_tmp, cov_tmp, 1)

    mean_chain.append(mean_itr[0])

    # Update cov
    data_demean = data - mean_itr
    scale_tmp = psi + (data_demean.transpose()).dot(data_demean).sum(axis=0)

    cov_itr = sps.invwishart.rvs(N-1, scale_tmp)

    cov_chain.append(cov_itr)

mean_chain = np.array(mean_chain)
cov_chain = np.array(cov_chain)

# Make FIgures
dim = 1
sns.distplot(mean_chain[show_num: , dim], hist=True, kde=False)
plt.plot([mean_vec[dim], mean_vec[dim]], [0, (iter_num-show_num)*0.2], linewidth=2, color='red')

index = (1,1)
sns.distplot(cov_chain[show_num:, index[0], index[1]], hist=True, kde=False)
plt.plot([cov_mat[index[0], index[1]], cov_mat[index[0], index[1]]], [0, (iter_num-show_num)*0.2], linewidth=2, color='red')

Parece que las medias están bien, pero sobreestima la matriz de covarianza cuando los valores son pequeños (arriba: una de las medias, abajo: una de las covarianzas, las líneas rojas son los valores verdaderos).

Actualizado: Uso de Normal-Inverse-Wishart

# Prepare priors
## mean
mu_0 = np.repeat(0, D)
k0 = 0.1
cov_0 = np.diag(np.repeat(0.5, D))

## cov
v0 = D + 1.5
psi = (v0 - D - 1) * np.identity(D) # I am not sure what is proper

# Initialization
mean_itr = npr.uniform(0, max_mean*2, D)
cov_itr = sps.invwishart.rvs(v0, psi)

# Iteration
mean_chain = []
cov_chain = []
for i in range(iter_num):
    # Update mean
    data_mean = data.mean(axis=0)
    mean_tmp = (k0*mu_0 + N * data_mean ) / (k0 + N)
    k = k0 + N
    mean_itr = npr.multivariate_normal(mean_tmp, cov_itr/k, 1)

    mean_chain.append(mean_itr[0])

    # Update cov
    data_demean = data - mean_itr
    scale_tmp = psi + (data_demean.transpose()).dot(data_demean) + (k0*N)/(k0+N) *
 (data_mean-mu_0).transpose().dot(data_mean-mu_0)
    v = v0 + N

    cov_itr = sps.invwishart.rvs(v, scale_tmp)

    cov_chain.append(cov_itr)

mean_chain = np.array(mean_chain)
cov_chain = np.array(cov_chain)

Answer 1

Al observar el muestreo de la $\Sigma$ Creo que se ha olvidado de reemplazar la media $\mu$ con la muestra $\mu(t)$ de su distribución posterior de la media.

Así que $\eta$ debe determinarse como

$\eta=\Phi+\sum_n^N(x_n-\mu(t))(x_n-\mu(t))^T$

El muestreo de Gibbs es un procedimiento iterativo: cuando se muestrea a partir de la distribución posterior de una sola variable, es necesario sustituir las demás variables de las que depende la variable muestreada por su valor muestreado anteriormente.

véase https://en.wikipedia.org/wiki/Gibbs_sampling por favor