Paraméterbecslés Maximum likelihood módszerrel — példa számítás

Statisztikai elemzéséknél gyakran kell paraméterbecslést végeznünk. Ez a következő probléma: van egy ismert típusú populációnk, viszont nem ismerjük azt a paramétert, ami a konkrét populációra jellemző. Ebben a posztban megpróbálom bemutatni, hogy határozhatjuk meg ezt a nem ismert paramétert a maximum likelihood módszerrel.1 Most érdekes elméleti kérdések helyett csak egy számítási sorvezető következik. Azért írtam meg ezt a bejegyzést, mert nem nagyon találtam olyan magyar leírást, ami tisztán csak a számítás kivitelezésére koncentrál. Úgy vélem mindenki példákon keresztül tanul a legjobban, megpróbálok bemutatni egy számítási példát, aminek a lépéseit könnyű átültetni bármilyen egyparaméterű eloszlás becslésére.


A fenti leírás alapján talán nem egészen tiszta: mi az a paraméter, szóval nézünk néhány példát. Első példa a klasszikus érmefeldobásos: van egy érménk és feldobjuk n alkalommal és megfigyeljük, hány alkalommal lett fej az eredmény. Mi a paraméter? Ugye a kísérlet leírásából tudjuk, hogy egy Bernoulli-eloszlás, tehát ez a populáció típusa. Mitől függ ennél a típusnál a konkrét eloszlása az eredménynek? Ez a paraméter, amit ebben az esetben p-nek, neveznek és lényegében annak a valószínűleg, hogy fej (ha azt jelöljük 1-es értékkel) az eredmény. Nézzünk egy másik példát: ez esetben Normális eloszlást követ a populációnk. Mitől függ ennél a típusnál a konkrét eloszlása az eredménynek? Mint az tudjuk: az átlagtól és a szórástól. Vagyis itt két paraméter van. És mi azzal a populációval számolunk, amit a következő sűrűségfüggvény ír le:

(1)  f_{\theta }(x) = \frac{x}{\theta^2} \cdot exp\left(-\frac{x^2}{2 \theta^2} \right)

Itt egyetlen paraméter van a théta (θ). Ezzel a példával fogom bemutatni a maximum likelihood számítást. Azért választottam ezt az eloszlást, mert nem nevezetes típus. Így mindent magunknak kell kiszámolnunk, a sűrűségfüggvényből kiindulva. Ez pedig rákényszerít minket arra, hogy általános formában lássuk a maximum likelihood számítást. Annyi extra információnk van még, hogy x és a théta egyaránt pozitív.

A populáció általános jellemzése

Most, hogy tudjuk a populációnk milyen típusba tartozik, nézzünk meg néhány alap- tulajdonságát ennek az eloszlásnak. Azt, hogy mi az átlag, variancia és az eloszlásfüggvény.

Átlag

Ez egy folytonos eloszlás, szóval az átlagát könnyen számíthatjuk a sűrűségfüggvény x-el való szorzatának integráljából :

(2)  E[X] = \int_0^{\infty}\! x \cdot \frac{x}{\theta^2} \cdot exp\left(-\frac{x^2}{2 \theta^2} \right)    \,  \mathrm{d}x

Ha az integrálást elvégezzük akkor az eredmény:

(3)  E[X] = \sqrt{\frac{\pi}{2}} \cdot\theta

Variancia

Most nézzünk mi a variancia. Ehhez pedig a jól ismert (4) formulát használjuk.

(4)  Var[X] =E[X^2]-(E[X])^2

Persze ehhez ki kell számolnunk mi a E[X²]-et. De ez nem probléma, megtehetjük az átlag számításához hasonló módon:

(5)  E[X^2] = \int_0^{\infty}\! x^2 \cdot \frac{x}{\theta^2} \cdot  exp\left(-\frac{x^2}{2 \theta^2} \right)    \,  \mathrm{d}x = 2\cdot\theta^2

Amit ha visszahelyettesítünk a (4)-be akkor:

(6)  Var[X] =  2\cdot\theta^2 - ( \sqrt{\frac{\pi}{2}} \cdot\theta )^2 = (2- \frac{\pi}{2} ) \cdot \theta^2

Eloszlásfüggvény

Ez pedig ugye:

(7)  F_{\theta}(x) = \int_0^{x} \frac{x}{\theta^2} exp\left(-\frac{x^2}{2 \theta^2} \right) \text{ dx} = 1 - exp\left(-\frac{x^2}{2 \theta^2}\right)

Minta vétel

Ideje, hogy készítsünk néhány teszt adatot. Bár magához a Maximum likelihood számításhoz nem kell, de jól fog jönni, hogy a végén ellenőrizzük az eredményt. Ehhez Mintavétel szimulálása posztban részletezett megvalósítást fogjuk használni. A példa kedvéért a valós Théta értéke 4 lesz. Ez az amit nem ismerünk és szeretnénk majd megbecsülni a mintavételből:

import numpy as np
import matplotlib.pyplot as plt

# a sűrűségfüggvény 
def pdf(theta, x):
    return( x/ ((np.power(theta,2)) * np.power(np.e, np.power(x,2)/(2*np.power(theta,2))) ) )
 
# az eloszlásfüggvény
def cdf(theta,x):
    return 1-np.power(np.e, np.multiply(np.divide(np.power(x, 2), 2*np.power(theta, 2)), -1))

# valós théta amit nem ismerünk
tcs = 4
# átlag
atlag = np.sqrt( np.divide(np.pi, 2) )*tcs
# variancia
var = (2-np.divide(np.pi, 2) )*np.power(tcs, 2)
# medián
median = tcs*np.sqrt(2*np.log(2))

# sűrűségfüggvény
bins_tavolsag = 10
bins = [ x/bins_tavolsag for x in range(15*bins_tavolsag) ]
px = pdf(tcs, bins)
# sűrűségfüggvény ábra
plt.clf()
plt.plot(bins, px, label="Sűrűségfüggvény" )
plt.axvline(atlag, color="red", label="Átlag: "+str(np.round(atlag,2)))
plt.xlabel("x")
plt.ylabel("p(x)")
plt.legend()
plt.show()

# eloszlásfüggvény
cx = cdf(tcs, bins)
# eloszlásfüggvény ábra
plt.clf()
plt.plot(bins, cx, color="green", label="Eloszlásfüggvény")
plt.axvline(atlag, color="red", label="Átlag: "+str(np.round(atlag,2)))
plt.axvline(median, label="Median: "+str(np.round(median,2)))
plt.xlabel("x")
plt.ylabel("Φ(x)")
plt.legend()
plt.show()

Tehát a sűrűségfüggvény és az eloszlásfüggvény, ha a paraméterünk 4:

Mint a fenti ábrákból is látható, az a típusú eloszlás nem szimmetrikus az átlagra, ezt mondjuk a (3)-ből is gondolhattuk.

Készítsünk egy 100 elemű mintát, és a rend kedvért nézzük meg, hogy illeszkedik a populáció tényleges sűrűségfüggvényre:

# mennyi mintát veszünk
n = 100
# a minták
minta = icdf(tcs, n)

# sűrűségfüggvény pontjai
bins_tavolsag = 10
bins = [ x/bins_tavolsag for x in range(15*bins_tavolsag) ]
px = pdf(tcs, bins)

# mintavétel összehasonlítása a pdf-el
plt.hist(minta, bins=bins, density=True, label="Minta, n="+str(n))
plt.plot(bins, px , label="Sűrűségfüggvény" )
plt.axvline(atlag, color="red", label="Minta átlag: "+str(np.round(np.mean(minta),2)))
plt.xlabel("x")
plt.ylabel("p(x)")
plt.legend()
plt.show()

A továbbiakban majd ebből a mintavételből fogjuk megbecsülni a thétát.

Maximum log likelihood becslés

És végre elérkeztünk a Maximum loglikelihood számításhoz. A likelihood nem más mint a sűrűségfüggvény n-szeres szorzata:

(8)  L(x_i, x_2 \dots x_n; \theta) =  \prod_{i=1}^{n} \frac{x_i}{\theta^2} \cdot exp\left(-\frac{x_i^2}{2 \theta^2} \right)

Az egyszerű likelihood helyett viszont ennek log-át használjuk, abból ez egyszerű okból, hogy a logaritmus sok hasznos tulajdonsággal rendelkezik. Pl: a szorzatot összeadássá lehet alakítani.

(9)  l(x_i, x_2 \dots x_n; \theta) =  ln\left(\prod_{i=1}^{n}  \frac{x_i}{\theta^2} \cdot exp\left(-\frac{x_i^2}{2 \theta^2}  \right)\right)

Most egyszerűsítsünk egy kicsit:

(10)  l(x_i, x_2 \dots x_n; \theta) =  -2\cdot n \cdot ln(\theta) + ln(\prod_{i=1}^{n} x_i) -\frac{\sum_{i=1}^{n}x_i^2}{2 \theta^2}

Ebből dobhatjuk a középső elemet, mivel nem függ theta-tól, így:

(11)  l(x_i, x_2 \dots x_n; \theta) =  -2\cdot n \cdot  ln(\theta) -\frac{\sum_{i=1}^{n}x_i^2}{2  \theta^2}

Az általunk keresett théta lesz az, ami ezt a loglikelihood függvényt maximalizálja. Tehát egy maximum-minimum problémával nézünk szembe. Ami pedig deriválást jelent. Szóval deriváljunk thétára (mivel azt szeretnénk maximalizálni):

(12)  \frac{d}{d\theta} \left(-2\cdot n \cdot ln(\theta) -\frac{\sum_{i=1}^{n}x_i^2}{2 \theta^2} \right) = \frac{\sum_{i=1}^{n}x_i^2 - 2\cdot n \cdot \theta^2}{\theta^3}

Majd állítsuk 0-ra és oldjuk meg thétára:

(13)  \frac{\sum_{i=1}^{n}x_i^2 - 2\cdot n \cdot \theta^2}{\theta^3} = 0

(14)  \hat{\theta} = \sqrt{ \frac{\sum_{i=1}^{n}x_i^2}{2n} }

Lényegében ennyi. Ez a maximum likelihood alapján a legjobb becslés a théta értékére. Mit is mond ez? Azt, hogy a legvalószínűbb théta a mintavételünk négyzetének átlagának felének a gyöke (Jézus segíts). Akkor nézzük is meg:

# mennyi mintát veszünk
n = 100
# a minták
minta = icdf(tcs, n)

def mle(x):
    return np.sqrt( np.sum(np.power(x, 2))/(2*len(x)) )

th_mle = mle(minta) 
print(th_mle)

Ami nekem a fenti mintában: 4.370462486244649. Vagyis az alapján a 100 minta alapján amit vettünk 4.3-ra becsüljük a thétát. Ez jó? Attól függ honnan nézzük. Ugye nem a valós 4, de közelebb áll hozzá mintha mondjuk ha 13-at becsülnénk. Itt álljunk meg egy szóra. Maga a paraméter becslés csak egy számot dob nekünk a folyamat végén. De azt nem mondja meg, hogy menyire vagyunk biztosak ebben a számban. Eleve folytonos eloszlásnál problematikus egy szám valószínűsége, mivel a konkrét számok valószínűsége 0. Még a valós értéké is. Tehát semmiféleképpen ne álljunk meg itt. Mi se fogunk, hanem megnézzük mennyire vagyunk biztosak ebben az értékben, ehhez pedig a Fisher információt fogjuk használni.

Fisher információ

A Fisher információ — példa számítás bejegyzés alapján ugye azt tudjuk, hogy a Fisher információn erre az eloszlásra:

(15)  F = \frac{4}{\theta^2 }

Ennek az információnak a reciproka pedig megadja nekünk a Maximum likelihood asymptotic eloszlás varianciáját. Vagyis:

(16)  Var(\theta) = \frac {\theta^2 }{4}

Most már számíthatjuk a Confidence interval-t a becslésünkre:

(17)  I = \hat{\theta} \pm \frac{  q(\alpha)  \cdot  \sqrt{  \frac {   \theta^2 }{4} }  }{\sqrt{n}}

Ugye ezt nem tudjuk egy az egyben használni, mivel a nem ismert valós théta értéktől függ, de ezt a thétát helyettesíthetjük a becsült théta-val, szóval:

(18)  I_{plugin}  = \hat{\theta} \pm \frac{  q(\alpha)  \cdot  \sqrt{  \frac {   \hat{\theta}^2   }{4} }  }{\sqrt{n}}

Ami Pythonban, 95% valószínűségre:

# variancia
var = np.power(th_mle,2)/(4)
q95 = 1.96
# távolság a középértéktől
diff = ( q95*np.sqrt(var)/np.sqrt(n) )

# határértékek
A = th_mle - diff
B = th_mle + diff

Ez pedig a fenti mintára a következő határértékeket adta:

(19)  I_{plugin}  = [3.94, 4.8]

A varianciára pedig 0.0477-et kaptam. Mit is jelent ez. Azt, hogy ha végtelen sokszor megismételném a fentiekhez hasonló mintavételt minden egyes alkalommal 100 darab mintával, akkor az valós théta az esetek 95%-ban az így számított határértékek között lenne. Persze ez nem garancia arra, hogy ebben a mintában is ott van (mi ugye a szimuláció miatt tudjuk, hogy ez igaz, de csak azért), de azért a határértékek közötti távolság támpont számunkra.

Végezetül csak a rend kedvéért, ha nem is végtelen esetben de m alkalommal ismételjük meg a mintavételt és nézzük meg az így kapott eloszlás varianciája nagyságrendileg hasonlít-e a Fisher információból számoltra:

# az egyes becslések
def th_becsles(n):    
    minta = icdf(tcs, n)
    return mle(minta)

# ismétlések száma
m = 10000
# mennyi minta egy mintavételnél
n = 100
# becsült MLE
result = []
for x in range(m):
    result.append(th_becsles(n))
print(np.mean(result), np.var(result))

Ami 0.0401740677671672-et adott nekem a varianciára. Ez nagyságrendileg megegyezik a Fisher információban számolttal, tehát boldog vagyok. És ahogy a centrális határeloszlás-tétel alapján elvárható a mintáink által számolt Maximum Likelihood egy szép normál függvénynek néz ki:

Remélem világosan látszanak az általános számítási lépések. Nem nagyon tértem ki az elméleti kérdésekre, mit miért csinálunk, de a “példa számítások” célja inkább, a gyakorlati megvalósítás szemléltetése.


Lábjegyzet

  1. Több módszer is van, a maximum likelihood az egyik legelterjedtebb de nem az egyetlen.
Hírdetés

Paraméterbecslés Maximum likelihood módszerrel — példa számítás” bejegyzéshez 7 hozzászólá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