Jak obliczyć skumulowany rozkład normalny?

100

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

martineau
źródło

Odpowiedzi:

125

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)
Alex Reynolds
źródło
9
Jako parametry można również określić średnią (loc) i wariancję (skalę). np. d = norm (loc = 10,0, skala = 2,0); d.cdf (12.0); Szczegóły tutaj: docs.scipy.org/doc/scipy-0.14.0/reference/generated/…
Irvan
6
@Irvan, parametr skali to w rzeczywistości odchylenie standardowe, a nie wariancja.
qkhhly
2
Dlaczego Scipy nazywa je jako loci scale? Skorzystałem z, help(norm.ppf)ale co do cholery są loci scale- potrzebuję pomocy w celu uzyskania pomocy ..
javadba
2
@javadba - lokalizacja i skala to bardziej ogólne terminy w statystykach, które służą do parametryzacji szerokiego zakresu dystrybucji. W przypadku rozkładu normalnego pokrywają się one ze średnią i sd, ale nie w przypadku innych rozkładów.
Michael Ohlrogge,
1
@MichaelOhlrogge. Dzięki! Oto strona NIST-u wyjaśniająca dalsze itl.nist.gov/div898/handbook/eda/section3/eda364.htm
javadba
40

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?

WTIFS
źródło
3
To było dokładnie to, czego szukałem. Jeśli ktoś inny niż ja zastanawia się, jak można to wykorzystać do obliczenia „procentu danych mieszczących się w standardowym rozkładzie”, cóż: 1 - (1 - phi (1)) * 2 = 0,6827 („68% danych w 1 standardzie odchylenie ”)
Hannes Landeholm
1
Dla ogólnego rozkładu normalnego byłoby to def phi(x, mu, sigma): return (1 + erf((x - mu) / sigma / sqrt(2))) / 2.
Bernhard Barker
18

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))
Nieznany
źródło
3
Ponieważ biblioteka std implementuje math.erf (), nie ma potrzeby oddzielnej implementacji.
Marc
nie mogłem znaleźć odpowiedzi, skąd się biorą te liczby?
TmSmth
17

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
Xavier Guihot
źródło
15

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
Cerin
źródło
9

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 .

Salvador Dali
źródło
2

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
David Miller
źródło
-9

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

platipodium
źródło
6
Pytanie dotyczy Pythona, a nie NetLogo. Tej odpowiedzi nie powinno tu być. Nie edytuj pytania, aby zmienić jego znaczenie.
interjay
Zdaję sobie sprawę, że nie jest to preferowany sposób, ale myślę, że jest to najbardziej pomocne, ponieważ ludzie są kierowani na tę stronę przez google (obecnie ...)
platipodium