Что эквивалентно функции ненормативной статистики Руби в Haskell? - PullRequest
Новая
       20

Что эквивалентно функции ненормативной статистики Руби в Haskell?

7 голосов
/ 25 мая 2011

Как видно здесь: 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 ]

5 голосов
/ 25 мая 2011

Это довольно просто перевести: 

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

Далее мы определяем w1 в коде рубина.Но мы переопределим это несколькими строками позже, что не очень рубиново.Значение, которое мы храним в w1 в первый раз, сразу используется в определении w3, так почему бы не пропустить его сохранение в w1?Нам даже не нужно делать шаг qn > 0.5 and w1 = 1.0 - w1, потому что мы используем произведение w1 * (1.0 - w1) в определении w3.

Поэтому мы пропускаем все это и переходим прямо к определению 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) в одно значение.Мы можем перефразировать это более функционально (но все еще в рубине), используя 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.

Теперь мы можем перенести это напрямую на 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.

Наконец, мы решаем, возвращать ли плюс или минус sqrt (w1 * w3), в зависимости от того, больше или меньше qn (мы уже знаем, что оно не равно).Поэтому вместо вычисления квадратного корня в двух разных местах, как в коде ruby, я рассчитал его один раз и умножил на +1 или -1 в соответствии со знаком qn - 0.5 (signum просто возвращаетзнак значения ).

5 голосов
/ 25 мая 2011

Покопавшись в Hackage, есть несколько библиотек для статистики:

  • hmatrix-gsl-stats - чистая привязка к GSL
  • hstatistics - интерфейс еще более высокого уровня с GSL
  • hstats - общие статистические методы
  • статистика - чаще всегостатистические методы
  • statistics-linreg - линейная регрессия между двумя выборками, основанная на другом пакете статистики.

Требуется версия pnormaldist, который "Возвращает P-значение normaldist (x)".

Возможно, что-то там обеспечивает то, что вам нужно

3 голосов
/ 26 мая 2011

Требуемая функция теперь доступна в пакете erf на hackage.Это называется invnormcdf.

1 голос
/ 27 августа 2014

вот доверительный интервал моего Вильсона для параметра Бернулли в 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;
}
0 голосов
/ 30 августа 2017

код 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 голосов
/ 25 мая 2011

Краткий обзор хакерства ничего не показал, поэтому я предлагаю вам перевести код рубина на Haskell.Это достаточно просто.

...