2017年4月13日木曜日

table関数の出力結果をmatrixに変換して、corresp関数で対応分析を行う(R言語)

 ↓この本に載っていた例で、

Rによるデータサイエンス データ解析の基礎から最新手法まで
Rによるデータサイエンス データ解析の基礎から最新手法まで

library(MASS)
caith

       fair red medium dark black
blue    326  38    241  110     3
light   688 116    584  188     4
medium  343  84    909  412    26
dark     98  48    403  681    85

↑このようなデータ(縦に並んでいるのが目の色、横に並んでいるのが髪の色)に対して、↓こんな感じで対応分析を行う、という例が載っていました。

caith.ca <- corresp(caith, nf=4)
biplot(caith.ca)

Rでのcaithデータの対応分析結果

と、ここまでが前置き。

今、手元にあるのが↓こんなデータだった、としましょう。

gender <- c("M","M","M","M","F","F","F","F","F","F")
blood  <- c("A","B","B","O","A","B","A","O","A","AB")

gender.blood <- data.frame(gender, blood)

gender.blood

   gender blood
1       M     A
2       M     B
3       M     B
4       M     O
5       F     A
6       F     B
7       F     A
8       F     O
9       F     A
10      F    AB

このデータに対して対応分析を行いたい、としましょう。

先ほどのデータと違うのは、まだ集計されていないデータであるということ。correspの入力はクロス集計して分割表になったものでしたので、その形式にしてやらなければなりません。

こういうときにはtable関数ですね

gender.blood.tbl <- table(gender.blood)
gender.blood.tbl

      blood
gender A AB B O
     F 3  1 1 1
     M 1  0 2 1

うん、集計されました。でも、これをそのまま corresp(gender.blood.tbl[,]) とやっちゃうと、「invalid table specification」というエラーが出ちゃうんですよね。どうやら、table型は受け付けてくれないようです。

そういうときは、↓こんな風に書いてやると、テーブルがmatrixになるみたいで、corresp関数が入力を受け付けてくれます。

gender.blood.tbl.ca <- corresp(gender.blood.tbl[,], nf=2)
biplot(gender.blood.tbl.ca)

データはデタラメなので、表示結果についてはスルーしてください。




2017年3月29日水曜日

Rで、データフレームの中身を一括で因子型に変換する方法

例えば、Rで↓こんなデータを扱っているとします。

# サンプルデータを作る
Q1 <- c(1, 1, 2, 2)
Q2 <- c(1, 2, 1, 2)
Q3 <- c(1, 2, 3, 1)
df <- data.frame(Q1, Q2 ,Q3)

# 中身の確認
df

  Q1 Q2 Q3
1  1  1  1
2  1  2  2
3  2  1  3
4  2  2  1

読み込んだときの都合か何かで、データはinteger型とかnumeric型になっていると。

でも、実は質問 Q1、Q2、Q3に、1:はい、2:いいえ、3:どちらともいえない、とかで答えたもので、因子型として扱いたい。

多重対応分析のmcaとかを使おうとすると、

mca(df = df) でエラー: all variables must be factors

が出ちゃうとか、そんなシチュエーション。

因子型に変換したいときには、as.factor関数ですが、これはデータフレームに対しては使えない。

こんな時は、あの一家。そう、applyファミリーの登場です。

lapplyを使って、1列ごとにas.factor関数を適用、リストとして返ってきたものを、またデータフレームに戻してやるという流れです。

df.fctr <- data.frame( lapply(df, as.factor) )

# 型の確認
df.fctr$Q1 

[1] 1 1 2 2
Levels: 1 2

無事、因子型になりました。




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個で発生させた正規分布
うん、正規分布っぽいですね。

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

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