logo标签列表

Haskell中与Ruby的pnormaldist统计函数等效的是什么？

rubyhaskellstatisticshackage
8

如下所示：http://www.evanmiller.org/how-not-to-sort-by-average-rating.html

这里是Ruby代码本身，是在Statistics2库中实现的：

# inverse of normal distribution ([2])
# Pr( (-\infty, x] ) = qn -> x
def pnormaldist(qn)
  b = [1.570796288, 0.03706987906, -0.8364353589e-3,
       -0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
       -0.104527497e-5, 0.8360937017e-7, -0.3231081277e-8,
       0.3657763036e-10, 0.6936233982e-12]

  if(qn < 0.0 || 1.0 < qn)
    $stderr.printf("Error : qn <= 0 or qn >= 1  in pnorm()!\n")
    return 0.0;
  end
  qn == 0.5 and return 0.0

  w1 = qn
  qn > 0.5 and w1 = 1.0 - w1
  w3 = -Math.log(4.0 * w1 * (1.0 - w1))
  w1 = b[0]
  1.upto 10 do |i|
    w1 += b[i] * w3**i;
  end
  qn > 0.5 and return Math.sqrt(w1 * w3)
  -Math.sqrt(w1 * w3)
end
6个回答
6

这个翻译比较简单：

module PNormalDist where

pnormaldist :: (Ord a, Floating a) => a -> Either String a
pnormaldist qn
  | qn < 0 || 1 < qn = Left "Error: qn must be in [0,1]"
  | qn == 0.5        = Right 0.0
  | otherwise        = Right $
      let w3 = negate . log $ 4 * qn * (1 - qn)
          b = [ 1.570796288, 0.03706987906, -0.8364353589e-3, 
                -0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5, 
                -0.104527497e-5, 0.8360937017e-7, -0.3231081277e-8, 
                0.3657763036e-10, 0.6936233982e-12]
          w1 = sum . zipWith (*) b $ iterate (*w3) 1
      in (signum $ qn - 0.5) * sqrt (w1 * w3)

首先，让我们看看 Ruby - 它返回一个值，但有时会打印错误消息（当给出不正确的参数时）。这不是很 Haskell 风格，因此让我们的返回值为 Either String a - 当给出不正确的参数时，我们将返回一个带有错误消息的 Left String，否则返回 Right a。
现在我们来检查顶部的两种情况：
  • qn < 0 || 1 < qn = Left "Error: qn must be in [0,1]" - 这是错误条件，当 qn 超出范围时。
  • qn == 0.5 = Right 0.0 - 这是 Ruby 检查 qn == 0.5 and return * 0.0
接下来，在 Ruby 代码中定义了 w1。但是我们稍后重新定义它，这不是很符合 Ruby 的风格。我们第一次存储在 w1 中的值立即在定义 w3 时使用，所以为什么不跳过将其存储在 w1 中的步骤呢？我们甚至不需要执行 qn > 0.5 and w1 = 1.0 - w1 步骤，因为我们在定义 w3 时使用了乘积 w1 * (1.0 - w1)。
因此，我们跳过所有这些步骤，直接移动到定义 w3 = negate . log $ 4 * qn * (1 - qn)。
接下来是定义 b，这是从 Ruby 代码中直接提取的（Ruby 中数组文字的语法是 Haskell 中列表的语法）。
这是最棘手的部分 - 定义 w3 的最终值。Ruby 代码所做的是
w1 = b[0]
1.upto 10 do |i|
  w1 += b[i] * w3**i;
end

所谓的折叠是将一组值（存储在Ruby数组中）缩减为单个值的过程。我们可以更加功能化地重新表述这个过程（但仍然使用Ruby），使用Array＃reduce：

w1 = b.zip(0..10).reduce(0) do |accum, (bval,i)|
  accum + bval * w3^i
end

注意我如何使用恒等式b [0] == b [0] * w3 ^ 0 将b [0]推入循环中。

现在我们可以直接将其移植到Haskell，但这有点丑陋。

w1 = foldl 0 (\accum (bval,i) -> accum + bval * w3**i) $ zip b [0..10]

相反的，我将其分成几个步骤 - 首先，我们并不真正需要i，我们只需要w3的幂次（从w3^0 == 1开始），因此让我们使用iterate (*w3) 1计算这些幂次。
然后，我们最终只需要它们的乘积，而不是将它们与b的元素进行配对，因此我们可以使用zipWith (*) b将它们与每个对的乘积一起配对。
现在我们的折叠函数非常简单 - 我们只需要对产品求和，我们可以使用sum来完成。
最后，根据qn是否大于0.5（我们已经知道它不相等），我们决定返回加号或减号sqrt (w1 * w3)。因此，与ruby代码中在两个不同位置计算平方根不同，我只计算了一次，并根据qn - 0.5的符号（signum 仅返回值的符号）乘以+1或-1。
5
在Hackage上搜索，有许多与统计相关的库：

您需要一个pnormaldist的版本，它“返回normaldist（x）的P值”。

也许其中的某些内容提供了您所需的东西？

我真的对统计一窍不通 :P。你知道这些函数中哪一个等同于pnormaldist吗？ - Chris Bolton
我认为这些函数都不完全符合你的需求。如果我没记错的话，你需要 erf 函数的反函数。 - augustss
3

您需要的功能现在已经可以在Hackage上的erf软件包中获得。它被称为invnormcdf。

1

这是我的Wilson得分置信区间在Node.js中用于伯努利参数的应用。

wilson.normaldist = function(qn) {
    var b = [1.570796288, 0.03706987906, -0.0008364353589, -0.0002250947176, 0.000006841218299, 0.000005824238515, -0.00000104527497, 0.00000008360937017, -0.000000003231081277,
        0.00000000003657763036, 0.0000000000006936233982
    ];
    if (qn < 0.0 || 1.0 < qn) return 0;
    if (qn == 0.5) return 0;
    var w1 = qn;
    if (qn > 0.5) w1 = 1.0 - w1;
    var w3 = -Math.log(4.0 * w1 * (1.0 - w1));
    w1 = b[0];

    function loop(i) {
        w1 += b[i] * Math.pow(w3, i);
        if (i < b.length - 1) loop(++i);
    };
    loop(1);
    if (qn > 0.5) return Math.sqrt(w1 * w3);
    else return -Math.sqrt(w1 * w3);
}

wilson.rank = function(up_votes, down_votes) {
    var confidence = 0.95;
    var pos = up_votes;
    var n = up_votes + down_votes;
    if (n == 0) return 0;
    var z = this.normaldist(1 - (1 - confidence) / 2);
    var phat = 1.0 * pos / n;
    return ((phat + z * z / (2 * n) - z * Math.sqrt((phat * (1 - phat) + z * z / (4 * n)) / n)) / (1 + z * z / n)) * 10000;
}
1

Ruby代码没有文档说明；没有规定这个函数应该做什么。任何人怎么知道它是否正确地执行了预期的操作？

我不会盲目地将这个算术从一个实现中复制并粘贴到另一个实现中（就像Ruby包的作者所做的那样）。

在注释中给出了引用作为（[2]），但是这是悬空的。我们在本地C代码的注释块中找到它，在_statistics2.c文件中。

/*
   statistics2.c
   distributions of statistics2
   by Shin-ichiro HARA
   2003.09.25
   Ref:
   [1] http://www.matsusaka-u.ac.jp/~okumura/algo/
   [2] http://www5.airnet.ne.jp/tomy/cpro/sslib11.htm
*/

仅引用从哪里抄袭系数的C源代码，而不是公式的原始来源，这种工作非常草率。

[1]链接已经失效了，服务器未找到。幸运的是，我们需要的是[2]。这是一个包含各种函数的C代码的日语页面。参考文献已给出。我们需要的是pnorm。在表格中，该算法归因于“戸田の近似式”，意思是“戸田的逼近”。

戸田是日本常见的姓氏；需要更多的侦探工作来找出这个人是谁。

经过很多努力，我们终于找到了：论文（日语）：标准正态分布百分点的极小极大逼近（1993），作者为戸田秀雄和小野晴美。

该算法归因于戸田（我假设是同一位与论文合著的戸田），日期为1967年第19页。

看起来相当晦涩；在 Ruby 包中使用它的可能原因是发现了国内学术界的源代码，其中引用了国内学者的姓名。

0

在 Hackage 上简单查看并没有发现任何相关内容，所以我建议你将 Ruby 代码翻译成 Haskell。这很简单。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接