Skip to content

Instantly share code, notes, and snippets.

@jg-you
Last active August 29, 2015 14:20
Show Gist options
  • Save jg-you/a1973f16cc02f7e47ee9 to your computer and use it in GitHub Desktop.
Save jg-you/a1973f16cc02f7e47ee9 to your computer and use it in GitHub Desktop.
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
Copy link
Author

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