「データの見えざる手」のU分布を、Rでシミュレート（改）

↓以前、こちらの記事を書いたのですが、

「データの見えざる手」のU分布を、Rでシミュレート

よくよく見てみると、書籍と軸の取り方が違ったりして、つっこみどころ満載だったので、悔い改めて、ちゃんとやることにしました。

書籍では、横軸がマスに入っている個数、縦軸が累積確率になっていました。

あと、初期値もちゃんとランダムで設定するようにしました。

それと、対数プロットする際に、軸の目盛ラベルの付けやすさから、ggplotを使ってみました。

library(ggplot2)
library(scales)

n <- 72000 # 点の個数
m <- 900   # マスの個数
masu <- numeric(m) # 空のマスを用意

# 点をランダムにマスに配置
indices <- sample(1:900, 72000, replace=T)
for(i in indices){
  masu[i] <- masu[i] + 1
}

for( i in 1:100000000 ) {
  s <- sample(1:m, 2) # ランダムに2つのマスを選ぶ
  if( masu[s[1]] != 0 ) { # 無い袖は振れないケースへの対処
    masu[s[1]] <- masu[s[1]] - 1 # 1つ目のマスから取って、
    masu[s[2]] <- masu[s[2]] + 1 # 2つ目のマスへ入れる
  }
}

tbl <- table(masu)    # 個数を集計
df <- data.frame(tbl) # データフレームにする

# 列名を分かりやすくする
colnames(df) <- c("num_of_dots", "freq")

# 点の数が因子型なので、整数型に変えておく
df$num_of_dots <- as.integer(df$num_of_dots)

# 累積確率を計算
df$cum_prob <- rev(cumsum(rev(df$freq))/m)

# プロット
ggplot(df, aes(x=num_of_dots, y=cum_prob)) +
  geom_point() +
  scale_y_continuous(
    trans = log2_trans(),
    breaks = c(1, 1/2, 1/4, 1/8, 1/16, 1/32, 1/64),
    labels = c("1", "1/2", "1/4", "1/8", "1/16", "1/32", "1/64")
  )

「# 累積確率を計算」のところがごちゃごちゃしてますが、個数が多い方からの累積なので、一旦 rev で降順にして、cumsumで累積度数を計算して、m(=900)で割って累積確率にして、再び rev で順番を戻しています。

1億回のループの結果は・・・、ジャン↓

シミュレーションで生成したU分布の曲線

書籍の図と似たものになりました。よかった、よかった。