Tương đương với chức năng thống kê pnormaldist của Ruby trong Haskell là gì?

Question

Như đã thấy ở đây: http://www.evanmiller.org/how-not-to-sort-by-average-rating.html

Dưới đây là các mã Ruby bản thân, thực hiện trong thư viện 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

này là khá dễ dàng để dịch:

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) 

Trước hết, chúng ta hãy nhìn vào ruby ​​- nó trả về một giá trị, nhưng đôi khi nó in một thông báo lỗi (khi đưa ra một lập luận không đúng cách). Đây không phải là rất kiêu ngạo, do đó, chúng ta hãy có giá trị trả về của chúng tôi là Either String a - nơi chúng tôi sẽ trả về một Left String với thông báo lỗi nếu đưa ra một đối số không đúng và Right a nếu không.

Bây giờ chúng tôi kiểm tra hai trường hợp ở đầu trang:

• qn < 0 || 1 < qn = Left "Error: qn must be in [0,1]" - đây là điều kiện lỗi, khi qn là ra khỏi phạm vi.
• qn == 0.5 = Right 0.0 - đây là việc kiểm tra ruby ​​qn == 0.5 and return * 0.0

Tiếp theo, chúng ta định nghĩa w1 trong mã ruby. Nhưng chúng tôi xác định lại nó một vài dòng sau đó, mà không phải là rất rubyish. Giá trị mà chúng tôi lưu trữ trong w1 lần đầu tiên được sử dụng ngay lập tức theo định nghĩa của w3, vậy tại sao chúng ta không bỏ qua lưu trữ nó trong w1? Chúng tôi thậm chí không cần thực hiện bước qn > 0.5 and w1 = 1.0 - w1, vì chúng tôi sử dụng sản phẩm w1 * (1.0 - w1) trong định nghĩa của w3.

Vì vậy, chúng tôi bỏ qua tất cả điều đó và chuyển thẳng đến định nghĩa w3 = negate . log $ 4 * qn * (1 - qn).

Tiếp theo là định nghĩa của b, là thang máy thẳng từ mã ruby ​​(cú pháp của ruby ​​cho mảng chữ là cú pháp của haskell cho danh sách).

Đây là bit phức tạp nhất - xác định giá trị cuối cùng của w3. Mã ruby ​​làm gì trong

w1 = b[0] 
1.upto 10 do |i| 
    w1 += b[i] * w3**i; 
end 

Được gọi là nếp gấp - giảm một tập hợp các giá trị (được lưu trữ trong một mảng ruby) thành một giá trị duy nhất. Chúng ta có thể xác định lại chức năng này hơn (nhưng vẫn trong ruby) sử dụng Array#reduce:

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

Lưu ý làm thế nào tôi đẩy b[0] vào vòng lặp, bằng cách sử dụng danh tính b[0] == b[0] * w3^0.

Bây giờ chúng ta có thể cổng này trực tiếp đến Haskell, nhưng đó là một chút xấu xí

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

Thay vào đó, tôi chia ra thành nhiều bước - trước hết, chúng ta không thực sự cần i, chúng ta chỉ cần quyền hạn của w3 (bắt đầu từ w3^0 == 1), do đó, hãy tính toán những người có iterate (*w3) 1.

Sau đó, thay vì nén chúng thành từng cặp với các phần tử của b, cuối cùng chúng ta chỉ cần các sản phẩm của chúng, vì vậy chúng tôi có thể nén chúng thành các sản phẩm của mỗi cặp sử dụng zipWith (*) b.

Bây giờ chức năng gấp của chúng tôi thực sự dễ dàng - chúng tôi chỉ cần tổng hợp các sản phẩm mà chúng tôi có thể thực hiện bằng cách sử dụng sum.

Cuối cùng, chúng tôi quyết định có trả lại dấu cộng hay trừ sqrt (w1 * w3), theo số liệu qn lớn hơn hoặc nhỏ hơn 0,5 (chúng tôi đã biết không bằng). Vì vậy, thay vì tính căn bậc hai ở hai vị trí riêng biệt như trong mã ruby, Tôi đã tính nó một lần và nhân nó theo +1 hoặc -1 theo ký hiệu qn - 0.5 (signumjust returns the sign of a value).

Answer 2

Đào xung quanh trên Hackage, có một số thư viện cho các thống kê:

• hmatrix-gsl-stats - liên kết thuần túy với GSL
• hstatistics - một giao diện mức thậm chí cao hơn để GSL
• hstats - phương pháp thống kê chung
• statistics - phương pháp thống kê phổ biến hơn
• statistics-linreg - một hồi quy tuyến tính giữa hai mẫu, dựa trên gói thống kê khác.

Bạn muốn có phiên bản pnormaldist, "Trả về giá trị P của normaldist (x)".

• Statistics.Distribution.Normal, từ gói thống kê, cung cấp nhiều chức năng để thao tác phân phối bình thường.
• Statistics.Test.NonParametric chứa một số điều cần làm với P-values.

Có lẽ thứ gì đó cung cấp những gì bạn cần?

Answer 3

Một cái nhìn ngắn gọn về hackage không tiết lộ bất cứ điều gì, vì vậy tôi khuyên bạn nên dịch mã ruby ​​thành Haskell. Nó đủ đơn giản.

Answer 4

Chức năng bạn muốn hiện có sẵn trong gói phần mềm erf về hackage. Nó được gọi là invnormcdf.

Answer 5

đây là khoảng tin cậy điểm Wilson của tôi cho một tham số Bernoulli trong 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

Mã Ruby là không có giấy tờ; không có đặc điểm kỹ thuật nào về chức năng này được yêu cầu làm. Làm thế nào để bất cứ ai biết liệu nó có đúng không?

Tôi sẽ không chỉ sao chép và dán một cách mù quáng số học này từ lần triển khai này sang lần thực hiện khác (như tác giả của gói Ruby đã làm).

Trích dẫn được đưa ra là ([2]) trong nhận xét, nhưng điều này đang lúng túng. Chúng tôi tìm thấy nó trong khối bình luận của mã C bản địa trong tập tin _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 
*/ 

Rất cẩu thả để chỉ trích dẫn mã nguồn C từ nơi hệ số đã được nôi, thay vì nguồn gốc của công thức.

Liên kết [1] không hoạt động nữa; không tìm thấy máy chủ. May mắn thay, cái chúng tôi muốn là [2]. Đây là một trang bằng tiếng Nhật với một số mã C cho các chức năng khác nhau. Tài liệu tham khảo được đưa ra. Người chúng tôi muốn là pnorm. Trong bảng, thuật toán được gán cho 戸 田 の 近似 式 có nghĩa là "Toda's Approximation".

Toda là họ phổ biến ở Nhật Bản; công việc thám tử nhiều hơn là cần thiết để tìm hiểu xem đây là ai.

Sau nhiều nỗ lực, ở đây chúng tôi đi: giấy (Nhật Bản): The Minimax Approximation for Percentage Points of the Standard Normal Distribution (1993) bởi Hideo Toda và Harumi Ono.

Thuật toán là do Toda (tôi giả sử một trong những tương tự mà là đồng tác giả của bài báo), năm 1967 trên P. 19.

Có vẻ như khá mơ hồ; lý do cơ bản để sử dụng nó trong gói Ruby là nó được tìm thấy trong mã nguồn của nguồn gốc trong nước trích dẫn tên của một học giả trong nước.