Két kevert Gaussian eloszlás paraméterbecslése — Momentumok módszerével

A mai bejegyzés egy statisztikatörténeti szempontból fontos cikkel fog foglalkozni. Mi történik, ha a megfigyelésünk nem egy populációból, hanem két populáció keverékéből származik. Ez a kérdés foglalkoztatta Karl Pearson-t a XIX. század végén. Megoldása áttörést hozott a statisztikában és elterjesztette a Momentumok módszerét.


1892-ben Walter Frank Raphael Weldon Cambridge-i biológus és felesége húsvéti vakációra utaztak Málta és Nápoly környékére. Az út során valami fura okból azzal mulattak, hogy megmérték az útjukba akadó parti tarisznyarákok (Carcinus maenas) 23 különböző tulajdonságát. Sok ilyen volt, összesen 1000 rákból álló mintát gyűjtöttek. Visszatérve Cambridge Weldon nekiállt feldolgozni az adatokat. Az derült ki, hogy 22 tulajdonsága a rákoknak Gaussian eloszlást követ, viszont volt egy (a Carapace és a testhossz aránya) ami nem. Ez így nézet ki:

Az egyértelmű volt, hogy nem Gaussian eloszlás. A bal oldali farok sokkal vastagabb. Weldon elképzelése az volt, hogy talán két normál eloszlású populáció keveredett a mintákban. Ez lehetséges, mivel több helyszínen is jártak az utazásuk során. A kérdés az volt, hogy ha ez igaz, akkor mik ezeknek az alpopulációknak a tulajdonságai. Vagyis paraméterbecslést szeretett volna végezni a biológus. Mivel 1892-ben még bőven a számítógépek megjelenése előtt volt, ez a feladat meghaladta Weldon matematikai és statisztikai tudását, szóval segítségre volt szüksége. Így jött képbe Pearson. Ő nagyon érdekesnek találta a kihívást és meglátta benne a lehetőséget, hogy egy sor kutatási területen később alkalmazható módszert dolgozzon ki.

Pearson a momentumok módszerével oldotta meg a dolgot. A megoldás jelentős áttörés volt a statisztika alkalmazásában, mivel lehetővé tette komplikált modellek paraméterbecslését viszonylag egyszerű lépésekben. Itt a hangsúly azon van, hogy a módszer eszközt adott a Waldon-hoz hasonló kutatók kezébe. Ezzel a statisztika alkalmazásának új lehetőségei nyíltak meg.

1894-ben jelent meg Pearson cikke, a Contributions to the Mathematical Theory of Evolution, amiben részletezte a megoldást. Ennek első részét fogjuk ebben a bejegyzésben átnézni.

Momentumok

Kezdjük először is azzal: mi a momentumok módszere. Ahogy korábban már volt róla szó:

Az alap ötlet lényegében a többismeretlenes egyenletrendszerek megoldásán alapul. Mint tudjuk ezekben az esetekben annyi egyenletre van szükségünk a megoldáshoz amennyi ismeretlenünk van. Csebisev lényegében átültette ezt a gondolatot a paraméterbecslés feladatára. Semmi mást nem kell csinálnunk a paraméter megbecslésére, mint annyi egyenletet előállítanunk, amennyi paramétert keresünk. Ehhez pedig egyszerűen elég hatványoznunk az elvárt értéket, ezek lesznek az úgynevezett Momentumok.

Most nézzük meg meg a Gaussian eloszlás első hat momentumát, azért ennyit, mert ennyi kell a fenti probléma megoldásához. Az első maga az átlag:

(1) M_{g_1} = E[x] = \int_{-\infty}^{\infty} x \cdot \frac{1}{\sqrt{2 \pi \sigma^2}} \cdot \exp{\frac{-(x-\mu)^2}{2 \sigma^2}} = \mu

Az integrálások elvégzése után a Momentumok az (1) nagyon hasonlóan:

M_{g_1} \mu
M_{g_2} \sigma^2+\mu^2
M_{g_3} 3 \mu \sigma^2+\mu^3
M_{g_4} 3\sigma^4+6\mu^2\sigma^2+\mu^4
M_{g_5} 15b\sigma^4+10\mu^3\sigma^2+\mu^5
M_{g_6} 15\sigma^6+45\mu^2\sigma^4+15\mu^4\sigma^2+\mu^6
1. táblázat — Gaussian eloszlás momentumai

Ezután a rövid bevezetés után nézzük a konkrét problémát.

Teszt adatok és az empirikus momentumok

Most, hogy tisztában vagyunk a momentumainkkal, készítsünk néhány mesterséges megfigyelést, amin bemutathatjuk a megoldást. Tegyük fel, hogy van két Gaussian eloszlású alpopulációnk. Az egyik átlaga 20, szórása 3. A másik paraméterei 13 és 4.5. A két alpopuláció mérete nem ugyanakkora, az első 50%-al nagyobb mint a második. Vegyünk mintát mindkét csoportból és keverjük össze a mintákat (a bejegyzésben használt teljes Python kód megtekinthető itt ):


import numpy as np

# mintavetelek szama
n = 1000
# how many come from the first population
p = 0.6

# mintak az elso populaciobol
elso_populacio = np.random.normal(20, 3, int(n*p))
# mintak a masodik populaciobol
masodik_populacio = np.random.normal(13, 4.5, int(n*(1-p)))

# osszekeverjuk a megfigyeleseket
observation = np.concatenate((elso_populacio, masodik_populacio))
indexes = [ i for i in range(len(observation))]
np.random.shuffle( indexes ) 
observation = observation[indexes]

Most, hogy meg vannak a megfigyeléseink nézzük meg, mik az empirikus Momentumaink. Jelöljük ezeket \mu'_i-vel. Ezek nem mások mint a megfigyelt Momentum értékek:


ns = [ i+1 for i in range(len(observation))]
# megfigyelt momentumok
mu1_vesszo = np.divide(np.cumsum(observation), ns)[-1]
mu2_vesszo = np.divide(np.cumsum( [ o**2 for o in observation ] ), ns)[-1]
mu3_vesszo = np.divide(np.cumsum( [ o**3 for o in observation ] ), ns)[-1]
mu4_vesszo = np.divide(np.cumsum( [ o**4 for o in observation ] ), ns)[-1]
mu5_vesszo = np.divide(np.cumsum( [ o**5 for o in observation ] ), ns)[-1]
mu6_vesszo = np.divide(np.cumsum( [ o**6 for o in observation ] ), ns)[-1]

Ami nekem ebben a mintában:

\mu'_1 17.0186351153081
\mu'_2 315.1198527459101
\mu'_3 6150.547485597151
\mu'_4 124794.48506271579
\mu'_5 2611413.5974093084
\mu'_6 56078524.60021198
2. táblázat — Empirikus momentumok

Megoldás

Ez eddig elég könnyű volt. Most jön az a rész ami miatt a cikk áttörő volt. A fentiek alapján most már világosnak kell lennie előttünk, hogy azért kézzel megoldani a (4) egyenletrendszereket igen nagy kihívás. Pearson ezt megtette, amivel lényegében egy könnyen követhető sorvezetőt adott mások kezébe.

Első lépésben normalizáljuk a megfigyelt adatokat. Ehhez használjuk a következő ábrát:

1. ábra

Ahol:

  • y’y’ — egy tetszöleges origó, esetünkben a 0 pont lesz.
  • yy — a kevert populáció átlaga
  • d — a yy távolsága a y’y’-től,

Könnyű belátni, hogy ekkor a normalizált Momentumok így alakulnak:

(2) \mu_n = \int_{-\infty}^{\infty} (x-d)^n \cdot y \text{ dx}

Ahol y egy tetszőleges eloszlás sűrűségfüggvénye.

Legyen d = h \cdot q . Ahol h egy tetszőleges hosszúság. Igazából nekünk a h-ra nem lesz szükségünk mert analitikusan fogjuk megoldani a feladatot, de Pearson a grafikus megoldás lehetősége miatt alkalmazta. Ezért, hogy a cikkel összhangban maradjak én is megteszem. A h használatával tetszőlegesen lehet nagyítani és kicsinyíteni az ábrákat. Ne felejtsük el, hogy a cikk írásának idejében a kézzel kockás papírra rajzolt ábrák korában vagyunk. Esetünkben h = 1 , így d = q .

A fentieknek megfelelően a közös populáció momentumai így fognak alakulni:

\mu_1 \mu'_1-q
\mu_2 \mu_2'-2q\mu'_1+q^2
\mu_3 \mu'_3-3q\mu'_2+3q^2\mu'_1-q^3
\mu_4 \mu'_4-4q\mu'_3+6q^2\mu'_2-4q^3\mu'_1+q^4
\mu_5 \mu'_5-5q\mu'_4+10q^2\mu'_3 - 10q^3\mu'_2+ 5q^5\mu'_1- q^5
3. táblázat — Normalizált momentumok számítása

Számítsuk akkor ki ezt a mi mintánkra:


q = mu1_vesszo
mu1 = mu1_vesszo-q
mu2 = mu2_vesszo-2*q*mu1_vesszo+q**2
mu3 = mu3_vesszo-3*q*mu2_vesszo+3*q**2*mu1_vesszo-q**3
mu4 = mu4_vesszo-4*q*mu3_vesszo+6*q**2*mu2_vesszo-4*q**3*mu1_vesszo+q**4
mu5 = mu5_vesszo-5*q*mu4_vesszo+10*q**2*mu3_vesszo-10*q**3*mu2_vesszo + 5 *q**4*mu1_vesszo-q**5

Ami eredménye:

\mu_1 0
\mu_2 25.485911557912175
\mu_3 -79.83316464535255
\mu_4 2051.761115700443
\mu_5 -15854.546063608956
4. táblázat — Normalizált Momentumok

Vezessünk be most két új jelölést: \alpha-t és c-t. Az alfa legyen a mintáink száma, míg a c az egyes alpopulációkból származó minták száma. Értelemszerűen alfa ismert míg a c-et nem ismerjük. Könnyű belátni, hogy:

(3) \alpha = c_1 + c_2

Most jön a trükk. Az alpopulációk átlagát nem az origótól mért távolságban akarjuk meghatározni, hanem a kevert populáció átlagától vett távolságban. Ehhez használhatjuk az 1. ábrát, csak annyit kell csinálnunk, hogy a y’y’ legyen most az első momentum a 0 helyett. Próbáljuk ennek megfelelően definiálni az alpopulációt:

2. ábra

Ahol a b az alpopuláció átlagának távolsága a kevert populációtól.

Most írjuk át ennek megfelelően a Gaussian eloszlás Momentumait, amit az 1. táblázatban számítottuk. Az átírásnál vegyük figyelembe a mintavételek számát is:

M_{1} b \cdot c
M_{2} (\sigma^2+b^2) \cdot c
M_{3} (3 b \sigma^2+b^3) \cdot c
M_{4} (3\sigma^4+6b^2\sigma^2+b^4) \cdot c
M_{5} (15b\sigma^4+10b^3\sigma^2+b^5) \cdot c
5. táblázat — Gaussian eloszlás momentumai

Pearson ezután bevezetett néhány új jelölést:

(4) z = \frac{c}{\alpha}

(5) u = \frac{\sigma}{b}

(6) \gamma = \frac{b}{h}

A z nem más mint, hogy mekkora részét képezi ez egyes alpopuláció a kevert megfigyelésnek. Esetünkben mivel a h az 1, a \gamma megegyezik a b-vel. Most írjuk át ennek megfelelően a fenti momentumokat:

M_{1} \gamma z \alpha h
M_{2} \gamma^2 z (1+u^2) \alpha h^2
M_{3} \gamma^3 z (1+3u^2) \alpha h^3
M_{4} \gamma^4 z (1+6u^2+3u^4) \alpha h^4
M_{5} \gamma^5 z (1+10u^2+15u^4) \alpha h^5
6. táblázat — Gaussian eloszlás momentumai Pearson átírása

Eddig ugye csak egy alpopulációval foglalkoztunk, de a probléma kulcsa, hogy a megfigyelt populáció két alpopuláció keveréke. Nem nehéz látni, hogy a megfigyelt sűrűségfüggvény a két alpopuláció sűrűségfüggvényének összege lesz. Így a normalizált kevert populáció Momentumai felírhatók a részek összegeként:

1 z_1+z_2
0\gamma_1 z_1+\gamma_2 z_2
\mu_2 \gamma_1^2 z_1 (1+u_1^2) + \gamma_2^2 z_2 (1+u_2^2)
\mu_3  \gamma_1^3 z_1 (1+3u_1^2) + \gamma_2^3 z_2 (1+3u_2^2)
\mu_4  \gamma_1^4 z_1 (1+6u_1^2+3u_1^4) + \gamma_2^4 z_2 (1+6u_2^2+3u_2^4)
\mu_5 \gamma_1^5 z_1 (1+10u_1^2+15u_1^4) + \gamma_2^5 z_2 (1+10u_2^2+15u_2^4)
7. táblázat — Kevert normalizált eloszlás Momentumai

Ezzel meg is vannak azok az egyenletek, amik segítségével Pearson megoldotta a paraméterbecslést.

A 7. táblázat első két sora segítségével ki tudjuk fejezni a z-ket:

(7) z_1 = -\frac{\gamma_2}{\gamma_1-\gamma_2}

(8) z_2 = -\frac{\gamma_1}{\gamma_1-\gamma_2}

Ezután nekiállt megoldani az egyenletrendszert. A megoldáshoz bevezetett néhány új jelölést:

(9) \lambda_4 = 9\mu_2^2 - 3 \mu_4

(10) \lambda_5 = 30\mu_2 \mu_3 - 3 \mu_5

Mint látható ebben a kettőben nincs semmi ismeretlen és a normalizált kevert eloszlásból számolható. Alkalmazásuk célja csak a számítás egyszerűsítése. Ha már itt vagyunk mi is nézzük meg az értéküket a mintánkra:


lambda4 = 9*mu2**2-3*mu4
lambda5 = 30*mu2*mu3-3*mu5

Szintén bevezetett három másik jelölést amiben már van ismeretlen:

(11) p_1 = \gamma_1+\gamma_2

(12) p_2 = \gamma_1 \gamma_2

(13) p_3 = p_1 p_2

Ezek elsőre elég légből kapott definícióknak tűnnek. Az igazság, hogy Pearson említi is, hogy sok időt vett igénybe a megoldás. Úgy gondolom, hogy anélkül, hogy mi magunk is végigjárnánk az egyenletrendszer megoldásának útját, nem nagyon fogjuk látni, hogy mi az a pont amikor ezek előjöttek. Elégedjünk meg azzal, hogy egy ponton Pearson bevezette őket, hogy jobban kordában tudja tartani a változók tengerét és ez működött. A megoldás során kiderült, hogy a p_3 felírható csak a p_2 segítségével:

(14) p3 = \frac{2\mu_3^3-2\mu_3\lambda_4 p_2 - \lambda_5 p_2^2 - 8 \mu_3 p_2^3 }{4\mu_3^2 - \lambda_4 p2 + 2 p_2^3 }

Ez meg ugye azt jelenti, hogy ha ismerjük a p_2 akkor tudjuk számolni belőle a p_1-et. Ezzel meg megvannak a \gamma-k, ezzel egyutt a b-k.

A p-k és a \gamma segítségével sikerült kifejeznie a \sigma-at:

(15) \sigma_1 = \mu_2-\frac{\mu_3}{3\gamma_2} - \frac{p_1\gamma_1}{3} + p_2

(16) \sigma_1 = \mu_2-\frac{\mu_3}{3\gamma_1} - \frac{p_1\gamma_2}{3} + p_2

Pearson végül sikerült megoldani a p_2-t, egy kilencedik fokú polinommal:

(17) p_2^9 + (-28/24\lambda_4)\cdot p_2^7 + (36/24\mu_3^2)\cdot p_2^6 + (-\mu_3\lambda_5-10/24\lambda_4^2)\cdot p_2^5 + (-(148\mu_3^2\lambda_4+2\lambda_5^2)/24 )\cdot p_2^4 + ((288\mu_3^4-12\lambda_4\lambda_5\mu_3-\lambda_4^3)/24 )\cdot p_2^3 + ( (24\mu_3^3 \lambda_5-7\mu_3^2\lambda_4^2)/24 )\cdot p_2^2 + (32\mu_3^4\lambda_4)/24)\cdot p_2 + \mu_3^6 = 0

Vegyük észre, hogy minden a fenti egyenletben csak a p_2 amit nem tudunk a mintából számolni. Vagyis ha megoldjuk, akkor megoldottuk az egyenletrendszert. A számítástechnikai probléma, hogy Abel–Ruffini bizonyítása után tudjuk, hogy negyediknél nagyobb polinomot nem lehet általánosan megoldani. De meg lehet konkrétan. Tegyük meg. Először nézzük meg a polinom számait (a szokásos jelölés szerint a_i-vel jelölve ezeket):


a2 = -28/24*lambda4
a3 = 36/24*mu3**2
a4 = -mu3*lambda5-10/24*lambda4**2
a5 = -(148*mu3**2*lambda4+2*lambda5**2)/24
a6 = (288*mu3**4-12*lambda4*lambda5*mu3-lambda4**3)/24
a7 = (24*mu3**3*lambda5-7*mu3**2*lambda4**2)/24
a8 = (32*mu3**4*lambda4)/24
a9 = -(24*mu3**6)/24

Majd oldjuk meg x-re:

from sympy import symbols, Eq, solve
x = symbols('x')
eq = Eq(x**9 + a2*x**7 + a3*x**6 + a4*x**5 + a5*x**4 + a6*x**3 + a7*x**2 + a8*x + a9, 0 )
sol = solve(eq)

Ugye ennek több megoldása lesz. Számunkra csak azok a megoldások érdekesek amik valós számok. Ezen kívül a p_2-nek vannak természetes korlátai a megfigyelésünk alapján. Mik ezek? A (12) alapján ezek a két alpolulació b értékének szorzatának minimuma és maximuma. Ugye ezek akkor vannak a legtávolabb egymástól, ha a két szélső megfigyelések. A határértékek ennek a távolságnak a mínusz és plusz előjeles értékei. Hogy miért kell plusz? Mert előfordulhat, hogy a két alpopuláció ugyanazon az oldalán van a kevert populáció átlagának. Szóval szűrjük az eredményt:


sol = [ complex(s) for s in sol  ]
# szűrő a p_2 tartományhoz 
realfilter = ( max(observation)-mu1_vesszo) * (mu1_vesszo-min(observation))
p2 = np.array([ float(s.real) for s in sol if s.imag == 0j and s.real >= -realfilter and s.real <= realfilter ])

Most már, hogy megvan a becslésünk a p_2-re, kiszámíthatjuk az összes minket érdeklő paramétert. Először a p értékeket a (13) és a (11) alapján:


p3 = np.divide( 
         np.subtract(
             np.subtract(
                 np.subtract( 
                     2*mu3**3, 
                     2*mu3*lambda4*p2
                 ), 
                 lambda5*np.power(p2, 2)), 
             8*mu3*np.power(p2, 3) 
         ), 
         np.add(
             np.subtract(
                 4*mu3**2, 
                 lambda4*p2
             ) ,
             2*np.power(p2, 3)
         ) 
     )
p1 = np.divide(p3, p2)

Majd a végiglépkedünk a megoldásokon és megbecsüljük a paramétereket a. A \gamma-át a (11) és (12) alapján felírhatjuk mint egy másodfokú polinomként:

(18) \gamma^2 + p_1 \gamma - p_2 = 0

A z értékekében a (7) és a (8) segít, a \sigma -ákat pedig a (15) és (16) segítségével számíthatjuk:


from math import isinf, isnan
models = {}
for idx in range(len(p2)):
    # gammak megoldasa a jól ismert masodfokú megoldóképlettel
    gamma = np.array( [ (p1[idx]-np.sqrt(p1[idx]**2-4*p2[idx]))/2, (p1[idx]+np.sqrt(p1[idx]**2-4*p2[idx]))/2 ] ) 
    # a gammak dobása amiknek nincs valós megoldása
    if isnan(gamma[0]) or isnan(gamma[1]):
        continue
    mu = mu1_vesszo+gamma
    z = np.array([ -gamma[1]/(gamma[0]-gamma[1]), gamma[0]/(gamma[0]-gamma[1]) ])
    if isinf(z[0]):
        continue

    sigma2 = np.array( [
        mu2 - mu3/(3*gamma[1]) - p1[idx]*gamma[0]/3 + p2[idx],
        mu2 - mu3/(3*gamma[0]) - p1[idx]*gamma[1]/3 + p2[idx]
    ])
    
    
    # 6. momentum
    check = float(np.sum(np.multiply(z, ( 15*sigma2**3+45*mu**2*sigma2**2+15*mu**4*sigma2+mu**6 )  )))
    
    models[idx] = { '6th moment loss': np.abs(mu6_vesszo-check), 'mean': mu, 'std': np.sqrt(sigma2), 'p': z }
    
print(models)

Ezzel lényegében készen is lennénk a paraméterbecslésekkel. De van itt egy probléma. Több lehetséges modellünk van. Nekem ebben az esetben kettő:

1. modell2.modell
\mu[13.71336346, 20.0742309 ][14.18106269, 20.20129223]
\sigma[4.7551137 , 2.78850582][4.89055873, 2.65094848]
z[0.46198655, 0.53801345][0.5092325, 0.4907675]
8. táblázat — Eredmények

Mint látható az első modell esetén a két alpopuláció átlaga 13.71 és 20.07, a második modell alapján pedig: 14.18 és 20.20.

Melyik a jobb? Talán feltűnt, hogy én a bejegyzés elején az első hat Momentumot számoltam ki, viszont a paraméterbecslésre csak 5 kellet. Ez nem volt véletlen. Pearson azt javasolja, hogy azt a modellt válaszuk aminek a hatodik Momentuma közelebb van az megfigyelt közös eloszlás hatodik Momentumához. Ennek számítása a fenti kódban is látható a 6th moment loss számításánál. Szóval akkor válaszuk ki a legjobb modellt:


# a legjobb modell indexe
bestidx = list(models.keys())[0]
# megkeressük a legjobb modellt
for k in models.keys():
    if models[k]['6th moment loss'] < models[bestidx]['6th moment loss']:
        bestidx = k

Ha ábrázoljuk ezután az eredményt akkor az én szimulált mintavételemre a következő:

3. ábra

Ez elég jónak néz ki. Viszont nem feltétlen igaz. Alapvetően két fajta módon állíthatunk elő normál eloszlásokból ilyen ferde eloszlásokat. Az egyik a fentebb is látható eset, amikor is összeadjuk a két alpopulációt. De lehet egy másik eset is. Igen. Amikor kivonjuk egymásból őket. Ilyen eset lehet például amikor egy populáció egyedszámát vizsgáljuk. Itt a születés és a halálozás két ellentétes irányú hatást fejt ki. Ilyenkor azt a modellt kell választanunk, ami ennek jobban megfelel. Ezeket az eseteket könnyen felismerhetjük: az alpopuláciok átlaga ugyanazon az oldalon helyezkedik el a kevert populáció átlagához képest. Ennek az eredménye így néz ki:

Azt is meg kell jegyeznünk, hogy két Gaussian eloszlás keverékeként tudunk modellezni egy adatsort az nem bizonyítja, hogy a populációnak tényleg két alpoplulációja van. Ez csak azt mondja, hogy lehet, hogy így van. De az is lehet, hogy nincsenek semmilyen alpopuláciok, csak egyszerűen Weibull eloszlást követ a populáció sűrűségfüggvénye.

Végkövetkeztetés

Pearson cikke után a Momentumok módszerének alkalmazása rohamos terjedésnek indult. Viszont nem véletlen, hogy napjainkban már nem annyira népszerű. A sokadfokú polinom miatt eléggé instabil. Napjainkban inkább az Elvárás Maximáló algoritmust vagy Markov-lánc Monte Carlo-t használnak inkább.

Irodalom

Hírdetés

Vélemény, hozzászólás?

Adatok megadása vagy bejelentkezés valamelyik ikonnal:

WordPress.com Logo

Hozzászólhat a WordPress.com felhasználói fiók használatával. Kilépés /  Módosítás )

Twitter kép

Hozzászólhat a Twitter felhasználói fiók használatával. Kilépés /  Módosítás )

Facebook kép

Hozzászólhat a Facebook felhasználói fiók használatával. Kilépés /  Módosítás )

Kapcsolódás: %s