iamsyu

## 1.
#    x 2 3 4 ... 7  8 ... 12
# p(x) 1 2 3 ... 6  5 ...  1 / 36
generating_f1 = function(N){
  X = rep(0, N)
  p = c(1:6, 5:1) / sum(c(1:6, 5:1))
  p_cum = cumsum(p)
  values = 2:12
  u = runif(N)
  X[u <= p_cum[1]] = values[1]
  for (i in 2:length(p_cum))
    X[u <= p_cum[i] & u >= p_cum[i-1]] = values[i]
  X
}
(XX = generating_f1(1000))
(freq_table = table(XX))
(output = data.frame(points = names(freq_table),
    simu_prob = as.vector(freq_table)/1000,
    real_prob = c(1:6, 5:1) / sum(c(1:6, 5:1))))


## 2.
# a.
# f(x) / g(x) = 20*x*(1-x)^3
# d/dx (f/g) = 20*[(1-x)^3 - 3*x*(1-x)^2] = 0 => x = 1/4
# c = f(1/4)/g(1/4) = 20*(1/4)*(1-1/4)^3 = 135/64
# accept if U < f(x) /(g(x)*c) = 256/27*x*(1-x)^3
generating_f2 = function(N){
  X = rep(NA, N)
  f = function(x) 256/27*x*(1-x)^3
  while(any(is.na(X)))
  {
    locToGen = which(is.na(X))
    u1 = runif(length(locToGen))
    u2 = runif(length(locToGen))
    X[locToGen[u2 < f(u1)]] = u1[u2 < f(u1)]
  }
  X
}
qqplot(generating_f2(1000), rbeta(1000, 2, 4))


## 3.
#    x    0    1    2    3    4
# p(x)  0.1  0.2  0.2  0.2  0.3
generating_f3 = function(N){
  X = rep(0, N)
  p = c(1,2,2,2,3) / sum(c(1,2,2,2,3))
  p_cum = cumsum(p)
  values = 0:4
  u = runif(N)
  X[u <= p_cum[1]] = values[1]
  for (i in 2:length(p_cum))
    X[u <= p_cum[i] & u >= p_cum[i-1]] = values[i]
  X
}
(XX = generating_f3(1000))
(freq_table = table(XX))
(output = data.frame(points = names(freq_table),
    simu_prob = as.vector(freq_table)/1000,
    real_prob = c(1,2,2,2,3) / sum(c(1,2,2,2,3))))

generating_f4 = function(N) sample(0:4, N, TRUE,
  c(1,2,2,2,3) / sum(c(1,2,2,2,3)))
(XX = generating_f4(1000))
(freq_table = table(XX))
(output = data.frame(points = names(freq_table),
    simu_prob = as.vector(freq_table)/1000,
    real_prob = c(1,2,2,2,3) / sum(c(1,2,2,2,3))))


## 4.
# a.
# f(x) / g(x) = 3*x^2
# d/dx (f/g) = 6*x = 0 => x = 0
# d^2/dx^2 (f/g) = 6 > 0 => f(0) / g(0) is minimum
# d/dx (f/g) > 0 for 0 <= x <= 1 => increasing
# c = f(1)/g(1) = 3
# accept if U < f(x) /(g(x)*c) = x^2
generating_f5 = function(N){
  X = rep(NA, N)
  f = function(x) x^2
  while(any(is.na(X)))
  {
    locToGen = which(is.na(X))
    u1 = runif(length(locToGen))
    u2 = runif(length(locToGen))
    X[locToGen[u2 < f(u1)]] = u1[u2 < f(u1)]
  }
  X
}
qqplot(generating_f5(1000), rbeta(1000, 3, 1))