00_distributions.r

# Preparing
install.packages("Rlab")
install.packages("ggplot2")
install.packages("actuar")
install.packages("plyr")
require(plyr)
require(actuar)
require(ggplot2)
require(Rlab)

01_bernoulli.r

###############################################################
# ベルヌーイ分布
###############################################################
# ベルヌーイ分布からのサンプリングを実行

# パラメーター
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)

02_binomial.r

###############################################################
# 二項分布
###############################################################
# ベルヌーイ分布に従う変数を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.")

03-02.standard_norm.r

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.") 

03_normal.r

###############################################################
# 正規分布
###############################################################
# 二項分布の正規近似 (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.") 

04_chi_square_df1.r

###############################################################
# カイ二乗分布 (自由度: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"
  )

05_chi_square_df3.r

###############################################################
# カイ二乗分布 (自由度: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")

06-02_poisson.r

###############################################################
# ポアソン分布
###############################################################
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.")

06_poison.r

###############################################################
# ポアソン分布
###############################################################

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"
  )

07_exponential.r

###############################################################
# 指数分布
###############################################################

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.")

08_gamma.r

###############################################################
# ガンマ分布
###############################################################

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")

09_inv_gamma.r

###############################################################
# 逆ガンマ分布
###############################################################
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")

10-2.uniform.r

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

10_std_uniform.r

###############################################################
# 一様分布
###############################################################
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))

11_beta.r

###############################################################
# ベータ分布
###############################################################
# 独立に一様分布 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")

12_hypergeometric.r

###############################################################
# 超幾何分布
###############################################################
# 非復元抽出
# 全部で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))

# １セット実行した結果を返す
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")

13_t_dist.r

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)