Monte Carlo szimulációk Pythonban I. – Pi becslés és a születésnap-paradoxon5 perc olvasás

Monte Carlo szimulációnak azokat a számítógépes algoritmusokat nevezzük, amelyek valószínűségi változók ismételt generálásán alapulnak. Ennek a módszernek rengeteg felhasználási területe van, például a matematikában, fizikában és pénzügyekben. A módszer a nevét a Monacoban található szerencsejáték-központról kapta, hiszen a véletlennek és a valószínűségeknek központi szerepe van a módszerben, hasonlóan a szerencsejátékokhoz.

Ilyen szimulációk segítségével megérthetőek, vagy bebizonyíthatóak egyes elvont matematikai problémák. Továbbá, ideális gyakorlás lehet egy ilyen program írása egy kezdő programozónak, adatelemzőnek. Ebben a cikkben kettő egyszerű szimulációt fogok szemléltetni.

Az első a Pi becslése. Gimnáziumban megtanultuk, hogy a Pi értéke 3.14159…., és ezt elfogadtuk. Arra is emlékszünk talán, hogy az ókori görögök bebizonyították, hogy miért pont ennyi, de nem volt módszerünk arra, hogy mi magunk is megbizonyosudjunk róla. Ha viszont rendelkezésünkre áll egy számítógép, amin tudunk futtatni pár egyszerű Python kódot, most ezt megtehetjük.

Először képzeljünk el egy négyzetet, amelynek oldalának hossza 1 egység.  A négyzet területe tehát 1. Ezután képzeljünk el a négyzeten belül egy negyedkört, aminek sugarának hossza szintén 1 egység. Ennek a körcikknek a területe tehát .

Következő lépésként kezdjünk el véletlenszerűen pontokat dobálni a négyzetbe. Számoljuk meg, hogy az összes pontból hány esik a negyedkörön belülre. A Pi ennek az aránynak a négyszereseként becsülhető, hiszen:

Pythonban ez a következő kóddal kivitelezhető:

import numpy as np
import matplotlib.pyplot as plt

samples=2000 #A generált véletlen pontok száma

'''
A Numpy csomag segítségével generálhatunk 0 és 1 közötti számokat
egyenletes eloszlásból. Ezekre tekintsünk mint koordináta párokra.
'''
data=np.random.rand(samples,2) 

'''
A negyedkörön belül eső pontokat úgy határozhatjuk meg, hogy ezeknek
a pontoknak a két koordinátájának négyzetösszegének a négyzetgyökének
kisebbnek kell lennie, mint 1, hiszen a kör képlete x^2+y^2=1, ahol
x és y a koordináták.
'''
circle=data[np.sqrt(data[:,0]**2+data[:,1]**2)<1]

plt.figure(figsize=(5,5))
plt.scatter(data[:,0],data[:,1],color='#429fc4')
plt.scatter(circle[:,0],circle[:,1],color='#d85468')
plt.show()

Ezek után az arány kiszámításával meg is becsülhetjük a pi-t:

pi_hat=float(len(circle))/samples*4

Az így kapott eredményem 3,164. Ha növeljük a generált pontok számát (samples), azzal növelhető a becslés pontossága. 100000 pont esetén a becslés 3,14112, ami a valódi pi értékhez képest mindössze 0,0001%-kal tér el.

A második szimulációs példa a születésnap-paradoxon elfogadásához nyújt segítséget. A születésnap probléma arra a kérdésre keresi a választ, hogy egy szobában, ahol n ember tartózkodik, milyen valószínűséggel lesz két ember, aki egy napon született. Biztosan akkor lehetséges ez, hogyha a szobában 367 ember tartózkodik, hiszen a Február 29-et beleértve 366 féle születésnap lehetséges. Azonban 70 ember esetén annak a valószínűsége, hogy van két ember, aki egyazon napon született 99%, 23 ember esetén pedig 50%, azt feltételezve, hogy az egyes szülinapok előfordulási esélye azonos. A valószínűség úgy számolható ki, hogy egyből kivonjuk annak a valószínűségét, hogy mindenki másik napon született, tehát

ahol:

Ezt még így is viszonylag nehéz elhinni, ezért próbáljunk meg leszimulálni 10000 darab olyan szobát, ahol 23 ember van és nézzük meg, hogy a szobák hány százalékában fordul elő, hogy két embernak ugyanakkor van a szülinapja. A szülinapok helyett elegendő 0 és 367 közötti számokat generálni. Ehhez a random packaget használhatjuk.

 

import random

found_dupe=0 # Azoknak a szobáknak a száma, ahol van azonos szülinap
samples=10000 # Generált szobák száma
for i in range(samples):

    birthdays=[]
    n=23
    dupes=0 #Hány azonos szülinap van
    for j in range(n):
        bday=random.randrange(367)
        if bday in birthdays:
            dupes+=1
        birthdays.append(bday)

    if dupes>=1:
        found_dupe+=1

probability=found_dupe/float(samples)

Az általam kapott eredmény 0,5051, azaz tényleg az esetek 50%-ában van egyezés, ha 23 ember van a szobában. Ha le akarjuk tesztelni 70 főre is, a kódban csak az n-t kell átírni 70-re. Az így kapott eredmény nekem 0,9994 volt. Ezzel a kétségeink a paradoxon igazával kapcsolatban remélhetőleg elmúltak.

Leave a Reply

Your email address will not be published. Required fields are marked *