[SciPy] Inverse der Chi Quadrat Verteilung

[AxXel001]

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

[MagBen]

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.

[AxXel001]

Danke für die Information!
Seit meinem Post heute morgen habe ich mir bereits eine Lösung zusammengehackt, ziemlich dirty aber funktioniert gut
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.

[snafu]

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]

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

[snafu]

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]

@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