Rejection sampling, part II

class: center, middle, inverse, title-slide

# Rejection sampling, part II
### Niels Richard Hansen
### September 15, 2020

---

## von Mises rejection sampling

* `$Y \sim \mathrm{unif}(-\pi, \pi)$`  (proposal)

* `$U \sim \mathrm{unif}(0, 1)$`

* Accept if `$U \leq e^{\kappa(\cos(Y) - 1)}$`

---
## Vectorized von Mises simulation

```r
vMsim_vec_ran <- function(m, kappa) {
  y <- runif(m, -pi, pi)
  u <- runif(m)
  accept <- u <= exp(kappa * (cos(y) - 1))
  y[accept] 
}
```

While the implementation above is simple and vectorized it returns a random 
number of samples. To return a fixed number of samples we can write a wrapper 
of the call to the function.

In fact, we can write a generic function factory that can convert any generator
of a random size output into a generator of a fixed size output.

---
## Vectorization function factory

```r
vec_sim <- function(generator) {
  alpha <- 1 
  function(n, ...) {
    y <- list()
    j <- 1
    l <- 0  # The number of accepted samples 
    while(l < n) {
      m <- ceiling((n - l) / alpha)  
*     y[[j]] <- generator(m, ...)
      l <- l + length(y[[j]])
      if (j == 1)
*       alpha <<- (l + 1) / (m + 1)  # Estimate of alpha
      j <- j + 1
    }
    unlist(y)[1:n]
  }
}
```

```r
vMsim_vec <- vec_sim(vMsim_vec_ran)
```

---
## Profiling

```r
profvis(vMsim_vec_ran(5000000, 5))  
```

---
## Profiling

```r
profvis(vMsim_vec(1000000, 5))  
```

---
## Gamma distribution

Can we find a suitable envelope?

Perhaps, but here rejection sampling of a non-standard distribution will 
be used in combination with a simple transformation.

Let `$t(y) = a(1 + by)^3$` for `$y \in (-b^{-1}, \infty)$`, then `$t(Y) \sim \Gamma(r,1)$` if
`$Y$` has density 
`$$f(y) \propto t(y)^{r-1}t'(y) e^{-t(y)} = e^{(r-1)\log t(y) + \log t'(y) - t(y)}.$$`

[Proof: Use univariate density transformation theorem.]

The density `$f$` will be the *target density* for a rejection sampler.

---
## Gamma distribution

With 
`$$f(y) \propto e^{(r-1)\log t(y) + \log t'(y) - t(y)},$$`
`$a = r - 1/3$` and `$b = 1/(3 \sqrt{a})$` 
`$$f(y) \propto e^{a \log t(y)/a - t(y) + a \log a} \propto \underbrace{e^{a \log t(y)/a - t(y) + a}}_{q(y)}.$$`

An analysis of `$w(y) := - y^2/2 - \log q(y)$` shows that it is convex on `$(-b^{-1}, \infty)$` 
and it attains its minimum in 0 with `$w(0) = 0$`, whence 
`$$q(y) \leq e^{-y^2/2}.$$`

---
## gamma rejection sampling

* `$Y \sim \mathcal{N}(0,1)$`  (proposal)

* `$U \sim \mathrm{unif}(0, 1)$`

* Accept if `$U \leq q(Y) e^{Y^2/2}$`

* If accept, return `$t(Y)$`

---
## Gamma distribution implementation

```r
# r >= 1 
tfun <- function(y, a) {
  b <- 1 / (3 * sqrt(a))
  (y > -1/b) * a * (1 + b * y)^3  # 0 when y <= -1/b
}
```

```r
qfun <- function(y, r) {
  a <- r - 1/3
  tval <- tfun(y, a)
  exp(a * log(tval / a) - tval + a)
}
```

---
## Gamma distribution implementation

```r
gammasim_ran <- function(m, r) {
  y <- rnorm(m)
  u <- runif(m)
  accept <- u <= qfun(y, r) * exp(y^2/2) 
  tfun(y[accept], r - 1/3)
}
```

```r
gammasim_vec <- vec_sim(gammasim_ran)
```

---
## Gamma distribution test

```r
tmp <- gammasim_vec(10000, 8)
hist(tmp, freq = FALSE)
curve(dgamma(x, 8), col = "blue", lwd = 2, add = TRUE)
```

---
## Gamma distribution test

```r
tmp <- gammasim_vec(10000, 3)
hist(tmp, freq = FALSE)
curve(dgamma(x, 3), col = "blue", lwd = 2, add = TRUE)
```

---
## Rejection probabilities

```r
tmp <- gammasim_vec(1e5, 16)
1 - environment(gammasim_vec)$alpha
```

```
[1] 0.001870583
```

```r
tmp <- gammasim_vec(1e5, 8)
1 - environment(gammasim_vec)$alpha
```

```
[1] 0.003413548
```

```r
tmp <- gammasim_vec(1e5, 4)
1 - environment(gammasim_vec)$alpha
```

```
[1] 0.008151957
```

```r
tmp <- gammasim_vec(1e5, 1)
1 - environment(gammasim_vec)$alpha
```

```
[1] 0.04782639
```

---
## Adaptive envelopes

When `$f$` is *log-concave* on `$I$` we can construct bounds of the form 
`$$f(x) \leq e^{V(x)}$$`

where 
`$$V(x) = \sum_{i=1}^m  (a_i x + b_i) 1_{I_i}(x),$$`
for intervals `$I_i$` forming a partition of `$I$`.

Typically, `$a_i x + b_i$` is tangent to the graph of `$\log(f)$` at 
`$x_i \in I_i = (z_{i-1}, z_i]$` for 
`$$z_0 < x_1 < z_1 < x_2 < \ldots < z_{m-1} < x_m < z_m.$$`

---
## Beta distribution

```r
Betasim_ran <- function(m, x1, x2, alpha, beta) {
  lf <- function(x) (alpha - 1) * log(x) + (beta - 1) * log(1 - x) 
  lf_deriv <- function(x) (alpha - 1)/x - (beta - 1)/(1 - x)
  a1 <- lf_deriv(x1)
  a2 <- lf_deriv(x2)
  if(a1 == 0 || a2 == 0 || a1 - a2 == 0) 
    stop("\nThe implementation requires a_1 and a_2 different and both different from zero. Choose different values of x_1 and x_2.")
  b1 <- lf(x1) - a1 * x1
  b2 <- lf(x2) - a2 * x2
  z1 <- (b2 - b1) / (a1 - a2)
  Q1 <- exp(b1) * (exp(a1 * z1) - 1) / a1 
  c = Q1 + exp(b2) * (exp(a2 * 1) - exp(a2 * z1)) / a2
  
  y <- ratio <- numeric(m)
  uy <- c * runif(m)
  u <- runif(m)
  i <- uy < Q1
  y[i] <- z <- log(a1 * exp(-b1) * uy[i] + 1) / a1
  y[!i] <- zz <- log(a2 * exp(-b2) * (uy[!i] - Q1) + exp(a2 * z1)) / a2
  ratio[i] <- exp(lf(z) - a1 * z - b1)
  ratio[!i] <- exp(lf(zz) - a2 * zz - b2)
  accept <- u <= ratio
  y[accept]
}
```

---
## Beta distribution test

```r
Betasim_vec <- vec_sim(Betasim_ran)
```