NetLogo 计算 beta 分布/功能答案

【问题标题】：NetLogo computing beta distribution/functionNetLogo 计算 beta 分布/功能
【发布时间】：2015-09-18 17:35:42
【问题描述】：

我正在使用 NetLogo 中的 stats 扩展来计算 beta 函数：

set z (stats:beta (H + 1) (T + 1))

其中 H 和 T 是抛硬币时正面和反面的数量。 stats扩展的使用不是必须的，我在使用beta函数的阶乘表达式时也遇到了同样的问题。

问题是当H + T > 168，NetLogo 报告z = 0 并且有一些程序我不能执行[特别是beta 分发]

有什么方法可以近似 Netlogo 中的 beta 函数（或分布），使其不会遇到这个问题？

【问题讨论】：

如果您在这里没有得到任何乐趣，您可以尝试在统计论坛中提问（以稍微不同的方式关注近似方面而不是 NetLogo）。他们可能更适合提供有关替代计算方法的建议。
谢谢 JenB，如果我没有得到回复，我会在几天后这样做。实际上，我认为核心问题是 NetLogo 无法计算大于零和小于 -aprox- 2 E-51 的数字。它只报告这些值为零。
请注意，NetLogo 使用 IEEE 754 双精度浮点，与大多数编程语言遵循的浮点标准相同，所以我认为您的问题并不是 NetLogo 特有的。跨度>
P.S.我不认为 2E-51 是正确的界限。 IEEE 754 降至 4.9E-324。
我认为问题在于计算大阶乘。 β 函数 B(H,T) 的分母是阶乘 (H+T-1)。但是如果 H+T-1 > 170，那么我们有一些 > E325。事实上，当我尝试计算阶乘（171）时，Netlogo 报告“Infinity”。

标签： netlogo

【解决方案1】：

基于 Chris 的回答，stats 扩展 确实具有 Gamma 函数的日志，stats:loggamma。它处理一个大大超过 171 的参数，因此 a 和 b 也可以增加到超过 172。

observer> show stats:beta 85 86
observer: 1.2864854397253604E-52
observer> show exp (stats:loggamma 85 + stats:loggamma 86 - stats:loggamma (85 + 86))
observer: 1.2864854397251408E-52
observer> show stats:beta 86 86
observer: 0
observer> show exp (stats:loggamma 86 + stats:loggamma 86 - stats:loggamma (86 + 86))
observer: 6.394810665301235E-53
observer> show exp (stats:loggamma 200 + stats:loggamma 200 - stats:loggamma (200 + 200))
observer: 9.713217247613997E-122

新版本的统计扩展 (v1.4.0) 已发布，其中包含使用日志的“bigBeta”功能。

【讨论】：

【解决方案2】：

我终于设法解决了我自己的问题。我希望这对其他人也有帮助。主要问题是计算大阶乘 - 因为它出现在 cmets 中。特别是171！实在是太大了，无法计算。因此，最佳解决方案是尽量减少计算这些数字的需要。

对于离散情况，beta 函数如下：

B(alpha,beta) = (alpha - 1)!(beta - 1)! / (alpha + beta - 1)！

错误是尝试先计算分母，然后是分母，然后是除法；因为分母变得太大。我找到了两种解决方案，其中一种比另一种更好，但要复杂得多——但实际上我认为它接近最优。

次优但更短是基于使用最大提名者简化分母的想法：

to-report betafunc [alpha beta]
  ifelse (alpha >= beta)
  [let v1 (factorial (beta - 1) (0))
   let v2 (factorial (alpha + beta - 1) (alpha - 1))
   report (v1 / v2)]
  [let v1 (factorial (alpha - 1) (0))
   let v2 (factorial (alpha + beta - 1) (beta - 1))
  report (v1 / v2)]
end

在哪里

to-report factorial [n m]
 if (n = m) [report 1]
 report (n * factorial (n - 1) (m))
end

这个问题是分母仍然增长非常快。第二种解决方案的想法是通过简化阶乘的每个步骤来进行。

例如，如果我们有 (5!8!)/(12)!我们可以将其表示为 (5.8/12).(4.7/11).(3.6/10) ... (1.4/8).(1.3/7)...(1.1/5)(1.1/4)...(1.1/ 1)

通过逐步简化序列的每个项，我们确保事情不会增长得那么快。这是我的解决方案：

to-report stepwisefactorial3 [n1 n2 d]
  if (n1 = 0) [report stepwisefactorial2 (n2) (d)]
  if (n2 = 0) [report stepwisefactorial2 (n1) (d)]
  if (d = 0)  [report (stepwisefactorial1 (n1) * stepwisefactorial1 (n2))]
  report ((((n1 * n2) / d)) * stepwisefactorial3 (n1 - 1) (n2 - 1) (d - 1))
end

to-report stepwisefactorial2 [n d]
  if (n = 0) [report (1 / (stepwisefactorial1 (d)))]
  if (d = 0) [report stepwisefactorial1 (n)]
  report ((n / d)* stepwisefactorial2 (n - 1) (d - 1))
end

to-report stepwisefactorial1 [d]
  if d = 0 [ report 1 ]
  report d * stepwisefactorial1 (d - 1)
end

如果 alpha + beta > 345，事情仍然会爆炸；但这是进步。

【讨论】：

不错的功能。请参阅sl8r000.github.io/ab_testing_statistics//… 进行讨论。从vosesoftware.com/ModelRiskHelp/index.htm#Distributions/… 开始，对于您关心的大多数情况，正常的近似值是可以的。您可以测试 A、B 是否 > 50，然后测试 B 附近的 A 的正态近似值，如果它们相差太大，则使用阶乘简化方法。请注意，不同的 A、B 是反正有更多的抵消，所以更小简化。
酷。感谢您将此留给未来的搜索者。

【解决方案3】：

（这是一个旁注，但问题已经回答了，无论如何。）

我想知道如何使用日志。

Beta(x,y) = Gamma(x) * Gamma(y) / Gamma(x + y)
= exp ( ln(Gamma(x)) + ln(Gamma(y)) - ln(Gamma(x + y)) )

稍微玩一下就会发现 Stats 扩展程序足够聪明，可以用这种方式计算 Beta，而它真正归结为 Gamma 中的溢出。

observer> show stats:gamma 171
observer: 7.257415615308E306
observer> show stats:gamma 172
observer: Infinity

(我在这里抄袭 http://www.stata.com/support/faqs/statistics/calculating-beta-function/.)

【讨论】：