机器算法验证 - Normal-Inverse-Wishart 的后协方差未正确收敛 - 吾爱随笔录

我正在尝试为 numpy/scipy 中具有未知均值和协方差的多元正态分布实现简单的正态-逆-Wishart 共轭先验分布，以便它可以采用数据向量并构造后验。我正在使用 Wikipedia 为 NIW 指定的更新方程：http ://en.wikipedia.org/wiki/Conjugate_prior

我的分发类如下：

import numpy as np
from scipy.stats import chi2

class NormalInverseWishartDistribution(object):
    def __init__(self, mu, lmbda, nu, psi):
        self.mu = mu
        self.lmbda = float(lmbda)
        self.nu = nu
        self.psi = psi
        self.inv_psi = np.linalg.inv(psi)

    def sample(self):
        sigma = np.linalg.inv(self.wishartrand())
        return (np.random.multivariate_normal(self.mu, sigma / self.lmbda), sigma)

    def wishartrand(self):
        dim = self.inv_psi.shape[0]
        chol = np.linalg.cholesky(self.inv_psi)
        foo = np.zeros((dim,dim))

        for i in range(dim):
            for j in range(i+1):
                if i == j:
                    foo[i,j] = np.sqrt(chi2.rvs(self.nu-(i+1)+1))
                else:
                    foo[i,j]  = np.random.normal(0,1)
        return np.dot(chol, np.dot(foo, np.dot(foo.T, chol.T)))

    def posterior(self, data):
        n = len(data)
        mean_data = np.mean(data, axis=0)
        sum_squares = np.sum([np.array(np.matrix(x - mean_data).T * np.matrix(x - mean_data)) for x in data], axis=0)
        mu_n = (self.lmbda * self.mu + n * mean_data) / (self.lmbda + n)
        lmbda_n = self.lmbda + n
        nu_n = self.nu + n
        psi_n = self.psi + sum_squares + self.lmbda * n / float(self.lmbda + n) * np.array(np.matrix(mean_data - self.mu).T * np.matrix(mean_data - self.mu))
        return NormalInverseWishartDistribution(mu_n, lmbda_n, nu_n, psi_n)

我正在运行一个简单的健全性检查，看看后验是否收敛到真实分布：

x = NormalInverseWishartDistribution(np.array([0,0])-3,1,3,np.eye(2))
samples = [x.sample() for _ in range(100)]
data = [np.random.multivariate_normal(mu,cov) for mu,cov in samples]
y = NormalInverseWishartDistribution(np.array([0,0]),1,3,np.eye(2))
z = y.posterior(data)

print 'mu_n: {0}'.format(z.mu)

print 'psi_n: {0}'.format(z.psi)

均值适当收敛，但比例矩阵似乎收敛到沿对角线不正确的大值，而不是 1 的真实值。

据我所知，我正在完全复制更新规则。我在这里实施了一些不恰当的事情吗？

编辑：看起来实际上后验正在收敛，但示例程序正在返回有偏见的样本。我的采样方法有问题吗？

Edit2：我已经确认在 R 中的 MCMCpack riwish 函数中发生了同样的现象：

> library(MCMCpack)
> samples <- replicate(100000, riwish(3, matrix(c(1,0,0,1),2,2)))
> mean(samples[1,1,])
[1] 4.889211

这让我相信我一定是误解了一些东西。从维基百科页面（http://en.wikipedia.org/wiki/Inverse-Wishart_distribution），我们有：

$\newcommand{\E}{\mathrm{E}}$ $\E[{\Sigma}] = \frac{\Psi}{\nu - p - 1}$

但是，在我的测试用例中，，所以。因此，如果我进行多次采样，我不应该平均 $\nu = p + 1$ $\E[{\Sigma}] = \Psi$ $\Sigma$ $\Psi$

编辑3：意识到我的问题是简单的代数疏忽：我需要设置。问题解决了 $\nu = p+2$