## Function to sample Chinese restaurant process

crp = function(num.customers, alpha) {
  table <- c(1)
  next.table <- 2
  for (i in 1:num.customers) {
    if (runif(1,0,1) < alpha / (alpha + i)) {
      table <- c(table, next.table)
      next.table <- next.table+1
    } else {
      select.table <- table[sample(1:length(table), 1)]
      table <- c(table, select.table)
    }
  }
  table
}

## Visualize CRP
crp(100, 10)
plot(
  table( crp(10000, 2) )
  ,xlab="Table Index", ylab="Frequency"
)

## Code to Sample Stick breaking weights

stickBreakingModel = function(num.vals, alpha) {
  betas = rbeta(num.vals, 1, alpha)
  stick.to.right = c(1, cumprod(1 - betas))[1:num.vals]
  weights = stick.to.right * betas
  weights
}

## Visualize the weights
plot( stickBreakingModel(100,1), pch=16, cex=.5 )
plot(
  table( stickBreakingModel(100,5) )
  ,xlab="Table Index", ylab="Frequency"
)
barplot(stickBreakingModel(100,1))



alpha <- 2
G_0 <- function(n) rnorm(n, 0, 10)
n <- 100
nu <- rbeta(n, 1, alpha)
p <- numeric(n)
p[1] <- b[1]
p[2:n] <- sapply(2:n, function(i) nu[i] * prod(1 - nu[1:(i-1)]))
y <- G_0(n)
theta <- sample(y, prob = p, replace = TRUE)
plot(
  table( theta )
  ,xlab="Table Index", ylab="Frequency"
)