[SciPy] Inverse der Chi Quadrat Verteilung

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
AxXel001
User
Beiträge: 29
Registriert: Sonntag 7. Juni 2015, 22:22

Hey Leute!
Da ich statistisch nicht allzu begabt bin und auch das SciPy Framework erst seit knapp 2 Stunden kenne, denke ich es ist einfach hier eine frage zu stellen, als sich durch die ewig langen Docs zu kämpfen. Ich weiß, dass das gegen die Philosophie des Forums verstößt aber mal ehrlich, SciPy ist echt riesig und ich hab ja bereits etwas gefunden.

Mein Problem: Ich möchte Quantile der Chi Quadrat Verteilung berechnen. Da ein Quantil Q(p) angibt, wie hoch die Wahrscheinilchkeit ist, dass P(X<=p), muss ich also G(x) berechnen, wenn G die Inverse Verteilungsfunktion der Verteilung ist.
Wikipedia liefert hübsche lange Formeln dafür, die ich sicherlich nicht in Python abtippen und auf kleine Rundungsfehler und gute Performance hoffen will. Darum versuche ich, das mit SciPy zu lösen.

Scheinbar kann man diese Aufgabe mit der ppf (Percent Point Function) einer Verteilung berechnen:
http://docs.scipy.org/doc/scipy/referen ... inuous.ppf

Allerdings ist mir noch schleierhaft, wie ich die Funktion nun genau nutzen muss.

Liebe Grüße und vielen Dank für jede Hilfe,
Alex
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Das wird so nicht funktionieren.
Ich würde Dir empfehlen zunächst mit Numpy und Matplotlib anzufangen. Numpy-Arrays brauchst Du für deine Daten und Matplotlib brauchst Du um Dir Ergebnisse anschauen zu können. Danach kannst Du anfangen Statistik mit Scipy zu machen, hier ist z.B. ein guter Startpunkt http://docs.scipy.org/doc/scipy/referen ... stats.html.
Danach bist Du vielleicht immer noch nicht in der Lage Dein konkretes Problem zu lösen, aber Du kannst hier Code posten und konkret fragen warum Dein Code nicht das tut, was Du gerne hättest und auf sowas bekommst Du hier dann meistens Antworten, die Dir weiterhelfen.
a fool with a tool is still a fool, www.magben.de, YouTube
AxXel001
User
Beiträge: 29
Registriert: Sonntag 7. Juni 2015, 22:22

Danke für die Information!
Seit meinem Post heute morgen habe ich mir bereits eine Lösung zusammengehackt, ziemlich dirty aber funktioniert gut :D
Ich hab quasi die Chi Quadrat Verteilungsfunktion selbst implementiert (Fehlerfunktion, Gammafunktion, etc. liefert ja das math Modul) und nähere mich per Intervallschachtelung an die Lösung der Gleichung

F(x) = c

mit c konstant an, bestimmte also das x. Glücklicherweise muss dieses auch nicht sonderlich genau sein, 2 Nachkommastellen reichen da aus, geht aber auch für 10 Nachkommastellen beinahe genau so schnell.
Benutzeravatar
snafu
User
Beiträge: 6731
Registriert: Donnerstag 21. Februar 2008, 17:31
Wohnort: Gelsenkirchen

Es wäre nett, wenn du den Code hier in den Thread posten könntest. Dies hilft Leuten weiter, die ein ähnliches Problem haben und z.B. durch Googlen auf diesen Thread gestoßen sind. Außerdem ist die Wahrscheinlichkeit hoch, dass dein "Anfänger-Code" noch verbessert werden kann. Das hilft dann auch dir weiter.
AxXel001
User
Beiträge: 29
Registriert: Sonntag 7. Juni 2015, 22:22

Joa, gute Idee! Da habters

Code: Alles auswählen

def chi_square_dist(alpha, dof):
	"""Calculates F(x,n) where x=alpha and n=dof (degrees of freedom)"""
	if dof % 2 == 0:
		resl = 1
		s = 0
		k = 0
		while k <= 0.5*dof - 1:
			tmp = (1.0/gamma(k+1)) * pow(0.5*alpha,k)
			s += tmp
			k += 1
		resl -= exp(-0.5*alpha) * s
	else:
		resl = erf(sqrt(0.5*alpha))
		s = 0
		k = 0
		while k <= (dof//2)-1:
			tmp = (1.0/gamma(k+1.5)) * pow(0.5*alpha, k+0.5)
			s += tmp
			k += 1
		resl -= exp(-0.5*alpha) * s
	return resl

def get_quantile(resl, dof):
	"""Calculates the (1-resl)-quantile of chi square distribution
	with 'dof' degrees of freedom"""
	alpha_low = 1.0
	alpha_high = 1.0
	while chi_square_dist(alpha_high, dof) < resl:
		alpha_high *= 2
	avrg = 0.5*(alpha_low + alpha_high)
	val = chi_square_dist(avrg, dof)
	# Precision is set inside code... Ugly!
	while abs(val - resl) > 0.0001:
		if val > resl:
			alpha_high = avrg
		else:
			alpha_low = avrg
		avrg = 0.5*(alpha_low + alpha_high)
		val = chi_square_dist(avrg, dof)
	return avrg

if __name__ == "__main__":
    alpha = 0.05 # a standard value
    dof = 5
    print get_quantile(alpha, dof) # >> ~= 11.07
    # see https://de.wikipedia.org/wiki/Chi-Quadrat-Test#Tabelle_der_Quantile_der_Chi-Quadrat-Verteilung
    # for more german information
Benutzeravatar
snafu
User
Beiträge: 6731
Registriert: Donnerstag 21. Februar 2008, 17:31
Wohnort: Gelsenkirchen

Aus den beiden `while`-Schleifen in `chi_square_dist()` sollte man besser `for`-Schleifen machen:

Code: Alles auswählen

while k <= 0.5*dof - 1:
   ...
# ...wird zu:
for k in xrange(0.5 * dof):
   ...

while k <= (dof//2)-1:
   ...
# ...wird zu:
for k in xrange(dof // 2):
   ...
Das jeweilige ``k = 0`` und ``k += 1`` fällt dann natürlich weg.
Sirius3
User
Beiträge: 17710
Registriert: Sonntag 21. Oktober 2012, 17:20

@AxXel001: die beiden if-Zweige von chi_square_dist sind fast identisch, also kann man sie auch zusammenfassen:

Code: Alles auswählen

def chi_square_dist(alpha, dof):
    """Calculates F(x,n) where x=alpha and n=dof (degrees of freedom)"""
    if dof % 2 == 0:
        result = 1
        b = 0
    else:
        result = erf(sqrt(0.5 * alpha))
        b = 0.5
    s = sum((0.5 * alpha) ** (k + b) / gamma(k + 1 + b) for k in xrange(dof // 2))
    result -= exp(-0.5 * alpha) * s
    return result
Antworten