¿Cuál es el equivalente de la función de estadísticas pnormaldist de Ruby en Haskell?

statistics hackage

Question

Como se ve aquí: http://www.evanmiller.org/how-not-to-sort-by-average-rating.html

Aquí está el código de Ruby, implementado en la biblioteca de 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

Buscando en Hackage, hay una serie de bibliotecas para estadísticas:

• hmatrix-gsl-stats - una unión pura a GSL
• hstatistics - una interfaz de nivel aún más alto para GSL
• Hstats - métodos estadísticos comunes
• estadísticas - métodos estadísticos más comunes
• statistics-linreg - una regresión lineal entre dos muestras, basada en el otro paquete de estadísticas.

Desea una versión de pnormaldist , que "Devuelve el valor P de normaldist (x)".

• Statistics.Distribution.Normal , del paquete de estadísticas, proporciona muchas funciones para manipular distribuciones normales.
• Statistics.Test.NonParametric contiene varias cosas que hacer con P-values .

Tal vez algo allí proporciona lo que necesita?

Answer 2

Una breve mirada al hackage no reveló nada, así que sugiero que traduzcas el código Ruby a Haskell. Es bastante simple.

Answer 3

Esto es bastante sencillo de traducir:

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)

Primero, veamos el rubí: devuelve un valor, pero a veces imprime un mensaje de error (cuando se le da un argumento impropio). Esto no es muy atractivo, así que tengamos que nuestro valor de retorno sea Either String a - donde le devolveremos Left String con un mensaje de error si se le da un argumento incorrecto, y a Right a contrario.

Ahora comprobamos los dos casos en la parte superior:

• qn < 0 || 1 < qn = Left "Error: qn must be in [0,1]" qn < 0 || 1 < qn = Left "Error: qn must be in [0,1]" - esta es la condición de error, cuando qn está fuera de rango.
• qn == 0.5 = Right 0.0 - este es el cheque de rubí qn == 0.5 and return * 0.0

A continuación, definimos w1 en el código ruby. Pero lo redefinimos unas líneas más tarde, que no es muy rubí. El valor que almacenamos en w1 la primera vez se usa inmediatamente en la definición de w3 , entonces ¿por qué no salteamos el almacenamiento en w1 ? Ni siquiera necesitamos hacer el paso qn > 0.5 and w1 = 1.0 - w1 , porque usamos el producto w1 * (1.0 - w1) en la definición de w3.

Así que omitimos todo eso, y pasamos directamente a la definición w3 = negate . log $ 4 * qn * (1 - qn) w3 = negate . log $ 4 * qn * (1 - qn) .

La siguiente es la definición de b , que es un levantamiento recto del código ruby ​​(la sintaxis de ruby ​​para una matriz literal es la sintaxis de Haskell para una lista).

Aquí está el truco: definir el valor máximo de w3 . Lo que hace el código Ruby en

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

Es lo que se denomina pliegue, lo que reduce un conjunto de valores (almacenados en una matriz de ruby) en un solo valor. Podemos replantear esto de manera más funcional (pero aún en rubí) usando Array#reduce :

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

Observe cómo presioné b[0] en el ciclo, usando la identidad b[0] == b[0] * w3^0 .

Ahora podríamos portar esto directamente a Haskell, pero es un poco feo

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

En cambio, lo dividí en varios pasos: primero, realmente no necesitamos i , solo necesitamos los poderes de w3 (comenzando en w3^0 == 1 ), así que calculemos aquellos con iterate (*w3) 1 .

Entonces, en lugar de comprimirlos en pares con los elementos de b, finalmente solo necesitamos sus productos, de modo que podemos comprimirlos en los productos de cada par usando zipWith (*) b .

Ahora nuestra función de doblado es realmente fácil, solo tenemos que resumir los productos, lo que podemos hacer usando la sum .

Por último, decidimos si devolver más o menos sqrt (w1 * w3) , de acuerdo con si qn es mayor o menor que 0.5 (ya sabemos que no es igual). Entonces, en lugar de calcular la raíz cuadrada en dos ubicaciones separadas como en el código ruby, lo calculé una vez y lo multipliqué por +1 o -1 acuerdo con el signo de qn - 0.5 ( signum simplemente devuelve el signo de un valor ).

Answer 4

La función que desea ahora está disponible en el paquete erf en hackage. Se llama invnormcdf .

Answer 5

aquí está el intervalo de confianza de puntuación de mi Wilson para un parámetro de Bernoulli en 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

El código de Ruby no está documentado; no hay especificaciones de lo que se supone que debe hacer esta función. ¿Cómo sabe alguien si hace correctamente lo que se pretende?

No solo copiaría y pegaría a ciegas esta aritmética de una implementación a otra (como lo hizo el autor del paquete Ruby).

Una cita se da como ([2]) en un comentario, pero esto está colgando. Lo encontramos en el bloque de comentarios del código C nativo en el archivo _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
*/

Un trabajo muy descuidado para citar solo el código fuente C desde donde se cifraron los coeficientes, en lugar de la fuente original de la fórmula.

El enlace [1] ya no funciona; Servidor no encontrado. Afortunadamente, el que queremos es [2] . Esta es una página en japonés con un código C para varias funciones. Las referencias son dadas El que queremos es pnorm . En la tabla, el algoritmo se atribuye a 戸 田 近似 式 式, que significa "Aproximación de Toda".

Toda es un apellido común en Japón; Se requiere más trabajo de detective para averiguar quién es.

Después de mucho esfuerzo, aquí vamos: papel (japonés): la aproximación de Minimax para puntos porcentuales de la Distribución Normal Estándar (1993) de Hideo Toda y Harumi Ono.

El algoritmo se atribuye a Toda (estoy asumiendo el mismo que es el coautor del artículo), fechado en 1967 en la página 19.

Parece bastante oscuro; la razón probable para usarlo en el paquete de Ruby es que se encontró en el código fuente de origen nacional citando el nombre de un académico nacional.