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

@jg-you jg-you 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.