Verteilung einer quadratischen Form für Einheitsvektoren
Andreas Tell
ich habe ein mathematisches Problem, das ich trotz Recherche bisher
nicht lösen konnte:
Sei M die Matrix einer symmetrischen quadratischen Form im R^n mit spec
(M) \subset [0,1].
Und sei v eine Zufallsvariable gleichförmig verteilter
Einheitsvektoren im R^n.
Was ist die Verteilung von v^t M v?
Zunächst ist M o.B.d.A diagonal. Mein Ansatz war bisher, v durch einen
Vektor unabhängiger normalverteilter Zufallsvariablen zu ersetzen um
die Kugelsymmetrie zu erhalten. Die quadratische Form erzeugt daraus
dann eine verallgemeinerte Chi^2 Verteilung. Da ich allerdings nur an
den Einheitsvektoren interessiert bin, muss ich noch durch die Norm
von v teilen. Das Ergebnis ist also eine Ratio-Verteilung zweier Chi^2
Verteilungen, allerdings keine F-Verteilung, da die beiden Chi^2
Verteilungen nicht unabhängig sind.
Ich benötige zunächst nur die Klasse der resultierenden Distribution
um sie an meine Daten zu fitten. Falls aber jemand sogar die Parameter
aus den Eigenwerten von M ableiten kann, wäre ich um so glücklicher.
Kann mir jemand hier weiterhelfen oder mich auf entsprechende
Literatur verweisen? Meine Versuche, die Ratio-Verteilung selbst zu
ermitteln sind leider fehlgeschlagen.
Freundliche Grüße,
Andreas Tell
Bastian Erdnuess
> Hallo,
>
> ich habe ein mathematisches Problem, das ich trotz Recherche bisher
> nicht l�sen konnte:
>
> Sei M die Matrix einer symmetrischen quadratischen Form im R^n mit spec
> (M) \subset [0,1].
> Und sei v eine Zufallsvariable gleichf�rmig verteilter
> Einheitsvektoren im R^n.
> Was ist die Verteilung von v^t M v?
Wie sieht es denn f�r den Fall n = 2 aus? Wenn M = diag(a,b) ist, mit
0 <= a,b <= 1 und v = (cos T, sin T), dann ist
X := v^t M v = a cos^2 T + b sin^2 T. Damit X monoton steigend ist und
wegen der Symmetrie von v kann man OE a <= b und T gleichverteilt in
[0,p/2] annehmen.
F�r die Dichte fX von X gilt nun mit (FX, FT Verteilung von X, resp. T)
fX(x) = d/dx FX(x) = d/dx P(X <= x) = d/dx P(T <= X^(-1)(x))
= d/dx FT(X^(-1)(x)) = fT(X^(-1)(x)) d/dx X^(-1)(x)
= 2/pi / X'(X^(-1)(x))
da die Dichte von T fT(t) = 2/pi ist.
Es sind weiter
X' = 2 (b-a) cos T sin T ,
cos^2 T = (b-X)/(b-a) , (aus X = a cos^2 T + b sin^2 T)
sin^2 T = (X-a)/(b-a) , also
X' = 2 sqrt( (b-X)(X-a) ) , und damit
fX(x) = 1/pi / sqrt( (b-x)(x-a) ) f�r a <= x <= b .
Die Verteilung kommt mir jetzt aber leider gar nicht bekannt vor.
Evtl. kannst du die Dichte f�r den Fall n = 3 nachrechnen und erkennst
das Bildungsgesetz f�r fX, das sich dann mit etwas Gl�ck per Induktion
beweisen l�sst.
Bastian
Andreas Tell
> Wie sieht es denn f r den Fall n = 2 aus?
...
> fX(x) = 1/pi / sqrt( (b-x)(x-a) ) f r a <= x <= b .
> Evtl. kannst du die Dichte f r den Fall n = 3 nachrechnen und erkennst
> das Bildungsgesetz f r fX, das sich dann mit etwas Gl ck per Induktion
> beweisen l sst.
Vielen Dank Bastian für Deine Ausführliche Herleitung. Leider war n=2
auch der einzige Fall, für den ich eine Lösung gefunden habe und dann
den direkten Ansatz verworfen habe. Für n=3 wird das ganze schon so
schwierig, vor allem wegen der Parametrisierung der Gleichverteilung
auf S^2 und der fehlenden Umkehrbarkeit von X, dass ich nicht mehr bis
zum Ergebnis vordringen konnte.
Daher bin ich auf den oben beschriebenen Ansatz gekommen. Meine
Versuche, auch ohne Kenntnis von n>2 eine Verallgemeinerung der
Distribution zu finden die zu meinen Daten passt erwiesen sich auch
als Sackgasse.
Vielleicht fallen Dir ja noch mögliche Formen für n > 2 ein? Jede Idee
wäre Hilfreich.
Danke,
Andreas
Bastian Erdnuess
> On 13 Dez., 20:29, earth...@web.de (Bastian Erdnuess) wrote:
>
> > Wie sieht es denn f r den Fall n = 2 aus?
> ...
> > fX(x) = 1/pi / sqrt( (b-x)(x-a) ) f r a <= x <= b .
>
> > Evtl. kannst du die Dichte f r den Fall n = 3 nachrechnen und erkennst
> > das Bildungsgesetz f r fX, das sich dann mit etwas Gl ck per Induktion
> > beweisen l sst.
>
> Vielen Dank Bastian f�r Deine Ausf�hrliche Herleitung. Leider war n=2
> auch der einzige Fall, f�r den ich eine L�sung gefunden habe und dann
> den direkten Ansatz verworfen habe. F�r n=3 wird das ganze schon so
> schwierig, vor allem wegen der Parametrisierung der Gleichverteilung
> auf S^2 und der fehlenden Umkehrbarkeit von X, dass ich nicht mehr bis
> zum Ergebnis vordringen konnte.
>
> Daher bin ich auf den oben beschriebenen Ansatz gekommen. Meine
> Versuche, auch ohne Kenntnis von n>2 eine Verallgemeinerung der
> Distribution zu finden die zu meinen Daten passt erwiesen sich auch
> als Sackgasse.
>
> Vielleicht fallen Dir ja noch m�gliche Formen f�r n > 2 ein? Jede Idee
> w�re Hilfreich.
Wenn man S^n mit
vn = ( cos x1 cos x2 ... cos xn,
sin x1 cos x2 ... cos xn,
sin x2 ... cos xn,
...,
sin x[n-1] cos xn,
sin xn ) e IR^(n+1)
in (x1, ..., xn) parametrisiert und Mn = diag(p0, p1, ..., pn) setzt mit
0 <= p0 <= p1 <= ... <= pn <= 1, findet man leicht die Rekursion
Xn := vn^t Mn vn = X[n-1] cos^2 xn + pn sin^2 xn , X0 := p0
f�r die gesuchten ZVs Xn.
Um auf S^n die Gleichverteilung hinzukriegen, kann man die xn selber als
ZVs auffassen mit Dichte f[xn](t) = cn cos^(n-1) t f�r n >= 1, mit einer
geeigneten Normierungskonstante cn.
Dazu �berlegt man sich am besten, dass man sich S^n "von Norden nach
S�ben"-Zerlegt als eine Sammlung von lauter skalierten S^(n-1)
vorstellen kann (vgl. Breitenkreise auf S^2). Beim Abstand sin xn vom
"�quator" hat die eingespannte n-1 Sph�re einen Radius von cos xn und
entsprechend eine Oberfl�che proportional zu cos^(n-1) xn.
Mit entsprechender Gewichtung darf dann ein Zuf�llspunkt auch nur in
diese n-1 Sph�re fallen, damit sich insgesamt eine Gleichverteilung
ergeben kann.
Mit g[n](x) := Xn|X[n-1]|xn=x = X[n-1] cos^2 x + pn sin^2 x (lies "Xn
gegeben X[n-1] an der Stelle xn=x", g[n](x) ist dabei selbs wieder eine
ZV "vom Typ von X[n-1]") findet man analog zum 2D-Fall (wenn man sich
wieder auf 0 <= x <= pi/2 einsch�nkt)
f[Xn](x) = d/dx F[Xn](x) = d/dx P(Xn <= x)
= d/dx P(X[n-1] cos^2 xn + pn sin^2 xn <= x)
= d/dx E[P(X[n-1] cos^2 xn + pn sin^2 xn <= x | X[n-1])]
= d/dx E[P(g[n](xn) <= x)] = d/dx E[P(xn <= g[n]^-1(x))]
= d/dx E[ F[xn](g[n]^-1(x)) ] = E[ d/dx F[xn](g[n]^-1(x)) ]
= E[ f[xn](g[n]^-1(x)) / g[n]'(g[n]^-1(x)) ]
= cn/2 E[sqrt( ((pn-x)/(pn-X[n-1]))^(n-1) / ((pn-x)(x-X[n-1])) )]
= cn/2 (pn-x)^(n/2-1) E[1/sqrt( (pn-X[n-1])^(n-1) (x-X[n-1]) )]
da aus g[n](x) = X[n-1] cos^2 x + pn sin^2 x folgen
cos^2 x = (pn - g[n](x)) / (pn - X[n-1])
sin^2 x = (g[n](x) - X[n-1]) / (pn - X[n-1])
g[n]'(x) = 2 (pn - X[n-1]) cos x sin x
= 2 sqrt( (pn - g[n](x)) (g[n](x) - X[n-1]) )
g[n]'(g[n]^-1(x)) = 2 sqrt( (pn - x) (x - X[n-1]) )
cos g[n]^-1(x) = sqrt( (pn - x) / (pn - X[n-1]) )
f[xn](g[n]^-1(x)) = cn sqrt( (pn - x) / (pn - X[n-1]) )^(n-1) .
Mit cn = 1 / int[0..pi/2] cos^(n-1) t dt = ((n-1)/2 �ber 1/2) ergibt
sich also die Funktional-Rekursion f�r die Dichten fn := f[Xn] von Xn
fn(x) = ((n-1)/2 �ber 1/2) / 2 (pn - x)^(n/2-1)
* int[0..pi/2] f[n-1](t) / sqrt( (pn - t)^(n-1) (x - t) ) dt
f0(x) = delta(p0 - x) , (delta = Dirac-Distribution) .
Zum Test f�r n = 1 (nanu, das war urspr�nglich mal n = 2) ergibt sich
f1(x) = (0 �ber 1/2) / 2 (p1 - x)^(-1/2) / sqrt( (x-p0) (p1-p0)^0 )
= 1/pi / sqrt( (p1 - x) (x - p0) ) . Passt!
F�r n = 2 kann man dann auch noch ziemlich direkt eine Formel mit einem
elliptischen Integral hinschreiben. Die kann man vielleicht mal mit
simulierten Werten vergleichen.
Aufpassen solltest du evtl. noch mit den Integrationsgrenzen, um die hab
ich mir n�mlich �berhaupt keine Gedanken gemacht; und auch bei anderen
Stellen, wie beim Vertauschen von Erwartungswert und Differentation hab
ich mir keine Gedanken um die Vorraussetzungen f�r Konvergenz gemacht.
Bastian
Andreas Tell
> fn(x) = ((n-1)/2 über 1/2) / 2 (pn - x)^(n/2-1)
> * int[0..pi/2] f[n-1](t) / sqrt( (pn - t)^(n-1) (x - t) ) dt
> f0(x) = delta(p0 - x) , (delta = Dirac-Distribution) .
>
> Zum Test für n = 1 (nanu, das war ursprünglich mal n = 2) ergibt sich
>
> f1(x) = (0 über 1/2) / 2 (p1 - x)^(-1/2) / sqrt( (x-p0) (p1-p0)^0 )
> = 1/pi / sqrt( (p1 - x) (x - p0) ) . Passt!
>
Vielen Dank Bastian! Du hast mir sehr geholfen. Man sieht, dass Du
schon einiges an Erfahrung mitbringst. Ich bin leider zu diesem
Problem gekommen wie die Jungfrau zum Kind und hatte davor nie viel
mit Statistik zu tun. Aus Deiner Lösung habe ich viel gelernt.
Liebe Grüße,
Andreas
Bastian Erdnuess
> On 14 Dez., 06:12, earth...@web.de (Bastian Erdnuess) wrote:
>
> > fn(x) = ((n-1)/2 �ber 1/2) / 2 (pn - x)^(n/2-1)
> > * int[0..pi/2] f[n-1](t) / sqrt( (pn - t)^(n-1) (x - t) ) dt
> > f0(x) = delta(p0 - x) , (delta = Dirac-Distribution) .
> >
> > Zum Test f�r n = 1 (nanu, das war urspr�nglich mal n = 2) ergibt sich
> >
> > f1(x) = (0 �ber 1/2) / 2 (p1 - x)^(-1/2) / sqrt( (x-p0) (p1-p0)^0 )
> > = 1/pi / sqrt( (p1 - x) (x - p0) ) . Passt!
> >
>
>
> Vielen Dank Bastian! Du hast mir sehr geholfen. Man sieht, dass Du
> schon einiges an Erfahrung mitbringst. Ich bin leider zu diesem
> Problem gekommen wie die Jungfrau zum Kind und hatte davor nie viel
> mit Statistik zu tun. Aus Deiner L�sung habe ich viel gelernt.
Vielen Dank f�r die Blumen.
Eine Anmerkung vielleicht noch. Du musst evtl. mit dem Fall p[n-1] = pn
ein bisschen aufpassen, da dann g[n] nicht invertierbar ist, falls
X[n-1] = pn ist. I. d. R. hat das Ereignis 0-Wk't und spielt daher
m�glicherweise keine gro�e Rolle, aber wenn am Anfang p0 = p1 = ... = pn
< p[n+1] ist, ergibt sich mit Xn = pn = p0 ein andere Rekursionsanfang.
Viele Gr��e,
Bastian
Bastian Erdnuess
Ich bin da heute nochmal etwas sorgf�ltiger dr�ber gegangen, und habe
auch auf die Grenzen der Integrale geachtet. Dabei habe einen Fehler in
meiner letzten Rechnung entdeckt.
> f[Xn](x) = d/dx F[Xn](x) = d/dx P(Xn <= x)
> = d/dx P(X[n-1] cos^2 xn + pn sin^2 xn <= x)
> = d/dx E[P(X[n-1] cos^2 xn + pn sin^2 xn <= x | X[n-1])]
F�r P(...) = 0 f�r x < X[n-1], daher ist hier
> = d/dx E[P(g[n](xn) <= x)] = d/dx E[P(xn <= g[n]^-1(x))]
das Integral vom Erwartungswert oben bei x abgeschnitten und bei
> = d/dx E[ F[xn](g[n]^-1(x)) ] = E[ d/dx F[xn](g[n]^-1(x)) ]
kann das Differenzial nicht so einfach in den Erwartungswert gezogen
werden. Stattdessen kommt noch ein Summand
F[xn](g[n]^-1(x))
dazu, falls x < p[n-1] ist, und der st�rt im Folgenden erst mal
gewaltig.
Sorry.
Bastian