Монте Карло методе

Монте Карло методе, или Монте Карло експерименти, представљају широку класу рачунарских алгоритама који се ослањају на поновљено случајно узорковање ради добијања нумеричких резултата. Основни концепт је коришћење случајности за решавање проблема који би у принципу могли бити детерминисани.

Монте Карло методе се углавном користе у три различите класе проблема: * оптимизација, * нумеричка интеграција, * генерисања узорака из расподела вероватноће.

Монте Карло методе се често примењују коришћењем рачунарских симулација и могу пружити приближна решења проблема који су иначе нерешиви или превише сложени за математичку анализу. Монте Карло методе такође имају нека ограничења и изазове, као што су компромис између тачности и рачунарског трошка, проблем димензионалности, поузданост генератора случајних бројева и верификација и валидација резултата.

У наставку су два уобичајена примера, Закон великих бројева (ЗВБ) и Централна гранична теорема (ЦГТ), који користе Монте Карло симулациону методу.

Централна гранична теорема (ЦГТ)

Теорема: Нека је \((X_n)\) низ независних случајних величина са истом расподелом, математичким очекивањем \(EX_1=\mu\) и коначном дисперзијом \(DX_1 = \sigma^2\). Ако је \(S_n=\sum\limits_{k=1}^n X_k\), онда за свако \(x\in\mathbb{R}\) при \(n\to\infty\) важи \[P\left\{\frac{S_n - n\mu}{\sigma\sqrt{n}}\le x\right\} \to \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x e^{-\frac{t^2}{2}}dt,\] односно \[\frac{S_n - ES_n}{\sqrt{DS_n}} \to X^*, \quad X^*\in \mathcal{N}(0,1).\]

# Central Limit Theorem Simulation
# n: sample size of each sample
# dist: underlying distribution where the sample is drawn

simulation <- function(n, dist=NULL, param1=NULL, param2=NULL) {
r <- 10000 # Number of replications/samples 

# This produces a matrix of observations with
# n columns and r rows. Each row is one sample:

my.samples <- switch(dist,
"E" = matrix(rexp(n*r,param1),r),
"N" = matrix(rnorm(n*r,param1,param2),r),
"U" = matrix(runif(n*r,param1,param2),r),
"P" = matrix(rpois(n*r,param1),r),
"C" = matrix(rcauchy(n*r,param1,param2),r),
"B" = matrix(rbinom(n*r,param1,param2),r),
"G" = matrix(rgamma(n*r,param1,param2),r),
"X" = matrix(rchisq(n*r,param1),r),
"T" = matrix(rt(n*r,param1),r))

all.sample.sums <- apply(my.samples,1,sum)
all.sample.means <- apply(my.samples,1,mean)
all.sample.vars <- apply(my.samples,1,var)

par(mfrow=c(2,2))
hist(my.samples[1,], col="gray", main="Distribution of One Sample", xlab="")
hist(all.sample.sums, col="gray", main="Sampling Distribution of the Sum", xlab="")
hist(all.sample.means,col="gray", main="Sampling Distribution of the Mean", xlab="")
hist(all.sample.vars,col="gray", main="Sampling Distribution of the Variance", xlab="")
}

simulation(n=1000, dist="P", param1=1)

simulation(n=10000, dist="E", param1=1)

Закон великих бројева (ЗВБ)

Теорема: Нека је \((X_n)\) низ случајних величина дефинисаних на истом простору вероватноћа. Нека је \(S_n = \sum\limits_{k=1}^n X_k\). Ако \(\frac{S_n-ES_n}{n}\) конвергира скоро сигурно (у вероватноћи) ка нули када \(n\to\infty\), онда важи јаки (слаби) закон великих бројева за низ \((X_n)\).

set.seed(1212)
n = 50000

x = sample(0:1, n, repl = TRUE)
s = cumsum(x)

r = s/(1:n)

plot(r, ylim = c(0.4, 0.6), type = "l", 
     xlab = "tosses", ylab = "probability of heads up", main = "LLN simulation graph")
lines(c(0,n), c(0.5, 0.5))

round(cbind(x, s, r), 5)[1:10, ]

##       x s       r
##  [1,] 1 1 1.00000
##  [2,] 1 2 1.00000
##  [3,] 0 2 0.66667
##  [4,] 0 2 0.50000
##  [5,] 1 3 0.60000
##  [6,] 1 4 0.66667
##  [7,] 1 5 0.71429
##  [8,] 1 6 0.75000
##  [9,] 0 6 0.66667
## [10,] 1 7 0.70000

r[n]

## [1] 0.5013

Апроксимација броја \(\pi\)

Посматрајмо квадрат \(\mathcal{S}\subset \mathbb{R}^2\) димензије \(2\times 2\) у који је уписана кружница \(\mathcal{D}\) полупречника 1.

Тада је \[\frac{\iint_{\mathcal{D}}dx_1dx_2}{\iint_{\mathcal{S}}dx_1dx_2} = \frac{\pi}{4}.\]

Претходне интеграле можемо оценити симулацијама јер је \[\frac{\iint_{\mathcal{D}}dx_1dx_2}{\iint_{\mathcal{S}}dx_1dx_2} = \iint_{\mathcal{S}} \mathbb{I}((x_1,x_2)\in\mathcal{D}) \frac{1}{4}dx_1dx_2 = \mathbb{E}[\phi(X_1,X_2)] = \theta.\]

Да бисмо генерисали униформну расподелу на квадрату \(\mathcal{S} = (-1,1)\times (-1,1)\), можемо посматрати \[X_1 = 2U_1 - 1, \, X_2 = 2U_2-1,\] где је \(U_1,U_2\in\mathcal{U}(0,1)\).

# Number of points to generate
n <- 10000

# Generate random points
set.seed(123)
xs <- runif(n, min=-1, max=1)
ys <- runif(n, min=-1, max=1)

# Calculate distance from (0,0) and check if inside the unit circle
in.circle <- xs^2 + ys^2 <= 1^2

# Estimate π
mc.pi <- (sum(in.circle)/n)*4
mc.pi

## [1] 3.1576

# Plot the points
plot(xs,ys,pch='.',col=ifelse(in.circle,"blue","grey")
     ,xlab='',ylab='',asp=1,
     main=paste("MC Approximation of Pi =",mc.pi))

Модел динамике популације

Основни модел динамике популације: на почетку прве године, бројност неке житарице износи 1000 јединки и расте просечним фактором од 1.1 годишње (процеси репродукције и смрти резултују стопом раста од \(10\%\)) пре жетве. Међутим, стопа раста варира насумично тако да сваке године, фактор раста од 1.1 има варијабилност која се уводи малим променама у процесима преживљавања и репродукције. Моделујте те варијације као лог-нормалне случајне променљиве. Након производње, \(8\%\) популације се уклања жетвом сваке године. Симулирајте бројност популације на крају сваке године током наредних 100 година.

nt = 100       # number of years
N = NULL       # container for abundance
N[1] = 1000    # first end-of-year abundance

for (t in 2:nt) {
  # N this year is N last year * growth *
    # randomness * fraction that survive harvest
  N[t] = (N[t-1] * 1.1 * rlnorm(1, 0, 0.1)) * (1 - 0.08)
}

plot(N, type = "l", pch = 15, xlab = "Year", ylab = "Abundance")

Често, симулације издвајамо у посебене функције како бисмо их репродуковали одређени број пута.

pop_sim = function(nt, grow, sd_grow, U, plot = F) {
  N = NULL # empty flexible vector container
  N[1] = 1000
  for (t in 2:nt) {
    N[t] = (N[t-1] * grow * rlnorm(1, 0, sd_grow)) * (1 - U)
  }
  
  if (plot) {
    plot(N, type = "l", pch = 15, xlab = "Year", ylab = "Abundance")
  }
  
  N
}

pop_sim(100, 1.1, 0.1, 0.08, T)

##   [1] 1000.0000  815.4564  815.4167  817.0242  865.3543 1004.3206 1113.3513
##   [8] 1043.4588 1117.4382 1139.5047 1212.2803 1179.8882 1116.3126 1251.9174
##  [15] 1281.0640 1387.0378 1316.4175 1010.1124  975.4666 1001.4428 1239.0185
##  [22] 1250.0057 1007.1856  971.3801  890.4844  947.6862  934.8841 1008.8605
##  [29] 1117.3076 1030.6148 1159.2436 1091.4611 1084.5315 1058.8098  868.2929
##  [36]  847.8565  697.4140  607.6462  627.3836  656.5371  711.3591  673.7285
##  [43]  716.5796  813.3450  786.1690  841.8792  721.0427  702.2719  710.0972
##  [50]  756.2345  697.2046  687.6382  704.1375  780.0096  854.7118  850.8419
##  [57]  886.5550  923.4473  864.1649 1109.7360 1146.7963 1217.7398 1078.7320
##  [64] 1066.1903 1059.8952 1162.7587 1235.9274 1225.3796 1301.3378 1367.9983
##  [71] 1536.6519 1742.0344 1854.3053 1672.5129 1641.0145 1600.9557 1888.2666
##  [78] 1870.9247 1774.0627 2187.6086 2590.0406 2591.4637 2759.9756 2491.9945
##  [85] 2673.8482 2557.7268 2551.3342 2374.5391 1942.3399 1709.0153 1881.5158
##  [92] 1663.7164 1585.9097 1458.6386 1454.9041 1562.3943 1528.3993 1405.2286
##  [99] 1509.1027 1544.7643

out = replicate(n = 1000, expr = pop_sim(100, 1.1, 0.1, 0.08, F))

Одредимо основне дескриптивне статистике на основу симулација.

N_mean = apply(out, 1, mean)
N_mean[1:10]

##  [1] 1000.000 1024.193 1042.864 1055.457 1072.159 1090.415 1109.529 1130.774
##  [9] 1153.229 1170.192

N_quants = apply(out, 1, function(x) quantile(x, c(0.1, 0.9)))
N_quants[,1:10]

##     [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
## 10% 1000  895.7908  863.7127  841.4886  816.9554  810.9626  801.2245  785.6669
## 90% 1000 1158.0943 1235.2097 1289.0765 1351.1922 1406.4950 1452.1362 1523.6120
##          [,9]     [,10]
## 10%  770.3563  776.1377
## 90% 1584.2758 1642.5945

plot(N_mean, type = "l", ylim = c(0, 10000))
lines(N_quants[1,], lty = 2)
lines(N_quants[2,], lty = 2)

Да бисмо одредили фреквенције појављивања јединки можемо користити функцију table. Нпр. занима нас колико итерација је у десетој години било мање од 1000 јединки.

out10 = ifelse(out[10,] < 1000, "less10", "greater10")
table(out10)

## out10
## greater10    less10 
##       636       364

Случајно лутање

Једнодимензионо случајно лутање: Честица се креће по правој, полазећи из тачке 0 и прави јединични корак у позитивном смеру са вероватночом \(p\), а јединични корак у негативном смеру са вероватночом \(1-p\). Кораци у лутању су дати низом независних случајних величина \(\{X_k, \, k\in\mathbb{N}\}\), за који важи \[P\{X_k = 1\} = p, \,\, P\{X_k = -1\} = 1- p.\] Случајни низ \(\{S_n, n\in\mathbb{N}_0\}\), чији је општи члан дефинисан са: \[S_0 = 0, \,\, S_n = X_1+\dots+X_n, \, ,\ge 1,\] представља положај честице после \(n\) корака.

set.seed(1)

n <- 1000
x <- cumsum(sample(c(-1, 1), n, TRUE))
plot(x, type = "l")

Нека је \(T\) случајна величина која представља тренутак када се честица први пут врати у нулу, тј.

\[T = \min\{n\in\mathbb{N}_+:\, X_n = 0\}.\] Приметимо да је повратак у нулу могућ само након парног броја корака, те случајна величина \(T\) узима вредности \(\{2,4,\dots\} \cup\{+\infty\}\).

Функција расподеле случајне величине \(T_{2n}\) дата је са \[\mathbb{P}\{T = 2 n\} = \binom{2 n}{n} \frac{1}{2 n - 1} p^n (1 - p)^n, \quad n \in \mathbb{N}_+.\]

Пример: Ана и Тамара бацају фер новчић. Ана добија бод ако падне глава, а Тамара добија бод ако падне писмо. Одредити вероватночу да ће имати исти број поена први пут након \(n\) бацања новчића, за свако \(n\in\{2,4,6,8,10\}\). (решење: \(f(2) = \frac{1}{2}\), \(f(4) = \frac{1}{8}\), \(f(6) = \frac{1}{16}\), \(f(8) = \frac{5}{128}\), \(f(10) = \frac{7}{512}\))

simulate_equal_scores <- function(n, trials = 100000) {
  count <- 0 
  
  for (trial in 1:trials) {
    tosses <- sample(c("H", "T"), size = n, replace = TRUE)
    
    heads <- cumsum(tosses == "H")
    tails <- cumsum(tosses == "T")
    
    if (heads[n] == tails[n] && all(heads[-n] != tails[-n])) {
      count <- count + 1
    }
  }
  
  return(count / trials)
}

n_values <- c(2, 4, 6, 8, 10)

set.seed(123) 
simulated_probabilities <- sapply(n_values, simulate_equal_scores)

names(simulated_probabilities) <- n_values
print(simulated_probabilities)

##       2       4       6       8      10 
## 0.50196 0.12507 0.06229 0.03848 0.02801

Stock Market

Акцијама компаније CountBayes тргује се на берзи под ознаком BAYZ. У просеку акција добија 1.001 пута већу вредност од своје почетне цене током једног дана, али то може варирати са стандардном девијацијом од 0.005 на било који дан (волатилност). Можемо симулирати један узорак путање за BAYZ акције тако што ћемо израчунати кумулативни производ из нормалне расподеле са средњом вредношћу од 1.001 и стандардном девијацијом од 0.005. Под претпоставком да BAYZ акције имају почетну цену 20 долара по акцији, илуструјемо пример путање акција за 200 дана.

days <- 200
changes <- rnorm(200, mean = 1.001, sd = 0.005)
plot(cumprod(c(20,changes)), type='l', ylab="Price", xlab="day", main="BAYZ closing price (sample path)")

Али ово је само једна могућа реализација кретања акција! Ако размишљате о инвестирању у BAYZ акције, желите да знате које су могуће цене акције на крају 200. дана. Да бисте проценили ризик ове акције, морате знати који су интервали поверења горње и доње границе будуће цене.

Да бисмо решили ово, симулираћемо 100.000 различитих могућих путања које акција може узети и затим ћемо анализирати расподелу цена за последњи дан.

runs <- 100000

generate.path <- function(){
  days <- 200
  changes <- rnorm(200,mean=1.001,sd=0.005)
  sample.path <- cumprod(c(20,changes))
  closing.price <- sample.path[days+1] # +1 because we add the opening price
  return(closing.price)
}

mc.closing <- replicate(runs,generate.path())

quantile(mc.closing, 0.95)

##      95% 
## 27.38742

quantile(mc.closing, 0.05)

##       5% 
## 21.69453

Напомена: Ово је поједностављен модел кретања акција на берзи. Чак и модели који се генерално сматрају слабим моделима за цене акција користе бар логнормалну расподелу… У реалним сценаријима, квантитативне финансије интензивно користе Монте Карло симулације.

Games of chance

У игри “Games of chance” играч осваја 1 бод ако се стрелица заустави на жутој боји, губи 1 бод ако се заустави на црвеној и добија 2 бода ако се заустави на плавој боји.

Лако налазимо очекивани добитак играча: \[E(X) = \frac{1}{2}\cdot 1 + \frac{1}{4}\cdot (-1) + \frac{1}{4}\cdot 2 = 0.75.\]

Колика је вероватноћа да играч има мање од 0 поена након што казаљку заврти 10 пута?

runs <- 100000

#simulates on game of 10 spins, returns whether the sum of all the spins is < 1
play.game <- function(){
  results <- sample(c(1,1,-1,2), 10, replace=T)
  return(sum(results) < 0)
}

mc.prob <- sum(replicate(runs,play.game()))/runs
mc.prob

## [1] 0.01388

ViS лото игра

Да би освојио одређену лутрију, играч мора саставити реч “ViS”. \(60\%\) листића садржи слово “V”, \(30\%\) садржи слово “i”, а \(10\%\) садржи слово “S”. Одредити просечан број листића које особа мора купити да би освојила награду. (30 пута поновити експеримент у циљу доношења закључка)

Свако слово има дату вероватноћу да се нађе на листићу и то редом:

• “V” -> 0.6,
• “i” -> 0.3,
• “S” -> 0.1.

Успоставимо везу између датих вероватноћа и (псеудо)случајних бројева. Сваком од слова доделићемо бројеве који ће бити резултат случајног генератора и то:

• “V” -> 0,1,2,3,4,5,
• “i” -> 6,7,8,
• “S” -> 9.

Дефинишимо функцију која симулира једну реализацију лото игре.

# Create a function to get result
Random_Number_List <- c() # Create an empty vector to append the result

Lottery <- function(){     
  V <- FALSE # Set default False for all outcome
  i <- FALSE # Set default False for all outcome
  S <- FALSE # Set default False for all outcome
  Total_Pick <- 0 # Number of Pick is 0
  
  while(!V | !i | !S){ # Only stop when temporary values of V, i, S holing TRUE at the same time
    Total_Pick <- Total_Pick + 1 # Count the number of Pick
    Random_Number <- sample(0:9,1) # Get random numbers from 0 to 9 with 1 unit space
    
    Random_Number_List <- append(Random_Number_List, Random_Number) #  Record the random number
    Random_Number_List <- paste(Random_Number_List, collapse = ", ") # Convert to character to illustrate the consequence of number
    
    if(Random_Number<=5){V <- TRUE} # Assign value to random number for V
    else if(Random_Number<=8){i <- TRUE} # Assign value to random number for i
    else {S <- TRUE} # Assign a value to random number for S
    
    Result <- data.frame(Random_Number_List=Random_Number_List,Total_Pick=Total_Pick) # Data frame the result
    }
  return(Result)
}

Lottery()

##          Random_Number_List Total_Pick
## 1 5, 4, 8, 4, 8, 4, 0, 0, 9          9

df <- data.frame(Random_Number_List = character(),
                 Total_Pick = integer(),
                 stringsAsFactors = FALSE) # Empty data frame to record result of each trial
Experiment <- function(n){                  # Function with input n
  trial <- c(1:n)                           # Vector of trial from 1 to n
    for( k in trial){                     
      if(k<=n){                              # Repeat calling function n times
        Random_Number_List <- c()           # Clear the list from last record
        Result <- Lottery()                 # Run Lottery function
      }
    
    df <- rbind(df,Result)                 # Append result for each trial
    Sum_Trial <- sum(df$Total_Pick)       #  Sum of picks of total trial
    Average_Pick <- Sum_Trial/n           # Average pick
    }

  return( list( data = df,                   # Create a list of return
              Sum_Trial = Sum_Trial,
              Average_Pick = Average_Pick))
}

exp1 <- Experiment(30)
exp1$Average_Pick

## [1] 11.6

exp2 <- Experiment(1000)
exp2$Average_Pick

## [1] 11.009

exp3 <- Experiment(3000)
exp3$Average_Pick

## [1] 11.134

Game Mathematician Test - “Pick A Prize” Game

Game rules: Players are presented with num_options hidden options. One of the hidden options is “collect”, and the others are credit prizes. Players will pick (and that way, reveal) options one by one until they pick the “collect” option. When they pick the “collect” option, all revealed credit prizes are awarded to the player and the game ends. Every credit prize is predetermined before the first pick and is chosen independently with a new random draw. Whenever a player picks a credit prize, and there are more hidden credit prizes on the screen (at least 2 hidden options left), there is a chance for a “bonus” feature. Othervise, the “bonus” feature is not available. The “bonus” feature reveals one more hidden credit prize automatically. After the “bonus” credit prize, the players continue picking, and the “bonus” feature is still available in the following picks.

Inputs you will need:

• num_options – number of hidden options. Example: num_options = 15
• award_values – possible award values. Example: award_values = [10, 20, 50, 100]
• award_values_probabilities – probabilities for choosing a value from award_values. Example: [0.5, 0.25, 0.2, 0.05]
• prob_bonus – the odds of activating a bonus award after a credit prize pick. Example: prob_bonus = 0.5

1. Implement a function single_game() that simulates a single game. This function should return the number of revealed awards and the awarded credit prize in a single game play.

2. Call a function simulate(num_games) which plays single_game(), num_games amount of times. It should return the average number of revealed awards and the average awarded credit prize. What are the results for the example inputs?

single_game <- function(num_options, award_values, award_values_probabilities, prob_bonus) {
  
  revealed_awards <- 0
  total_credit_award <- 0
  
  hidden_options <- c("collect", sample(award_values, num_options - 1, replace = TRUE, prob = award_values_probabilities))
  hidden_options <- sample(hidden_options) # Shuffle
  
  while (length(hidden_options) > 0) {
    
    # Player picks an option
    index <- sample(length(hidden_options), 1)
    pick <- hidden_options[index]
    hidden_options <- hidden_options[-index]
    
    revealed_awards <- revealed_awards + 1
    
    # If option is "collect", the game ends
    if (pick == "collect") {
      break
    }
    
    total_credit_award <- total_credit_award + as.numeric(pick)
    
    # Check if bonus can activate
    if (length(hidden_options) >= 2 && runif(1) < prob_bonus) {
      
      # Bonus reveals one more hidden credit prize, it cannot reveal collect
      hidden_options_bonus <- hidden_options[hidden_options != "collect"]
      
      if (length(hidden_options_bonus) > 0) {
        index <- sample(length(hidden_options_bonus), 1)
        bonus_pick <- hidden_options_bonus[index]

        revealed_awards <- revealed_awards + 1
        total_credit_award <- total_credit_award + as.numeric(bonus_pick)
      }
    }
  }
  
  return(list(revealed_awards = revealed_awards, total_credit_award = total_credit_award))
}

single_game(15, c(10, 20, 50, 100), c(0.5, 0.25, 0.2, 0.05), 0.5)

## $revealed_awards
## [1] 4
## 
## $total_credit_award
## [1] 120

simulate <- function(num_games) {
  
  num_options <- 15
  award_values <- c(10, 20, 50, 100)
  award_values_probabilities <- c(0.5, 0.25, 0.2, 0.05)
  prob_bonus <- 0.5
  
  total_revealed_awards <- 0
  total_credit_prizes <- 0
  
  # Simulate single game num_games times
  for (i in 1:num_games) {
    game_result <- single_game(num_options, award_values, award_values_probabilities, prob_bonus)
    total_revealed_awards <- total_revealed_awards + game_result$revealed_awards
    total_credit_prizes <- total_credit_prizes + game_result$total_credit_award
  }
  
  avg_revealed_awards <- total_revealed_awards / num_games
  avg_credit_prize <- total_credit_prizes / num_games
  
  return(list(avg_revealed_awards = avg_revealed_awards, avg_credit_prize = avg_credit_prize))
}

simulate(10000)

## $avg_revealed_awards
## [1] 11.5929
## 
## $avg_credit_prize
## [1] 265.592