Jūs esate čia

Apie Monte Carlo metodą

Ekonometriniuose modeliuose neišvengiamai susiduriame su atsitiktinumu (pvz., akivaizdu, kad tiksliai prognozuoti būsimo ketvirčio Lietuvos bendrajį vidaus produktą (BVP) neįmanoma). Mokykloje moksleiviai yra šiek tiek supažindinami su tikimybių teorijos pagrindais, o tiksliau, su vadinamuoju klasikiniu tikimybės apibrėžimu: jei visi tikimybinio eksperimento rezultatai arba elementarieji įvykiai iš aibės $\{w_1, w_2, ... w_N\}$  turi vienodus šansus, yra vienodai galimi arba tam tikra prasme simetriški, tai atsitiktinio įvykio $A$  tikimybė $P(A)=N(A)/N$ ; čia $N(A)$  yra įvykį $A$  sudarančių elementarių įvykių skaičius.

1 pavyzdys.    Nutarta, kad komitetą sudarys 5 asmenys, kuriuos reikia išrinkti iš 10 moterų ir 5 vyrų.
(a)    Kokia tikimybė, kad komitetą sudarys 3 moterys ir 2 vyrai?
(b)    Kokia tikimybė, kad komitetą sudarys tik moterys?

Sprendimas.

(a) Tradicinis sprendimas. Šis eksperimentas (atsitiktinai atrenkame 5 elementus iš aibės $\{v_1,...,v_5,m_1,...,m_{10})\}$) gali baigtis bet kuriuo iš  $N={15\choose 5 } = \frac{15!}{5!(15-5)!}=3003$

elementariųjų įvykių (pvz., tokiu:  $(m_1,m_5,m_8,v_2,v_4)$- kadangi tvarka nėra svarbi, komitetą užrašome abėcėlės ir indeksų didėjimo tvarka).  Tris moteris iš dešimties galime atrinkti ${10\choose 3}$ būdais (tai vadinamųjų derinių iš 10 po 3 skaičius), o likusius du vyrus - $5\choose 2$ būdais, taigi  $P(a)=\binom{10}{3}\binom{5}{2}/\binom{15}{5}=400/1001\approx 0,4$.

(a) Sprendimas su R.  Aukščiau pateikti skaičiavimai nėra sudėtingi, tačiau kartais tiek tikimybių teorijoje tiek ekonometrijoje susiduriame su labai komplikuotomis užduotimis. Pademonstruosime kaip šiuos skaičiavimus galima atlikti su ekonometristų  mėgiama programa R (tai galinga ir nemokama programa, kurią galima atsisiųsti internetu (su Googlu paieškokite cran)). Ją instaliavę, nukopijuokite žemiau pateiktą kodą į R‘o konsolę.

kand.m=paste("m",1:10)    # paste = suklijuoti
kand.m
kand.v=paste("v",1:5)
kand.m
kand=c(kand.m,kand.v)     # c = apjungti
kand                      # kandidatų sąrašas
choose(15,5)              # derinių iš 15 po 5 skaičius
head(t(combn(kand,5)),20) # pirmieji 20 derinių
dim(t(combn(kand,5)))     # t  reiškia  „transponuok“ (sukeisk eilutes ir
                          # stulpelius vietomis); dim = dimension
Pa=choose(10,3)*choose(5,2)/choose(15,5)
Pa                        # = 0.3996004
400/1001                  # toks pat rezultatas: 0.3996004
 
Šioje programėlėje naudojome derinių skaičių skaičiuojančią funkciją choose(n,k) (ji dar vadinama binominių koeficientų funkcija, nes pagal Newtono binomo formulę $(a+b)^n=\sum_{\substack{1 \le k \le n}} \binom{n}{k} {a^{k}}{b^{n-k}}=\sum_{\substack{1 \le k \le n}}$ choose(n,k) ${a^{k}}{b^{n-k}}$. Su R'u nesunku išbrėžti ir jos grafiką:
 
n=15
choose(n,0:n)         # k kinta nuo 0 iki n (=15)
#  1   15  105  455 1365 3003 5005 6435 6435 5005 3003 1365  455  105   15    1
plot(choose(n,0:n),type="h")

(a) Sprendimas Monte Carlo metodu. Klasikiniu atveju, tikimybę apibrėžus formule $P(A)=N(A)/N$, galima sukurti visai turiningą jos teoriją. Kita vertus, nėra visai aišku, kaip derėtų suprasti ar interpretuoti šį skaičių. Vienas iš pasiūlymų yra tikimybę aiškinti kaip santykinio dažnio ribą: jei eksperimentą pakartotume n kartų, tai tikimybė $P(A)$ lygi santykio $n(A)/n$ ribai (čia $n(A)$ nurodo kelis kartus, atlikus $n$ eksperimentus, pasirodė įvykis $A$). Monte Carlo metodas siūlo eksperimentus kartoti „daug“ kartų (1000, 10000 ir pan.), gautas santykis $n(A)/n$ ir bus apytiksliai lygus tikimybei.

set.seed(12345)
NN=5000
kand.v=paste("v",1:5)
kand.m=paste("m",1:10)
kand=c(kand.v,kand.m)
kand
sample(kand,5)        # negrąžintinės dydžio 5 atrankos iš kand pavyzdys

Dabar šį atrankos eksperimentą pakartosime „daug“ (NN=5000) kartų.

aa=numeric(NN)       
aa
 
for(i in 1:NN)
{
kom=sample(kand,5)   
aa[i]=ifelse(length(grep("v",kom))==2 & length(grep("m",kom))==3, 1, 0)
}
aa
 
kom            # [1] "v 3"  "m 3"  "m 4"  "m 10" "m 9"
grep("v",kom)  # [1] 1
tik.a=sum(aa)/NN
tik.a          # atsakymas: 0.4002
 
par(mfrow=c(1,2))                          # bus du grafikai vienoje eilutėje
plot(cumsum(aa)/(1:NN),type="l")           # santykinio dažnio grafikas
?cumsum
abline(400/1001,0,col=2)                   # col=2  ~ color=”red”
?abline
plot(cumsum(aa)[100:NN]/(100:NN),type="l") # vaizdesni santykinio dažnio
                                           # grafikas nuo 100-jo sekos nario
abline(400/1001,0,col=2)

 

1 užduotis. Trimis aptartais būdais atsakykite į (b) klausimą. 

2 pavyzdys. Vadinamasis Polya (tark Pojà, Pòjos, …, su Pojà)  urnų modelis yra apibrėžiamas taip: iš pradžių urnoje yra 1 baltas ir 1 juodas rutulys; pirmuoju žingsniu atsitiktinai traukiame vieną rutulį, po to jį gražiname atgal ir urną papildome dar vienu tos pačios spalvos rutuliu; kitais žingsniais ta pati procedūra kartojama. Kam lygi tikimybė $P_{n}(i), 1 \le i \le{n-1}$, kad kai urnoje bus $n$ rutulių, $i$ iš jų bus balti.

Polya modelis yra tam tikra prasme atvirkščias negrąžintinei atrankai (ten, ištraukus baltą rutulį, šansai kitą kartą ištraukti baltą sumažėja, o čia padidėja). Tradicinis sprendimas šį kartą sudėtingas, o atsakymas gana netikėtas: $P_{n}(i) \equiv 1/(n-1)$, t.y., lygiai tikėtina, kad urnoje galų gale bus 1, 2, ... ar $n-1$ baltas rutulys. Mes šį uždavinį išspręsime Monte Carlo metodu (išspręsime dalinai kadangi i) tik su viena $n$ reikšme ($n=11$ ) ir ii) atsakymas tik „beveik“ garantuoja, kad $1/(n-1)$ yra teisingas atsakymas). 

Pradėsime pavyzdžiu: devynis kartus trauksime rutulį (po to urnoje bus 11 rutulių).

set.seed(1)                      # Fiksuojame "atsitiktinumo" pradžią, todėl visi
                                 # kartojantys šį sprendimą gaus tą patį atsakymą
n=11
urnoje=c("b","j")
urnoje
 
for(i in 1:(n-2))                # 9 kartus modeliuosime rutulio traukimą
{
res=sample(urnoje,1)
res
if(res=="b") (urnoje = c(urnoje,"b")) else (urnoje = c(urnoje,"j"))
print(urnoje)                    # pamatysime kas yra urnoje po kiekvieno traukimo
}
 
[1] "b" "j" "b"
[1] "b" "j" "b" "j"
[1] "b" "j" "b" "j" "b"
[1] "b" "j" "b" "j" "b" "b"
[1] "b" "j" "b" "j" "b" "b" "j"
[1] "b" "j" "b" "j" "b" "b" "j" "j"
[1] "b" "j" "b" "j" "b" "b" "j" "j" "j"
[1] "b" "j" "b" "j" "b" "b" "j" "j" "j" "b"
[1] "b" "j" "b" "j" "b" "b" "j" "j" "j" "b" "j"  # 5 balti rutulai urnoje
 

Dabar mūsų eksperimentą pakartosime „daug“ (10000) kartų ir nustatysime keli eksperimentai baigėsi vienu baltu rutuliu urnoje, dviem, ..., dešimčia; santykinis dažnis bus maždaug lygus tikimybei.

set.seed(2)
NN=10000
n=11
ss=numeric(NN)
for(ii in 1:NN)
{
urnoje=c("b","j")
for(i in 1:(n-2))
{
res=sample(urnoje,1)
res
if(res=="b") (urnoje = c(urnoje,"b")) else (urnoje = c(urnoje,"j"))
}
ss[ii]=sum(urnoje=="b")    # baltų rutulių skaičius po 9 traukimų
}
table(ss)/NN               # santykinis įvairių baigčių dažnis
plot(table(ss)/NN)
ss
     1      2      3      4      5      6      7      8      9     10
0.0993 0.0989 0.1022 0.0979 0.1021 0.1032 0.0997 0.0981 0.0973 0.1013

(taigi visos tikimybės artimos 1/10).

3 pavyzdys. Pademonstruosime kaip, gana netikėtai, MC metodas gali būti panaudotas skaičiui $\pi  (=3.14159)$ apskaičiuoti (ar prisimenate kaip apibrėžiamas skaičius $\pi$?). Su R‘u generuosime seką $\vec{\alpha_1},\vec{\alpha_2},... $ dvimačių atsitiktinių vektorių arba, kitaip sakant, seką plokštumos taškų $\vec{\alpha_i}=(\alpha_i^{(1)},\alpha_i^{(2)})$, turinčių tolygų skirstinį kvadrate su viršūnėmis taškuose (-1,-1), (1,-1), (1,1) ir (-1,1). Vieną tokio atsitiktinio vektoriaus kopiją (t.y., dvi vektoriaus koordinates, t.y., du skaičius) galima generuoti su runif(2,-1,1) (funkcija runif generuoja random uniform skaičius), o išbrėžti penkiasdešimties tokių taškų grafiką galima taip:

xx <- c(-1,1,1,-1,-1)
yy <- c(-1,-1,1,1,-1)
par(mfrow=c(1,3))   # brėš tris grafikus vienoje eilutėje         
plot(xx,yy,type="l", main="50 taškų kvadrate")  
set.seed(12)
runif(2,-1,1)       # [1] -0.8612782  0.6355504
points(runif(50,-1,1),runif(50,-1,1),pch="*")

(pakartokite šią programą kelis kartus, gausite vis kitus taškus).

Priminsime, kad tolygiojo skirstinio kvadrate $S$ atveju tikimybė taškui pakliūti į aibę $A$  apibrėžiama taip: $P(\vec \alpha \in A)=L(A)/L(S)=L(A)/4$ ; čia $L(A)$  yra aibės $A$  Lebesgue (tark Lebègas, Lebego, su Lebeg̀u) matas, t.y., tiesiog jos plotas. Iš čia išplaukia, kad tikimybė pakliūti į vienetinį skritulį $C$  lygi  $\pi/4$. Jei atliksime daug (10000 ar pan.) MC eksperimentų, tai santykinis dažnis taškų pakliuvusių į $C$  bus maždaug lygus $\pi/4$ . Su žemiau pateiktu R'o kodu apytikslai apskaičiuosime $\pi$: 

xx <- c(-1,1,1,-1,-1)
yy <- c(-1,-1,1,1,-1)
plot(xx,yy,type="l", main="500 taškų")  # brėžiame kvadratą S
x <- cos(seq(0,2*pi,length=100))
y <- sin(seq(0,2*pi,length=100))  # apskritimas polinėse koordinatėse
polygon(x,y,col=3)                # nuspalviname apskritimą
points(runif(500,-1,1),runif(500,-1,1),pch="*") # žr. paveikslą žemiau
set.seed(10)
xxx <- runif(100000,-1,1) # atliekame 100000 MC eksperimentų, t.y.,
yyy <- runif(100000,-1,1) # metame 100000 taškų į S
s <- numeric(500)
print(s)                  # vektorius sudarytas iš 500 nulių
# ciklas: vertinsime santykinį dažnį grupėse, sudarytose iš 200*i taškų
for(i in 1:500) {s[i] <- 4*sum(ifelse(xxx[1:(200*i)]^2+yyy[1:(200*i)]^2<=1,1,0))/(200*i)}   
print(s)                  # santykiniai dažniai * 4 ($\approx \pi$ )
plot(1:500,s,type="l")
abline(pi,0)
s[500]                    # 500*200=100000 įvertis pagal visus taškus  
# [1] 3.14252
 
            Tolygaus skirstinio iliustracija (kairėje), pirmieji 500 taškų (viduryje); santykinis dažnis artėja prie $\pi$  (dešinėje)

 

Aukščiau pateiktas kodas atrodo sudėtingas dėl grafikų brėžimo. Pati MC procedūrą yra visai paprasta:

set.seed(10)
NN=100000
xxx <- runif(NN,-1,1)
yyy <- runif(NN,-1,1)
pii=4*sum(ifelse(xxx^2+yyy^2<=1,1,0))/NN
pii        # 3.14252

Atkreipsime jūsų dėmesį į tai, kad, nors „metėme” labai daug (100000) taškų, tačiau įverčio tikslumas nėra labai didelis (tai būdinga MC metodui). Antra vertus, jis nepakeičiamas, kai tikimybių skaičiavimai sudėtingi. Beje, pirmuosius 2000 skaičiaus $\pi$ skaitmenis galite pamatyti su

library(UsingR)
pi2000  

2 užduotis. Naudodami MC metodą apskaičiuokite plotą po viršutine sinusoidės banga. O tiems, kurie žino integralą, priminsiu - šis plotas lygus $\int_{0}^{\pi} \sin x dx$, t.y., ???

curve(sin, 0, 2*pi)
abline(0,0)

Na o dabar – su interneto naršykle nueikite į Googlą, surinkite Monte Carlo method ir pasiruoškite ilgoms vakaronėms…

 

Komentuoti

Basic HTML

  • Web puslapių adresai ir el. pašto adresai automatiškai tampa nuorodomis.
  • Tags allowed: a, em, strong, u, s, cite, code, blockquote, ol, ul, li, dl, dt, dd, pre, p, br
  • Mathematics inside the configured delimiters is rendered by MathJax. The default math delimiters are $$...$$ and \[...\] for displayed mathematics, and $...$ and \(...\) for in-line mathematics.
  • Syntax highlight code surrounded by the <pre class="brush: lang">...</pre> tags, where lang is one of the following language brushes: php, python, r, sass, sql, vb.
  • Linijos ir paragrafai atskiriami automatiškai

Markdown

  • You can enable syntax highlighting of source code with the following tags: <code>, <blockcode>, <c>, <cpp>, <drupal5>, <drupal6>, <java>, <javascript>, <php>, <python>, <r>, <ruby>. The supported tag styles are: <foo>, [foo].
  • Quick Tips:
    • Two or more spaces at a line's end = Line break
    • Double returns = Paragraph
    • *Single asterisks* or _single underscores_ = Emphasis
    • **Double** or __double__ = Strong
    • This is [a link](http://the.link.example.com "The optional title text")
    For complete details on the Markdown syntax, see the Markdown documentation and Markdown Extra documentation for tables, footnotes, and more.
  • Typographic refinements will be added.
  • Mathematics inside the configured delimiters is rendered by MathJax. The default math delimiters are $$...$$ and \[...\] for displayed mathematics, and $...$ and \(...\) for in-line mathematics.

Plain text

  • HTML žymės neleidžiamos.
  • Web puslapių adresai ir el. pašto adresai automatiškai tampa nuorodomis.
  • Linijos ir paragrafai atskiriami automatiškai
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Target Image

Theme by Danetsoft and Danang Probo Sayekti inspired by Maksimer