Instantly share code, notes, and snippets.

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

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)
Owner Author

### jg-you commented Apr 27, 2015

 Produces a rough equivalent of this Wikipedia figure.
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.