dirichlet_process.R

library(magrittr)

F <- function(y, theta) {
  dnorm(y, mean = theta, sd = 1)
}

G0 <- function(theta) {
  dnorm(theta, mean = 0, sd = 1)
}

sample_from_distribution <- function(min, max, distribution, threshold) {
  while (TRUE) {
    x <- runif(n = 1, min = min, max = max)
    dice <- runif(n = 1, min = 0, max = 1)
    if (dice < distribution(x) / threshold) {
      return(x)
    }
  }
}

alpha <- 1
num_iterations <- 4096
min_theta <- -10
max_theta <-  10

y_array <- c(
  rnorm(64, mean = -7.5, sd = 1),
  rnorm(64, mean =    0, sd = 1),
  rnorm(64, mean =  7.5, sd = 1)
)
# y_array <- rnorm(64, mean = 1, sd = 0.1)
num_samples <- length(y_array)
theta_array <- rep(0, num_samples)

theta_samples <- matrix(nrow = num_iterations, ncol = num_samples)
  
for (iteration in 1:num_iterations) {
  for (y_index in 1:num_samples) {
    y_i <- y_array[y_index]
    
    # theta-wise array
    q_over_b <- F(y_i, theta_array[-y_index])
    # constant
    r_over_b <- 0.5 * alpha
    
    b <- (sum(q_over_b) + r_over_b)^(-1)
    r <- r_over_b * b
    q <- q_over_b * b
    
    dice <- runif(n = 1, min = 0, max = 1)
    if (dice < r) {
      # sample a new point from G0 * F
      distribution_function <- function(theta) G0(theta) * F(y_i, theta)
      # totally empirical
      threshold <- max(distribution_function(c(0, 10))) * 3
      new_theta <- sample_from_distribution(
        min = min_theta,
        max = max_theta,
        distribution = distribution_function,
        threshold = threshold)
    } else {
      # sample from existing theta with probability proportonial to q
      new_theta <- sample(size = 1, x = theta_array[-y_index], prob = q)
    }
    theta_array[y_index] <- new_theta
    theta_samples[iteration, y_index] <- new_theta
  }
}
    
null_axis <- rep(0, y_array %>% length)
plot(y_array, null_axis)
hist(y_array)
hist(theta_samples, probability = TRUE)

# histogram of unique elements
unique_elements <- apply(theta_samples, MARGIN = 1, FUN = function(x) x %>% unique %>% length)
hist(unique_elements)