本笔记基于Bayesian Data Analysis Third Edition， Andrew Gelman et. al. 学习编写，由于是英文教材，可能在学习过程中一些翻译或者内容有误，如有问题或错误，可以发送至邮箱[email protected]
假设我们现在有一个向量 θ = ( θ 1 , θ 2 ) \theta = (\theta_1, \theta_2) θ = ( θ 1 , θ 2 ) θ 1 \theta_1 θ 1 θ 2 \theta_2 θ 2
y ∣ μ , σ 2 ∼ N ( μ , σ 2 ) , μ = θ 1 , σ 2 = θ 2 y|\mu,\sigma^2 \sim \operatorname{N}(\mu,\sigma^2),\quad \mu = \theta_1, \quad \sigma^2 = \theta_2 y ∣ μ , σ 2 ∼ N ( μ , σ 2 ) , μ = θ 1 , σ 2 = θ 2 其中，均值和方差均未知。我们现在只对均值感兴趣，所以我们像要求的是均值的后验分布，即 p ( θ 1 ∣ y ) p(\theta_1|y) p ( θ 1 ∣ y )
p ( θ 1 , θ 2 ∣ y ) ∝ p ( y ∣ θ 1 , θ 2 ) p ( θ 1 , θ 2 ) p(\theta_1,\theta_2|y) \propto p(y|\theta_1,\theta_2)p(\theta_1,\theta_2) p ( θ 1 , θ 2 ∣ y ) ∝ p ( y ∣ θ 1 , θ 2 ) p ( θ 1 , θ 2 ) 如果我们只对 θ 1 \theta_1 θ 1 θ 2 \theta_2 θ 2
p ( θ 1 ∣ y ) = ∫ p ( θ 1 , θ 2 ∣ y ) d θ 2 p(\theta_1|y) = \int p(\theta_1,\theta_2|y)d\theta_2 p ( θ 1 ∣ y ) = ∫ p ( θ 1 , θ 2 ∣ y ) d θ 2 又或者，我们可以因子化后：
p ( θ 1 ∣ y ) = ∫ p ( θ 1 ∣ θ 2 , y ) p ( θ 2 ∣ y ) d θ 2 p(\theta_1|y) = \int p(\theta_1|\theta_2,y)p(\theta_2|y)d\theta_2 p ( θ 1 ∣ y ) = ∫ p ( θ 1 ∣ θ 2 , y ) p ( θ 2 ∣ y ) d θ 2 在上面的等式里，p ( θ 2 ∣ y ) p(\theta_2|y) p ( θ 2 ∣ y ) θ 2 \theta_2 θ 2 θ 1 \theta_1 θ 1 边缘和条件模拟 去进行计算，首先我们从 θ 2 \theta_2 θ 2 θ 2 \theta_2 θ 2 θ 2 \theta_2 θ 2 θ 1 \theta_1 θ 1
首先我们来考虑多参数模型中的第一个，正态分布。正态分布里有两个参数，一个是均值，一个是标准差。我们现在都不知道，所以我们首先使用非信息先验。我们在单参数模型（三） 提到过，均值是位置参数，标准差是缩放参数，因此构建非信息先验 时，我们会使用 ( μ , log σ ) (\mu, \log\sigma) ( μ , log σ ) 均匀分布 ，但我们需要转化为 p ( μ , σ 2 ) p(\mu,\sigma^2) p ( μ , σ 2 )
p ( μ , σ 2 ) = p ( μ , log σ ) ∣ J ∣ = p ( μ , log σ ) ∣ ∂ μ ∂ μ ∂ μ ∂ σ 2 ∂ log σ ∂ μ ∂ log σ ∂ σ 2 ∣ = p ( μ , log σ ) ∣ 1 0 0 1 2 σ 2 ∣ = p ( μ , log σ ) 1 2 σ 2 ∝ σ − 2 \begin{aligned}
p(\mu,\sigma^2) &= p(\mu,\log\sigma)|J| \\
&= p(\mu, \log\sigma) \left|
\begin{matrix}
\frac{\partial \mu}{\partial \mu} & \frac{\partial \mu}{\partial \sigma^2} \\
\frac{\partial \log\sigma}{\partial \mu} & \frac{\partial \log\sigma}{\partial \sigma^2}
\end{matrix}
\right| \\
&= p(\mu, \log\sigma) \left|
\begin{matrix}
1 & 0 \\
0 & \frac{1}{2\sigma^2}
\end{matrix}
\right| \\
&= p(\mu, \log\sigma)\frac{1}{2\sigma^2} \\
&\propto \sigma^{-2}
\end{aligned} p ( μ , σ 2 ) = p ( μ , log σ ) ∣ J ∣ = p ( μ , log σ ) ∂ μ ∂ μ ∂ μ ∂ l o g σ ∂ σ 2 ∂ μ ∂ σ 2 ∂ l o g σ = p ( μ , log σ ) 1 0 0 2 σ 2 1 = p ( μ , log σ ) 2 σ 2 1 ∝ σ − 2 由于是均匀分布，且 μ \mu μ log σ \log\sigma log σ p ( μ , log σ ) = p ( μ ) p ( l o g σ ) = C o n s t a n t p(\mu, \log\sigma) = p(\mu)p(log\sigma) = Constant p ( μ , log σ ) = p ( μ ) p ( l o g σ ) = C o n s t an t
似然函数为：
p ( y ∣ μ , σ 2 ) ∝ ∏ i = 1 n σ − 1 exp ( − 1 2 σ 2 ( y i − μ ) 2 ) = σ − n exp ( − 1 2 σ 2 ∑ i = 1 n ( y i − μ ) 2 ) \begin{aligned}
p(y|\mu,\sigma^2) &\propto \prod^n_{i=1} \sigma^{-1} \exp\left(-\frac{1}{2\sigma^2} (y_i - \mu)^2\right) \\
&= \sigma^{-n} \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - \mu)^2\right)
\end{aligned} p ( y ∣ μ , σ 2 ) ∝ i = 1 ∏ n σ − 1 exp ( − 2 σ 2 1 ( y i − μ ) 2 ) = σ − n exp ( − 2 σ 2 1 i = 1 ∑ n ( y i − μ ) 2 ) 有了这个以后，我们就可以写出联合后验分布 的计算
p ( μ , σ 2 ∣ y ) ∝ σ − n − 2 exp ( − 1 2 σ 2 ∑ i = 1 n ( y i − μ ) 2 ) = σ − n − 2 exp ( − 1 2 σ 2 [ ∑ i = 1 n ( y i − y ˉ ) 2 + n ( y ˉ − μ ) 2 ] ) = σ − n − 2 exp ( − 1 2 σ 2 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) \begin{aligned}
p(\mu, \sigma^2|y) &\propto \sigma^{-n-2} \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - \mu)^2\right) \\
&= \sigma^{-n-2} \exp\left(-\frac{1}{2\sigma^2}\left[\sum_{i=1}^n (y_i - \bar{y})^2 + n(\bar{y} - \mu)^2\right]\right) \\
&= \sigma^{-n-2} \exp\left(-\frac{1}{2\sigma^2}\big[(n-1)s^2 + n(\bar{y} - \mu)^2\big]\right)
\end{aligned} p ( μ , σ 2 ∣ y ) ∝ σ − n − 2 exp ( − 2 σ 2 1 i = 1 ∑ n ( y i − μ ) 2 ) = σ − n − 2 exp ( − 2 σ 2 1 [ i = 1 ∑ n ( y i − y ˉ ) 2 + n ( y ˉ − μ ) 2 ] ) = σ − n − 2 exp ( − 2 σ 2 1 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) 其中，
s 2 = 1 n − 1 ∑ i = 1 n ( y i − y ˉ ) 2 s^2 = \frac{1}{n-1} \sum_{i=1}^{n} (y_i - \bar{y})^2 s 2 = n − 1 1 i = 1 ∑ n ( y i − y ˉ ) 2 其是y i y_i y i y ˉ \bar y y ˉ s 2 s^2 s 2
然后，根据我们之前说的，我们写出条件后验分布 ，p ( μ ∣ σ 2 , y ) p(\mu|\sigma^2,y) p ( μ ∣ σ 2 , y )
我们在这里将 σ 2 \sigma^2 σ 2
μ ∣ σ 2 , y ∼ N ( y ˉ , σ 2 / n ) \mu|\sigma^2,y \sim \operatorname{N}(\bar y, \sigma^2/n) μ ∣ σ 2 , y ∼ N ( y ˉ , σ 2 / n ) 然后我们写出边缘后验分布 ，p ( σ 2 ∣ y ) p(\sigma^2|y) p ( σ 2 ∣ y )
所以：
p ( σ 2 ∣ y ) ∝ ∫ σ − n − 2 exp ( − 1 2 σ 2 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) d μ = σ − n − 2 exp ( − 1 2 σ 2 [ ( n − 1 ) s 2 ] ) ∫ exp ( − n 2 σ 2 ( y ˉ − μ ) 2 ) d μ = σ − n − 2 exp ( − 1 2 σ 2 ( n − 1 ) s 2 ) 2 π σ 2 / n ∝ σ − n − 2 exp ( − 1 2 σ 2 ( n − 1 ) s 2 ) = ( σ 2 ) − ( n + 1 ) / 2 exp ( − 1 2 σ 2 ( n − 1 ) s 2 ) \begin{aligned}
p(\sigma^2|y) &\propto \int \sigma^{-n-2} \exp\left(-\frac{1}{2\sigma^2}\big[(n-1)s^2 + n(\bar{y} - \mu)^2\big]\right) d\mu \\
&= \sigma^{-n-2}\exp\left(-\frac{1}{2\sigma^2}\big[(n-1)s^2\big]\right)\int\exp\left(-\frac{n}{2\sigma^2}(\bar{y} - \mu)^2\right) d\mu \\
&= \sigma^{-n-2}\exp\left(-\frac{1}{2\sigma^2}(n-1)s^2\right) \sqrt{2\pi\sigma^2/n} \\
&\propto \sigma^{-n-2}\exp\left(-\frac{1}{2\sigma^2}(n-1)s^2\right) \\
&=(\sigma^2)^{-(n+1)/2}\exp\left(-\frac{1}{2\sigma^2}(n-1)s^2\right)
\end{aligned} p ( σ 2 ∣ y ) ∝ ∫ σ − n − 2 exp ( − 2 σ 2 1 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) d μ = σ − n − 2 exp ( − 2 σ 2 1 [ ( n − 1 ) s 2 ] ) ∫ exp ( − 2 σ 2 n ( y ˉ − μ ) 2 ) d μ = σ − n − 2 exp ( − 2 σ 2 1 ( n − 1 ) s 2 ) 2 π σ 2 / n ∝ σ − n − 2 exp ( − 2 σ 2 1 ( n − 1 ) s 2 ) = ( σ 2 ) − ( n + 1 ) /2 exp ( − 2 σ 2 1 ( n − 1 ) s 2 ) 其中，第三个等号是由高斯积分得到，即：
∫ − ∞ ∞ e − a ( x + b ) 2 d x = π a \int_{-\infty}^{\infty} e^{-a(x+b)^2} dx = \sqrt{\frac{\pi}{a}} ∫ − ∞ ∞ e − a ( x + b ) 2 d x = a π 此时，符合 Scaled inverse-χ 2 \chi^2 χ 2
σ 2 ∣ y ∼ Inv- χ 2 ( n − 1 , s 2 ) \sigma^2|y \sim \operatorname{Inv-}\chi^2 (n-1,s^2) σ 2 ∣ y ∼ Inv- χ 2 ( n − 1 , s 2 ) 因此，我们就可以先从 σ 2 ∣ y \sigma^2|y σ 2 ∣ y σ \sigma σ μ ∣ σ 2 , y \mu|\sigma^2,y μ ∣ σ 2 , y μ \mu μ
如果，我们直接计算边缘分布
p ( μ ∣ y ) = ∫ 0 + ∞ p ( μ ∣ σ 2 , y ) d σ 2 = ∫ 0 + ∞ σ − n − 2 exp ( − 1 2 σ 2 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) d σ 2 \begin{aligned}
p(\mu|y) &= \int_0^{+\infin} p(\mu|\sigma^2,y) d\sigma^2 \\
&= \int_0^{+\infin}\sigma^{-n-2} \exp\left(-\frac{1}{2\sigma^2}\big[(n-1)s^2 + n(\bar{y} - \mu)^2\big]\right)d\sigma^2
\end{aligned} p ( μ ∣ y ) = ∫ 0 + ∞ p ( μ ∣ σ 2 , y ) d σ 2 = ∫ 0 + ∞ σ − n − 2 exp ( − 2 σ 2 1 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) d σ 2 我们把与 σ 2 \sigma^2 σ 2 A = ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 A = (n-1)s^2 + n(\bar{y} - \mu)^2 A = ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 z = A 2 σ 2 z = \frac{A}{2\sigma^2}\\ z = 2 σ 2 A
于是乎：
d σ 2 = d z ⋅ − A 2 z 2 d\sigma^2 = dz \cdot -\frac{A}{2z^2} d σ 2 = d z ⋅ − 2 z 2 A 换元代入：
p ( μ ∣ y ) = ∫ 0 + ∞ ( σ 2 ) − n 2 − 1 exp ( − A 2 σ 2 ) d σ 2 = ∫ 0 + ∞ ( A 2 z ) − n 2 − 1 exp ( − z ) ( A 2 z 2 ) d z ∝ A − n / 2 ∫ 0 ∞ z ( n − 2 ) / 2 exp ( − z ) d z ∝ [ ( n − 1 ) s 2 + n ( μ − y ˉ ) 2 ] − n / 2 ∝ [ 1 + n ( μ − y ˉ ) 2 ( n − 1 ) s 2 ] − n / 2 \begin{aligned}
p(\mu|y) &= \int_0^{+\infin} (\sigma^2)^{-\frac{n}{2}-1}\exp(-\frac{A}{2\sigma^2})d\sigma^2\\
&= \int_0^{+\infin} (\frac{A}{2z})^{-\frac{n}{2}-1}\exp(-z)(\frac{A}{2z^2})dz \\
&\propto A^{-n/2} \int_{0}^{\infty} z^{(n-2)/2} \exp(-z) dz \\
&\propto \Big[(n-1)s^{2} + n(\mu - \bar{y})^{2}\Big]^{-n/2} \\
&\propto \Big[1 + \frac{n(\mu - \bar{y})^{2}}{(n-1)s^{2}}\Big]^{-n/2}
\end{aligned} p ( μ ∣ y ) = ∫ 0 + ∞ ( σ 2 ) − 2 n − 1 exp ( − 2 σ 2 A ) d σ 2 = ∫ 0 + ∞ ( 2 z A ) − 2 n − 1 exp ( − z ) ( 2 z 2 A ) d z ∝ A − n /2 ∫ 0 ∞ z ( n − 2 ) /2 exp ( − z ) d z ∝ [ ( n − 1 ) s 2 + n ( μ − y ˉ ) 2 ] − n /2 ∝ [ 1 + ( n − 1 ) s 2 n ( μ − y ˉ ) 2 ] − n /2 其中，第二个等号是由于变换变量后积分上下限变换的负号与被积项变量转换负号相抵消，第四个等比例符号是由于积分内的项目与 μ \mu μ
因此我们意识到，这个分布服从 t t t
μ ∣ y ∼ t n − 1 ( y ˉ , s 2 / n ) \mu|y \sim t_{n-1}(\bar y,s^2/n) μ ∣ y ∼ t n − 1 ( y ˉ , s 2 / n )
t t t
记作 θ ∼ t ν ( μ , σ 2 ) \theta \sim t_\nu(\mu,\sigma^2) θ ∼ t ν ( μ , σ 2 ) p ( θ ) = t ν ( θ ∣ μ , σ 2 ) p(\theta) = t_{\nu}(\theta|\mu,\sigma^2) p ( θ ) = t ν ( θ ∣ μ , σ 2 )
其中 ν \nu ν μ \mu μ σ 2 \sigma^2 σ 2
密度函数：
p ( θ ) = Γ ( ( ν + 1 ) / 2 ) Γ ( ν / 2 ) ν π σ ( 1 + 1 ν ( θ − μ σ ) 2 ) − ( ν + 1 ) / 2 p(\theta) = \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2) \sqrt{\nu \pi}\sigma} \left(1 + \frac{1}{\nu} \Big(\frac{\theta - \mu}{\sigma}\Big)^{2}\right)^{-(\nu+1)/2} p ( θ ) = Γ ( ν /2 ) ν π σ Γ (( ν + 1 ) /2 ) ( 1 + ν 1 ( σ θ − μ ) 2 ) − ( ν + 1 ) /2 期望：E ( θ ) = μ , for ν > 1 E(\theta) = \mu,\quad \text{for } \nu > 1 E ( θ ) = μ , for ν > 1
方差：var ( θ ) = ν ν − 2 σ 2 , for ν > 2 \operatorname{var}(\theta) = \frac{\nu}{\nu-2} \sigma^2,\quad \text{for } \nu > 2 var ( θ ) = ν − 2 ν σ 2 , for ν > 2
众数：mode ( θ ) = μ \operatorname{mode}(\theta) = \mu mode ( θ ) = μ
我们转化为标准 t t t
μ − y ˉ s / n ∣ y ∼ t n − 1 \frac{\mu - \bar{y}}{s/\sqrt{n}}\bigg| y \sim t_{n-1} s / n μ − y ˉ y ∼ t n − 1 其中 t n − 1 t_{n-1} t n − 1 t t t
这种边缘后验分布为与抽样理论的比较提供了另一个有趣的视角，在抽样分布 p ( y ∣ μ , σ 2 ) p(y|\mu,\sigma^2) p ( y ∣ μ , σ 2 )
y ˉ − μ s / n ∣ y ∼ t n − 1 \frac{\bar{y}-\mu}{s/\sqrt{n}}\bigg| y \sim t_{n-1} s / n y ˉ − μ y ∼ t n − 1 枢轴量的抽样分布 ( y ˉ − μ ) / ( s / n ) (\bar{y}-\mu)/(s/\sqrt{n}) ( y ˉ − μ ) / ( s / n ) σ 2 \sigma^2 σ 2
后验预测分布可以写作：
p ( y ~ ∣ y ) = ∬ p ( y ~ ∣ μ , σ 2 , y ) p ( μ , σ 2 ∣ y ) d μ d σ 2 p(\tilde{y}|y) = \iint p(\tilde{y}|\mu,\sigma^2, y)p(\mu,\sigma^2|y)d\mu d\sigma^2 p ( y ~ ∣ y ) = ∬ p ( y ~ ∣ μ , σ 2 , y ) p ( μ , σ 2 ∣ y ) d μ d σ 2 其中，第一项是对于未来观测值的给定 μ , σ 2 \mu,\sigma^2 μ , σ 2 μ , σ 2 \mu,\sigma^2 μ , σ 2 y ~ \tilde y y ~ y y y
p ( y ~ ∣ μ , σ 2 , y ) = p ( y ~ ∣ μ , σ 2 ) = N ( y ~ ∣ μ , σ 2 ) p(\tilde y|\mu,\sigma^2,y)=p(\tilde y|\mu,\sigma^2)=\operatorname{N}(\tilde y|\mu,\sigma^2) p ( y ~ ∣ μ , σ 2 , y ) = p ( y ~ ∣ μ , σ 2 ) = N ( y ~ ∣ μ , σ 2 ) 而第二项为参数的联合后验分布。为了从后验预测分布里抽样，我们第一步先从联合后验分布 p ( μ , σ 2 ∣ y ) p(\mu,\sigma^2|y) p ( μ , σ 2 ∣ y ) ( μ , σ 2 ) (\mu,\sigma^2) ( μ , σ 2 )
y ~ ∣ μ , σ 2 ∼ N ( μ , σ 2 ) \tilde y|\mu,\sigma^2 \sim \operatorname{N}(\mu,\sigma^2) y ~ ∣ μ , σ 2 ∼ N ( μ , σ 2 ) 实际上，后验预测分布也同样满足 t t t
y ~ ∣ y ∼ t n − 1 ( y ˉ , ( 1 + 1 n ) 1 / 2 s ) \tilde y |y\sim t_{n-1}\left(\bar y, (1+\frac{1}{n})^{1/2}s\right) y ~ ∣ y ∼ t n − 1 ( y ˉ , ( 1 + n 1 ) 1/2 s ) 其推导方式同之前一样（见后），都是把干扰参数积分掉。不过这里更容易理解的方法是：先固定 σ 2 \sigma^2 σ 2 μ \mu μ σ 2 \sigma^2 σ 2
首先将联合后验分布因子化：
p ( μ , σ 2 ∣ y ) = p ( μ ∣ σ 2 , y ) p ( σ 2 ∣ y ) p(\mu,\sigma^2|y)=p(\mu|\sigma^2,y)p(\sigma^2|y) p ( μ , σ 2 ∣ y ) = p ( μ ∣ σ 2 , y ) p ( σ 2 ∣ y ) 于是后验预测分布可以写成：
p ( y ~ ∣ y ) = ∬ p ( y ~ ∣ μ , σ 2 , y ) p ( μ ∣ σ 2 , y ) p ( σ 2 ∣ y ) d μ d σ 2 = ∫ [ ∫ p ( y ~ ∣ μ , σ 2 , y ) p ( μ ∣ σ 2 , y ) d μ ] p ( σ 2 ∣ y ) d σ 2 \begin{aligned}
p(\tilde y|y)
&= \iint p(\tilde y|\mu,\sigma^2,y)p(\mu|\sigma^2,y)p(\sigma^2|y)d\mu d\sigma^2 \\
&= \int \left[\int p(\tilde y|\mu,\sigma^2,y)p(\mu|\sigma^2,y)d\mu\right]p(\sigma^2|y)d\sigma^2
\end{aligned} p ( y ~ ∣ y ) = ∬ p ( y ~ ∣ μ , σ 2 , y ) p ( μ ∣ σ 2 , y ) p ( σ 2 ∣ y ) d μ d σ 2 = ∫ [ ∫ p ( y ~ ∣ μ , σ 2 , y ) p ( μ ∣ σ 2 , y ) d μ ] p ( σ 2 ∣ y ) d σ 2 中括号里的部分为
p ( y ~ ∣ σ 2 , y ) = ∫ p ( y ~ ∣ μ , σ 2 , y ) p ( μ ∣ σ 2 , y ) d μ p(\tilde y|\sigma^2,y)
=\int p(\tilde y|\mu,\sigma^2,y)p(\mu|\sigma^2,y)d\mu p ( y ~ ∣ σ 2 , y ) = ∫ p ( y ~ ∣ μ , σ 2 , y ) p ( μ ∣ σ 2 , y ) d μ 它的含义是：如果先假设 σ 2 \sigma^2 σ 2 y ~ \tilde y y ~
均值 μ \mu μ
即使 μ \mu μ y ~ \tilde y y ~ μ \mu μ
我们已经知道：
μ ∣ σ 2 , y ∼ N ( y ˉ , σ 2 n ) \mu|\sigma^2,y\sim \operatorname{N}\left(\bar y,\frac{\sigma^2}{n}\right) μ ∣ σ 2 , y ∼ N ( y ˉ , n σ 2 ) 同时：
y ~ ∣ μ , σ 2 , y ∼ N ( μ , σ 2 ) \tilde y|\mu,\sigma^2,y\sim \operatorname{N}(\mu,\sigma^2) y ~ ∣ μ , σ 2 , y ∼ N ( μ , σ 2 ) 因此可以把未来观测值写成：
y ~ = μ + ϵ , ϵ ∣ σ 2 ∼ N ( 0 , σ 2 ) \tilde y=\mu+\epsilon,\qquad \epsilon|\sigma^2\sim \operatorname{N}(0,\sigma^2) y ~ = μ + ϵ , ϵ ∣ σ 2 ∼ N ( 0 , σ 2 ) 在给定 σ 2 \sigma^2 σ 2 y y y
μ = y ˉ + η , η ∣ σ 2 , y ∼ N ( 0 , σ 2 n ) \mu=\bar y+\eta,\qquad \eta|\sigma^2,y\sim \operatorname{N}\left(0,\frac{\sigma^2}{n}\right) μ = y ˉ + η , η ∣ σ 2 , y ∼ N ( 0 , n σ 2 ) 所以：
y ~ = y ˉ + η + ϵ \tilde y=\bar y+\eta+\epsilon y ~ = y ˉ + η + ϵ 其中 η \eta η μ \mu μ ϵ \epsilon ϵ
y ~ ∣ σ 2 , y ∼ N ( y ˉ , σ 2 n + σ 2 ) = N ( y ˉ , ( 1 + 1 n ) σ 2 ) \tilde y|\sigma^2,y
\sim \operatorname{N}\left(\bar y,\frac{\sigma^2}{n}+\sigma^2\right)
=\operatorname{N}\left(\bar y,\left(1+\frac{1}{n}\right)\sigma^2\right) y ~ ∣ σ 2 , y ∼ N ( y ˉ , n σ 2 + σ 2 ) = N ( y ˉ , ( 1 + n 1 ) σ 2 ) 其和 μ ∣ σ 2 , y \mu|\sigma^2,y μ ∣ σ 2 , y
μ ∣ σ 2 , y ∼ N ( y ˉ , σ 2 n ) \mu|\sigma^2,y\sim \operatorname{N}\left(\bar y,\frac{\sigma^2}{n}\right) μ ∣ σ 2 , y ∼ N ( y ˉ , n σ 2 ) 接下来再对 σ 2 \sigma^2 σ 2
σ 2 ∣ y ∼ Inv- χ 2 ( n − 1 , s 2 ) \sigma^2|y\sim \operatorname{Inv-}\chi^2(n-1,s^2) σ 2 ∣ y ∼ Inv- χ 2 ( n − 1 , s 2 ) 所以：
p ( y ~ ∣ y ) = ∫ p ( y ~ ∣ σ 2 , y ) p ( σ 2 ∣ y ) d σ 2 ∝ ∫ 0 ∞ ( σ 2 ) − 1 / 2 exp [ − ( y ~ − y ˉ ) 2 2 ( 1 + 1 n ) σ 2 ] ( σ 2 ) − ( n + 1 ) / 2 exp [ − ( n − 1 ) s 2 2 σ 2 ] d σ 2 = ∫ 0 ∞ ( σ 2 ) − n / 2 − 1 exp { − 1 2 σ 2 [ ( n − 1 ) s 2 + ( y ~ − y ˉ ) 2 1 + 1 n ] } d σ 2 \begin{aligned}
p(\tilde y|y)
&= \int p(\tilde y|\sigma^2,y)p(\sigma^2|y)d\sigma^2 \\
&\propto \int_0^{\infty}
(\sigma^2)^{-1/2}
\exp\left[-\frac{(\tilde y-\bar y)^2}{2(1+\frac{1}{n})\sigma^2}\right]
(\sigma^2)^{-(n+1)/2}
\exp\left[-\frac{(n-1)s^2}{2\sigma^2}\right]
d\sigma^2 \\
&= \int_0^{\infty}
(\sigma^2)^{-n/2-1}
\exp\left\{
-\frac{1}{2\sigma^2}
\left[(n-1)s^2+\frac{(\tilde y-\bar y)^2}{1+\frac{1}{n}}\right]
\right\}
d\sigma^2
\end{aligned} p ( y ~ ∣ y ) = ∫ p ( y ~ ∣ σ 2 , y ) p ( σ 2 ∣ y ) d σ 2 ∝ ∫ 0 ∞ ( σ 2 ) − 1/2 exp [ − 2 ( 1 + n 1 ) σ 2 ( y ~ − y ˉ ) 2 ] ( σ 2 ) − ( n + 1 ) /2 exp [ − 2 σ 2 ( n − 1 ) s 2 ] d σ 2 = ∫ 0 ∞ ( σ 2 ) − n /2 − 1 exp { − 2 σ 2 1 [ ( n − 1 ) s 2 + 1 + n 1 ( y ~ − y ˉ ) 2 ] } d σ 2 令
A = ( n − 1 ) s 2 + ( y ~ − y ˉ ) 2 1 + 1 n A=(n-1)s^2+\frac{(\tilde y-\bar y)^2}{1+\frac{1}{n}} A = ( n − 1 ) s 2 + 1 + n 1 ( y ~ − y ˉ ) 2 则上面的积分与之前推导 p ( μ ∣ y ) p(\mu|y) p ( μ ∣ y )
∫ 0 ∞ ( σ 2 ) − n / 2 − 1 exp ( − A 2 σ 2 ) d σ 2 ∝ A − n / 2 \int_0^{\infty}(\sigma^2)^{-n/2-1}\exp\left(-\frac{A}{2\sigma^2}\right)d\sigma^2
\propto A^{-n/2} ∫ 0 ∞ ( σ 2 ) − n /2 − 1 exp ( − 2 σ 2 A ) d σ 2 ∝ A − n /2 因此：
p ( y ~ ∣ y ) ∝ [ ( n − 1 ) s 2 + ( y ~ − y ˉ ) 2 1 + 1 n ] − n / 2 ∝ [ 1 + 1 n − 1 ( y ~ − y ˉ s 1 + 1 n ) 2 ] − n / 2 \begin{aligned}
p(\tilde y|y)
&\propto
\left[(n-1)s^2+\frac{(\tilde y-\bar y)^2}{1+\frac{1}{n}}\right]^{-n/2} \\
&\propto
\left[
1+\frac{1}{n-1}
\left(
\frac{\tilde y-\bar y}{s\sqrt{1+\frac{1}{n}}}
\right)^2
\right]^{-n/2}
\end{aligned} p ( y ~ ∣ y ) ∝ [ ( n − 1 ) s 2 + 1 + n 1 ( y ~ − y ˉ ) 2 ] − n /2 ∝ 1 + n − 1 1 s 1 + n 1 y ~ − y ˉ 2 − n /2 这正是自由度为 n − 1 n-1 n − 1 y ˉ \bar y y ˉ s 1 + 1 n s\sqrt{1+\frac{1}{n}} s 1 + n 1 t t t
y ~ ∣ y ∼ t n − 1 ( y ˉ , s 1 + 1 n ) \tilde y|y\sim t_{n-1}\left(\bar y,s\sqrt{1+\frac{1}{n}}\right) y ~ ∣ y ∼ t n − 1 ( y ˉ , s 1 + n 1 ) 因此，后验预测分布比 μ \mu μ
μ ∣ y ∼ t n − 1 ( y ˉ , s n ) , y ~ ∣ y ∼ t n − 1 ( y ˉ , s 1 + 1 n ) \mu|y\sim t_{n-1}\left(\bar y,\frac{s}{\sqrt n}\right),
\qquad
\tilde y|y\sim t_{n-1}\left(\bar y,s\sqrt{1+\frac{1}{n}}\right) μ ∣ y ∼ t n − 1 ( y ˉ , n s ) , y ~ ∣ y ∼ t n − 1 ( y ˉ , s 1 + n 1 ) 这是因为预测一个新的观测值时，不仅要考虑均值 μ \mu μ
在将模型推向更一般的情况时，我们为二参数单变量正态抽样模型假设一个共轭先验分布会比非信息先验分布更值得考虑。对于
p ( μ , σ 2 ∣ y ) ∝ σ − n − 2 exp ( − 1 2 σ 2 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) p(\mu, \sigma^2|y) \propto \sigma^{-n-2} \exp\left(-\frac{1}{2\sigma^2}\big[(n-1)s^2 + n(\bar{y} - \mu)^2\big]\right) p ( μ , σ 2 ∣ y ) ∝ σ − n − 2 exp ( − 2 σ 2 1 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) 这种形式的，共轭先验密度必须具有乘积形式 p ( σ 2 ) p ( μ ∣ σ 2 ) p(\sigma^2)p(\mu|\sigma^2) p ( σ 2 ) p ( μ ∣ σ 2 ) σ 2 \sigma^2 σ 2 χ 2 \chi^2 χ 2 μ ∣ σ 2 \mu|\sigma^2 μ ∣ σ 2 μ \mu μ t t t
μ ∣ σ 2 ∼ N ( μ 0 , σ 2 / κ 0 ) σ 2 ∼ I n v − χ 2 ( ν 0 , σ 0 2 ) \mu|\sigma^{2} \sim \mathrm{N}(\mu_{0}, \sigma^{2}/\kappa_{0}) \\ \sigma^{2} \sim \mathrm{Inv-}\chi^{2}(\nu_{0}, \sigma_{0}^{2}) μ ∣ σ 2 ∼ N ( μ 0 , σ 2 / κ 0 ) σ 2 ∼ Inv − χ 2 ( ν 0 , σ 0 2 ) 此时，对应的先验分布为：
p ( μ , σ 2 ) ∝ p ( μ ∣ σ 2 ) p ( σ 2 ) ∝ κ 0 2 π σ exp ( − κ 0 2 σ 2 ( μ − μ 0 ) 2 ) ( σ 2 ) − ( ν 0 / 2 + 1 ) exp ( − 1 2 σ 2 ν 0 σ 0 2 ) ∝ σ − 1 ( σ 2 ) − ( ν 0 / 2 + 1 ) exp ( − 1 2 σ 2 [ ν 0 σ 0 2 + κ 0 ( μ 0 − μ ) 2 ] ) \begin{aligned}
p(\mu,\sigma^2) &\propto p(\mu|\sigma^2)p(\sigma^2) \\
&\propto \frac{\sqrt{\kappa_0}}{\sqrt{2\pi}\sigma}\exp\left({-\frac{\kappa_0}{2\sigma^2}(\mu-\mu_0)^2}\right) (\sigma^2)^{-(\nu_0/2+1)}\exp\left({-\frac{1}{2\sigma^2}\nu_0\sigma_0^2}\right) \\
&\propto \sigma^{-1}(\sigma^{2})^{- (\nu_{0}/2+1)} \exp\left(-\frac{1}{2\sigma^{2}}\left[ \nu_{0}\sigma_{0}^{2} + \kappa_{0}(\mu_{0}-\mu)^{2} \right]\right)
\end{aligned} p ( μ , σ 2 ) ∝ p ( μ ∣ σ 2 ) p ( σ 2 ) ∝ 2 π σ κ 0 exp ( − 2 σ 2 κ 0 ( μ − μ 0 ) 2 ) ( σ 2 ) − ( ν 0 /2 + 1 ) exp ( − 2 σ 2 1 ν 0 σ 0 2 ) ∝ σ − 1 ( σ 2 ) − ( ν 0 /2 + 1 ) exp ( − 2 σ 2 1 [ ν 0 σ 0 2 + κ 0 ( μ 0 − μ ) 2 ] ) 对于这个共轭先验，我们记作
N − I n v − χ 2 ( μ , σ 2 ∣ μ 0 , σ 2 / κ 0 ; ν 0 , σ 0 2 ) \mathrm{N-Inv-}\chi^2(\mu,\sigma^2|\mu_0,\sigma^2/\kappa_0;\nu_0,\sigma_0^2) N − Inv − χ 2 ( μ , σ 2 ∣ μ 0 , σ 2 / κ 0 ; ν 0 , σ 0 2 ) 这个分布里有四个未知参数: ( μ 0 , κ 0 , ν 0 , σ 0 2 ) (\mu_0,\kappa_0,\nu_0,\sigma_0^2) ( μ 0 , κ 0 , ν 0 , σ 0 2 ) μ \mu μ σ 2 \sigma^2 σ 2
μ ∣ σ 2 \mu|\sigma^2 μ ∣ σ 2 σ 2 \sigma^2 σ 2 μ \mu μ σ 2 \sigma^2 σ 2 σ 2 \sigma^2 σ 2 μ \mu μ μ \mu μ y y y κ 0 \kappa_0 κ 0
我们有了先验之后，我们乘以似然得到后验分布：
p ( μ , σ 2 ∣ y ) ∝ p ( μ , σ 2 ) p ( y ∣ μ , σ 2 ) ∝ σ − 1 ( σ 2 ) − ( ν 0 / 2 + 1 ) exp ( − 1 2 σ 2 [ ν 0 σ 0 2 + κ 0 ( μ − μ 0 ) 2 ] ) × ( σ 2 ) − n / 2 exp ( − 1 2 σ 2 [ ( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) \begin{aligned}
p(\mu,\sigma^2|y) &\propto p(\mu,\sigma^2)p(y|\mu,\sigma^2)\\
&\propto \sigma^{-1}(\sigma^{2})^{-(\nu_{0}/2+1)} \exp\left(-\frac{1}{2\sigma^{2}}[\nu_{0}\sigma_{0}^{2}+\kappa_{0}(\mu-\mu_{0})^{2}]\right)\times (\sigma^{2})^{-n/2}\exp\left(-\frac{1}{2\sigma^{2}}[(n-1)s^{2}+n(\bar{y}-\mu)^{2}]\right)
\end{aligned} p ( μ , σ 2 ∣ y ) ∝ p ( μ , σ 2 ) p ( y ∣ μ , σ 2 ) ∝ σ − 1 ( σ 2 ) − ( ν 0 /2 + 1 ) exp ( − 2 σ 2 1 [ ν 0 σ 0 2 + κ 0 ( μ − μ 0 ) 2 ] ) × ( σ 2 ) − n /2 exp ( − 2 σ 2 1 [( n − 1 ) s 2 + n ( y ˉ − μ ) 2 ] ) 通过代数转变：
ν n = ν 0 + n κ n = κ 0 + n μ n = κ 0 κ 0 + n μ 0 + n κ 0 + n y ˉ ν n σ n 2 = ν 0 σ 0 2 + ( n − 1 ) s 2 + κ 0 n κ 0 + n ( y ˉ − μ 0 ) 2 \nu_n = \nu_0+n \\
\kappa_n = \kappa_0 +n \\
\mu_n = \frac{\kappa_0}{\kappa_0+n} \mu_0+ \frac{n}{\kappa_0+n} \bar{y} \\
\nu_n\sigma_n^2 = \nu_0\sigma_0^2 + (n-1)s^2 + \frac{\kappa_0n}{\kappa_0+n}(\bar{y}-\mu_0)^2 ν n = ν 0 + n κ n = κ 0 + n μ n = κ 0 + n κ 0 μ 0 + κ 0 + n n y ˉ ν n σ n 2 = ν 0 σ 0 2 + ( n − 1 ) s 2 + κ 0 + n κ 0 n ( y ˉ − μ 0 ) 2 于是乎，就可以写作：
p ( μ , σ 2 ∣ y ) ∝ N − I n v − χ 2 ( μ , σ 2 ∣ μ n , σ n 2 / κ n ; ν n , σ n 2 ) p(\mu,\sigma^2|y) \propto \mathrm{N-Inv-}\chi^{2}(\mu,\sigma^{2}|\mu_{n},\sigma^{2}_{n}/\kappa_{n}; \nu_{n},\sigma^{2}_{n}) p ( μ , σ 2 ∣ y ) ∝ N − Inv − χ 2 ( μ , σ 2 ∣ μ n , σ n 2 / κ n ; ν n , σ n 2 ) 这时候，我们一样的，可以写出条件后验分布 p ( μ ∣ σ 2 , y ) p(\mu|\sigma^2,y) p ( μ ∣ σ 2 , y )
μ ∣ σ 2 , y ∼ N ( μ n , σ 2 / κ n ) = N ( ( κ 0 / σ 2 ) μ 0 + ( n / σ 2 ) y ˉ κ 0 / σ 2 + n / σ 2 , 1 κ 0 / σ 2 + n / σ 2 ) \begin{aligned}
\mu |\sigma^2, y &\sim N(\mu_{n}, \sigma^2/\kappa_{n}) \\&= N\left( \frac{(\kappa_{0}/\sigma^{2}) \mu_{0} + (n/\sigma^{2}) \bar{y}}{\kappa_{0}/\sigma^{2} + n/\sigma^{2}}, \frac{1}{\kappa_{0}/\sigma^{2} + n/\sigma^{2}} \right)
\end{aligned} μ ∣ σ 2 , y ∼ N ( μ n , σ 2 / κ n ) = N ( κ 0 / σ 2 + n / σ 2 ( κ 0 / σ 2 ) μ 0 + ( n / σ 2 ) y ˉ , κ 0 / σ 2 + n / σ 2 1 ) 我们知道，σ 2 \sigma^2 σ 2 χ 2 \chi^2 χ 2
σ 2 ∣ y ∼ I n v − χ 2 ( ν n , σ n 2 ) \sigma^{2}|y \sim \mathrm{Inv-}\chi^{2}(\nu_{n}, \sigma_{n}^{2}) σ 2 ∣ y ∼ Inv − χ 2 ( ν n , σ n 2 ) 随后抽样的方法如前一样，我们先从标准差的边缘后验分布抽样得到 σ 2 \sigma^2 σ 2 σ 2 \sigma^2 σ 2 μ \mu μ
p ( μ ∣ y ) ∝ ( 1 + κ n ( μ − μ n ) 2 ν n σ n 2 ) − ( ν n + 1 ) / 2 = t ν n ( μ ∣ μ n , σ n 2 / κ n ) \begin{aligned}
p(\mu|y) &\propto \left( 1 + \frac{\kappa_n(\mu - \mu_n)^2}{\nu_n \sigma_n^2} \right)^{-(\nu_n+1)/2} \\
&= t_{\nu_n}(\mu|\mu_n, \sigma_n^2/\kappa_n)
\end{aligned} p ( μ ∣ y ) ∝ ( 1 + ν n σ n 2 κ n ( μ − μ n ) 2 ) − ( ν n + 1 ) /2 = t ν n ( μ ∣ μ n , σ n 2 / κ n ) 我们知道当结果只有两个时，y i = ( 0 , 1 ) y_i = (0,1) y i = ( 0 , 1 )
p ( y ∣ θ ) ∝ ∏ j = 1 k θ j y j p(y|\theta) \propto \prod_{j=1}^{k} \theta^{y_j}_{j} p ( y ∣ θ ) ∝ j = 1 ∏ k θ j y j 其中， ∑ j = 1 k θ j = 1 \sum^k_{j=1}\theta_j = 1 ∑ j = 1 k θ j = 1 ∑ j = 1 k y j = n \sum^k_{j=1}y_j = n ∑ j = 1 k y j = n
那么，此时共轭先验分布就是广义的Beta分布，又称为狄利克雷(Dirichlet)分布。
是的，就是那个德国数学家狄利克雷，其提出了狄利克雷函数，验证了连续不可导的函数，为数学分析里典型的病态函数。
Dirichlet Distribution:
记作：θ ∼ Direchlet ( α 1 , … , α k ) \theta \sim \operatorname{Direchlet}(\alpha_1,\dots,\alpha_k) θ ∼ Direchlet ( α 1 , … , α k )
密度函数：
p ( θ ) = Γ ( α 1 + ⋯ + α k ) Γ ( α 1 ) ⋯ Γ ( α k ) θ 1 α 1 − 1 … θ k α k − 1 , θ 1 … θ k > 0 , ∑ j = 1 k θ j = 1 p(\theta) = \frac{\Gamma(\alpha_1 + \cdots + \alpha_k)}{\Gamma(\alpha_1) \cdots \Gamma(\alpha_k)} \theta_1^{\alpha_1 - 1} \dots \theta_k^{\alpha_k - 1}, \quad \theta_1 \dots \theta_k >0,\quad \sum^k_{j=1}\theta_j = 1 p ( θ ) = Γ ( α 1 ) ⋯ Γ ( α k ) Γ ( α 1 + ⋯ + α k ) θ 1 α 1 − 1 … θ k α k − 1 , θ 1 … θ k > 0 , j = 1 ∑ k θ j = 1 期望：E ( θ j ) = α j α 0 E(\theta_j) = \frac{\alpha_j}{\alpha_0} \\ E ( θ j ) = α 0 α j
方差：var ( θ j ) = α j ( α 0 − α j ) α 0 2 ( α 0 + 1 ) \operatorname{var}(\theta_j) = \frac{\alpha_j(\alpha_0-\alpha_j)}{\alpha_0^2(\alpha_0+1)} \\ var ( θ j ) = α 0 2 ( α 0 + 1 ) α j ( α 0 − α j )
cov ( θ i , θ j ) = − α i α j α 0 2 ( α 0 + 1 ) \operatorname{cov}(\theta_i,\theta_j)= -\frac{\alpha_i\alpha_j}{\alpha_0^2(\alpha_0+1)}\\ cov ( θ i , θ j ) = − α 0 2 ( α 0 + 1 ) α i α j
众数：mode ( θ j ) = α j − 1 α j − k \operatorname{mode}(\theta_j) = \frac{\alpha_j-1}{\alpha_j-k} \\ mode ( θ j ) = α j − k α j − 1
先验分布写作：
p ( θ ∣ α ) ∝ ∏ j = 1 k θ j α j − 1 p(\theta|\alpha) \propto \prod_{j=1}^{k} \theta_j^{\alpha_j-1} p ( θ ∣ α ) ∝ j = 1 ∏ k θ j α j − 1 这里需要主要的是，这里所有的 θ \theta θ α \alpha α ∑ j = 1 k ( α j − 1 ) \sum_{j=1}^{k} (\alpha_j - 1) ∑ j = 1 k ( α j − 1 ) j j j α \alpha α j j j α j = 1 \alpha_j=1 α j = 1 α j = 0 \alpha_j=0 α j = 0 k k k y y y
在这一节里，这个在很多贝叶斯的模型里是十分重要的，在原著中标名为“ central activity of much applied statistical work“
因为现在是多元正态模型，因此我们有了很多的变量，这些变量都符合正态分布，因此，他们之间就有了不同的均值和不同的方差。
本文将讨论的基本模型涉及一个具有 d 个分量的可观测向量 y，该向量服从多元正态分布。
y ∣ μ , Σ ∼ N ( μ , Σ ) y | \mu , \Sigma \sim N( \mu , \Sigma ) y ∣ μ , Σ ∼ N ( μ , Σ ) 这里的 μ \mu μ d d d Σ \Sigma Σ d × d d\times d d × d
对于单一观测值，其似然可以写作：
p ( y ∣ μ , Σ ) ∝ ∣ Σ ∣ − 1 / 2 exp ( − 1 2 ( y − μ ) T Σ − 1 ( y − μ ) ) p(y|μ,\Sigma) ∝ |\Sigma|^{-1/2} \exp\left(-\tfrac{1}{2}(y-μ)^T \Sigma^{-1} (y-μ)\right) p ( y ∣ μ , Σ ) ∝ ∣Σ ∣ − 1/2 exp ( − 2 1 ( y − μ ) T Σ − 1 ( y − μ ) ) 更一般的，对于多个观测值，其观测值之间是互相独立的，则有
p ( y 1 , … , y n ∣ μ , Σ ) ∝ ∣ Σ ∣ − n / 2 exp ( − 1 2 ∑ i = 1 n ( y i − μ ) T Σ − 1 ( y i − μ ) ) = ∣ Σ ∣ − n / 2 exp ( − 1 2 t r ( Σ − 1 S 0 ) ) \begin{aligned}
p(y_1,\dots,y_n\mid \mu, \Sigma) &\propto |\Sigma|^{-n/2}\exp\left(-\tfrac{1}{2}\sum_{i=1}^{n}(y_i-\mu)^T\Sigma^{-1}(y_i-\mu)\right) \\
&= |\Sigma|^{-n/2}\exp\left(-\tfrac{1}{2}\mathrm{tr}\left(\Sigma^{-1}S_0\right)\right)
\end{aligned} p ( y 1 , … , y n ∣ μ , Σ ) ∝ ∣Σ ∣ − n /2 exp ( − 2 1 i = 1 ∑ n ( y i − μ ) T Σ − 1 ( y i − μ ) ) = ∣Σ ∣ − n /2 exp ( − 2 1 tr ( Σ − 1 S 0 ) ) 其中
S 0 = ∑ i = 1 n ( y i − μ ) ( y i − μ ) T S_0 = \sum^n_{i=1}(y_i-\mu)(y_i-\mu)^T S 0 = i = 1 ∑ n ( y i − μ ) ( y i − μ ) T 因为我们现在已知方差，我们对 μ \mu μ μ \mu μ
于是乎，我们写出先验：
μ ∼ N ( μ 0 , Λ 0 ) \mu \sim \operatorname{N}(\mu_0, \Lambda_0) μ ∼ N ( μ 0 , Λ 0 ) 于是乎，我们就可以写出后验分布：
p ( μ ∣ y , Σ ) ∝ exp ( − 1 2 ( ( μ − μ 0 ) T Λ 0 − 1 ( μ − μ 0 ) + ∑ i = 1 n ( y i − μ ) T Σ − 1 ( y i − μ ) ) ) p(\mu|y,\Sigma) \propto \exp\left(-\tfrac{1}{2}\left(\left(\mu - \mu_0\right)^T \Lambda_0^{-1} \left(\mu - \mu_0\right) + \sum_{i=1}^{n} (y_i - \mu)^T \Sigma^{-1} (y_i - \mu)\right)\right) p ( μ ∣ y , Σ ) ∝ exp ( − 2 1 ( ( μ − μ 0 ) T Λ 0 − 1 ( μ − μ 0 ) + i = 1 ∑ n ( y i − μ ) T Σ − 1 ( y i − μ ) ) ) 由于：
( μ − μ 0 ) T Λ 0 − 1 ( μ − μ 0 ) = μ T Λ 0 − 1 μ − 2 μ T Λ 0 − 1 μ 0 + const , ∑ i = 1 n ( y i − μ ) T Σ − 1 ( y i − μ ) = n μ T Σ μ − 2 μ T Σ − 1 ∑ i = 1 n y i + const (\mu-\mu_0)^T \Lambda_0^{-1}(\mu-\mu_0) = \mu^T \Lambda_0^{-1}\mu-2\mu^T\Lambda_0^{-1}\mu_0+\text{const}, \\
\sum_{i=1}^n(y_i-\mu)^T\Sigma^{-1}(y_i-\mu)=n\mu^T\Sigma\mu-2\mu^T\Sigma^{-1}\sum_{i=1}^n y_i+\text{const} ( μ − μ 0 ) T Λ 0 − 1 ( μ − μ 0 ) = μ T Λ 0 − 1 μ − 2 μ T Λ 0 − 1 μ 0 + const , i = 1 ∑ n ( y i − μ ) T Σ − 1 ( y i − μ ) = n μ T Σ μ − 2 μ T Σ − 1 i = 1 ∑ n y i + const 又因为 (\sum_i y_i=n\bar y)，两项相加得到
μ T ( Λ 0 − 1 + n Σ − 1 ) μ − 2 μ T ( Λ 0 − 1 μ 0 + n Σ − 1 y ˉ ) + const . \mu^T(\Lambda_0^{-1}+n\Sigma^{-1})\mu -2\mu^T(\Lambda_0^{-1}\mu_0+n\Sigma^{-1}\bar y) +\text{const}. μ T ( Λ 0 − 1 + n Σ − 1 ) μ − 2 μ T ( Λ 0 − 1 μ 0 + n Σ − 1 y ˉ ) + const . 令：
Λ n − 1 = Λ 0 − 1 + n Σ − 1 μ n = ( Λ 0 − 1 + n Σ − 1 ) − 1 ( Λ 0 − 1 μ 0 + n Σ − 1 y ˉ ) . \Lambda_n^{-1}=\Lambda_0^{-1}+n\Sigma^{-1} \\
\mu_n=(\Lambda_0^{-1}+n\Sigma^{-1})^{-1}(\Lambda_0^{-1}\mu_0+n\Sigma^{-1}\bar y). Λ n − 1 = Λ 0 − 1 + n Σ − 1 μ n = ( Λ 0 − 1 + n Σ − 1 ) − 1 ( Λ 0 − 1 μ 0 + n Σ − 1 y ˉ ) . 于是乎，线性项满足：
Λ 0 − 1 μ 0 + n Σ − 1 y ˉ = Λ n − 1 μ n \Lambda_0^{-1}\mu_0+n\Sigma^{-1}\bar y=\Lambda_n^{-1}\mu_n Λ 0 − 1 μ 0 + n Σ − 1 y ˉ = Λ n − 1 μ n 代回原式，
μ T Λ n − 1 μ − 2 μ T Λ n − 1 μ n = ( μ − μ n ) T Λ n − 1 ( μ − μ n ) − μ n T Λ n − 1 μ n \begin{aligned}
&\mu^T\Lambda_n^{-1}\mu-2\mu^T\Lambda_n^{-1}\mu_n \\&\quad =(\mu-\mu_n)^T\Lambda_n^{-1}(\mu-\mu_n)-\mu_n^T\Lambda_n^{-1}\mu_n
\end{aligned} μ T Λ n − 1 μ − 2 μ T Λ n − 1 μ n = ( μ − μ n ) T Λ n − 1 ( μ − μ n ) − μ n T Λ n − 1 μ n 最后一项不含 (\mu)，在 (\propto) 中被吸收到归一化常数。因此
p ( μ ∣ y , Σ ) ∝ exp [ − 1 2 ( μ − μ n ) T Λ n − 1 ( μ − μ n ) ] p(\mu\mid y,\Sigma)\propto\exp\left[-\frac12(\mu-\mu_n)^T\Lambda_n^{-1}(\mu-\mu_n)\right] p ( μ ∣ y , Σ ) ∝ exp [ − 2 1 ( μ − μ n ) T Λ n − 1 ( μ − μ n ) ] 于是乎有：
p ( μ ∣ y , Σ ) ∝ exp ( − 1 2 ( μ − μ n ) T Λ n − 1 ( μ − μ n ) ) = N ( μ ∣ μ n , Λ n ) \begin{aligned}
p(\mu|y,\Sigma) &\propto \exp\left(-\tfrac{1}{2}(\mu-\mu_n)^T\Lambda_n^{-1}(\mu-\mu_n)\right) \\
&= N(\mu|\mu_n,\Lambda_n)
\end{aligned} p ( μ ∣ y , Σ ) ∝ exp ( − 2 1 ( μ − μ n ) T Λ n − 1 ( μ − μ n ) ) = N ( μ ∣ μ n , Λ n ) 这跟我们之前看到的单变量的正态模型很像，后验均值是数据和先验均值的加权平均值，权重由数据和先验精度矩阵给出。后验精度是先验精度和数据精度之和。
我们很容易知道，参数子集 μ ( 1 ) \mu^{(1)} μ ( 1 ) 边缘后验分布 也是多元正态模型，均值向量等于后验均值向量的相应子向量 μ n \mu_n μ n Λ n \Lambda_n Λ n μ ( 2 ) \mu^{(2)} μ ( 2 ) μ ( 1 ) \mu^{(1)} μ ( 1 )
μ ( 1 ) ∣ μ ( 2 ) , y ∼ N ( μ n ( 1 ) + β 1 ∣ 2 ( μ ( 2 ) − μ n ( 2 ) ) , Λ 1 ∣ 2 ) \mu^{(1)} \mid \mu^{(2)}, y \sim \mathrm{N}\left( \mu_n^{(1)} + \beta^{1|2}\left( \mu^{(2)} - \mu_n^{(2)} \right), \Lambda^{1|2} \right) μ ( 1 ) ∣ μ ( 2 ) , y ∼ N ( μ n ( 1 ) + β 1∣2 ( μ ( 2 ) − μ n ( 2 ) ) , Λ 1∣2 ) 对于回归系数 β 1 ∣ 2 \beta^{1\mid2} β 1 ∣ 2 Λ 1 ∣ 2 \Lambda^{1\mid2} Λ 1 ∣ 2
β 1 ∣ 2 = Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ 1 ∣ 2 = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) \beta^{1\mid2} = \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\\
\Lambda^{1\mid2} = \Lambda_n^{(11)} - \Lambda_n^{(12)} \left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)} β 1 ∣ 2 = Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ 1 ∣ 2 = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 )
推导 :
将后验均值和后验方差矩阵按照 μ ( 1 ) \mu^{(1)} μ ( 1 ) μ ( 2 ) \mu^{(2)} μ ( 2 )
μ = ( μ ( 1 ) μ ( 2 ) ) , μ n = ( μ n ( 1 ) μ n ( 2 ) ) , Λ n = ( Λ n ( 11 ) Λ n ( 12 ) Λ n ( 21 ) Λ n ( 22 ) ) \mu =
\begin{pmatrix}
\mu^{(1)}\\
\mu^{(2)}
\end{pmatrix},\qquad
\mu_n =
\begin{pmatrix}
\mu_n^{(1)}\\
\mu_n^{(2)}
\end{pmatrix},\qquad
\Lambda_n =
\begin{pmatrix}
\Lambda_n^{(11)} & \Lambda_n^{(12)}\\
\Lambda_n^{(21)} & \Lambda_n^{(22)}
\end{pmatrix} μ = ( μ ( 1 ) μ ( 2 ) ) , μ n = ( μ n ( 1 ) μ n ( 2 ) ) , Λ n = ( Λ n ( 11 ) Λ n ( 21 ) Λ n ( 12 ) Λ n ( 22 ) ) 由于
μ ∣ y , Σ ∼ N ( μ n , Λ n ) \mu\mid y,\Sigma \sim \mathrm{N}(\mu_n,\Lambda_n) μ ∣ y , Σ ∼ N ( μ n , Λ n ) 所以
E ( μ ( 1 ) ∣ y ) = μ n ( 1 ) , E ( μ ( 2 ) ∣ y ) = μ n ( 2 ) \operatorname{E}(\mu^{(1)}\mid y)=\mu_n^{(1)},\qquad
\operatorname{E}(\mu^{(2)}\mid y)=\mu_n^{(2)} E ( μ ( 1 ) ∣ y ) = μ n ( 1 ) , E ( μ ( 2 ) ∣ y ) = μ n ( 2 ) 并且
Var ( μ ( 1 ) ∣ y ) = Λ n ( 11 ) , Var ( μ ( 2 ) ∣ y ) = Λ n ( 22 ) , Cov ( μ ( 1 ) , μ ( 2 ) ∣ y ) = Λ n ( 12 ) \operatorname{Var}(\mu^{(1)}\mid y)=\Lambda_n^{(11)},\qquad
\operatorname{Var}(\mu^{(2)}\mid y)=\Lambda_n^{(22)},\qquad
\operatorname{Cov}(\mu^{(1)},\mu^{(2)}\mid y)=\Lambda_n^{(12)} Var ( μ ( 1 ) ∣ y ) = Λ n ( 11 ) , Var ( μ ( 2 ) ∣ y ) = Λ n ( 22 ) , Cov ( μ ( 1 ) , μ ( 2 ) ∣ y ) = Λ n ( 12 ) 为了得到 μ ( 1 ) \mu^{(1)} μ ( 1 ) μ ( 2 ) \mu^{(2)} μ ( 2 )
ε = μ ( 1 ) − μ n ( 1 ) − B ( μ ( 2 ) − μ n ( 2 ) ) \varepsilon=
\mu^{(1)}-\mu_n^{(1)}-
B\left(\mu^{(2)}-\mu_n^{(2)}\right) ε = μ ( 1 ) − μ n ( 1 ) − B ( μ ( 2 ) − μ n ( 2 ) ) 我们希望这个残差与 μ ( 2 ) \mu^{(2)} μ ( 2 )
Cov ( ε , μ ( 2 ) ∣ y ) = Cov ( μ ( 1 ) , μ ( 2 ) ∣ y ) − B Var ( μ ( 2 ) ∣ y ) = Λ n ( 12 ) − B Λ n ( 22 ) \begin{aligned}
\operatorname{Cov}(\varepsilon,\mu^{(2)}\mid y)
&=
\operatorname{Cov}(\mu^{(1)},\mu^{(2)}\mid y)-
B\operatorname{Var}(\mu^{(2)}\mid y)\\
&=
\Lambda_n^{(12)}-B\Lambda_n^{(22)}
\end{aligned} Cov ( ε , μ ( 2 ) ∣ y ) = Cov ( μ ( 1 ) , μ ( 2 ) ∣ y ) − B Var ( μ ( 2 ) ∣ y ) = Λ n ( 12 ) − B Λ n ( 22 ) 令其等于 0 0 0
B Λ n ( 22 ) = Λ n ( 12 ) B\Lambda_n^{(22)}=\Lambda_n^{(12)} B Λ n ( 22 ) = Λ n ( 12 ) 因此
B = Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 B=\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1} B = Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 即回归系数：
β 1 ∣ 2 = Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 \beta^{1\mid2}=
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1} β 1 ∣ 2 = Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 因为 ( μ ( 1 ) , μ ( 2 ) ) (\mu^{(1)},\mu^{(2)}) ( μ ( 1 ) , μ ( 2 ) ) B B B ε \varepsilon ε μ ( 2 ) \mu^{(2)} μ ( 2 ) μ ( 1 ) \mu^{(1)} μ ( 1 )
μ ( 1 ) = μ n ( 1 ) + β 1 ∣ 2 ( μ ( 2 ) − μ n ( 2 ) ) + ε \mu^{(1)}=
\mu_n^{(1)}+
\beta^{1\mid2}\left(\mu^{(2)}-\mu_n^{(2)}\right)+
\varepsilon μ ( 1 ) = μ n ( 1 ) + β 1 ∣ 2 ( μ ( 2 ) − μ n ( 2 ) ) + ε 接下来只需要计算残差 ε \varepsilon ε
ε = μ ( 1 ) − μ n ( 1 ) − β 1 ∣ 2 ( μ ( 2 ) − μ n ( 2 ) ) \varepsilon=
\mu^{(1)}-\mu_n^{(1)}-
\beta^{1\mid2}\left(\mu^{(2)}-\mu_n^{(2)}\right) ε = μ ( 1 ) − μ n ( 1 ) − β 1 ∣ 2 ( μ ( 2 ) − μ n ( 2 ) ) 得
Var ( ε ∣ y ) = Var ( μ ( 1 ) ∣ y ) − β 1 ∣ 2 Cov ( μ ( 2 ) , μ ( 1 ) ∣ y ) − Cov ( μ ( 1 ) , μ ( 2 ) ∣ y ) ( β 1 ∣ 2 ) T + β 1 ∣ 2 Var ( μ ( 2 ) ∣ y ) ( β 1 ∣ 2 ) T = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) + Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 22 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) \begin{aligned}
\operatorname{Var}(\varepsilon\mid y)
&=
\operatorname{Var}(\mu^{(1)}\mid y)-
\beta^{1\mid2}\operatorname{Cov}(\mu^{(2)},\mu^{(1)}\mid y)\\
&\quad-
\operatorname{Cov}(\mu^{(1)},\mu^{(2)}\mid y)(\beta^{1\mid2})^T+
\beta^{1\mid2}\operatorname{Var}(\mu^{(2)}\mid y)(\beta^{1\mid2})^T\\
&=
\Lambda_n^{(11)}-
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}-
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}\\
&\quad+
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(22)}
\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}\\
&=
\Lambda_n^{(11)}-
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}
\end{aligned} Var ( ε ∣ y ) = Var ( μ ( 1 ) ∣ y ) − β 1 ∣ 2 Cov ( μ ( 2 ) , μ ( 1 ) ∣ y ) − Cov ( μ ( 1 ) , μ ( 2 ) ∣ y ) ( β 1 ∣ 2 ) T + β 1 ∣ 2 Var ( μ ( 2 ) ∣ y ) ( β 1 ∣ 2 ) T = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) + Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 22 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) 因此条件方差矩阵为：
Λ 1 ∣ 2 = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) \Lambda^{1\mid2}=
\Lambda_n^{(11)}-
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)} Λ 1 ∣ 2 = Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) 最后，由于 ε \varepsilon ε μ ( 2 ) \mu^{(2)} μ ( 2 ) μ ( 2 ) \mu^{(2)} μ ( 2 ) ε \varepsilon ε
μ ( 1 ) ∣ μ ( 2 ) , y ∼ N ( μ n ( 1 ) + Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 ( μ ( 2 ) − μ n ( 2 ) ) , Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) ) \mu^{(1)}\mid \mu^{(2)},y
\sim
\mathrm{N}\left(
\mu_n^{(1)}+
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}
\left(\mu^{(2)}-\mu_n^{(2)}\right),
\Lambda_n^{(11)}-
\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}
\right) μ ( 1 ) ∣ μ ( 2 ) , y ∼ N ( μ n ( 1 ) + Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 ( μ ( 2 ) − μ n ( 2 ) ) , Λ n ( 11 ) − Λ n ( 12 ) ( Λ n ( 22 ) ) − 1 Λ n ( 21 ) )
我们一样能写出后验分布的表达式：
p ( y ~ , μ ∣ y ) = p ( y ~ ∣ y ) p ( μ ∣ y ) = N ( y ~ ∣ μ , Σ ) N ( μ ∣ μ n , Λ n ) \begin{aligned}
p(\tilde y,\mu|y) &= p(\tilde y|y) p(\mu|y) \\
&= \mathrm{N}(\tilde y|\mu,\Sigma) \mathrm{N}(\mu|\mu_n,\Lambda_n)
\end{aligned} p ( y ~ , μ ∣ y ) = p ( y ~ ∣ y ) p ( μ ∣ y ) = N ( y ~ ∣ μ , Σ ) N ( μ ∣ μ n , Λ n ) 于是乎，我们也可以写出 y ~ \tilde y y ~
E ( y ~ ∣ y ) = E ( E ( y ~ ∣ μ , y ) ∣ y ) = E ( μ ∣ y ) = μ n v a r ( y ~ ∣ y ) = E ( v a r ( y ~ ∣ μ , y ) ∣ y ) + v a r ( E ( y ~ ∣ μ , y ) ∣ y ) = E ( Σ ∣ y ) + v a r ( μ ∣ y ) = Σ + Λ n E(\tilde{y}|y) = E(E(\tilde{y}|\mu, y)|y) = E(\mu|y) = \mu_n \\
\mathrm{var}(\tilde{y}|y) = E(\mathrm{var}(\tilde{y}|\mu, y)|y) + \mathrm{var}(E(\tilde{y}|\mu, y)|y) = E(\Sigma|y) + var(\mu|y) = \Sigma + \Lambda_n E ( y ~ ∣ y ) = E ( E ( y ~ ∣ μ , y ) ∣ y ) = E ( μ ∣ y ) = μ n var ( y ~ ∣ y ) = E ( var ( y ~ ∣ μ , y ) ∣ y ) + var ( E ( y ~ ∣ μ , y ) ∣ y ) = E ( Σ∣ y ) + v a r ( μ ∣ y ) = Σ + Λ n 在这里，我们对 μ \mu μ p ( μ ) ∝ C o n s t a n t p(\mu) \propto \mathrm{Constant} p ( μ ) ∝ Constant ∣ Λ 0 − 1 ∣ → 0 \left|\Lambda_0^{-1}\right| \to 0 Λ 0 − 1 → 0 n ≥ d n\ge d n ≥ d S 0 S_0 S 0
在之前的单元正态模型中，我们使用了 Scaled Inverse-χ 2 \chi^2 χ 2 χ 2 \chi^2 χ 2 χ 2 \chi^2 χ 2 Inverse-Wishart 分布 。
Inverse-Wishart Distribution
记作：W ∼ I n v − W i s h a r t ν ( S − 1 ) W \sim \mathrm{Inv-Wishart}_{\nu}(S^{-1}) W ∼ Inv − Wishart ν ( S − 1 )
其中， ν \nu ν S S S k × k k \times k k × k
密度方程：
p ( W ) = ( 2 ν k / 2 π k ( k − 1 ) / 4 ∏ i = 1 k Γ ( ν + 1 − i 2 ) ) − 1 × ∣ S ∣ ν / 2 ∣ W ∣ − ( ν + k + 1 ) / 2 × exp ( − 1 2 t r ( S W − 1 ) ) , W pos. definite p(W) = \left( 2^{\nu k/2}\pi^{k(k-1)/4} \prod_{i=1}^{k} \Gamma\left( \frac{\nu+1-i}{2} \right) \right)^{-1} \times |S|^{\nu/2}|W|^{- (\nu+k+1)/2} \times \exp\left(-\tfrac{1}{2} \mathrm{tr}(SW^{-1})\right), \; W \text{ pos. definite} p ( W ) = ( 2 ν k /2 π k ( k − 1 ) /4 i = 1 ∏ k Γ ( 2 ν + 1 − i ) ) − 1 × ∣ S ∣ ν /2 ∣ W ∣ − ( ν + k + 1 ) /2 × exp ( − 2 1 tr ( S W − 1 ) ) , W pos. definite 期望：E ( W ) = ( ν − k − 1 ) − 1 S E(W) = (\nu - k - 1)^{-1} S E ( W ) = ( ν − k − 1 ) − 1 S
威沙特分布看着确实很吓人
写出我们的共轭先验：
Σ ∼ Inv-Wishart ν 0 ( Λ 0 − 1 ) μ ∣ Σ ∼ N ( μ 0 , Σ / κ 0 ) \Sigma \sim \text{Inv-Wishart}_{\nu_0}(\Lambda_0^{-1})\\
\mu|\Sigma \sim \mathrm{N}(\mu_0, \Sigma/\kappa_0) Σ ∼ Inv-Wishart ν 0 ( Λ 0 − 1 ) μ ∣Σ ∼ N ( μ 0 , Σ/ κ 0 ) 于是乎，我们就可以写出联合先验分布 p ( μ , Σ ) p(\mu,\Sigma) p ( μ , Σ )
p ( μ , Σ ) ∝ ∣ Σ ∣ − ( ( ν 0 + d ) / 2 + 1 ) exp ( − 1 2 t r ( Λ 0 Σ − 1 ) − κ 0 2 ( μ − μ 0 ) T Σ − 1 ( μ − μ 0 ) ) p(\mu, \Sigma) \propto |\Sigma|^{-((\nu_0 + d)/2 + 1)} \exp\left(-\frac{1}{2} \mathrm{tr}(\Lambda_0 \Sigma^{-1}) - \frac{\kappa_0}{2}(\mu - \mu_0)^T \Sigma^{-1} (\mu - \mu_0)\right) p ( μ , Σ ) ∝ ∣Σ ∣ − (( ν 0 + d ) /2 + 1 ) exp ( − 2 1 tr ( Λ 0 Σ − 1 ) − 2 κ 0 ( μ − μ 0 ) T Σ − 1 ( μ − μ 0 ) ) 类似的，我们记作： μ , Σ ∼ N − I n v − W i s h a r t ( μ 0 , Λ 0 / κ 0 ; ν 0 , Λ 0 ) \mu,\Sigma \sim \mathrm{N-Inv-Wishart}(\mu_0, \Lambda_0/\kappa_0;\nu_0,\Lambda_0) μ , Σ ∼ N − Inv − Wishart ( μ 0 , Λ 0 / κ 0 ; ν 0 , Λ 0 )
同样的，这里我们也有四个未知参数，我们就可以写出后验分布的表达式：
类似的，我们最终得到的后验分布，也会是与先验同样的分布，我们写作：
μ , Σ ∣ y ∼ N − I n v − W i s h a r t ( μ n , Λ n / κ n ; ν n , Λ n ) \mu,\Sigma |y \sim \mathrm{N-Inv-Wishart}(\mu_n, \Lambda_n/\kappa_n;\nu_n,\Lambda_n) μ , Σ∣ y ∼ N − Inv − Wishart ( μ n , Λ n / κ n ; ν n , Λ n ) 与前面提到的未知均值和方差的一元正态模型类似，我们可以写出他们的表达式：
μ n = κ 0 κ 0 + n μ 0 + n κ 0 + n y ˉ κ n = κ 0 + n ν n = ν 0 + n Λ n = Λ 0 + S + κ 0 n κ 0 + n ( y ˉ − μ 0 ) ( y ˉ − μ 0 ) T \begin{aligned}
\mu_n &= \frac{\kappa_0}{\kappa_0 + n}\mu_0 + \frac{n}{\kappa_0 + n}\bar{y} \\
\kappa_n &= \kappa_0 + n \\
\nu_n &= \nu_0 + n \\
\Lambda_n &= \Lambda_0 + S + \frac{\kappa_0 n}{\kappa_0 + n}(\bar{y} - \mu_0)(\bar{y} - \mu_0)^T
\end{aligned} μ n κ n ν n Λ n = κ 0 + n κ 0 μ 0 + κ 0 + n n y ˉ = κ 0 + n = ν 0 + n = Λ 0 + S + κ 0 + n κ 0 n ( y ˉ − μ 0 ) ( y ˉ − μ 0 ) T 其中： S S S
S = ∑ i = 1 n ( y i − y ˉ ) ( y i − y ˉ ) T S = \sum_{i=1}^{n} (y_i - \bar{y})(y_i - \bar{y})^T S = i = 1 ∑ n ( y i − y ˉ ) ( y i − y ˉ ) T 于是乎，在未知均值和未知方差矩阵的情况下，我们仍然可以通过共轭结构得到一套比较完整的后验推断结果。
在联合后验分布中，μ \mu μ
μ ∣ Σ , y ∼ N ( μ n , Σ / κ n ) \mu|\Sigma,y \sim \mathrm{N}(\mu_n,\Sigma/\kappa_n) μ ∣Σ , y ∼ N ( μ n , Σ/ κ n ) 而 Σ ∣ y \Sigma|y Σ∣ y μ \mu μ Σ \Sigma Σ Σ \Sigma Σ μ \mu μ t t t
μ ∣ y ∼ t ν n − d + 1 ( μ n , Λ n κ n ( ν n − d + 1 ) ) \mu|y \sim t_{\nu_n-d+1}\left(\mu_n,\frac{\Lambda_n}{\kappa_n(\nu_n-d+1)}\right) μ ∣ y ∼ t ν n − d + 1 ( μ n , κ n ( ν n − d + 1 ) Λ n ) 这里的自由度是 ν n − d + 1 \nu_n-d+1 ν n − d + 1 μ \mu μ t t t 矩阵维度 d d d
对于一个新的观测值 y ~ \tilde{y} y ~
y ~ ∣ μ , Σ , y ∼ N ( μ , Σ ) \tilde{y}|\mu,\Sigma,y \sim \mathrm{N}(\mu,\Sigma) y ~ ∣ μ , Σ , y ∼ N ( μ , Σ ) 在后验预测分布中，y ~ \tilde{y} y ~ μ , Σ \mu,\Sigma μ , Σ μ \mu μ Σ \Sigma Σ y ~ \tilde{y} y ~ t t t
y ~ ∣ y ∼ t ν n − d + 1 ( μ n , ( κ n + 1 ) Λ n κ n ( ν n − d + 1 ) ) \tilde{y}|y \sim t_{\nu_n-d+1}\left(\mu_n,\frac{(\kappa_n+1)\Lambda_n}{\kappa_n(\nu_n-d+1)}\right) y ~ ∣ y ∼ t ν n − d + 1 ( μ n , κ n ( ν n − d + 1 ) ( κ n + 1 ) Λ n ) 与 μ ∣ y \mu|y μ ∣ y κ n + 1 \kappa_n+1 κ n + 1 y ~ \tilde{y} y ~ μ \mu μ
Var ( y ~ ∣ y ) = Var ( μ ∣ y ) + E ( Σ ∣ y ) \operatorname{Var}(\tilde{y}|y) =
\operatorname{Var}(\mu|y) +
\operatorname{E}(\Sigma|y) Var ( y ~ ∣ y ) = Var ( μ ∣ y ) + E ( Σ∣ y ) 它通常会比 μ ∣ y \mu|y μ ∣ y
由于后验分布保持了N-inverse-Wishart 的共轭形式，所以从联合后验分布中抽样非常直接。正如一元正态模型中我们的处理方式一样，只不过，我们现在抽的不再是单一的方差和均值标量：
先抽取方差矩阵：
Σ ∣ y ∼ I n v - W i s h a r t ν n ( Λ n − 1 ) \Sigma|y \sim \mathrm{Inv\text{-}Wishart}_{\nu_n}(\Lambda_n^{-1}) Σ∣ y ∼ Inv - Wishart ν n ( Λ n − 1 )
在已经抽到的 Σ \Sigma Σ
μ ∣ Σ , y ∼ N ( μ n , Σ / κ n ) \mu|\Sigma,y \sim \mathrm{N}(\mu_n,\Sigma/\kappa_n) μ ∣Σ , y ∼ N ( μ n , Σ/ κ n )
如果还想从后验预测分布中抽取一个新的观测值，则在每一组已经抽到的 ( μ , Σ ) (\mu,\Sigma) ( μ , Σ )
y ~ ∣ μ , Σ , y ∼ N ( μ , Σ ) \tilde{y}|\mu,\Sigma,y \sim \mathrm{N}(\mu,\Sigma) y ~ ∣ μ , Σ , y ∼ N ( μ , Σ ) 因此，一次完整的后验预测模拟可以写成：
Σ ( s ) ∣ y → μ ( s ) ∣ Σ ( s ) , y → y ~ ( s ) ∣ μ ( s ) , Σ ( s ) , y \Sigma^{(s)}|y \rightarrow \mu^{(s)}|\Sigma^{(s)},y \rightarrow \tilde{y}^{(s)}|\mu^{(s)},\Sigma^{(s)},y Σ ( s ) ∣ y → μ ( s ) ∣ Σ ( s ) , y → y ~ ( s ) ∣ μ ( s ) , Σ ( s ) , y 上面使用的是共轭先验：
Σ ∼ I n v - W i s h a r t ν 0 ( Λ 0 − 1 ) , μ ∣ Σ ∼ N ( μ 0 , Σ / κ 0 ) \Sigma \sim \mathrm{Inv\text{-}Wishart}_{\nu_0}(\Lambda_0^{-1}),\qquad
\mu|\Sigma \sim \mathrm{N}(\mu_0,\Sigma/\kappa_0) Σ ∼ Inv - Wishart ν 0 ( Λ 0 − 1 ) , μ ∣Σ ∼ N ( μ 0 , Σ/ κ 0 ) 但在实际建模时，我们也常常希望使用更弱的信息先验 。多元正态模型中常见的非信息先验 主要有下面几种。
第一种做法是用 d + 1 d+1 d + 1
Σ ∼ I n v - W i s h a r t d + 1 ( I ) \Sigma \sim \mathrm{Inv\text{-}Wishart}_{d+1}(I) Σ ∼ Inv - Wishart d + 1 ( I ) 这个设定有一个比较好的性质：Σ \Sigma Σ
第二种做法是使用多元 Jeffreys 先验 (有 d − 1 d-1 d − 1
p ( μ , Σ ) ∝ ∣ Σ ∣ − ( d + 1 ) / 2 p(\mu,\Sigma)\propto |\Sigma|^{-(d+1)/2} p ( μ , Σ ) ∝ ∣Σ ∣ − ( d + 1 ) /2 它可以看作是共轭先验在下面极限情形下的结果：
κ 0 → 0 , ν 0 → − 1 , ∣ Λ 0 ∣ → 0 \kappa_0\rightarrow 0,\qquad
\nu_0\rightarrow -1,\qquad
|\Lambda_0|\rightarrow 0 κ 0 → 0 , ν 0 → − 1 , ∣ Λ 0 ∣ → 0 在这个先验下，如果后验分布是 proper 的，那么后验可以写成：
Σ ∣ y ∼ I n v - W i s h a r t n − 1 ( S − 1 ) μ ∣ Σ , y ∼ N ( y ˉ , Σ / n ) \Sigma|y \sim \mathrm{Inv\text{-}Wishart}_{n-1}(S^{-1}) \\
\mu|\Sigma,y \sim \mathrm{N}(\bar{y},\Sigma/n) Σ∣ y ∼ Inv - Wishart n − 1 ( S − 1 ) μ ∣Σ , y ∼ N ( y ˉ , Σ/ n ) 进一步把 Σ \Sigma Σ
μ ∣ y ∼ t n − d ( y ˉ , S n ( n − d ) ) \mu|y \sim t_{n-d}\left(\bar{y},\frac{S}{n(n-d)}\right) μ ∣ y ∼ t n − d ( y ˉ , n ( n − d ) S ) 直接给 Σ \Sigma Σ
这个模型把协方差矩阵写成：
Σ = Diag ( ξ ) Σ η Diag ( ξ ) \Sigma=\operatorname{Diag}(\xi)\Sigma_{\eta}\operatorname{Diag}(\xi) Σ = Diag ( ξ ) Σ η Diag ( ξ ) 其中，ξ \xi ξ Σ η \Sigma_{\eta} Σ η Σ η \Sigma_{\eta} Σ η
Σ η ∼ I n v - W i s h a r t d + 1 ( I ) \Sigma_{\eta}\sim \mathrm{Inv\text{-}Wishart}_{d+1}(I) Σ η ∼ Inv - Wishart d + 1 ( I ) 这样可以让相关系数具有比较弱的信息先验；同时再单独给尺度参数 ξ \xi ξ
这种写法的好处是：相关矩阵和方差大小可以被分开建模，不会因为选择了某个 inverse-Wishart 先验而同时对方差和相关性施加过强约束。在层级模型，尤其是随机截距、随机斜率 这类模型中，这种分离方差和相关结构的思路会非常有用。
写出似然
写出后验分布，如果先验分布是高度格式化的，将其纳入。否则使用弱信息先验或者非信息先验分布 （应注意后验分布是否合理）
粗略估计参数 θ \theta θ
抽样模拟，然后计算感兴趣的含 θ \theta θ
预测未知的 y ~ \tilde y y ~
