Finite Dirichlet process with N(0,1)

dirichlet_process.r

diri <- function(alpha) {
    # Sample from the Dirichlet distribution with parameter (vector) alpha
    k <- length (alpha)
    Z <- rep(0,k)
    for(i in 1:k) {
        Z[i]  <- rgamma(n=1,shape=alpha[i], rate =1)
        S  <- sum(Z)
        P  <- Z/S
    }
    return(P)
}

draw_and_plot <- function(alpha, binning) {
    # Sample from the Dirichlet distribution and plot the "histogram"
    x <- diri(alpha)
    centers<-rep(0,length(binning$breaks)-1)
    for(i in 1:length(centers)) {
        centers[i] <- (binning$breaks[i+1] + binning$breaks[i])/2
    }
    plot(centers,x,type='b')
}

# Partition the real line
partition <- seq(-4,4,by=0.05)

# Evaluate the probability that we draw from these intervals (numerically)
repetitions <- 10000000
data <- rnorm(repetitions,0,1)
binning <- hist(data, breaks=c(min(data),partition,max(data)),plot=FALSE)
p0 <- rep(0,length(binning$counts))
for(i in seq(1,length(p0), by=1)) {
    p0[i] <- binning$counts[i] / repetitions
}


# For multiple measures alpha or equivalently for multiple n0, draw and plot 3 times
par(mfrow=c(4,3))
n0 <- 1000
alpha <- p0 * n0
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
n0 <- 100
alpha <- p0 * n0
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
n0 <- 10
alpha <- p0 * n0
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
n0 <- 1
alpha <- p0 * n0
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)
draw_and_plot(alpha, binning)

Produces a rough equivalent of this Wikipedia figure.