可視化で理解する「負の二項分布」 http://d.hatena.ne.jp/hoxo_m/20151012/p1

negative_binomial_distribution1.R

library(animation)
library(ggplot2)
library(gridExtra)

xlim <- xlim(-0.5, 20.5)

saveGIF({
  g_nbiom <- qplot(x=NULL, y=NULL) + 
    geom_bar(stat="identity") + xlim + ylim(0, 1) +
    theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
  for(i in 1:16) {
    g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + 
      geom_bar(stat="identity") + 
      theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
    g_gamma <- ggplot(data=data.frame(X=c(0,20)), aes(x=X)) + 
      stat_function(fun=dgamma, args = list(shape=4, rate=0.3/(1-0.3)), geom="density", fill="black") + 
      geom_vline(aes(xintercept = gamma_values[i], color="red"), size=2) + 
      theme(axis.text.y=element_blank()) + xlim + xlab("ガンマ分布") + ylab("")
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
    if(poisson_values[i] < 20) {
      g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + 
        geom_bar(stat="identity") + geom_vline(aes(xintercept = poisson_values[i], color="red"), size=2) + 
        theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
    }
    table <- table(poisson_values[1:i])
    g_nbiom <- qplot(x=as.integer(names(table)), y=unclass(table)) + 
      geom_bar(stat="identity", width=0.9) + xlim + 
      theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
  }
  for(i in 17:60) {
    if(poisson_values[i] < 20) {
      g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + 
        geom_bar(stat="identity") + geom_vline(aes(xintercept = poisson_values[i], color="red"), size=2) + 
        theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
      table <- table(poisson_values[1:i])
      g_nbiom <- qplot(x=as.integer(names(table)), y=unclass(table)) + 
        geom_bar(stat="identity") + xlim + 
        theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
      g_gamma <- ggplot(data=data.frame(X=c(0,20)), aes(x=X)) + 
        stat_function(fun=dgamma, args = list(shape=4, rate=0.3/(1-0.3)), geom="density", fill="black") + 
        geom_vline(aes(xintercept = gamma_values[i], color="red"), size=2) + 
        theme(axis.text.y=element_blank()) + xlim + xlab("ガンマ分布") + ylab("")
    }
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
  }
  for(i in seq(100, 4000, by=80)) {
    g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + 
      geom_bar(stat="identity") + geom_vline(aes(xintercept = poisson_values[i], color="red"), size=2) + 
      theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
    table <- table(poisson_values[1:i])
    g_nbiom <- qplot(x=as.integer(names(table)), y=unclass(table)) + 
      geom_bar(stat="identity") + xlim + 
      theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
    g_gamma <- ggplot(data=data.frame(X=c(0,20)), aes(x=X)) + 
      stat_function(fun=dgamma, args = list(shape=4, rate=0.3/(1-0.3)), geom="density", fill="black") + 
      geom_vline(aes(xintercept = gamma_values[i], color="red"), size=2) + 
      theme(axis.text.y=element_blank()) + xlim + xlab("ガンマ分布") + ylab("")
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
  }
}, interval=rep(c(0.5, 0.2, 0.1, 3), c(32, 44, 48, 1)), movie.name="nbiom_animation3.gif")

negative_binomial_distribution2.R

library(animation)
library(ggplot2)

saveGIF({
  for(i in 1:100) {
    df <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:i]))
    g <- ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none", axis.text.y=element_blank()) + ggtitle(sprintf("n = %d", i))
    print(g)
  }
}, interval= c((1-(1:99/100)^2)/2, 3), movie.name="nbiom_animation2.gif")