Skip to content

Instantly share code, notes, and snippets.

@jg-you jg-you/dirichlet_process.r
Last active Aug 29, 2015

Embed
What would you like to do?
Finite Dirichlet process with N(0,1)
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)
@jg-you

This comment has been minimized.

Copy link
Owner Author

commented Apr 27, 2015

Produces a rough equivalent of this Wikipedia figure.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.