2017年2月15日水曜日

Rでヒストグラムの一部に色をつける(colオプション指定で可)

「R ヒストグラム 一部 色をつける」で検索してみると、hist関数でヒストグラムを描いた後に、polygon関数で色をつける、なんて方法がヒットしました。

polygon使えばなんでもできそうだけど、なんか、ちょっと違うよなあ、とか思ってしまいまして。

で、実はhist関数のcolオプションでも、できるんですよね。

colオプションに1つの値(スカラー)を指定すると、全体が一色で塗りつぶされてしまいますが、ここにベクトルを指定すると、それぞれの棒の色を指定することができます。

例えば、ヒストグラムに10個のビンがあって、それぞれを任意の色で塗りたい場合は、10個の要素を持つベクトルをcolオプションに指定すればOKです。

set.seed(0)      # 再現性のために
rd <- rnorm(100) # 100個の乱数
cols <- c("white", "white", "red"  , "white", "white",
          "blue" , "white", "white", "white", "white")
hist(rd, col=cols)

ヒストグラムの一部に色をつける

色を塗りたくない場合は(パワポなどの「塗りつぶしなし」みたいな感じ)、色名の代わりにNAを指定すればいいです。

cols <- c(NA    , NA, "red", NA, NA,
          "blue", NA, NA   , NA, NA)
hist(rd, col=cols)

最初の例と全く同じ見た目になると思いますが、add=T指定で重ねたときなんかに差がでますね。

階級がいっぱいあって、いちいち全部書き出すのが面倒なときは、下記のような感じで、塗りたいところだけを指定すればいいですね。

cols <- rep("white", 100) # ”white”を詰めた、長めのベクトルを作っておく
cols[3] <- "red"
cols[6] <- "blue"
hist(rd, col=cols)

下記のようにすれば、最頻値に対応する棒を赤く塗る、なんてこともできます。

最初のhistは絵を描くためではなく、階級がどのように分けられるか、それぞれの度数がいくつか、なんて情報を得るためにやっています。

h <- hist(rd)                           # 度数分布に関する情報を得るために
cols <- rep("white", length(h$counts))  # 今度はビンの個数が分かるぞ
cols[which.max(h$counts)] <- "red"      # 最頻値のビンの色を赤にする
hist(rd, col=cols)

ヒストグラムの最頻値に色をつける

あと、線の色を変えたい場合は、borderオプションですね。

やってみると・・・

cols <- rep("black", length(h$counts)) # 枠線のデフォルトはblackなので
cols[which.max(h$counts)] <- "red"     # 最頻値の枠線を赤にする
hist(rd, border=cols)

border指定は右側の線で上書きされてしまう・・・

あー、重なった線の部分が、右側の黒で上書きされてしまった。

こういうときは・・・、polygon関数でも使ってください




2017年2月1日水曜日

「データの見えざる手」の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分布の曲線

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




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

さて、本筋とは関係のないところでケチばっかりつけていた、↓前回と前々回の記事でしたが、

「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した

「データの見えざる手」の図が分かりにくかったので、Rで一次元プロット

今回は、著者の矢野和男さんの言うところの「U分布」なるものを、計算機シミュレーションで作り出してみましょう。

次に、このようにランダムに玉を分配した後で、マス目間で玉をやりとりさせてみよう。
ランダムにマス目を二つ選んで、一方から他方に玉を1個移す。そして、これを繰り返してみよう。もともと、ランダムに置いた玉なのだから、そこからランダムにマス目を選んで、玉を動かしても、結果は変わらない、と思うだろう。この問題を多くの人に出題してみたが、全員が「結果は変わらない」と答えた。

たしかに、直感的には、ランダムに配置後にランダムに交換しても、マクロな状況は変わらないような気がしますね。でも、そうじゃないのが興味深いところ。

書籍では初期値はランダムとありましたが、手を抜いて、1マス80個の「平等」状態からスタートさせてみました。(結果は同じになりますよね、たぶん)

指定回数の交換を行ったあと、それぞれのマスが持っている個数でソートして、少ない方が左になるようにプロットしています。

1万回、2万回、・・・と実行しながら、プロット結果を画像として出力していきます。jのところのループ回数を変えて、1億回まで実行してみました。


n <- 72000  # 点の個数
m <- 900    # マスの個数
masu <- rep(n/m, m) # 平等に配分

for( i in 1:9 ) {
  for( j in 1:10000 ) {
    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つ目のマスへ入れる
    }
  }
  s1 <- paste(i, "万回(リニア)", sep="")
  s2 <- paste(i, "万回(対数)", sep="")

  # リニアでプロット
  png(paste(s1,"png",sep="."))
  plot(masu[order(masu)], main=s1)
  dev.off()

  # y軸を片対数でプロット
  png(paste(s2,"png",sep="."))
  plot(log(masu[order(masu)],2), main=s2)
  dev.off()
}


まず、1万回↓

1万回 リニアプロット

全然、指数関数っぽくないですよね。まだ「交換」が十分に行われていないので、80個付近のマスが多いため、真ん中あたりが踊り場のようになっています。そんな中でも、0個近くの負け組と200個以上の勝ち組が現れてきています。

次、10万回↓

10万回 リニアプロット

あまり変わらないですね。(機械的にやったのでグラフのタイトル表記が変ですが、気にしないでね)

そして、100万回↓

100万回 リニアプロット

踊り場が無くなってきました。

さらに、1000万回↓

1000万回 リニアプロット

おお、指数関数っぽくなってきました。

最後に、1億回

1億回 リニアプロット

1000万回とあまり変わっていません、収束してきたのかな。


さて、曲線を見る限りは指数関数っぽいですが、対数をとってグラフが直線になっているのかどうか見てみましょう。

1万回 対数プロット
10万回 対数プロット
100万回 対数プロット
1000万回 対数プロット
1億回 対数プロット

う~ん、直線にはなりませんねえ。ということは、このシミュレーション結果は厳密には指数関数にはならないものなのでしょうか、それとも私の組んだロジックに不備があるのでしょうか。

でも、中間のあたりに関しては直線になっていると言えなくもないですね。