Jak obliczyć skumulowany rozkład normalny?

Szukam funkcji w Numpy lub Scipy (lub dowolnej rygorystycznej bibliotece Pythona), która da mi funkcję skumulowanego rozkładu normalnego w Pythonie.

Oto przykład:

>>> from scipy.stats import norm
>>> norm.cdf(1.96)
0.9750021048517795
>>> norm.cdf(-1.96)
0.024997895148220435

Innymi słowy, około 95% standardowego przedziału normalnego mieści się w granicach dwóch odchyleń standardowych, wyśrodkowanych na standardowej średniej zerowej.

Jeśli potrzebujesz odwrotnego CDF:

>>> norm.ppf(norm.cdf(1.96))
array(1.9599999999999991)

Może być za późno, aby odpowiedzieć na to pytanie, ale ponieważ Google wciąż prowadzi ludzi, postanowiłem napisać tutaj swoje rozwiązanie.

Oznacza to, że od Pythona 2.7 mathbiblioteka ma zintegrowaną funkcję błędumath.erf(x)

erf()Funkcja może być wykorzystywana do obliczania tradycyjne funkcje statystyczne takie jak skumulowanego rozkładu normalnego:

from math import *
def phi(x):
    #'Cumulative distribution function for the standard normal distribution'
    return (1.0 + erf(x / sqrt(2.0))) / 2.0

Odniesienie:

https://docs.python.org/2/library/math.html

https://docs.python.org/3/library/math.html

W jaki sposób są powiązane funkcja błędu i standardowa funkcja rozkładu normalnego?

Na podstawie: http://mail.python.org/pipermail/python-list/2000-June/039873.html

from math import *
def erfcc(x):
    """Complementary error function."""
    z = abs(x)
    t = 1. / (1. + 0.5*z)
    r = t * exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
        t*(.09678418+t*(-.18628806+t*(.27886807+
        t*(-1.13520398+t*(1.48851587+t*(-.82215223+
        t*.17087277)))))))))
    if (x >= 0.):
        return r
    else:
        return 2. - r

def ncdf(x):
    return 1. - 0.5*erfcc(x/(2**0.5))

Zaczynając Python 3.8, biblioteka standardowa udostępnia NormalDistobiekt jako część statisticsmodułu.

Można go użyć do uzyskania funkcji rozkładu skumulowanego ( cdf- prawdopodobieństwo, że próbka losowa X będzie mniejsza lub równa x) dla danej średniej ( mu) i odchylenia standardowego ( sigma):

from statistics import NormalDist

NormalDist(mu=0, sigma=1).cdf(1.96)
# 0.9750021048517796

Które można uprościć w przypadku standardowego rozkładu normalnego ( mu = 0i sigma = 1):

NormalDist().cdf(1.96)
# 0.9750021048517796

NormalDist().cdf(-1.96)
# 0.024997895148220428

Aby zbudować na przykładzie Unknown, odpowiednikiem Pythona funkcji normdist () zaimplementowanej w wielu bibliotekach byłoby:

def normcdf(x, mu, sigma):
    t = x-mu;
    y = 0.5*erfcc(-t/(sigma*sqrt(2.0)));
    if y>1.0:
        y = 1.0;
    return y

def normpdf(x, mu, sigma):
    u = (x-mu)/abs(sigma)
    y = (1/(sqrt(2*pi)*abs(sigma)))*exp(-u*u/2)
    return y

def normdist(x, mu, sigma, f):
    if f:
        y = normcdf(x,mu,sigma)
    else:
        y = normpdf(x,mu,sigma)
    return y

Odpowiedź Alexa pokazuje rozwiązanie dla standardowego rozkładu normalnego (średnia = 0, odchylenie standardowe = 1). Jeśli masz rozkład normalny z meani std(co jest sqr(var)) i chcesz obliczyć:

from scipy.stats import norm

# cdf(x < val)
print norm.cdf(val, m, s)

# cdf(x > val)
print 1 - norm.cdf(val, m, s)

# cdf(v1 < x < v2)
print norm.cdf(v2, m, s) - norm.cdf(v1, m, s)

Przeczytaj więcej o cdf tutaj i scipy implementacji normalnej dystrybucji z wieloma formułami tutaj .

Zrobione z góry:

from scipy.stats import norm
>>> norm.cdf(1.96)
0.9750021048517795
>>> norm.cdf(-1.96)
0.024997895148220435

W przypadku testu dwustronnego:

Import numpy as np
z = 1.96
p_value = 2 * norm.cdf(-np.abs(z))
0.04999579029644087

Tak proste:

import math
def my_cdf(x):
    return 0.5*(1+math.erf(x/math.sqrt(2)))

Znalazłem wzór na tej stronie https://www.danielsoper.com/statcalc/formulas.aspx?id=55

Ponieważ Google podaje tę odpowiedź dla wyszukiwania netlogo pdf , oto wersja netlogo powyższego kodu Pythona

    ;; Skumulowana funkcja gęstości rozkładu normalnego
    zgłosić normcdf [x mu sigma]
        niech tx - mu
        let y 0.5 * erfcc [- t / (sigma * sqrt 2.0)]
        if (y> 1.0) [set y 1.0]
        zgłoś y
    koniec

    ;; Funkcja gęstości prawdopodobieństwa rozkładu normalnego
    zgłosić normpdf [x mu sigma]
        niech u = (x - mu) / abs sigma
        niech y = 1 / (sqrt [2 * pi] * abs sigma) * exp (- u * u / 2.0)
        zgłoś y
    koniec

    ;; Uzupełniająca funkcja błędu
    zgłosić erfcc [x]
        niech z abs x
        niech t 1,0 / (1,0 + 0,5 * z)
        niech rt * exp (- z * z -1,26551223 + t * (1,00002368 + t * (0,37409196 +
            t * (0,09678418 + t * (-0,18628806 + t * (.27886807 +
            t * (-1,13520398 + t * (1,48851587 + t * (-0,82215223 +
            t * .17087277)))))))))
        ifelse (x> = 0) [raport r] [raport 2.0 - r]
    koniec