1. 均值和方差未知的多元正态分布的后验（Multivariate normal with unknown mean and variance）

\quad假设有N个观测值{xi|i=1,2,...,N}\{x_{i}|i=1,2,...,N\}，且服从均值为μ\mu方差为Σ\Sigma的多元正态分布，即:

xi∼N(μ,Σ)

\begin{equation} \begin{aligned} x_{i}\sim N(\mu,\Sigma) \end{aligned} \end{equation}
均值和方差都未知的情况下，多元正态分布的共轭先验是正态逆威沙特分布（Normal-Inverse-Wishart），即有：

(μ,Σ)Σμ|Σ∼NIW(μ0,κ0;ν0,Λ0)∼Inv−Wishart(ν0,Λ0)∼N(μ0,Σ/κ0)

\begin{equation} \begin{aligned} (\mu,\Sigma)&\sim NIW(\mu_{0},\kappa_{0};\nu_{0},\Lambda_{0})\\ \Sigma & \sim Inv-Wishart(\nu_{0},\Lambda_{0})\\ \mu|\Sigma & \sim N(\mu_{0},\Sigma /\kappa_{0}) \end{aligned} \end{equation}
其中逆威沙特分布的概率密度函数为如下形式：

p(Σ|Λ0,ν0)=|Λ0|ν0/2|Σ|−(ν0+k+1)/2exp(−tr(Λ0Σ−1)/2)2ν0k/2Γk(ν0/2)

\begin{equation} p(\Sigma|\Lambda_{0},\nu_{0})=\frac{|\Lambda_{0}|^{\nu_{0}/2}|\Sigma|^{-(\nu_{0}+k+1)/2}exp(-tr(\Lambda_{0} \Sigma^{-1})/2)}{2^{\nu_{0} k/2}\Gamma_{k}(\nu_{0}/2)} \end{equation}
Γk(⋅)\Gamma_{k}(\cdot)是多变量Gamma分布， ν0\nu_{0}和 Λ0\Lambda_{0}分别是逆威沙特分布的自由度和尺度矩阵， kk是数据的维度。
依据文献[1]，在观测到数据{xi|i=1,2,...,N}\{x_{i}|i=1,2,...,N\}后，均值 μ\mu和方差 Σ\Sigma的后验分布依然服从正态逆威沙特分布，如下所示：

(μ,Σ)∼NIW(μ′,κ′;ν′,Λ′)

\begin{equation} (\mu,\Sigma) \sim NIW(\mu',\kappa';\nu',\Lambda') \end{equation}
其中：

μ′κ′ν′Λ′=κ0κ0+nμ0+Nκ0+Nx¯=κ0+Nν0+N=Λ0+∑n=1N(xi−x¯)(xi−x¯)T+κ0Nκ0+N(x¯−μ0)(x¯−μ0)T

\begin{equation} \begin{aligned} \mu'&=\frac{\kappa_{0}}{\kappa_{0}+n}\mu_{0}+\frac{N}{\kappa_{0}+N}\bar{x}\\ \kappa' & = \kappa_{0}+N\\ \nu' & \nu_{0}+N\\ \Lambda'&=\Lambda_{0}+\sum_{n=1}^{N}(x_{i}-\bar{x})(x_{i}-\bar{x})^{T}+\frac{\kappa_{0}N}{\kappa_{0}+N}(\bar{x}-\mu_{0})(\bar{x}-\mu_{0})^{T} \end{aligned} \end{equation}
得到了后验分布的概率密度函数，我们就可以通过其采样多元正态分布的均值 μ\mu和方差 Σ\Sigma了。

2. 从后验分布中采样均值μ\mu和方差Σ\Sigma

均值μ\mu的采样需要依赖于Σ\Sigma，因此采样顺序一般为：首先采样Σ∼Inv−Wishart(ν′,Λ′)\Sigma\sim Inv-Wishart(\nu',\Lambda')，然后采样μ|Σ,x∼N(μ′,Σ/κ′)\mu|\Sigma,x\sim N(\mu',\Sigma/\kappa')。关于均值的采样，可以看这篇博客。下面介绍一下如何从逆威沙特分布中采样方差Σ\Sigma。首先介绍一下Odell&Feiveson于1966年提出的基本采样思路[2]，然后给出Java代码。

一、假设Vi(1≤i≤k)V_{i}(1\leq i\leq k)是独立的随机变量，并且采样自自由度为ν−i+1\nu - i + 1的卡方分布，所有有ν−k+1≤ν−i+1≤ν\nu-k+1\leq \nu-i + 1\leq \nu.
二、假设NijN_{ij}是独立的采样自均值为0方差为1的正态分布中的随机变量，且有1≤i≤j≤k1\leq i \leq j\leq k，NijN_{ij}独立于ViV_{i}.
三、定义随机变量bijb_{ij}，且1≤i,j≤k1\leq i,j\leq k，当1≤i≤j≤k1\leq i \leq j \leq k时，有bji=bijb_{ji}=b_{ij}，我们通过如下公式构造bijb_{ij}。

biibij=Vi+∑r=1i−1N2ri,1≤i≤k=NijVi−−√+∑r=1i−1NriNrj,i<j≤k

\begin{equation} \begin{aligned} b_{ii}& = V_{i}+\sum_{r=1}^{i-1}N_{ri}^{2}, 1\leq i\leq k\\ b_{ij}& = N_{ij}\sqrt{V_{i}}+\sum_{r=1}^{i-1}N_{ri}N_{rj}, i
四、对方阵Λ\Lambda进行Cholesky分解，即LLT=Λ−1LL^{T}=\Lambda^{-1}
五、构造矩阵R=LBLTR=LBL^{T}
六、则Σ′=R−1\Sigma'=R^{-1}为该逆威沙特分布的样本。
至于为什么这么做，大家可参考文献[3]或者[2]。上面的过程已经很清晰了，下面我们给出具体的Java代码，来源自GitHub，并且做了一点的修改（Note，下面的代码使用的依然是commons.math的3.0版本，事实上，现在已经更新到4.0版本的。）

import java.util.Arrays;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.apache.commons.math3.distribution.GammaDistribution;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.CholeskyDecomposition;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;/*** Inverse Wishart distribution implementation, to sample random covariances matrices for* multivariate gaussian distributions.* <p/>* The sampling method follows the procedure described by Odell & Feiveson, 1966 to get samples* from a Wishart distribution, and then computes the inverse of the obtained samples.** @author Marc Pujol <mpujol@iiia.csic.es>*/
public class InverseWishartDistribution {private static final Logger LOG = Logger.getLogger(InverseWishartDistribution.class.getName());private GammaDistribution[] gammas;private double df;private RealMatrix scaleMatrix;private CholeskyDecomposition cholesky;private RandomGenerator random;/*** Builds a new Inverse Wishart distribution with the given scale and degrees of freedom.** @param scaleMatrix scale matrix(Λ)* @param df degrees of freedom.*/public InverseWishartDistribution(RealMatrix scaleMatrix, double df) {if (!scaleMatrix.isSquare()) {throw new RuntimeException("scaleMatrix must be square.");}this.scaleMatrix = scaleMatrix;this.df = df;this.random = new Well19937c();initialize();}private void initialize() {final int dim = scaleMatrix.getColumnDimension();// Build gamma distributions for the diagonalgammas = new GammaDistribution[dim];for (int i = 0; i < dim; i++) {gammas[i] = new GammaDistribution((df-i+0.0)/2, 2);}// Build the cholesky decompositioncholesky = new CholeskyDecomposition(inverseMatrix(scaleMatrix));}/*** Reseeds the random generator.** @param seed new random seed.*/public void reseedRandomGenerator(long seed) {random.setSeed(seed);for (int i = 0, len = scaleMatrix.getColumnDimension(); i < len; i++) {gammas[i].reseedRandomGenerator(seed+i);}}/*** Returns the inverse matrix.* @return inverted matrix.*/public RealMatrix inverseMatrix(RealMatrix A) {RealMatrix result = new LUDecomposition(A).getSolver().getInverse();return result; }/*** Returns a sample matrix from this distribution.* @return sampled matrix.*/public RealMatrix sample() {for (int i=0; i<100; i++) {try {RealMatrix A = sampleWishart();RealMatrix result = inverseMatrix(A);LOG.log(Level.FINE, "Cov = {0}", result);return result;} catch (SingularMatrixException ex) {LOG.finer("Discarding singular matrix generated by the wishart distribution.");}}throw new RuntimeException("Unable to generate inverse wishart samples!");}private RealMatrix sampleWishart() {final int dim = scaleMatrix.getColumnDimension();// Build N_{ij}double[][] N = new double[dim][dim];for (int j = 0; j < dim; j++) {for (int i = 0; i < j; i++) {N[i][j] = random.nextGaussian();}}if (LOG.isLoggable(Level.FINEST)) {LOG.log(Level.FINEST, "N = {0}", Arrays.deepToString(N));}// Build V_jdouble[] V = new double[dim];for (int i = 0; i < dim; i++) {V[i] = gammas[i].sample();}if (LOG.isLoggable(Level.FINEST)) {LOG.log(Level.FINEST, "V = {0}", Arrays.toString(V));}// Build Bdouble[][] B = new double[dim][dim];// b_{11} = V_1 (first j, where sum = 0 because i == j and the inner//               loop is never entered).// b_{jj} = V_j + \sum_{i=1}^{j-1} N_{ij}^2, j = 2, 3, ..., pfor (int j = 0; j < dim; j++) {double sum = 0;for (int i = 0; i < j; i++) {sum += Math.pow(N[i][j], 2);}B[j][j] = V[j] + sum;}if (LOG.isLoggable(Level.FINEST)) {LOG.log(Level.FINEST, "B*_jj : = {0}", Arrays.deepToString(B));}// b_{1j} = N_{1j} * \sqrt V_1for (int j = 1; j < dim; j++) {B[0][j] = N[0][j] * Math.sqrt(V[0]);B[j][0] = B[0][j];}if (LOG.isLoggable(Level.FINEST)) {LOG.log(Level.FINEST, "B*_1j = {0}", Arrays.deepToString(B));}// b_{ij} = N_{ij} * \sqrt V_1 + \sum_{k=1}^{i-1} N_{ki}*N_{kj}for (int j = 1; j < dim; j++) {for (int i = 1; i < j; i++) {double sum = 0;for (int k = 0; k < i; k++) {sum += N[k][i] * N[k][j];}B[i][j] = N[i][j] * Math.sqrt(V[i]) + sum;B[j][i] = B[i][j];}}if (LOG.isLoggable(Level.FINEST)) {LOG.log(Level.FINEST, "B* = {0}", Arrays.deepToString(B));}RealMatrix BMat = new Array2DRowRealMatrix(B);RealMatrix A = cholesky.getL().multiply(BMat).multiply(cholesky.getLT());if (LOG.isLoggable(Level.FINER)) {LOG.log(Level.FINER, "A* = {0}", Arrays.deepToString(N));}return A;}}

其中因为commons.math中的卡方分布没有采样函数，所以我们借助于commons.math提供的Gamma分布进行采样，事实上，卡方分布和Gamma概率密度函数非常相近。上述采样的核心其实是先从威沙特分布中采样一个方阵，然后求其逆矩阵，则得到逆威沙特分布的一个样本。代码中inverseMatrix(scaleMatrix)是将逆威沙特分布的尺度矩阵求逆，这样就得到威沙特分布的尺度矩阵。此外近一段时间找资料的过程还发现了其一些代码，如下：

Java代码：链接，其介绍文档。链接，其介绍文档。
c#代码：链接，其对应的介绍。
Matlab：其中有一个iwishrnd方法，其介绍在这里。

[1] Gelman, A., Carlin et al., Bayesian data analysis. London: Chapman & Hall
[2] Stanley Sawyer, Wishart Distributions and Inverse-Wishart Sampling
[3] Odell, P.L., and A.H. Feiveson (1966) A numerical procedure to generate a sample covariance matrix

多元正态分布的后验采样(包含程序)相关推荐

多元正态分布的后验采样
均值和方差未知的多元正态分布的后验Multivariate normal with unknown mean and variance 从后验分布中采样均值mu和方差Sigma 1. 均值和方差未知的 ...
numpy笔记整理 multivariate_normal（多元正态分布采样）
1 基本用法 np.random.multivariate_normal(mean, cov, size=None, check_valid=None, tol=None) 根据均值和协方差矩阵的情况 ...
【ML学习笔记】17：多元正态分布下极大似然估计最小错误率贝叶斯决策
简述多元正态分布下的最小错误率贝叶斯如果特征的值向量服从d元正态分布,即其概率密度函数为: 即其分布可以由均值向量和对称的协方差矩阵唯一确定. 如果认为样本的特征向量在类内服从多元正态分布: 即对 ...
R语言学习——一元与多元正态分布检验(也可以用于其他分布的检验)
文章目录 1 一元正态的评估 1.1 图像法 1.1.1 直方图 1.1.2 Q-Q图 1.2 峰度和偏度 1.3 统计检验 1.3.1 Shapiro-Wilks检验 1.3.2 Kolmogoro ...
【应用多元统计分析】CH3 多元正态分布
目录一.多元正态分布的定义 1.定义 2.二元正态分布二.多元正态分布的性质 [property1*] [property2] [property3] [property4] [property5 ...
【C++】13 多元正态分布抽样
目的 2022/4/28 用C++实现多元正态分布抽样本地目录:E:\Master\study\Cpp\MultivariateNormalDistributionSamples 参考资料 [1]h ...
透彻理解多元正态分布
文章目录前言多元标准正态分布多元一般正态分布缩放与位移相同尺度缩放与位移不同尺度旋转变换多元正态分布的概率密度函数多元正态分布的性质本篇内容主要是对于基本书籍教材多元正态分布相关章节 ...
R语言：多元正态分布的检验
多元正态分布的检验多元正态分布 mshapiro.test {mvnormtest} mvn {MVN} 多元正态分布转自个人微信公众号[Memo_Cleon]的统计学习笔记:多元正态分布检验的R ...
java后验条件_JAVA并发实战学习笔记——3，4章~
JAVA并发实战学习笔记第三章对象的共享失效数据: java程序实际运行中会出现①程序执行顺序对打乱:②数据对其它线程不可见--两种情况上述两种情况导致在缺乏同步的程序中出现失效数据这一现象, ...

多元正态分布的后验采样(包含程序)

1. 均值和方差未知的多元正态分布的后验（Multivariate normal with unknown mean and variance）

2. 从后验分布中采样均值μ\mu和方差Σ\Sigma

多元正态分布的后验采样(包含程序)相关推荐

最新文章

热门文章