MONTE KARLO SIMULACIJE

Igrač baca strelicu na tablu na kojoj se nalaze tri koncentrična kruga poluprečnika 5, 10 i 15 cm. Najveći krug je upisan u kvadrat stranice 30 cm. Smatra se da je svaki dio table ravnopravan.

Napisati funkciju strelica() koja simulira opisanu igru. (Funkcija treba da vrati 1-ako je pogođen prvi krug, 2 za drugi i 3 za treći, -1 ako je strelica van krugova)

strelica<-function(){

# Postavljamo koordinatni pocetak u centar kruga i na slucajan nacin biramo x i 
# y koor- dinatu za pogodjenu tacku, (x,y) su iz [-15,15]x[-15,15]
  
  poz<-runif(2, -15, 15)
  
# racunamo rastojanje tacke od centra
  
  r<-norm(poz, type="2")  # ekvivalentno: sqrt(sum(poz^2))

  if(r<=5)              return(1)
    else if(r>5 & r<=10)  return(2)
      else if(r>10 & r<=15) return(3)
        else                  return(-1)      

}

Ocijeniti vjerovatnocu da je da je pogodjen 1., odnosno 2. krug.

# Koristimo funkciju replicate cijim pozivom 
# 1000 puta pozivamo funkciju strelica() i rezultate tih poziva smjestamo 
# u jedan vektor

simulacije<-replicate(1000, strelica(), simplify = "array")


# Napomena: Ako fji replicate proslijedimo argument simplify= F, vratice listu

table(simulacije)

## simulacije
##  -1   1   2   3 
## 219 100 285 396

frekvencije<-table(simulacije)/1000
frekvencije

## simulacije
##    -1     1     2     3 
## 0.219 0.100 0.285 0.396

# Ako izracunamo vjerovatnoce teorijski, gledajuci kao gemetrijsku vjerovatnocu 
# po formuli: P(A)=m(A)/m(kvadrat) , dobicemo priblizno:

#  -1           1             2         3

#  0.21460    0.08727    0.26180    0.43633

Igrači A i B naizmjenično bacaju po tri strelice na tablu. Na početku igre imaju po 20 poena. Ako se pogodi prvi krug gubi se 1 poen, drugi 3, treci 5 i van kruga 10 poena. Igru započinje igrač A, a završava se kada jedan od igrača izgubi sve poene. Simulirati igru. Napisati funkciju pikado() koja vraća 1 ako je pobijedio igrač A, odnosno 2 ako je pobijedio igrač B.

 # Pravimo pomocnu fju koja racuna koliko treba oduzeti poena igracu nakon jedne
 # partije. Prosledjujemo joj vektor sa oznakama pogodjenih polja


poeni<-function(x){
  
  # racunamo koliko je puta pogodjeno svako polje
  
  krug.1<-sum(x==1)
  krug.2<-sum(x==2)
  krug.3<-sum(x==3)
  kvadrat<-sum(x==-1)
  
  return(krug.1 + 3*krug.2 + 5*krug.3 +10*kvadrat)
  
}



# Hocemo da napisemo funckiju pikado() koja vraca 1 ako je pobijedio igrac A, a
# 2 u suprotnom.


pikado<-function(){
  
  poeni.A<-20
  poeni.B<-20
  
# Napomena:  
# Indikator da je igrac A na redu: U aritmeticki izrazima T ce imati vr 1, a F 0
  
  A.na.potezu<-TRUE
  
  
  while(poeni.A*poeni.B>0){
    
  # Oduzimamo poeni igracu koji je na potezu
    
  poeni.A<-poeni.A - A.na.potezu*poeni(replicate(3,strelica()))
  poeni.B<-poeni.B - (!A.na.potezu)*poeni(replicate(3,strelica()))
    
  # Mijenjamo igraca
    
  A.na.potezu<-!A.na.potezu
    
    
  }
  
  # Ko prvi izgubi poene izgubio je igru.
       
  ifelse(poeni.A<=0 ,return (2), return(1))
      
  
}

# Ocijeniti vjerovatnocu da pobijedi igrac koji zapocinje igru, odnosno igrac A.
# Sta se ocekuje?


s<-replicate(1000,pikado())


mean(s==1)

## [1] 0.202

mean(s==2)

## [1] 0.798

MONTE KARLO METODE INTEGRACIJE

“Hit or miss”

Ideja ove metode je da se simulira veliki broj tačaka (nezavisno) iz dvodimenzione uniformne raspodjele na nekom pravougaoniku na kojem se nalazi grafik funkcije koju želimo da integralimo. Kako je integral jednak površini ispod grafika funkcije, njegovu ocjenu ćemo dobiti tako što procenat tačaka koje se nađu ispod grafika pomnožimo sa površinom cijelog pravougaonika.
Mi ovdje zapravo posmatramo srednju vrijednost nezavisnih indikatora, pa će po zakonu velikih brojeva ona konvergirati kao očekivanju, s.s. Preciznije: Posmatramo integral \(\int\limits_{a}^{b}f(x)dx\). Pretpostavimo da se maksimum te funkcije postiže u nekoj tački \(c \in [a,b]\). Generišemo \((U_1^{i},U_2^{i}),\ i=1,2,\dots,n\), gdje je \(U_1 \in \mathcal{U}[a,b]\ \text{ i } \mathcal{U_2} \in \mathcal{U}[0,f(c)]\). Opisana ocjena će zapravo biti \(\frac{1}{n}\sum\limits_{i=1}^{n}I\{U_2\leq f(U_1)\}\) i po zakonu velikih brojeva ona konvergira ka \(E(I\{U_2 \leq f(U_1)\})=\frac{1}{(b-a)f(c)}\int\limits_{a}^{b}f(x)dx\), pa pomnoženo sa površinom pravougaonika \([a,b]\times[0,f(c)]\) je baš jednako integralu koji ocjenjujemo.

Primjeri

Izračunati integral \(\int\limits_{0}^{2}e^{-\frac{x^2}{2}}dx\).

curve(exp(-(x ^ 2 / 2)), 0, 2)

# max je 1 pa tražimo slučajne tačke na pravougaoniku [0,2]x[0,1]

# neka je d=f(c), tj max 
mk.integral <- function(n, a, b, d, f) {
  u1 <- runif(n, a, b)
  u2 <- runif(n, 0, d)
  mean(u2 <= f(u1)) * (b - a) * d
  
}

mk.integral(10000, 0, 2, 1, function(x)
exp(-(x ^ 2 / 2)))

## [1] 1.205

integrate(function(x)
exp(-(x ^ 2 / 2)), 0, 2)

## 1.196288 with absolute error < 1.3e-14

# ili
sqrt(2 * pi) * (pnorm(2) - pnorm(0))

## [1] 1.196288

Aproksimacija broja \(\pi\).

curve(sqrt(1-x^2),0,1)

# Iskoristimo funckiju iz prethodnog primjera
4 * mk.integral(100000, 0, 1, 1, function(x)
  sqrt(1 - x ^ 2))

## [1] 3.13992

Monte Karlo integracija

Vratimo se na problem računanja očekivanja neke slučajne veličine \[ E(h(X))=\int\limits_{\mathcal{X}}h(x)f(x)dx \] Ako generišemo uzorak iz te raspodjele koju ima \(X\) i uzmemo uzoračku sredini, po zakonu velikih brojeva, ona će konvergirati skoro sigurno ka integralu \(E(h(X))\).

Za isti primjer \(\int\limits_{0}^{2}e^{-\frac{x^2}{2}}dx\) možemo uzeti da je gustina \(f(x)=1/2,\ \text{na } [0,2]\), a \(h(x)=2e^{-\frac{x^2}{2}}\), tj. treba da generišemo uzorak iz uniformne \(\mathcal{U}[0,2]\) raspodjele.

u <- runif(10000, 0, 2)
x <- 2 * exp(-(u ^ 2 / 2))
mean(x)

## [1] 1.197275

Računanje funkcije standardne normalne raspodjele u nekoj tački t.

\[ \Phi(t)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{t}e^{-\frac{x^2}{2}}dx=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}I\{x \leq t\}e^{-\frac{x^2}{2}}dx \] Generisaćemo brojeve iz normalne raspodjele i tražiti očekivanje indikatora iz izraza gore.

u<-rnorm(10^7)
g<-function(t){
  u<=t
}
mean(g(0))

## [1] 0.4999626

Izračunati \(P\{X>20\}\) ako je \(X \in \mathcal{N}(0,1)\). Kako je ovo vrijednost veoma mala, ako iskoristimo prethodnu metodu dobićemo uvijek vrijednost 0.

1-mean(g(20))

## [1] 0

Kako je normalna raspodjela simetrična, \(P\{X>20\}=P\{X<-20\}\).Zato ćemo posmatrati intgral oblika \(\int\limits_{-\infty}^{t}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}dx,\ t<0\) koji se smjenom \(x=\frac{1}{u}\) svodi na sledeće: \[ \int\limits_{\frac{1}{t}}^{0}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2u^2}}\frac{1}{u^2}du \] Nama u zadatku treba specijalno za \(t=-20\). Sada se ovo svodi na očekivanje \(g(U)\) gdje je je \(U \in \mathcal{U}[-\frac{1}{20},0]\)

u<-runif(10^5,-1/20,0)
g<-function(t){
  1/sqrt(2*pi)*exp(-1/(2*(t^2)))*1/(t^2)
}
mean(1/20*g(u))

## [1] 2.62105e-89

METODA INVERZNE TRANSFORMACIJE

Neka je X slučajna veličina sa neprekidnom i rastućom funkcijom raspodjele \(F\). Označimo sa \(F^{-1}\) inverz funkcije \(F\). Tada slučajnu veličinu X generišemo na sledeći način:

-Generišemo slučajan broj U iz \(U(0,1)\) raspodjele -Vratimo vrijednost \(F^{-1}(U)\)

Na ovaj način dobijamo vrijednosti iz raspodjele X jer pri ovim uslovima važi \(F(X) \in U(0,1)\), odosno \(F^{-1}(U)\) ima istu raspodjelu kao X.

Uniformna raspodjela na proizvoljnom intervalu (a,b)

\[ U \in \mathcal{U}(0,1) \rightarrow a+(b-a)U \in \mathcal{U}(a,b) \]

unif<-function(a,b){
U<-runif(1)
a+(b-a)*U
}

Eksponencijalna raspodjela

\[ U \in \mathcal{U}(0,1) \rightarrow -\frac{1}{\lambda}\log(1-U) \in \mathcal{E}(\lambda),\quad \lambda>0 \]

eksponencijalna<-function(lambda){

U <- runif(1)
-1/lambda * log(1-U)
# ili -1/lambda * log(U)
}

Hoćemo da provjerimo da je uzorak dobijen ovom metodom bas iz eksponencijalne raspodjele.

uzorak <- replicate(1000,eksponencijalna(3))

# srednja vrijednost (EX=1/3)
mean(uzorak)

## [1] 0.3381786

# disperzija (DX=1/9)
var(uzorak)

## [1] 0.1090279

# histogram vjerovatnoca

library ( lattice )
hist (uzorak, prob =TRUE , xlab ="",ylab ="",col ='coral1',border ='bisque', main="Exp(3)")
curve ( dexp (x , rate=3) ,add =TRUE , lwd =3, col ='cornflowerblue ')

# Ocjena vjerovatnoce P{X>2} za X iz Exp(3) raspodjele

mean(uzorak>2)

## [1] 0.003

# stvarno vrijednost 

1-pexp(2, rate=3)

## [1] 0.002478752

Primjer: gama raspodjela kao zbir nezavisnih eksponencijalnih

gama<-function(n,lambda){

r<-replicate(n,eksponencijalna(lambda))
sum(r)

}

uzorak <- replicate(1000,gama(10,0.5))
hist (uzorak, prob =TRUE , xlab ="",ylab ="",col ='coral1',border ='bisque', main="Gama(10,0.5)", ylim=c(0,0.07))
curve ( dgamma (x , rate=0.5, shape=10) ,add =TRUE , lwd =3, col ='cornflowerblue ')

mean(uzorak) # EX = 10/0.5

## [1] 19.86923

var(uzorak)  # DX = 10/0.25

## [1] 39.84313