# Stochasticke simulacne metody
# Radoslav Harman, KAMS FMFI UK
# Cvicenie 7 (Metropolisov-Hastingsov algoritmus) - strucny zaznam prikazov pre R


metro.disc <- function(p, x1, n, b=0, GraphA=FALSE, GraphH=FALSE) {

  # Metropolisov-Hastingsov algoritmus pre rozdelenia na konecnej mnozine
  # Independence sampler (efektivny len ak je cielove rozdelenie podobne kandidatskemu rozdeleniu)
  # p ... cielove rozdelenie ako vektor pravdepodobnosti (nemusi byt normalizovane)
  # x1 ... startovaci bod
  # n ... pozadovany pocet realizacii
  # b ... dlzka fazy burn-in
  # GraphA ... ci sa ma zobrazit ilustracny obrazok k priebehu algoritmu
  # GraphH ... ci ma zobrazit histogram vysledkov

  res <- rep(0, b + n) # vektor realizacii stavov
  can <- rep(0, b + n) # vektor kandidatov na dalsi krok (len na ilustracne ucely)
  v <- length(p)       # pocet stavov
  x <- x1              # momentalny stav

  for(i in 1:(b + n)) {
    res[i] <- x
    y <- ceiling(v * runif(1))  # Kandidatske rozdelenie nezavisi od x a je rovnomerne
    can[i] <- y
    alpha <- min(c(1, p[y] / p[x]))
    if(runif(1) < alpha) x <- y
  }

  can <- c(x1, can)
  resx <- res[(b + 1):(b + n)]
  resy <- can[(b + 1):(b + n)]

  if(GraphA) {
    plot(resy, pch=19, cex=1.5)
    points(resx, pch=19, col="red")
    if(GraphH) windows()
  }

  if(GraphH) hist(resx, breaks=(1:(v + 1)) - 0.5)

  return(list(x=resx, y=resy))
}

metro.disc(20:1, 1, 1000, 100, TRUE, TRUE)
res <- metro.disc(20:1, 1, 20000, 10000, FALSE, TRUE)
chisq.test(table(res$x), p=(20:1) / sum(20:1))


metro2d.cont <- function(f, x0, n, b=100) {

  # Normal random-walk Metropolis pre 2d hustoty
  # f je cielova hustota na R^2
  # x0 je startovaci bod
  # n je pozadovany pocet bodov
  # b je dlzka burn-in fazy

  x <- x0
  res <- matrix(ncol=2, nrow=n)
  for(i in 1:(b + n)) {
    y <- x + rnorm(2)
    alpha <- min(c(1, f(y) / f(x)))
    if(runif(1) < alpha) x <- y
    if(i > b) res[i - b,] <- x
  }

  return(res)
}

fn <- function(x) {

  # Hustota 2d normalneho rozdelenia s nulovou strednou hodnotou
  # a zadanou kovariancnou maticou (bez normalizacneho clena)

  Sigma <- matrix(c(2, 1, 1, 2), ncol=2)
  exp(-t(x) %*% solve(Sigma) %*% x / 2)
}
 
res <- metro2d.cont(fn, c(0, 0), 10000)
plot(res, pch=19, cex=0.5, asp=1)
 
fu <- function(x) {

  # Funkcia umerna hustote rovnomerneho rozdelenia na neznamej mnozine

  res <- 0; if(x[1]^2 + (x[2] - abs(x[1])^(2/3))^2 < 1) res <- 1
  return(res)
}
 
res <- metro2d.cont(fu, c(0, 0), 30000, b=10000)
plot(res, pch=19, cex=0.5, asp=1)