Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
可視化で理解する「負の二項分布」 http://d.hatena.ne.jp/hoxo_m/20151012/p1
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")
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")
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.