Skip to content

Instantly share code, notes, and snippets.

@matsuken92
Last active December 5, 2015 00:50
Show Gist options
  • Save matsuken92/671ae282429d6ad97aec to your computer and use it in GitHub Desktop.
Save matsuken92/671ae282429d6ad97aec to your computer and use it in GitHub Desktop.
# Preparing
install.packages("Rlab")
install.packages("ggplot2")
install.packages("actuar")
install.packages("plyr")
require(plyr)
require(actuar)
require(ggplot2)
require(Rlab)
###############################################################
# ベルヌーイ分布
###############################################################
# ベルヌーイ分布からのサンプリングを実行
# パラメーター
p = 0.7
trial_size = 10000
set.seed(71)
# ベルヌーイ分布に従う乱数生成
data <- rbern(trial_size, p)
# ベルヌーイ分布の確率分布
dens <- data.frame(y=c((1-p),p)*trial_size, x=c(0, 1))
# グラフ描画
ggplot() +
layer(data=data.frame(x=data), mapping=aes(x=x), geom="bar",
stat="bin", bandwidth=0.1
) + layer(data=dens, mapping=aes(x=x, y=y), geom="bar",
stat="identity", width=0.05, fill="#777799", alpha=0.7)
###############################################################
# 二項分布
###############################################################
# ベルヌーイ分布に従う変数をsample_size回行ったときの1の回数を
# trial_size分繰り返してベクトルに入れていく。
# パラメーター
p = 0.7
trial_size = 10000
sample_size = 30
set.seed(71)
# ベルヌーイ分布に従う乱数生成
gen_binom_var <- function() {
return(sum(rbern(sample_size, p)))
}
result <- rdply(trial_size, gen_binom_var())
# 二項分布の密度関数
dens <- data.frame(y=dbinom(seq(sample_size), sample_size, 0.7))*trial_size
# グラフ描画
ggplot() +
layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=1, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(sample_size)+.5, y=y),
geom="line", stat="identity", position="identity",colour="red"
) + ggtitle("Bernoulli to Binomial.")
set.seed(73)
# パラメーター
n <- 10000
p <- 0.7
trial_size = 30000
width=0.18
# ベルヌーイ分布に従う乱数生成
gen_binom_var <- function() {
return(sum(rbern(n, p)))
}
result <- rdply(trial_size, gen_binom_var())
m <- mean(result$V1)
sd <- sd(result$V1)
result <- (result - m)/sd
# 標準正規分布の密度関数
dens <- data.frame(y=dnorm(seq(-4,4,0.05), mean=0,
sd=1)*trial_size*width)
# グラフ描画
ggplot() +
layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(-4,4,0.05), y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Standard Normal.")
###############################################################
# 正規分布
###############################################################
# 二項分布の正規近似 (nを大きくする)
set.seed(72)
# パラメーター
n <- 10000
p <- 0.7
trial_size = 10000
width=10
# ベルヌーイ分布に従う乱数生成
gen_binom_var <- function() {
return(sum(rbern(n, p)))
}
result <- rdply(trial_size, gen_binom_var())
# 正規分布の密度関数
dens <- data.frame(y=dnorm(seq(6800,7200), mean=n*p,
sd=sqrt(n*p*(1-p)))*trial_size*width)
# グラフ描画
ggplot() +
layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(6800,7200), y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Normal.")
###############################################################
# カイ二乗分布 (自由度:1)
###############################################################
# 正規分布に従う乱数を二乗してヒストグラムにする。
# パラメーター
sample_size <- 1000
trial_size <- 100000
width <- 0.3
set.seed(71)
# ベルヌーイ分布に従う乱数生成
res <- c()
for (i in sequence(trial_size)) {
res <- rbind(res, sum(rbern(sample_size, p)))
}
res <- (res - mean(res))/sd(res) # 標準化
data <- data.frame(x=(res**2))
# カイ二乗分布の密度関数(自由度=1)
xx <- seq(0,20,0.1)
dens <- data.frame(y=dchisq(x=xx, df=1)*trial_size*width)
length(seq(sample_size))
# グラフ描画
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=xx, y=y),
geom="line", stat="identity", position="identity", colour="red"
)
###############################################################
# カイ二乗分布 (自由度:3)
###############################################################
# 正規分布に従う乱数を二乗してヒストグラムにする。
set.seed(71)
# パラメーター
p <- 0.7
n <- 1000
trial_size <- 100000
width <- 0.3
df <- 3
# ベルヌーイ分布に従う乱数生成(3まわし)
gen_binom_var <- function() {
return(sum(rbern(n, p)))
}
gen_chisq_var <- function() {
result <- rdply(trial_size, gen_binom_var())
return(((result$V1 - mean(result$V1))/sd(result$V1))**2)
}
# 自由度dfの分だけ生成する
result <- rlply(df, gen_chisq_var(),.progress = "text")
res <- result[[1]] + result[[2]] + result[[3]]
data <- data.frame(x=res)
# カイ二乗分布の密度関数(自由度=3)
xx <- seq(0,20,0.1)
dens <- data.frame(y=dchisq(x=xx, df=df)*trial_size*width)
# グラフ描画
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=xx, y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Chisquare")
###############################################################
# ポアソン分布
###############################################################
trial_size = 5000; width <- 1;
# もともとの問題設定
p = 0.7; n = 30;
np <- p*n
cat("p=",p, ", n=", n, "np=", np)
# n→∞、p→0、np=一定
n = 100000; p <- np/n
cat("p=",p, ", n=", n)
# ベルヌーイ分布に従う乱数生成
gen_binom_var <- function() {
return(sum(rbern(n, p)))
}
result <- rdply(trial_size, gen_binom_var())
# ポアソン分布の密度関数
dens <- data.frame(y=dpois(seq(20), np))*trial_size
# グラフ描画
ggplot() +
layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(20)+.5, y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Poisson.")
###############################################################
# ポアソン分布
###############################################################
trial_size = 100000
width <- 1
# もともとの問題設定
p = 0.7
sample_size = 30
np <- p*sample_size
cat("p=",p, ", n=", sample_size, "np=", np)
# n→∞、p→0、np=一定
sample_size = 10000
p <- np/sample_size
cat("p=",p, ", n=", sample_size)
# ベルヌーイ分布に従う乱数生成
gen_rand <- function(){
res <- c()
for (i in sequence(trial_size)) {
res <- rbind(res, sum(rbern(sample_size, p)))
}
return(res)
}
system.time(res <- gen_rand())
data <- data.frame(x=res)
# 正規分布の密度関数
dens <- data.frame(y=dpois(seq(30), np))*trial_size
# グラフ描画
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(30)+.5, y=y),
geom="line", stat="identity", position="identity", colour="blue"
)
########################
trial_size = 100000
width <- 1
# もともとの問題設定
p = 0.7
sample_size = 30
np <- p*sample_size
cat("p=",p, ", n=", sample_size, "np=", np)
# n→∞、p→0、np=一定
sample_size = 10000
p <- np/sample_size
cat("p=",p, ", n=", sample_size)
# ベルヌーイ分布に従う乱数生成
gen_rand <- function(){
res <- c()
for (i in sequence(trial_size)) {
res <- rbind(res, sum(rbern(sample_size, p)))
}
return(res)
}
system.time(res <- gen_rand())
data <- data.frame(x=res)
# ポアソン分布の密度関数
dens <- data.frame(y=dpois(seq(20), np))*trial_size
# グラフ描画
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(20)+.5, y=y),
geom="line", stat="identity", position="identity", colour="blue"
)
###############################################################
# 指数分布
###############################################################
trial_size = 7000; width <- .01;
# もともとの問題設定
p = 0.7; n = 10;
np <- p*n
cat("p=",p, ", n=", n, "np=", np)
p
# n→∞、p→0、np=一定
n = 10000; p <- np/n
cat("p=",p, ", n=", n)
# ベルヌーイ分布に従う乱数生成
gen_exp_var <- function() {
cnt <- 0
while (TRUE) {
cnt <- cnt + 1
if (rbern(1, p)==1){
return(cnt) # 1が出たら何回目かを返す
}
}
}
data <- data.frame(x=rdply(trial_size, gen_exp_var())/n)
names(data) <- c("n", "x")
# 指数分布の密度関数
dens <- data.frame(y=dexp(seq(0, 1.5, 0.1), np)*trial_size*width)
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(0, 1.5, 0.1), y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Exponential.")
###############################################################
# ガンマ分布
###############################################################
trial_size = 7000; width <- .035;
# もともとの問題設定
p = 0.7; n = 10;
np <- p*n
cat("p=",p, ", n=", n, "np=", np)
# n→∞、p→0、np=一定
n = 10000; p <- np/n
cat("p=",p, ", n=", n)
alpha <- 5
# ベルヌーイ分布に従う乱数生成
get_interval <- function(){
cnt <- 0
while (TRUE) {
cnt <- cnt + 1
if (rbern(1, p)==1){
return(cnt) # 1が出たら何回目かを返す
}
}
}
gen_exp_var <- function() {
data <- data.frame(x=rdply(trial_size, get_interval())/n)
names(data) <- c("n", "x")
return(data)
}
result <- rlply(alpha, gen_exp_var())
data <- data.frame(x=result[[1]]$x + result[[2]]$x + result[[3]]$x + result[[4]]$x + result[[5]]$x)
# ガンマ分布の密度関数
dens <- data.frame(y=dgamma(seq(0, 3,.01), shape=alpha, rate=np)*trial_size*width)
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(0,3,.01), y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Gamma")
###############################################################
# 逆ガンマ分布
###############################################################
trial_size = 7000; width <- .1;
# もともとの問題設定
p = 0.7; n = 10;
np <- p*n
cat("p=",p, ", n=", n, "np=", np)
# n→∞、p→0、np=一定
n = 10000; p <- np/n
cat("p=",p, ", n=", n)
alpha <- 5
# ベルヌーイ分布に従う乱数生成
get_interval <- function(){
cnt <- 0
while (TRUE) {
cnt <- cnt + 1
if (rbern(1, p)==1){
return(cnt) # 1が出たら何回目かを返す
}
}
}
gen_exp_var <- function() {
data <- data.frame(x=rdply(trial_size, get_interval())/n)
names(data) <- c("n", "x")
return(data)
}
result <- rlply(alpha, gen_exp_var())
data <- data.frame(x=1/(result[[1]]$x + result[[2]]$x +
result[[3]]$x + result[[4]]$x + result[[5]]$x))
# ガンマ分布の密度関数
dens <- data.frame(y=dinvgamma(seq(0, 23,.01), shape=5, rate=1/np)*trial_size*width)
ggplot() +
layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(0,23,.01), y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Inversegamma")
a <- 5; b <- 8;
width <- 0.05; p <- 0.5
sample_size <- 1000; trial_size <- 500000
gen_unif_rand <- function() {
# sample_size桁の2進少数をベルヌーイ分布に従う乱数から生成
return (sum(rbern(sample_size, p) * (rep(1/2, sample_size) ** seq(sample_size))))
}
gen_rand <- function(){
return( rdply(trial_size, gen_unif_rand()) )
}
system.time(res <- gen_rand())
res$V1 <- res$V1 * (b-a) + a
ggplot() +
layer(data=res, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + ggtitle("Bernoulli to Uniform") + xlim(4,9)
# 自己相関チェック
acf(res$V1)
small_data <- data.frame(x=res$V1[1:300])
ggplot(data=small_data) + geom_point(aes(x=seq(300), y=x))
###############################################################
# 一様分布
###############################################################
width <- 0.02; p <- 0.5
sample_size <- 1000; trial_size <- 100000
gen_unif_rand <- function() {
# sample_size桁の2進少数をベルヌーイ分布に従う乱数から生成
return (sum(rbern(sample_size, p) * (rep(1/2, sample_size) ** seq(sample_size))))
}
gen_rand <- function(){
return( rdply(trial_size, gen_unif_rand()) )
}
system.time(res <- gen_rand())
ggplot() +
layer(data=res, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + ggtitle("Bernoulli to Standard Uniform")
# 自己相関チェック
acf(res$V1)
small_data <- data.frame(x=res$V1[1:300])
ggplot(data=small_data) + geom_point(aes(x=seq(300), y=x))
###############################################################
# ベータ分布
###############################################################
# 独立に一様分布 U(0,1) に従う p+q-1 個の確率変数を大きさの順に並べ替えたとき,
# 小さい方から p 番め(大きい方からは q 番目)の確率変数 X の分布がベータ分布 B(p,q) となります。
# http://www.kwansei.ac.jp/hs/z90010/sugakuc/toukei/beta/beta.htm
install.packages("Rlab")
install.packages("ggplot2")
install.packages("actuar")
require(actuar)
require(ggplot2)
require(Rlab)
###############################################################
# 一様分布
###############################################################
width <- 0.03
p <- 0.5
digits_length <- 30
set_size <- 3
trial_size <- 30000
gen_unif_rand <- function() {
# digits_length桁の2進少数をベルヌーイ分布に従う乱数から生成
return (sum(rbern(digits_length, p) * (rep(1/2, digits_length) ** seq(digits_length))))
}
gen_rand <- function(){
return( rdply(set_size, gen_unif_rand())$V1 )
}
unif_dataset <- rlply(trial_size, gen_rand, .progress='text')
# 独立に一様分布 U(0,1) に従う p+q-1 個の確率変数を大きさの順に並べ替えたとき,
# 小さい方から p 番め(大きい方からは q 番目)の確率変数 X の分布がベータ分布 B(p,q) となります。
# 800個サンプルがあるので、400番目、後ろから399番目
p <- ceiling(set_size * 0.5)
q <- set_size - p + 1
cat("p:", p, ", q:",q)
get_nth_data <- function(a){
return(a[order(a)][p])
}
disp_data <- data.frame(lapply(unif_dataset, get_nth_data))
names(disp_data) <- seq(length(disp_data))
disp_data <- data.frame(t(disp_data))
names(disp_data) <- "V1"
x_range <- seq(0, 1, 0.001)
dens <- data.frame(y=dbeta(x_range, p, q)*trial_size*width)
ggplot() +
layer(data=disp_data, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=x_range, y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Bernoulli to Beta")
###############################################################
# 超幾何分布
###############################################################
# 非復元抽出
# 全部でN個から n個取り出す
# 1がM個、0がL個(N-M個)ある状況。
require(plyr)
N <- 50 # 全部でN個
M <- 12 # 1がM個
L <- N-M # 0がL個
n <- 15 # n個取り出す
# 試行するセット回数
trial_size = 10000
# 標本空間の生成
sample_space <- append(rep(1, M), rep(0, L))
# 1セット実行した結果を返す
gen_hypergeo_var <- function() {
return(sum(sample_space[sample.int(N)[seq(n)]]))
}
result <- rdply(trial_size, gen_hypergeo_var())
x_range <- seq(0,n)
dens <- data.frame(dhyper(x_range, M, L, n)) * trial_size
names(dens) <- c("p")
#### グラフ描画
ggplot() +
layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=x_range+.5, y=p),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Uniform to Hypergeometric")
mu <- 5
n <- 5
sigma <- 3
trial_size <- 10000
width <- 0.2
gen_t_val <- function(){
sample <- rnorm(n, mean=mu, sd=sigma)
t <- (mean(sample)-mu)/(sd(sample)/sqrt(n))
}
result <- rdply(trial_size, gen_t_val())
# ベルヌーイ分布の確率分布
dens <- data.frame(y=dt(seq(-15,15, width), n))*trial_size*width
# グラフ描画
ggplot() +
layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin",
binwidth=width, fill="#6666ee", color="gray"
) + layer(data=dens, mapping=aes(x=seq(-15,15, width), y=y),
geom="line", stat="identity", position="identity", colour="red"
) + ggtitle("Normal to t.") + xlim(-8, 8)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment