Acceptance-Rejection metod

Obavezno pročitati teorijski uvod u materijalima (ideja, algoritam, dokaz korektnosti metode i detaljno izvođenje primjera)

Kodovi (za primjere iz materijala)

# Napomena: U ovim primjerima generisacemo vektore iz trazene raspodjele A-R metodom a zatim uporediti histogram dobijenog uzorka sa teorijskom krivom gustine za odgovarajucu raspodjelu


# Beta raspodjela Beta(n+1,n+1)

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

## [1] 0.4984222

r <- replicate(1000, beta(3))

hist(r, probability = T, main = "")

curve(dbeta(x, 4, 4), add = T, lwd=2, col="blue")

# Normalna raspodjela N(0,1)

norm <- function() {
  while (1) {
    u1 <- runif(1)
    y <- -log(u1)
    u2 <- runif(1)
    
    if (u2 <= exp(-(y - 1) ^ 2 / 2)) {
      ifelse(runif(1) >= 0.5, return(y), return(-y))
      
    }
  }
}
r <- replicate(10000, norm())

hist(r, probability = T, main="")

curve(dnorm(x), add = T, col="orange", lwd=2)

# Drugi nacin, preko dvije eksponencijalne slucajne velicine


norm2 <- function() {
  while (1) {
    u1 <- runif(1)
    u2 <- runif(1)
    y1 <- -log(u1)
    y2 <- -log(u2)
    if (y2 >= (y1 - 1) ^ 2 / 2) {
      ifelse(runif(1) >= 0.5, return(y1), return(-y1))
      
    }
    
  }
  
  
  
}

r <- replicate(10000, norm2())

hist(r, probability = T, main = "")

curve(dnorm(x), add = T, lwd=2, col="red")

---

Još neke metode za generisanje normalne raspodjele

---

Box-Muller (više o ovoj metodi na sledećem času)

norm4 <- function() {
  u1 <- runif(1)
  u2 <- runif(1)
  z1 <- sqrt(-2 * log(u1)) * cos(2 * pi * u2)
  z2 <- sqrt(-2 * log(u1)) * sin(2 * pi * u2)
  return(c(z1, z2))
  
  
}

r<-replicate(10000,norm4())
hist(r[1,],probability = T, main = "")
curve(dnorm(x),add = T, col="blue", lwd=2)

hist(r[2,],probability = T, main="")
curve(dnorm(x),add = T, col="blue", lwd=2)

CGT i uniformna U(0,1)

Primjer

norm3 <- function(n) {
  u <- runif(n)
  x<-(sum(u)-n/2)/sqrt(n/12)
  return(x)
  
}

# Vec za n=12 se postizu dobri rezultati

r <- replicate(10000, norm3(12))

hist(r, probability = T, main = "")

curve(dnorm(x), add = T, lwd=2, col="green")