我想写一段代码，完成最高后验密度区间（Highest Posterior Density Interval）的动画演示。对于了解HPD的人来说，这个动画显得很无聊，不过也许会对初学者有用吧：）
HPD是贝叶斯统计中的一个基础概念，这个HPD区间跟置信区间有类似之处但有本质的不同，这里我就细说了。简要说来，假设现在我们手里已经有了参数的后验分布了，这个分布的形式是给了你一串x和y值，它们给出了后验密度的大小。现在找HPD区间的过程就是，从最高点的密度开始往下“扫描”，直到覆盖的分布区域达到给定的概率（如95%），此时把x的范围拿出来就是HPD的区间了。
我已经写了个框架，里面有个繁琐的地方我不想写了，就是当这个分布式双峰或者多峰的时候，这个区间可能会有多个，我懒得想怎么去找这些区间的“边界”了，也许游程rle()会有用，或者暴力循环。下面的演示有个大问题，就是阴影的标注是不对的，正确的范围应该是从曲线上一直画到y=0的地方（也就是我漏了个矩形的底），那里面才是概率（即积分）。除了解决多峰问题之外，还得加上概率的标注，即：某一块或多块阴影对应的概率是多少。本例最终会收入animation包中，当然我会注明作者。
library(animation)<br />
ani.options(interval = 0.2)<br />
HPD.ani = function(x = seq(-3, 3, length.out = 300), <br />
y = dnorm(x), delta = 0.01, p.col = "lightgray", ...) {<br />
<br />
if (min(y) < 0) <br />
warning("Probability density should be non-negative!")<br />
y0 = max(y)<br />
interval = ani.options("interval")<br />
while (y0 > 0) {<br />
y0 = y0 - delta<br />
plot.new()<br />
plot.window(range(x), c(0, max(y)))<br />
idx = y >= y0<br />
<br />
px = ifelse(idx, x, NA)<br />
<br />
polygon(c(px, rev(px)), c(y, rep(0, length(y))), border = NA, <br />
col = p.col)<br />
<br />
lines(x, y, ...)<br />
Sys.sleep(interval)<br />
}<br />
}<br />
HPD.ani()
以上代码也许有bug。
