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億回 対数プロット

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

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




2017年1月31日火曜日

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

↓この記事を書いていて思ったのですが、

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


一日の生活の900分は一次元的であるのに、それを30×30の二次元のマスで表現しているところが、そもそも分かりにくい。

人が理解するときのモデルとしても分かりにくいし、二次元になっているせいで、シミュレーションのスクリプトを書くときもいろいろと面倒な処理が必要になったりして(上記のリンクのinteractionのくだりとか)。

素直に、一次元的な図を載せた方が、読者の理解も進むのではと思って、Rで書いてみました。(点やマスの数は見た目がほどよくなるように減らしてあります)

n <- 80 # 点の数
m <- 10 # マスの数
x <- runif(n, min=0, max=m)
stripchart(x, pch=1, xlab="1分ごとのマス")
abline(v=0:m)
legend("topright", legend="手の動きのあった時点", pch=1, bg="white")




 こういう図の方が分かりやすいと思うけどなあ。




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

データの見えざる手」は読んでいて、引っ掛かりまくりでした。

読んだかた、「U分布」ってピンときました?

どこで躓いたかというと、こんな感じでコンピュータシミュレーションの結果が紹介されていたのですが、上の「正規分布(ポアソン分布)」と書かれている図↓の見た目が一様分布っぽいんですよね。

正規分布(ポアソン分布)とU分布 (「データの見えざる手」より)

本文を見てみると、
コンピュータシミュレーションでこれを実行するには、玉の位置をランダムに生成すればよい。横方向の位置(x)を決める1~30の乱数と縦方向の場所(y)を決める1~30の乱数を発生させ、(x,y)の位置に玉を置くのだ。

とあります。う~ん、xとyを一様分布に従って発生させて2次元にプロットしている、ってだけだよなー。

・・・(考え中)・・・

で、しばらく考えてみて、やっと分かりました。こうやって発生させたデータだと、マスの中の点の個数が正規分布になるのね。

実際に、一つずつやってみましょう。

n <- 1000 # 点の個数
m <- 10   # マスの区切りの数

# 座標は一様分布で従って発生させ、プロット
x <- runif(n, min=0, max=m)
y <- runif(n, min=0, max=m)
plot(x, y, pch="・")

# マスを書く
abline(h=0:m)
abline(v=0:m)
一様分布で発生させた点(マス内の個数は正規分布になる)

はい、本に載っているのと似たような絵になりました。

マスに入る個数を調べるには、引数を越えない整数を返すceilingが使えますね。

> interaction( ceiling(x), ceiling(y), sep="," )
   [1] 4,1   9,3   10,1  1,5   3,1   1,6   8,1   8,3   10,6  5,7   2,6
  [12] 2,5   2,7   7,6   6,10  4,9   3,3   1,1   2,4   9,6   10,10 2,3
  ...

マスが、(1,1), (1,2), ..., (1,10), (2,1), ..., (10,10)のように100個あって、1つ目の点は(4,1)のマスに、2つ目の点は(9.3)のマスに入っている、という感じです。

集計します。

> table( interaction( ceiling(x), ceiling(y), sep="," ) )

  1,1   2,1   3,1   4,1   5,1   6,1   7,1   8,1   9,1  10,1   1,2   2,2
   11    13    13    13     8     3    11    12    13     6    11    12
  3,2   4,2   5,2   6,2   7,2   8,2   9,2  10,2   1,3   2,3   3,3   4,3
    4     5     6     8    12    10     6    13     8     8    18     9
  ...

(1,1)のマスには11個の点、(2,1)のマスには13個の点が入っているようです。

ヒストグラムにしてみると、

マス100個、点1000個で発生させた正規分布
う~ん、正規分布かぁ?

点とマスが少な過ぎたようです。

書籍に書いてあった数字でやってみましょう。

n <- 72000 # 点の個数
m <- 30    # マスの区切りの数
x <- runif(n, min=0, max=m)
y <- runif(n, min=0, max=m)
t <- table( interaction( ceiling(x), ceiling(y), sep="," ) )
hist(t)

マス900個、点72000個で発生させた正規分布
うん、正規分布っぽいですね。

で、つまり手の動きが一様分布に従って発生するならば、単位時間あたりに動いた回数は正規分布に従うのだろうけど、実際はそうではなかった、ということみたいです。

考えている時間、キーボードを打つ時間、座っている時間、歩いている時間などがあることを考えると、手の動きが一様に発生するわけないってのは、なんだか当たり前のような気がします。




2016年12月8日木曜日

「エレガントな問題解決」演習問題 2.1.27(ひっかけ問題)(c)の解答

いつもは、↓別のブログの方に載せていたのですが、

「エレガントな問題解決」演習問題 2.1.25の解答(分母が3つの項の積になっている数列の和): 主張

今回はRを使うので、こちらのブログで。

「エレガントな問題解決」 P29 2.1.27(ひっかけ問題)(c)
偏りのないコインを投げて、表が3回以上出る確率が50%を超えるには、最低何回投げればよいだろうか。
自分が題意を理解できているのか、いまいち自信がありません。私は、ひっかかっているのか?

とにかく解答してみます。

コインを1回投げて、表が3回以上出ることは、もちろんありません。

なので、考慮に値するのは3回以上投げるケースから。

3回投げて、表が3回出る確率をRで求めてみると、

> (1/2)^3
[1] 0.125

なので、12.5%。これは50%を超えていませんね。

では、4回投げる場合。表が3回出る確率と、表が4回出る確率を足せばいいですね。

4回投げて、表が3回出る確率は、(Cは組み合わせの記号だと思ってください)

(1/2)^4 × 4C3

ですね。Rでは、組み合わせ数を求めるのにchooseコマンドが使えるので、

> (1/2)^4 * choose(4,3)
[1] 0.25

4回投げて、表が3回出る確率は、

> (1/2)^4 * choose(4,4)
[1] 0.0625

上記の2つを足しても、まだ50%は超えませんね。
いちいち足すのもまどろっこしいので、関数にしてみましょう。

#
# n回投げたときに、3回以上でる確率を返す関数
#

ProbMoreThan3 <- function(n){
  s <- 0    # 組み合わせ数の合計値格納用
  for(i in 3:n){
    s <- s + choose(n,i)
  }
  p <- 1/2^n * s
  return(p)
}


試しに呼んでみましょう。

> ProbMoreThan3(3)
[1] 0.125

> ProbMoreThan3(4)
[1] 0.3125

よさそうです。(エラー処理は入れていないので、2とか小数とか、そういうのを引数にするのは禁止ね)

じゃあ、いつ50%を超えるか、やっていきましょう。

> ProbMoreThan3(5)
[1] 0.5

> ProbMoreThan3(6)
[1] 0.65625

> ProbMoreThan3(7)
[1] 0.7734375

「確率が50%を超える」なので、正解は「6回」です。たぶん。

どのあたりが「ひっかけ」だったのでしょうか。

「3回以上」なのに、「ぴったし3回」の確率を求めしまうとか。

「50%を超える」なのに、50%になったとき(=5回)を答えてしまうとか。

それとも、何か別のひっかかる要素があって、私もまんまとそれにひっかかっているとか。
謎は深まります。

が、考えても分かる気がしないので、時間つぶしにグラフでも書いてみましょう。

n <- 3:20
p <- sapply(n, ProbMoreThan3)
plot(n,p)
abline(h=0.5, lty=2)


n回投げて3回以上表が出る確率

初めの方でぐっと上がって、それから1に漸近していくような感じですね。




2016年11月2日水曜日

Rで順列、組み合わせを計算する

最近知ったんですが、Rで順列(Permutation)や、組み合わせ(Combination)を計算するデフォルトの関数ってないんですね。

【2016年11月18日追記 開始】
勘違いしていました。訂正します。
組み合せの数を計算するchooseという関数があります。
例えば、4個から2個選ぶ組み合わせの数は

> choose(4,2)
[1] 6

みたいな感じで求めることができます。
【2016年11月18日追記 終了】

そういうパッケージはあるみたいですけど、わざわざパッケージを使うまでもないような気もしないでもないです。

高校数学を復習する意味でも、自分で実装してみましょう。

まず、順列。

n個の中からk個を取り出して、順番を意識して並べていくような場合の数ですね。

nから初めて、1ずつ減らしながら、r個分かければ計算できます。

    n × (n-1) × (n-2) × ... × {n-(r-2)} × {n-(r-1)}

上記をやりたい場合、受け取った引数を全部かけてくれるprodという関数が使えます。1ずつ減っていく引数は、 ":" (コロン)を使えばよさそうですね。

ということで、

    prod(n:(n-r+1))

とやればいいのですが、これだと r が 0 のときにうまくいきません。0個を選び出す順列の数は 1 なので、関数にしてしまって、↓こんな感じで、どうでしょうか?

# 順列を計算する関数定義
Permu <- function(n, r){
  ifelse( r==0, 1, prod(n:(n-r+1)) )
}


動作確認してみましょう。

> Permu(4, 0)
[1] 1
> Permu(4, 1)
[1] 4
> Permu(4, 2)
[1] 12
> Permu(4, 3)
[1] 24
> Permu(4, 4)

[1] 24

うん、よさそうです。

次は、組み合わせ。

nからr個選ぶ順列を、rの階乗で割れば(r個の順番を無視する)、OKですね。

階乗はfactorialという関数があるので、そのまま使えます。

こちらも関数定義にしてみましょう。

# 組み合わせを計算する関数定義
Combi <- function(n, r){
  Permu(n, r)/factorial(r)
}


こちらも動作確認してみましょう。

> Combi(4, 0)
[1] 1
> Combi(4, 1)
[1] 4
> Combi(4, 2)
[1] 6
> Combi(4, 3)
[1] 4
> Combi(4, 4)

[1] 1

うん、よさそうです。