Box-plot i autlajeri

Izgledaju ovako:

boxplot(discoveries, col="blue")

Probamo na nekim drugim podacima sa detaljnijim objašnjenjem.

Baza koja sadrži podatke o kvalitetu vazduha u Njujorku od maja do septembra 1973.

airquality[1:10, ]

##    Ozone Solar.R Wind Temp Month Day
## 1     41     190  7.4   67     5   1
## 2     36     118  8.0   72     5   2
## 3     12     149 12.6   74     5   3
## 4     18     313 11.5   62     5   4
## 5     NA      NA 14.3   56     5   5
## 6     28      NA 14.9   66     5   6
## 7     23     299  8.6   65     5   7
## 8     19      99 13.8   59     5   8
## 9      8      19 20.1   61     5   9
## 10    NA     194  8.6   69     5  10

Velicina baze

length(airquality[,1])

## [1] 153

I trazimo njene boks-plotove po svakoj promenljivoj

boxplot(airquality)

Posebno izdvojimo kolonu “Ozone”, jer se u njoj javili autlajeri

Ozon = airquality$Ozone

range(Ozon, na.rm=TRUE)

## [1]   1 168

median(Ozon, na.rm=TRUE)

## [1] 31.5

quantile(Ozon, na.rm=TRUE)

##     0%    25%    50%    75%   100% 
##   1.00  18.00  31.50  63.25 168.00

IQR(Ozon, na.rm=TRUE)

## [1] 45.25

Centralna debela crta u Boks-plotu predstavlja Medijanu posmatranog atributa. Granice pravougaonika su \[quantile(25\%)\] i \[quantile(75\%),\] dok preostale dve crte su \[quantile(25\%) - 1.5*IQR\] i \[quantile(75\%) + 1.5*IQR.\]

quantile(Ozon, na.rm=TRUE)[4] + 1.5*IQR(Ozon, na.rm=TRUE)

##     75% 
## 131.125

sum(Ozon > quantile(Ozon, na.rm=TRUE)[4] + 1.5*IQR(Ozon, na.rm=TRUE), na.rm=TRUE)

## [1] 2

Monte Karlo Simulacije

Igrač baca strelicu na tablu na kojoj se nalaze tri koncentrična kruga sa poluprečnicima 5, 10 i 15 cm. Najveći krug je upisan u kvadrat stranice 30 cm. Smatra se da je svaki deo 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 koordinate za pogodjenu tacku, (x,y) su iz [-15,15]x[-15,15]
  
  pozicija = runif(2, -15, 15)
  
# racunamo rastojanje tacke od centra
  
  rastojanje = sqrt(sum(pozicija^2))  # ekvivalentno: norm(pozicija, type="2")

  if(rastojanje<=5)              
    return(1)
    else if(rastojanje<=10)  
      return(2)
      else if(rastojanje<=15) 
        return(3)
        else          
          return(-1)      

}

Oceniti verovatnoću da je da je pogođen prvi(5cm), odnosno drugi(10cm) krug.

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

simulacije = replicate(n = 1000, expr = strelica(), simplify = "array")


# Napomena: Ako funkciji replicate() prosledimo argument simplify=FALSE, vratice listu

table(simulacije)

## simulacije
##  -1   1   2   3 
## 197  98 265 440

frekvencije = table(simulacije)/1000
frekvencije

## simulacije
##    -1     1     2     3 
## 0.197 0.098 0.265 0.440

# Ako izracunamo verovatnoce teorijski, gledajuci kao gemetrijsku verovatnocu 
# 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 naizmenič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, treći 5 i van kruga 10 poena. Igru započinje igrač A, a završava se kada jedan od igrača izgubi sve poene i samim tim izgubi igru. Simulirati igru. Napisati funkciju pikado() koja vraća 1 ako je pobedio igrač A, odnosno 2 ako je pobedio 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)
  
  # I izracunamo ukupni zbir poena
  return(krug.1 + 3*krug.2 + 5*krug.3 +10*kvadrat)
  
}



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


pikado <- function(){
  
  poeni.A = 20
  poeni.B = 20
  
# Napomena:  
# Indikator da je igrac A na redu: U aritmetickim izrazima TRUE ce imati vr 1, a FALSE 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()))
    
  # Promena igraca
    
  A.na.potezu = !A.na.potezu

  }
  
  # Ko prvi izgubi poene izgubio je igru.
       
  ifelse(poeni.A <= 0, return(2), return(1))
      
  
}


# Narednom linijom pozivamo veliki broj simulacija ove igre radi ocenjivanja verovatnoce pobede svakog igraca
s = replicate(10000, pikado())


mean(s == 1)

## [1] 0.2064

mean(s == 2)

## [1] 0.7936

Monte Karlo integracija

Integracija

Ideja ove metode je da se simulira veliki broj tačaka (nezavisno) iz dvodimenzione uniformne raspodele na nekom pravougaoniku na kojem se nalazi grafik funkcije koju želimo da integralimo. Kako je integral jednak površini ispod grafika funkcije, njegovu ocenu ćemo dobiti tako što procenat tačaka koje se nađu ispod grafika pomnožimo sa površinom celog pravougaonika. Zapravo posmatramo srednju vrednost nezavisnih indikatora, pa će po zakonu velikih brojeva ona skoro sigurno konvergirati ka očekivanju. Formalno posmatramo integral \(\int_a^bf(x)dx\). Pretpostavimo da se maksimum te funkcije postiže u nekoj tački \(c\in[a,b]\).

Generišemo \((U_i^1,U_i^2)\), \(i=1,2,\dots,n,\) gde je \(U_1\in\mathbb{U}[a,b]\) i \(U_2\in\mathbb{U}[0,f(c)]\). Opisana ocena će zapravo biti \[\dfrac{1}{n}\sum_{i=1}^n I\{U_2≤f(U_1)\}\] i po zakonu velikih brojeva ona konvergira ka \[\mathbb{E}(I\{U_2≤f(U_1)\})=\dfrac{1}{(b−a)f(c)}\int_a^bf(x)dx,\] pa pomnoženo sa površinom pravougaonika \([a,b]×[0,f(c)]\) je baš jednako integralu koji ocenjujemo.

Zadaci

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

# max je 1 pa tražimo slucajne tacke 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)
  return(mean(u2 <= f(u1)) * (b - a) * d)
  
}

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

## [1] 1.1932

Prava vrednost integrala:

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

## 1.196288 with absolute error < 1.3e-14

Isti integral pomocu znanja iz ViS:

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

## [1] 1.196288

Drugi način primene MK metode u integraciji je sledeći.

Setimo se kako se računa matematičko očekivanje neke slučajne veličine

\[\mathbb{E}(h(X)) = \int_{\mathcal{X}}h(x)f(x)dx. \] Ako generišemo uzorak iz raspodele koju ima \(X\) i nađemo uzoračku sredinu, onda prema zakonu velikih brojeva, ona će skoro sigurno konvergirati ka integralu \(\mathbb{E}(h(X)).\)

Za isti integral kao u zadatku 2. možemo uzeti da je gustina \(f(x)= \frac{1}{2}\) na \([0, 2]\) dok je \(h(x) = 2e^{-\frac{x^2}{2}}.\) Treba da generišemo uzorak iz uniformne \(\mathbb{U}[0, 2]\) raspodele.

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

## [1] 1.200828

Neka je data standardna normalna raspodela. Naći vrednost njene funkciju raspodele u proizvoljnoj tacki \(t\).

Resenje:

h <- function(t)
{
  u = rnorm(1000000)
  return(mean(u<= t))
}
# sad mozemo u proizvoljnoj tacki to da racunamo. Npr u 1
h(1)

## [1] 0.841256

# ili u 0
h(0)

## [1] 0.499818

Metoda inverzne transformacije

Neka je \(X\) slučajna veličina sa neprekidnom i strogo rastućom funkcijom raspodele \(F\). Sa \(F^{-1}\) ćemo označiti njenu inverznu funkciju. Tada \(X\) možemo generisati na sledeći način: 1) Generišemo broj \(u\) iz uniformne \(\mathbb{U}[0, 1]\) raspodele. 2) Vratimo vrednost \(F^{-1}(u).\)

Time smo dobili vrednost iz raspodele koju ima slučajna veličina \(X\), tj. \(F(X)\in \mathbb{U}[0, 1]\) dok \(F^{-1}(U)\) ima raspodelu kao i \(X\).

Naći inverznu transformaciju za proizvoljnu uniformnu \(\mathbb{U}[a, b]\) raspodelu.

(Rešenje:) kako je \(U\in\mathbb{U}[0, 1]\implies a+ (b-a)U\in \mathbb{U}[a, b]\) to onda funkcijom

proizvoljna_uniformna <- function(a, b)
{
  U = runif(1)
  
  return(a+(b-a)*U)
}

replicate(n = 10, proizvoljna_uniformna(a=2, b=7))

##  [1] 4.117767 3.014382 2.736245 4.888473 3.599582 3.417625 2.641420 3.624007
##  [9] 6.168204 6.531842

Naći inverznu transformaciju za proizvoljnu Eksponencijalnu \(\mathcal{E}(a), a > 0\) raspodelu.

(Rešenje:) kako je \(U\in\mathbb{U}[0, 1]\implies -\frac{1}{a}\log(1-U)\in \mathcal{E}(a)\) to onda funkcijom

proizvoljna_eksponencijalna <- function(a)
{
  U = runif(1)
  
  return(-1/a * log(1-U))

}

uzorak = replicate(n = 1000, proizvoljna_eksponencijalna(a=2))

Dobijamo traženi uzorak.

I možemo pomoću histograma da proverimo da li je uzorak zaista iz tražene raspodele

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

Naći inverznu transformaciju za proizvoljnu Gama \(Gamma(n,\lambda), \lambda> 0\) raspodelu.

(Rešenje:) kako je \(U\in\mathbb{U}[0, 1]\implies -\frac{1}{a}\log(1-U)\in \mathcal{E}(a)\), a znamo i da je Gama raspodela sa celobrojnim parametrom \(n\) zapravo zbir nezavisnih eksponencijalnih, tada funkcijom

# nije bas skroz proizvoljna, ipak imamo uslov ovde da je n - prirodan broj
proizvoljna_gama <- function(n, lambda)
{
  return(sum(replicate(n=n, proizvoljna_eksponencijalna(a=lambda))))
}

uzorak = replicate(n = 1000, proizvoljna_gama(n=13, lambda = 1/3))