Haskell中Ruby的pnormaldist統計函數的等價物是什麼？

Question

正如在這裏看到：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 

Answer 1

這是非常簡單的翻譯：

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) 

首先，讓我們來看看紅寶石 - 它返回一個值，但有時它會打印錯誤消息（給出不正確的參數時）。這不是非常危險，所以 讓我們返回值爲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 - 這是紅寶石檢查qn == 0.5 and return * 0.0

接下來，我們定義在Ruby代碼w1。但是我們在幾行後重新定義了它，這不是很紅寶石。我們在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 

是所謂的摺疊什麼 - 減少一組值（存儲在陣列紅寶石）爲單精度值。我們可以重申這一點使用Array#reduce更多功能（但仍紅寶石）：

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

注意我怎麼推b[0]進入死循環，使用標識b[0] == b[0] * w3^0。

現在我們可以口這直接哈斯克爾，但它是一個有點難看

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)。因此，不是像計算紅寶石代碼那樣在兩個分開的位置計算平方根，我計算了一次，並且根據qn - 0.5（signumjust returns the sign of a value）的符號將其乘以+1或-1。

Answer 2

上Hackage周圍挖，有統計的一些庫：

• hmatrix-gsl-stats - 與GSL的純結合
• hstatistics - 一個更高層次的接口GSL
• hstats - 共同統計方法
• statistics - 更常見的統計方法
• statistics-linreg - 兩個樣本之間的線性迴歸的基礎上，其他統計信息包。

你想要一個版本的pnormaldist，其中「返回normaldist（x）的P值」。

• Statistics.Distribution.Normal從統計數據包中提供了許多操作正態分佈的函數。
• Statistics.Test.NonParametric包含一些與P-values有關的事情。

也許這裏提供了你所需要的東西？

Answer 3

對hackage的簡要介紹沒有透露任何內容，所以我建議您將ruby代碼轉換爲Haskell。這很簡單。

Answer 4

你想要的功能現在可用於hafage上的erf包。它被稱爲invnormcdf。

Answer 5

這裏是在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; 
} 

Answer 6

的Ruby代碼是無證我威爾遜的得分置信區間;沒有規定這個功能應該做什麼。任何人如何知道它是否正確無誤地做了什麼？

我不會盲目地將這個算法從一個實現複製並粘貼到另一個實現（比如Ruby包的作者）。

在評論中引用爲([2])，但這是懸而未決的。我們在_statistics2.c文件的原始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。在表中，該算法歸因於戸田の近似式，意思是「戶田的近似」。

戶田是日本的常見姓氏，需要更多的偵探工作來找出這是誰。

經過一番努力，在這裏我們去：紙（日本）：The Minimax Approximation for Percentage Points of the Standard Normal Distribution (1993)由英朗戶和晴美小野。

的算法歸因於戶田（我假設同樣一個是該論文的合着者），可追溯至1967年19 P.

它似乎相當模糊的;在Ruby軟件包中使用它的可能的基本原理是它在國內的源代碼中被發現，引用國內學者的名字。