机器算法验证 - 寻找一个分布，也许不常见，与两个数据点和专家约束一致？ - 吾爱随笔录

寻找一个分布，也许不常见，与两个数据点和专家约束一致？

机器算法验证 r 分布可能性贝叶斯优化

2022-03-25 18:31:41

我试图说明贝叶斯元分析的先验分布。

我有关于随机变量的以下信息：

两个观察值：3.0、3.6
一位研究该变量的科学家告诉我，并且高达 6 的值具有非零概率。 $P(X<2)=P(X>8)=0$

我使用了以下优化方法（log-N =的模式： $e^{\mu-\sigma^2)}$

prior <- function(parms, x, alpha) {
  a <- abs(plnorm(x[1], parms[1], parms[2]) - (alpha/2))
  b <- abs(plnorm(x[2], parms[1], parms[2]) - (1-alpha/2))
  mode <- exp(parms[1] - parms[2]^2)
  c <- abs(mode-3.3)
  return(a + b + c)
}
v = nlm(prior,c(log(3.3),0.14),alpha=0.05,x=c(2.5,7.5))
x <- seq(1,10,0.1)
plot(x, dlnorm(x, v$estimate[1], v$estimate[2]))
abline(v=c(2.5,7.5), lty=2) #95%CI

替代文字

在图中，您可以看到 this 返回的分布，但我想找到更像我画的红线的东西。

这使用对数正态、伽马或正态提供相同形状的分布，并导致和的分布，即： $P(X=5)<0.05$ $P(X=6)< 0.01$

 plnorm(c(5,6), v$estimate[1],v$estimate[2])

任何人都可以提出替代方案吗？我宁愿坚持单一分布而不是混合分布。

谢谢！

4个回答

如果对我上面的评论给出答案，您希望限制分布的范围，为什么不简单地拟合一个重新调整到单位间隔的 Beta 分布呢？换句话说，如果您知道感兴趣的参数应该在之间，那么为什么不定义。我首先以零为中心，除以宽度，使 Y 的范围为 1，然后添加，使 Y 的范围为。（你可以想到任何一种方式：直接从或从 $[2 , 8]$ $Y = \frac{X - 5}{6} + \frac{1}{2} = \frac{X - 2}{6}$ $\frac{1}{2}$ $[0,1]$ $[2,8] \rightarrow [0,1]$ $[2,8] \rightarrow [-\frac{1}{2},\frac{1}{2}] \rightarrow [0,1]$ ，但我认为后者起初可能更容易）。

然后，使用两个数据点，您可以将 beta 后验与统一的 beta 先验拟合吗？

Kumaraswamy发行版怎么样，它有以下 pdf：

f (x; a, b) = a b x^{a - 1} {(1 - x^{a})}^{b - 1}

$f(x; a,b) = a b x^{a-1}{ (1-x^a)}^{b-1}$ 为了

a > 0

$a>0$ ,

b > 0

$b>0$ ,

0 < x < 1

$0 < x < 1$ . 可以重新调整此分布以获得所需的支持。

由于对数正态分布有两个参数，因此您不能将它令人满意地拟合到三个不适合它的约束。对于 2.5 和 7.5 的极端分位数，众数约为 4，对此您无能为力。由于和的误差规模a远b小于c，因此在优化过程中几乎会忽略其中之一。

为了更好地拟合，您可以选择三参数分布，例如广义 gamma 分布（在VGAM包中实现），或者向对数正态（或 gamma，...）分布添加一个移位参数。

最后一点，由于您正在寻找的分布显然不是对称的，因此两个给定观察值的平均值不是该模式的正确值。我将在 3.0 和 3.6 处最大化密度总和，同时将极端分位数保持在 2.5 和 7.5 - 如果您有三个参数，这是可能的。

您也可以尝试三角分布。为了适应这一点，您基本上指定了一个下限（这将是 X=2）、一个上限（这将是 X=8）和一个“最有可能”的值。维基百科页面http://en.wikipedia.org/wiki/Triangular_distribution有关于这个分布的更多信息。如果对“最有可能”的值没有太大的信心（看起来是，在观察任何数据之前），最好在其上放置一个非信息性的先验分布，然后使用这两个数据点来估计这个值。一个好的方法是 jeffrey 的先验，对于这个问题，它是 p(c) = 1/(pi*sqrt((c-2)*(c-8)))，其中“c”是“最可能的值” "（与维基百科符号一致）。

鉴于此，您可以通过分析或模拟计算出 c 的后验分布。似然的解析形式不是特别好，所以模拟似乎更有吸引力。这个例子特别适合拒绝抽样（有关拒绝抽样的一般描述，请参阅 wiki 页面），因为无论 c 的值如何，最大化的可能性都是 1/3^n，它提供了“上限”。因此，您从 jeffrey 的先验（称为 c_i）中生成一个“候选”，然后评估该候选 L(x1,..,xn|c_i) 的似然性，然后除以最大化似然，得到 (3^n )*L(x1,..,xn|c_i)。然后生成一个 U(0,1) 随机变量，如果 u 小于 (3^n)*L(x1,..,xn|c_i)，则接受 c_i 作为后验采样值，否则丢弃 c_i并重新开始。重复此过程，直到您有足够的可接受样本（100、500、1,000 或更多，具体取决于您想要的准确度）。然后只需对您感兴趣的任何 c 函数取样本平均值（新观察的可能性显然是您的应用程序的候选者）。

接受-拒绝的替代方法是使用可能性值作为权重（并且不生成 u），然后继续使用所有候选者进行加权平均，而不是接受候选者的未加权平均

其它你可能感兴趣的问题

上一篇prop.test 返回的耶茨对置信区间的连续性校正下一篇根据日期在 R 中对数据框进行子集