Mintavétel szimulálása — példa számítás

Gyakran előfordul, hogy szimulálni szeretnénk egy tetszőleges eloszlásból való véletlenszerű mintavételt. Az ismertebb eloszlásoknál ez nem probléma, általában van hozzá generátor függvény a Scipy-ban. A mai bejegyzésben azt fogjuk megnézni, mit tehetünk ha nem vagyunk ilyen szerencsések, és nekünk kell megírni a függvényt.


Vegyük azt a populációt 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)

Sajnos ez nem egy ismert eloszlásfüggvény, így nekünk kell megírni a véletlenszerű mintavételt hozzá. Rajta!

Ehhez az általánosított inverz eloszlásfüggvényre van szükségünk. Ami könnyen előállítható, ha az eloszlásfüggvényben az x-et kicseréljük, Q-ra (az általánosított inverz eloszlásfüggvény jele), és ezt az egyenlővé tesszük u-val, majd megoldjuk Q-ra. Mi az u? Egy egyenletes eloszlás 0 és 1 között.

Tehát a fentieket követve először az eloszlásfüggvényre van szükségünk, ami:

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

A későbbi személtetés kedvért nézzük ugyanezt Pythonban, a példa kedvéért a valós Théta értéke 4 lesz.1

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

Most a (2)-es ben számított eloszlásfüggvényt tegyük egyelővé u-val és cseréljük ki benne az x-et Q-ra:

(3)  1 - exp\left(-\frac{Q^2}{2 \theta^2}\right) = u

Amit, ha megoldunk Q-ra, akkor azt kapjuk, hogy:

(4)  Q = \sqrt{-2\cdot\theta^2\cdot ln(1-u) }

Ami Python kódban:

# Általánosított inverz eloszlásfüggvény
def icdf(theta, n):
    u = np.random.uniform(0, 1, n ) 
    return (np.sqrt(-2*np.power(theta, 2)*np.log(1-u)))

Most már minden rendelkezésre áll a mintavétel szimulálására. Tegyük is meg. Vegyünk 1000 mintát és hasonlítsuk össze a populáció valós sűrűségfüggvényével:

# mennyi mintát veszünk
n = 1000
# 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()

Az eredmény pedig a következő:

Ezzel készen is vagyunk, az eredmény pedig vizuálisan meggyőző.

 


Lábjegyzet

  1. Ha érdekel, hogy lehet a mintákból ezt a thétát megbecsülni, akkor olvasd el: Paraméterbecslés Maximum likelihood módszerrel
Hírdetés

Mintavétel szimulálása — példa számítás” bejegyzéshez egy 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