2011-05-24 14 views
7

come si vede qui: http://www.evanmiller.org/how-not-to-sort-by-average-rating.htmlQual è l'equivalente della funzione di statistica dei pormalisti di Ruby in Haskell?

Ecco il codice Ruby stessa, implementato nella biblioteca 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 

risposta

5

questo è abbastanza semplice tradurre:

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) 

Prima di tutto, diamo un'occhiata al rubino - restituisce un valore, ma a volte un messaggio d'errore (quando dato un argomento improprio). Questo non è molto haskellish, quindi il nostro valore di ritorno è Either String a - dove restituiremo un Left String con un messaggio di errore se viene fornito un argomento errato e un Right a altrimenti.

Ora controlliamo i due casi in alto:

  • qn < 0 || 1 < qn = Left "Error: qn must be in [0,1]" - questa è la condizione di errore, quando qn è fuori portata.
  • qn == 0.5 = Right 0.0 - questo è il controllo rubino qn == 0.5 and return * 0.0

Next up, definiamo w1 nel codice Ruby. Ma lo ridefiniamo poche righe dopo, che non è molto rubino. Il valore che memorizziamo nella per la prima volta è immediatamente utilizzato nella definizione di w3, quindi perché non ignorare la memorizzazione in w1? Non abbiamo nemmeno bisogno di fare il passo qn > 0.5 and w1 = 1.0 - w1, perché usiamo il prodotto w1 * (1.0 - w1) nella definizione di w3.

Quindi saltiamo tutto ciò e passiamo direttamente alla definizione w3 = negate . log $ 4 * qn * (1 - qn).

Il prossimo è la definizione di b, che è un passaggio diretto dal codice rubino (la sintassi di ruby ​​per un array letterale è la sintassi di haskell per un elenco).

Ecco il bit più difficile: definire il valore finale di w3. Quello che il codice Ruby fa in

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

E 'quello che viene chiamato una piega - la riduzione di un insieme di valori (memorizzati in una matrice rubino) in un singolo valore. Siamo in grado di riformulare questa più funzionale (ma ancora in Ruby) utilizzando Array#reduce:

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

Nota come ho spinto b[0] nel circuito, utilizzando l'identità b[0] == b[0] * w3^0.

Ora potremmo porta questo direttamente a Haskell, ma è un po 'brutto

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

Invece, ho rotto in su in diverse fasi - prima di tutto, non abbiamo davvero bisogno i, abbiamo solo bisogno di Potenze di w3 (a partire da w3^0 == 1), quindi calcoliamo quelli con iterate (*w3) 1.

Quindi, piuttosto che zippare quelli in coppie con gli elementi b, in definitiva abbiamo solo bisogno loro prodotti, in modo da poter comprimili in i prodotti di ciascuna coppia con zipWith (*) b.

Ora la nostra funzione di piegatura è davvero semplice: abbiamo solo bisogno di riassumere i prodotti, cosa che possiamo fare usando sum.

Infine, decidere se tornare più o meno sqrt (w1 * w3), a seconda che qn è maggiore o minore di 0,5 (che già sappiamo che non è uguale). Quindi, piuttosto che calcolare la radice quadrata in due posizioni separate come nel codice rubino, l'ho calcolato una volta e lo ho moltiplicato per +1 o -1 in base al segno di qn - 0.5 (signumjust returns the sign of a value).

5

Scavando in giro su Hackage, c'è un certo numero di librerie per le statistiche:

Si desidera una versione di pnormaldist, che "restituisce il valore P di normaldist (x)".

Forse qualcosa fornisce quello che ti serve?

+0

Non so davvero nulla delle statistiche: P. Sai quale di quelle funzioni è equivalente a un anonimo? –

+0

Non credo che nessuna di queste funzioni sia esattamente ciò di cui hai bisogno. Hai bisogno dell'inverso della funzione erf, se non sbaglio. – augustss

0

Una breve occhiata a hackage non ha rivelato nulla, quindi ti suggerisco di tradurre il codice rubino in Haskell. È abbastanza semplice

3

La funzione desiderata è ora disponibile nel pacchetto erf in hackage. Si chiama invnormcdf.

1

ecco il mio intervallo di confidenza del punteggio di Wilson per un parametro di Bernoulli in 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

Il codice Ruby è documentata; non esiste una specifica di ciò che questa funzione dovrebbe fare. Come fa qualcuno a sapere se fa tutto ciò che è giusto?

Non vorrei copiare e incollare ciecamente questa aritmetica da un'implementazione a un'altra (come ha fatto l'autore del pacchetto Ruby).

Una citazione è data come ([2]) in un commento, ma questo è penzolante. Lo troviamo nel blocco dei commenti del codice C nativo nel file _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 
*/ 

lavoro molto scrupoloso per citare solo il codice sorgente C da dove si cribbed i coefficienti, piuttosto che la fonte originale della formula.

Il collegamento [1] non funziona più; server non trovato. Fortunatamente, quello che vogliamo è [2]. Questa è una pagina in giapponese con qualche codice C per varie funzioni. I riferimenti sono dati. Quello che vogliamo è pnorm. Nella tabella, l'algoritmo è attribuito a 戸 田 の 近似 式 che significa "Approssimazione di Toda".

Toda è un cognome comune in Giappone; più lavoro investigativo è necessario per scoprire chi è questo.

Dopo tanto sforzo, eccoci: carta (giapponese): The Minimax Approximation for Percentage Points of the Standard Normal Distribution (1993) di Hideo Toda e Harumi Ono.

L'algoritmo è attribuita a Toda (sto assumendo la stessa che è co-autore della carta), datato 1967, P. 19.

Sembra piuttosto oscuro; la probabile ragione per usarlo nel pacchetto Ruby è che è stato trovato nel codice sorgente di origine nazionale che cita il nome di un accademico nazionale.

Problemi correlati