马尔可夫“折棍子”过程 Markovian Stick-breaking Process 在直方图平滑中的应用

用Dirichlet prior做Histogram Smoothing
用Markovian Stick-breaking Prior做Histogram Smoothing
- Dirichlet Generator
- Tridiagonal Generator

上一篇介绍了Markovian Stick-breaking Process的构造，这一篇介绍它在非参贝叶斯统计中的一个简单应用——Histogram Smoothing。

用Dirichlet prior做Histogram Smoothing

这个toy example来自WILLIAM LIPPITT AND SUNDER SETHURAMAN的manuscript（ON THE USE OF MARKOVIAN STICK-BREAKING PRIORS）。

假设一个鞋店老板想了解附近居民的鞋码分布，假设有ddd种尺码。在实际走访之前，因为他没有任何关于鞋码分布的先验信息，于是他决定使用non-informative measure：
μ=1⃗d\mu=\frac{\vec 1}{d}μ=d1

在他走访之后，记下了每种鞋码对应的人数，用f⃗\vec ff表示；如果先验为Dirichlet(θμ)Dirichlet(\theta \mu)Dirichlet(θμ)，则对于Multinomial的数据，后验为Dirichlet(θμ+f⃗)Dirichlet(\theta\mu+\vec f)Dirichlet(θμ+f)，后验均值为
θμ+f⃗θ+f⃗T1⃗\frac{\theta \mu+\vec f}{\theta+\vec f^T \vec 1}θ+fT1θμ+f

因此θ\thetaθ可以代表先验的强度，θ\thetaθ越大，先验对后验的影响就越大。

上图展示了Dirichlet prior的缺陷，在处理count data的时候，如果某些类别的数据样本很少（中图），那么在后验均值中，这些类别也不会被填充（右图），所以后验均值的直方图依旧无法展示数据的真实特征，也就是说Dirichlet prior并没有完全发掘出所有类别的先验信息。

用Markovian Stick-breaking Prior做Histogram Smoothing

Dirichlet Generator

要使用MSB(G)MSB(G)MSB(G)作为先验的话，需要先确定generator matrix GGG，因为Dirichlet prior是MSB(G)MSB(G)MSB(G)的一个特例，我们可以先看看Dirichlet prior Dirichlet(α)Dirichlet(\alpha)Dirichlet(α)对应的generator matrix。

用有向图代表Dirichlet prior不同category（在这个例子中就是鞋码）之间的dependence，但是Dirichlet prior中不同category是独立的，因此从节点iii到jjj的有向边的权重只取决于jjj，对于Dirichlet(α)Dirichlet(\alpha)Dirichlet(α)而言，这个权重与αj\alpha_jαj成正比（简单起见可以就用αj\alpha_jαj），因此Dirichlet prior有向图的伴随矩阵为

A=[α1α2α3⋯αdα1α2α3⋯αdα1α2α3⋯αd⋯⋯⋯⋯⋯α1α2α3⋯αd]A=\left[ \begin{matrix} \alpha_1 & \alpha_2 & \alpha_3 & \cdots & \alpha_d \\ \alpha_1 & \alpha_2 & \alpha_3 & \cdots & \alpha_d \\ \alpha_1 & \alpha_2 & \alpha_3& \cdots & \alpha_d \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \alpha_1 & \alpha_2 & \alpha_3 & \cdots & \alpha_d\end{matrix} \right]A=⎣⎢⎢⎢⎢⎡α1α1α1⋯α1α2α2α2⋯α2α3α3α3⋯α3⋯⋯⋯⋯⋯αdαdαd⋯αd⎦⎥⎥⎥⎥⎤

根据Generator的定义把伴随矩阵修正为generator：
G=[α1−∑i=1dαiα2α3⋯αdα1α2−∑i=1dαiα3⋯αdα1α2α3−∑i=1dαi⋯αd⋯⋯⋯⋯⋯α1α2α3⋯αd−∑i=1dαi]G=\left[ \begin{matrix} \alpha_1-\sum_{i=1}^d \alpha_i & \alpha_2 & \alpha_3 & \cdots & \alpha_d \\ \alpha_1 & \alpha_2-\sum_{i=1}^d \alpha_i & \alpha_3 & \cdots & \alpha_d \\ \alpha_1 & \alpha_2 & \alpha_3-\sum_{i=1}^d \alpha_i & \cdots & \alpha_d \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \alpha_1 & \alpha_2 & \alpha_3 & \cdots & \alpha_d -\sum_{i=1}^d \alpha_i \end{matrix} \right]G=⎣⎢⎢⎢⎢⎢⎡α1−∑i=1dαiα1α1⋯α1α2α2−∑i=1dαiα2⋯α2α3α3α3−∑i=1dαi⋯α3⋯⋯⋯⋯⋯αdαdαd⋯αd−∑i=1dαi⎦⎥⎥⎥⎥⎥⎤

这就是Dirichlet prior的generator matrix。

Tridiagonal Generator

根据Dirichlet generator的思想，用有向图表示不同category之间的dependence structure，再把有向图的伴随矩阵改写为generator的形式，这样我们就可以根据需求设计category之间的prior dependence structure。

在鞋码的这个例子中，我们可以假设鞋码其实是根据连续型随机变量——脚的长度——离散化之后的结果，因此鞋码之间的dependence structure应该是：

相邻鞋码之间可能存在dependence，用www来表示这种dependence的权重；
与相邻鞋码之外的其他鞋码之间无dependence

这两个结果非常直观，比如在买鞋的时候，平时穿40鞋的人最多再试一下39或者41，一般不会尝试38或者42；根据我们对鞋码之间dependence的认识，我们可以写出对应有向图的伴随矩阵：

A=[ww00wwww⋯00ww⋯0⋯⋯⋯⋯⋯w00⋯w]A=\left[ \begin{matrix} w & w & 0 & 0 & w \\ w & w & w & \cdots & 0 \\ 0 & w & w& \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ w & 0 & 0 & \cdots & w\end{matrix} \right]A=⎣⎢⎢⎢⎢⎡ww0⋯wwww⋯00ww⋯00⋯⋯⋯⋯w00⋯w⎦⎥⎥⎥⎥⎤

这里我做了一点标准化处理，假设category 1与category d相邻，改写为generator之后它就成了对角元相同的一个很像三对角矩阵的矩阵（因为有右上左下两个www，所以并不是三对角矩阵），这是就称这种generator为tridiagonal generator的原因：

G=[−2ww00ww−2ww⋯00w−2w⋯0⋯⋯⋯⋯⋯w00⋯−2w]G=\left[ \begin{matrix} -2w & w & 0 & 0 & w \\ w & -2w & w & \cdots & 0 \\ 0 & w & -2w& \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ w & 0 & 0 & \cdots & -2w\end{matrix} \right]G=⎣⎢⎢⎢⎢⎡−2ww0⋯ww−2ww⋯00w−2w⋯00⋯⋯⋯⋯w00⋯−2w⎦⎥⎥⎥⎥⎤

下面是原文的两个模拟结果：