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.¹ 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

Több módszer is van, a maximum likelihood az egyik legelterjedtebb de nem az egyetlen.