Моделовање дискретних случајних променљивих са коначно много вредности

Нека је \(X\) дискретна случајна величина са коначно много вредности и нека је познат закон расподеле:

\[P\{X=x_k\}=p_k, \quad k=1,2,...n.\] Нека је \(\gamma\) један (псеудо)случајни број.

Ако је:

• \(\gamma \leq p_1\), сматра се да се реализовала вредност \(x_1\),

• \(p_1 < \gamma \leq p_1+p_2\), сматра се да се реализовала вредност \(x_2\)

• \(\dots\)

• \(p_1+p_2+...+p_{n-1}< \gamma\), реализовала се вредност \(x_n\).

Значи да се за добијање једне реализације случајне величине користи један (псеудо)случајан број.

Моделовање случајних догађаја

На основу поступка за моделовање дискретних случајних величина могу се моделирати и реализације случајних догађаја. Ако је вероватноћа догађаја \(А\) једнака \(p\), тада је индикатор догађаја \(А\) случајна величина:

\[ I_A: \begin{pmatrix} 0 & 1 \\ 1-p & p \end{pmatrix} \]

Нека је \(\gamma\) један (псеудо)случајни број.

Ако је:

• \(\gamma \leq p\), реализовала се вредност \(1\) индикатора \(I_A\)

• \(p<\gamma\), реализовала се вредност 0 индикатора \(I_A\).

Моделовањем те случајне величине \(I_A\) добија се да се догађај \(А\) реализовао ако је добијена моделирана вредност \(1\), а да се догађај \(А\) није реализовао ако је добијена моделована вредност \(0\).

indikator <- function(p){
  gama<-runif(1)
  if(gama > p)
    x<-0
  if(gama <= p)
    x<-1
  return (x)
}

Моделовање биномне расподеле

Нека је вероватноћа реализације неког догађаја \(А\) у сваком експерименту једнака \(p\). Нека је случајна величина \(S_n\) једнака броју реализација догађаја \(А\) у \(n\) независних експеримената. Тада \(S_n\) има биномну расподелу \(\mathcal{B}(n,p)\). Њен закон расподеле је \[P\{S_n=k\}=\binom{n}{k}p^k(1-p)^{n-k}, \quad k=0,1,...,n.\]

На основу дефиниције, \(S_n\) се може представити као збир независних случајних величина \(I_k\) које све имају исту расподелу и представљају индикатор догађаја \(А\) у к–том експерименту, односно \[S_n = \sum_{k=0}^{n} I_k.\] Моделирају се вредности за \(I_k\), \(k=1,2,...,n\) и саберу се. Добијени број је реализована вредност случајне величине \(S_n\).

binomna <- function(n,p){
  x <- vector()
  for(i in 1:n) {
    x[i] <- indikator(p)
    }
  Sn <- sum(x)
  return (Sn)
}

Предност овог поступка је што се не морају рачунати вероватноће \(p_k\) из закона расподеле за \(S_n\), а недостатак је што је потребно \(n\) (псеудо) случајних бројева за добијање једне реализације случајне величине.

Моделовање дискретних случајних величина са пребројиво много вредности

Нека је \(X\) дискретна случајна величина са пребројиво много вредности и нека је познат закон расподеле: \[ X = \begin{pmatrix} x_1 & x_2 & \ldots & x_n & \ldots \\ p_1 & p_2 & \ldots & p_n & \ldots \end{pmatrix} \]

Нека је \(n\) природни број такав да је \(p_{n_0+1}+p_{n_0+2}+...< \delta\), где је \(\delta\) произвољно мали (унапред изабрани) позитиван реални број. Уместо случајне величине \(X\) посматра се засечена случајна величина:

\[ X^* = \begin{pmatrix} x_1 & x_2 & \ldots & x_{n_0} \\ p_1 & p_2 & \ldots & p_{n_0}^* \end{pmatrix} \]

где је \(p_{n_0}^* = 1- p_1-...-p_{n_0-1}\).

Вредност ове случајне величине \(X^*\) се моделира на начин који је већ описан за моделовање вредности дискретне случајне величине са коначно много вредности.

На тај начин се неће добити ни једна од вредности случајне величине \(X\) која је већа од \(x_n\).

Међутим, вероватноћа добијања било које вредности веће од \(x_n\) је мања од \(\delta\), а како је \(\delta\) унапред изабрани врло мали позитиван број, то су и вероватноће добијања вредности случајне величине \(X\) које су веће од \(x_n\) занемарљиво мале.

У зависности од задатка који се решава и броја вредности које је потребно моделирати бира се \(\delta\).

Моделовање геометријске расподеле

Нека је вероватноћа реализације неког догађаја \(А\) у сваком експерименту једнака \(p\) и нека се експерименти понављају (при истим условима) док се први пут не оствари догађај \(А\). Случајна величина \(X\) која представља број изведених експеримената, има геометријску расподелу. Њен закон расподеле је: \[ X = \begin{pmatrix} x_1 & x_2 & \ldots & x_n & \ldots\\ p & (1-p)p & \ldots & (1-p)^{n-1}p & \ldots \end{pmatrix} \]

Ова случајна величина се моделира тако што се, за изабрани позитивни број \(\delta\) и познату вероватноћу \(p\), одреди најмањи природни број који задовољава неједнакост: \[n_0>\frac{\ln \delta}{\ln (1-p)}.\] Формира се засечена случајна величина: \[ X^* = \begin{pmatrix} x_1 & x_2 & \ldots & x_{n_0} \\ p_1 & p_2 & \ldots & p_{n_0}^* \end{pmatrix} \]

где је \(p_{n_0}^* = 1-(p+(1-p)p+...+(1-p)^{n_0-2}p)\), а затим се моделира коришћењем поступка за моделовање дискретних случајних величина са коначно много вредности.

geometrijska <- function(p, delta){
  n0 <- trunc(log(delta)/log(1-p))+1
  X <- c(1:n0)
  P <- c()
  for(i in 1:n0)
    P[i]<-((1-p)^(i-1))*p
  gama<-runif(1)
  if(gama<=P[1])
    x<-X[1]
  pom<-P[1]
  for(i in 2:n0){
    if(gama>pom && gama<=pom+P[i])
      x<-X[i]
    pom<-pom+P[i]
    }
  return (x)
}

Моделовање Пуасонове расподеле

Случајна величина \(X\) има Пуасонову расподелу са параметром \(\lambda\), \(\lambda>0\) ако је њен закон расподеле облика \[p_k=P\{X=k\}=e^{-\lambda}\frac{\lambda^k}{k!}, \quad k \in \mathbf{N}.\]

Ова случајна величина се моделира тако што се, за изабрани позитивни број \(n_0\) и познати параметар \(\lambda\), одреди најмањи природни број који задовољава неједнакост: \[e^{-\lambda}\frac{\lambda^{n_0+1}}{(n_0+1)!}\frac{n_0+2}{n_0+2-\lambda}<\delta.\] Формира се засечена случајна величина: \[ X^* = \begin{pmatrix} x_1 & x_2 & \ldots & x_{n_0} \\ p_1 & p_2 & \ldots & p_{n_0}^* \end{pmatrix} \] где је \(p_k = e^{-\lambda}\frac{\lambda^k}{k!},\quad k=0,1,...,n_0-1, \quad p_{n_0}^*=1-\sum_{k=0}^{n_0-1} p_k\), а затим се моделира коришћењем поступка за моделовање дискретних случајних величина са коначно много вредности.

Моделовање непрекидних случајних величина - Метода инверзне функције

Теорема: Нека је дата случајна величина \(X\) чија је функција расподеле \(F\) строго монотона и непрекидна и нека је \(F^{-1}\) њена инверзна функција. Нека је \(Y\) случајна величина са униформном расподелом \(U(0, 1)\), тада случајна величина \(F^{-1}(Y)\) има функцију расподеле \(F\).

Дакле, реализована вредност \(x\) случајне величине \(X\) се добија помоћу једног (псеудо)случајног броја \(\lambda\), по формули \(x=F^{-1}(\lambda)\).

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

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

Моделовање експоненцијалне расподеле

Случајна величина \(X\) има експоненцијалну расподелу са параметром \(\lambda,\quad \lambda>0\) ако је њена функција расподеле облика: \[F(x)=1-e^{-\lambda x}, \quad x \geq 0.\] Ако је \(\gamma\) (псеудо)случајни број, тада је моделирана вредност случајне величине \(X\) једнака \(x =\frac{−1}{\lambda} \ln (1− \gamma)\).

Овај израз се може поједноставити коришћењем особине униформне расподеле: Ако случајна величина \(Y\) има униформну расподелу \(U(0, 1)\), тада и случајна величина \(1−Y\) има униформну расподелу \(U(0, 1)\).

На тај начин се добија једноставна формула за моделовање вредности експоненцијалне расподеле \[x=\frac{-1}{\lambda} \ln \gamma,\] где је \(\gamma\) (псеудо)случајни број.

r_exp <- function(n, lambda) {
  u <- runif(n)
  x <- - log(1 - u) / lambda
  return(x)
}

x <- r_exp(1000, 2)

par(mfrow = c(1,2))
hist(x, breaks = seq(0, 8, length.out = 30), prob = TRUE, xlim = c(0, 4))
curve(dexp(x, 2), 0, 4, col = "blue", lwd = 2, add = TRUE)
curve(pexp(x, 2), 0, 2, col = "blue", lwd = 2)
plot(ecdf(x), xlim = c(0, 2), add = TRUE, lwd = 2)

par(mfrow = c(1,1))

Моделовање Ерлангове (гама) расподеле

Случајна величина \(X\) има гама \(\gamma(n,1)\) расподелу ако је њена густина расподеле облика: \[f(x)= \frac{x^{n-1}e^{-x}}{\Gamma(n)}, \quad x >0.\]

Ова расподела се може представити као збир независних једнако расподељених случајних величина са експоненцијалном \(\epsilon(1)\) расподелом.

Ако имамо \(n\) (псеудо)случајних бројева \(\gamma_1,\gamma_2,...,\gamma_n\), тада се моделирана вредност случајне променљиве \(X\) добија по формули: \[x=- \ln (\gamma_1\cdot...\cdot\gamma_n).\]

gama_1 <- function(a){
  y<-runif(20)
  x <- -log(prod(y))
  return(x)
}

Erlang_random<-function(n,lambda){
  x <- 0
  for(i in 1:n){
    gama <- runif(1)
    x <- x + (-1/lambda)*log(1-gama)
  }
  return(x)
}
set.seed(1)
Erlang_random(10,2)

## [1] 5.778018

Моделовање Гумбелове расподеле

Случајна променљива \(X\) има Гумбелову (двоструку експоненцијалну) расподелу са параметрима \(u, u>0\) и \(a, a>0\) ако је њена функција расподеле \[F(x)=e^{-e^{-\frac{x-u}{a}}}, \quad x \in \mathbb{R}.\]

Може се применити метода инверзне функције.

Ако је \(\gamma\) (псеудо)случајни број, тада је моделована вредност случајне променљиве \(X\) једнака \(x=u-a \ln(-\ln \gamma)\).

gumbel <- function(a,b){
  u <- runif(1)
  x <- a-b*log(-log(u))
  return(x)
}

Моделовање Фрешеове расподеле

Фрешеова расподела са параметрима \(a \in \mathbb{R}\), \(b>0\) и \(k \geq 0\) одређена је функцијом расподеле \[F(x)=e^{-{(\frac{x-a}{b})}^k}, \quad x>a.\]

Примена методе инверзне функције даје \(x=a+b(- \ln \gamma )^{-\frac{1}{k}}\), где је \(\gamma\) (псеудо)случајни број.

freseova <- function(a,b,k){
  gama <- runif(1)
  x <- a+b*(-log(gama))^((-1)/k)
  return(x)
}

Kumaraswamy distribution

Функција расподеле је дата са \(F_X(x) = 1 - (1 - x^a)^b\), где је \(a,b\in\mathbb{R}\). Одавде је \[\begin{aligned} u &= F_X(x) \\ &= 1 - (1 - x^a)^b. \\ \end{aligned}\] \[\begin{aligned} &\implies 1 - u = (1 - x^a)^b \\ &\implies(1 - u)^{1/b} = 1 - x^a \\ &\implies1 - (1 - u)^{1/b} = x^a \\ &\implies x = (1 - (1 - u)^{1/b})^{1/a}. \end{aligned}\]

a <- 1/2
b <- 2

set.seed(1)
u <- runif(10000)

generated_vals <- (1 - (1 - u)^(1/b) ) ^(1/a)

library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra)

# get the exact values of the PDF
x <- seq(0.01, 1, length.out = 10000)
theoretical_pdf <- a * b * x^(a-1) * (1 - x^a)^(b-1)

# create list to hold our two plots
plts <- list()

# constuct histogram of generated values
plts[[1]] <- tibble(generated_vals = generated_vals) %>%
    ggplot() +
        geom_histogram(aes(x = generated_vals), bins = 30) +
        theme_bw() +
        labs(x = "x", y = "Frequency", title = "Generated Kumaraswamy density")

# construct line chart of exact pdf
plts[[2]] <- tibble(x = x, theoretical_pdf = theoretical_pdf) %>%
    ggplot() +
        geom_line(aes(x = x, y = theoretical_pdf, group = 1), colour = "blue") +
        theme_bw() +
        labs(x = "x", y = "f(x)", title = "Exact Kumaraswamy density")

grid.arrange(grobs = plts, nrow = 1)

Моделовање непрекидних случајних величина - Метода одбацивања функције (Нојманова метода)

Претпоставимо случајна величина \(X\) има густину расподеле \(f(x)\) која је ограничена на неком коначном интервалу.

Случајну променљиву \(X\) симулирамо на следећи начин:

• Генеришемо два случајна броја \(\gamma_1\) и \(\gamma_2\).

• Рачунамо \(x=A+(B-A)\gamma_1\) и \(y=M\gamma_2\).

Ако је \(y\le f(x)\), тада је \(x\) реализована вредност случајне величине \(X\).

Моделовање бета расподеле

Генерисање \(B(n+1,n+1)\) расподеле користећи методу одбацивања.

Функција густине је \[f(x)=\frac{x^n(1-x)^x}{B(n+1,n+1)}.\]

beta_sim <- function(n)
{
  while (1) {
    u2 <- runif(1)
    u1 <- runif(1)
    if (u2 <= u1^n*(1-u1)^n/beta(n+1,n+1))
      return(u1)
  }
}
# Ocekivanje za Beta(n+1,n+1) je 1/2 za svako n
mean(replicate(10000, beta_sim(3)))

## [1] 0.5028632

r <- replicate(1000, beta_sim(3))
hist(r, probability = T, main = "")
curve(dbeta(x, 4, 4), add = T, lwd=2, col="blue")

Моделовање нормалне расподеле

Постоје различите методе за моделовање нормалне расподеле, јер је овим методама посвећена посебна пажња, због значаја и честе примене нормалне расподеле.

Тврђење: Ако случајна величина \(X\) има нормалну расподелу \(N(m,\sigma^2)\), тада случајна величина \(Y = \frac{X-m}{\sigma}\) има стандардну нормалну расподелу \(N(0,1)\).

С обзиром на тврђење, довољно је навести поступке моделирања случајне величине која има стандардну нормалну расподелу \(N(0,1)\).

Моделовање нормалне расподеле на основу централне граничне теореме

Нека су дате независне случајне величине \(Y_1, Y_2, ..,Y_{n}\) са \(U(0,1)\) расподелом, тада случајна величина \(S_n=\sum\limits_{j=1}^{n} Y_j\) има математичко очекивање \(E(S_n)= \frac{n}{2}\) и дисперзију \(D(S_n)=\frac{n}{12}\).

Према централној граничној теореми: \(S_n^*=\frac{S_n-E(S_n)}{\sqrt{D(S_n)}} \sim N(0,1)\).

Конвергенција је брза и већ за \(n=12\) добијају се врло мала одступања, па ако су \(\gamma_1,...,\gamma_{12}\) (псеудо)случајни бројеви, може се сматрати да је \(\sum\limits_{i=1}^{12} \gamma_i - 6\) реализована вредност случајне величине са \(N(0,1)\) расподелом.

На овај начин за моделовање једне вредности случајне величине са стандардном нормалном расподелом \(N(0,1)\) је потребно \(12\) (псеудо) случајних бројева.

normalna_cgt_1 <- function (){
  gama<-c()
  for(i in 1:12)
    gama[i]<-runif(1)
  x<-sum(gama)-6
  return(x)
  }
normalna_cgt_2 <- function(n){
  Y<-c()
  for(i in 1:n)
    Y[i]<-runif(1)
  x<-sqrt(3/n)*sum(2*Y-1)
  return(x)
}

Моделовање нормалне расподеле Box-Muller-ов методом

Нормална нормирана расподела се може моделирати кориштењем Box-Muller-ове трансформације.

Тврђење: Ако су \(Y_1\) и \(Y_2\) независне случајне величине са униформном расподелом \(U(0,1)\), тада су случајне величине: \[Z_1=\sqrt{-2\ln Y_1}\cos(2\pi Y_2), \quad Z_2=\sqrt{-2\ln Y_1}\sin(2\pi Y_2)\] независне са \(N(0,1)\) расподелама.

Дакле, ако су \(\gamma_1\) и \(\gamma_2\) два (псеудо)случајна броја, помоћу њих се добијају реализације \(x_1\) и \(x_2\) две независне случајне величине са стандардном нормалном расподелом, по формулама: \[x_1=\sqrt{-2\ln \gamma_1}\cos(2\pi\gamma_2), \quad x_2=\sqrt{-2\ln \gamma_1}\sin(2\pi\gamma_2)\]

normalna_bm <- function (){
  x <- c()
  gama1<-runif(1)
  gama2<-runif(1)
  x1<-sqrt(-2*log(gama1))*cos(2*pi*gama2)
  x2<-sqrt(-2*log(gama1))*sin(2*pi*gama2)
  x[1] <- x1
  x[2] <- x2
  return(x)
}

Моделовање нормалне расподеле коришћењем поларних координата

Стандардна нормална расподела се може моделирати и преласком на поларне координате.

• Генеришемо два случајна броја \(u\) и \(v\) у интервалу \((-1,1)\) као \(u=2\gamma_1-1\) и \(v=2\gamma_2-1\)
• Рачунамо \(\omega=u^2+v^2\) и ако је \(\omega>1\) вратимо се на корак \(1\).
• Рачунамо \(x=uz\) и \(y=vz\), где је \(z=\sqrt{-2\frac{\ln \omega}{\omega}}\). Дакле, \(x\) и \(y\) су два независна случајна броја из стандардне нормалне расподеле.

Овај метод је бржи од претходног, јер смо избегли рачунање \(\sin\) и \(\cos\). Изрази \(\frac{u}{\sqrt{\omega}}\) и \(\frac{v}{\sqrt{\omega}}\) имају улогу \(\sin\) и \(\cos\) из претходне методе.

Моделовање нормалне расподеле коришћењем једне експоненцијалне расподеле

Нека су \(Y \sim U(0,1)\) и \(V \sim \epsilon(1)\) независне случајне величине. Ако је \(Y\leq e^{-\frac{(V-1)^2}{2}}\), тада се за реализацију случајне величине \(X\) са стандардном нормалном расподелом узима вредност случајне величине \(V\).

Помоћу \(\gamma_2\) се моделира вредност \(v=- \ln \gamma_2\) случајне величине \(V\). Ако је \(\gamma_1 \leq e^{-\frac{(v-1)^2}{2}}\), тада се сматра да је \(v\) реализована вредност случајне величине са стандардном нормалном расподелом. Ако неједнакост не важи, поступак се понавља са наредним паром (псеудо)случајних бројева.

Овај поступак за моделовање једне вредности стандарднe нормалне расподеле захтева коришћење бар два (псеудо)случајна броја \(\gamma_1\) и \(\gamma_2\).

normalna_1_eksp <- function(){
  gama1<-runif(1)
  gama2<-runif(1)
  v=-log(gama2)
  while (gama1 > exp(-((v-1)^2)/2)){
    gama1<-runif(1)
    gama2<-runif(1)
    v=-log(gama2)
    }
  return(v)
}

Моделовање нормалне расподеле коришћењем две експоненцијалне расподеле

Нека су случајне величине \(Y_1\) и \(Y_2\) независне са истом експоненцијалном расподелом \(\epsilon(1)\). Ако је \(Y_2 \geq \frac{(Y_1-1)^2}{2}\), тада се за реализацију случајне величине \(X\) са стандардном нормалном расподелом узима вредност случајне величине \(Y_1\).

Овај поступак за моделовање једне вредности стандарднe нормалнe расподелe захтева коришћење бар два (псеудо)случајна броја \(\gamma_1\) и \(\gamma_2\).

Помоћу \(\gamma_1\) и \(\gamma_2\) се моделираjу вредности \(y_1=− \ln \gamma_1\) и \(y_1=−\ln \gamma_1\) случајне величине \(Y_i\). Ако је \(y_2 \geq \frac{(y_1-1)^2}{2}\), тада се сматра да је \(y_1\) реализована вредност случајне величине \(X\) са стандардном нормалном расподелом. Ако неједнакост не важи, поступак се понавља са наредним паром (псеудо)случајних бројева.

normalna_2_eksp <- function(){
  gama1<-runif(1)
  gama2<-runif(1)
  v1<--log(gama1)
  v2<--log(gama2)
  while(v2<((v1-1)^2)/2){
    gama1<-runif(1)
    gama2<-runif(1)
    v1<--log(gama1)
    v2<--log(gama2)
    }
  return(v1)
}

Моделовање \(log\)–нормалне расподеле

Случајна променљива \(X\) има \(LN(m,\sigma^2)\) расподелу ако случајна променљива \(Y=\ln X\) има нормалну расподелу \(N(m,\sigma^2)\).

Моделовање вредности случајне променљиве \(X\) добијамо по формули \(x=e^y\), где је \(y\) моделована вредност случајне променљиве \(Y\).

log_normalna <- function(){
  y<-normalna_1_eksp() # z je modelovana vrednost normalne raspodele pomocu 1 eksponencijalne funkcije
  x<-exp(y)
  return(x)
}

Задатак за вежбу

1. Методом инверзне функције генерисати вредност случајне променљиве \(X\) чија је функција расподеле \(F(x)=1-e^{-\lambda\sqrt{x}},\; x \ge 0\).

2. Симулирајте \(B(2,4)\) расподелу користећи методу одбацивања.