2014年5月27日火曜日

Rの様々な種類のグラフでデータの分布を理解する(hist、boxplot、stripchart)

データがどんな感じで分布しているか見たいとき、とりあえず視覚化してみますよね。Rだと簡単だし。

で、どんな視覚化があるかというと、ヒストグラム(度数分布図)はお馴染みですが、その他にも、箱ひげ図や一次元散布図なんかがあります。

下記は、正規分布に従う100個のランダムデータを発生させ、同じデータをそれぞれのやり方で視覚化したものです。


x <- rnorm(100,50,10) # 正規分布に従うランダムデータ
par(mfcol=c(3,1))     # 出力粋を3行×1列に分割
hist(x,breaks=seq(20,80,5))               # ヒストグラム
boxplot(x,horizontal=T,ylim=c(20,80))     # 箱ひげ図
stripchart(x,pch=1,cex=1.5,xlim=c(20,80)) # 一次元散布図


ヒストグラムは学校でも習ったし、直感的にも分かりやすいですよね。度数が高さで表されているので、量的にも理解しやすい。

箱ひげ図の意味は、以下のようになっています。
 ひげの左端  ・・・ 最小値
  箱の左端   ・・・ 25%点
 中央の太線  ・・・ 中央値(50%点)
  箱の右端   ・・・ 75%点
 ひげの左端  ・・・ 最大値

慣れないとピンとこないかもしれません。

一次元散布図は、データをそのまま座標にプロットした感じですね。データ数が多くなってくると、データの粗密はなんとなく分かりますが、数はよく分からなくなってしまいます。

こうやって並べてみると、ヒストグラムでいいじゃんって気もしますが、そうでない場合もあります。複数のデータを比較したいような場合は下の2つの方が使いやすい気がします。

↓こんな感じで、標準偏差の異なる複数のデータを作ってみます。

> sd_05 <- rnorm(100,50, 5)
> sd_10 <- rnorm(100,50,10)
> sd_15 <- rnorm(100,50,15)
> sd_20 <- rnorm(100,50,20)
> sd_25 <- rnorm(100,50,25)
> sd_30 <- rnorm(100,50,30)
>
> x <- data.frame(sd_05,sd_10,sd_15,sd_20,sd_25,sd_30)
>
> head(x)

     sd_05    sd_10    sd_15    sd_20    sd_25    sd_30
1 46.86773 43.79633 56.14103 67.87347 76.86102 52.31909
2 50.91822 50.42116 75.33310 29.05404 97.39137 41.09394
3 45.82186 40.89078 73.79883 89.42675 34.92507 14.50273
4 57.97640 51.58029 45.03638 42.32736 40.22830 50.33878
5 51.64754 43.45415 15.72147 83.08291 39.59445 79.74803
6 45.89766 67.67287 87.46492 80.24425 40.60856 97.81902


複数のグラフを並べてみます。

↓ヒストグラム


par(mfcol=c(6,1))
for(i in 1:6){
  hist(x[,i],breaks=seq(-50,150,10),main=colnames(x)[i])
}

無理矢理に縦に並べてみましたが、やっぱり場所をとりますね。

↓下記みたいに、2列にすればレイアウト的にはすっきりしますが、比較しにくくなります。


par(mfcol=c(3,2))
for(i in 1:6){
  hist(x[,i],breaks=seq(-50,150,10),main=colnames(x)[i])
}


箱ひげ図だと、↓こんな感じになります。
boxplot(x,horizontal=T)



箱ひげ図のデフォルトは縦方向なので、「水平だぞ」というhorizontalオプションを付けています。

一次元散布図だと、↓こんな感じになります。
stripchart(x,pch=1)



プロットする記号は□(pch=0 デフォルト)よりは、○(pch=1)の方が見やすいかなと思います。

箱ひげ図も一次元散布図も、基本的には一次元的な表現方法なので場所をとりませんね。

以下は30個並べてみたものですが、なんとか見られますね。
n <- 30
x <- data.frame(matrix(rep(0,n*100),ncol=n))
for(i in 1:n){
  x[,i] <- rnorm(100,50,i)
}
boxplot(x,horizontal=T)


↓一次元散布図
n <- 30
x <- data.frame(matrix(rep(0,n*100),ncol=n))
for(i in 1:n){
  x[,i] <- rnorm(100,50,i)
}
stripchart(x,pch=1)


一次元散布図は密度が高いと、プロットがつぶれちゃうので、データ数が少ないとき向きのように思います。被験者が10人とか。↓こんな感じ。


x <- rnorm(10,50,10) # 正規分布に従うランダムデータ
par(mfcol=c(3,1))    # 出力粋を3行×1列に分割
hist(x,breaks=seq(20,80,10))              # ヒストグラム
boxplot(x,horizontal=T,ylim=c(20,80))     # 箱ひげ図
stripchart(x,pch=1,cex=1.5,xlim=c(20,80)) # 一次元散布図


これくらいのデータ数だと、個々のデータの値が見える一次元散布図のメリットが出てきますね。




Rのpolygon関数でグラフを塗りつぶす方法を理解する

グラフの一部を塗りつぶす方法として、segments関数で線分を敷き詰めるやり方を↓こちらで紹介しましたが、

segments関数でグラフを塗りつぶす - Rプログラミングの小ネタ

正直なところ、polygon関数で塗りつぶすやり方↓の方がスマートな気がします。

51. 低水準作図関数 - R-Source

でも、上記に掲載されているサンプルを見ても、すぐに理解できない人も多いのでは?ということで、分かりやすく解説する記事を書いてみようと思います。Rエキスパートの人は帰った、帰った。

例えば、サインカーブの下部を塗りつぶしたいとします。あえて分割数を少なくして、上記に紹介されている方法を適用してみると、↓こんな感じになります。
x座標を3点(polygon関数による描画)

plot(sin, 0, pi) # sinカーブをプロット(0からπまで)
n  <- 3                    # x軸上に取る点の数
xs <- seq(0, pi, length=n) # プロットするx座標を決める
ys <- sin(xs)              # 曲線上のy座標を求める
 
# x軸上→曲線上の順で、頂点を指定する
polygon( c(xs,rev(xs)), c(rep(0,n),rev(ys)), col="gray")


3点だけだと、全然曲線っぽくないですよね。

点の数を増やしていきます。(サンプルソースの中のnを変える)

x座標を4点(polygon関数による描画)


x座標を5点(polygon関数による描画)


x座標を10点(polygon関数による描画)
x座標を20点(polygon関数による描画)


↑20点も取れば、ほとんど曲線に見えますね。


x座標を100点(polygon関数による描画)


↑上記は100点。これ以上増やしても見栄えは変わらないと思います。

次は、polygonへの頂点指定の仕組みを見ていきましょう。

c関数で連結してあるので分かりにくいですが、最初にx軸上の点を、次に曲線上の点が指定されています。

polygon( c(xs,rev(xs)), c(rep(0,n),rev(ys)), col="gray")

青がx軸上の点の指定、赤がy軸上の点の指定となっています。

rep関数は、0をn個繰り返す(x軸上だからy座標は0なので)、rev関数は、ysベクトルの要素を反対に並べる(帰りがけにプロットしていくので)、という意味です。

点に番号を振って、分かりやすくしてみました。

頂点をナンバリングしたもの(polygon関数による描画)


plot(sin, 0, pi) # sinカーブをプロット(0からπまで)
n  <- 5                    # x軸上に取る点の数
xs <- seq(0, pi, length=n) # プロットするx座標を決める
ys <- sin(xs)              # 曲線上のy座標を求める
 
# x軸上→曲線上の順で、頂点を指定する
polygon( c(xs,rev(xs)), c(rep(0,n),rev(ys)), col="gray") 
 
text(xs, rep(0,n), 1:n)                    #行きがけ(x軸上)
text(rev(xs), rev(ys), (n+1):(2*n), pos=3) #帰りがけ(曲線上)


これらの番号が振られた点を頂点とするポリゴンを描画せよ、という命令なわけですね。




2014年5月26日月曜日

ブログにRスクリプトを載せるときにハイライト表示する

ブログ記事にjavaとかhtmlとかperlとかのソースコードを載せるときは、SyntaxHighlighterを使うといい感じになりますよね。予約語やコメントを強調表示してくれて見やすくなる

でも、SyntaxHighlighterはRには対応していない模様。

で、検索してみると、inside-RというサイトのPretty R syntax highlighterという機能で同じようなことができました。

べたにソースを貼り付けると↓こんな感じです。(コメントとソースが分かりやすいように色付けだけしてみました)

# 正規分布に従う乱数を1000個
x <- rnorm(1000)

# 表示領域を3行×2列に分割
par(mfrow=c(3,2))

# 階級幅の数を変えながら複数のヒストグラムを描画
for (n in c(5,10,15,20,25,30)){
  hist(x,breaks=seq(min(x),max(x),length=n+1),main=paste(n,"分割"))
}


Pretty R syntax highlighterを使ってみます。

ソースコードを貼り付けて「Go!」ボタンを押す





見た目とHTMLコードが表示される



あとは、HTMLをコピーして、ブログにHTMLとして貼り付ければOKです。

↓こんな感じにできあがります

# 正規分布に従う乱数を1000個
x <- rnorm(1000)
 
# 表示領域を3行×2列に分割
par(mfrow=c(3,2))
 
# 階級幅の数を変えながら複数のヒストグラムを描画
for (n in c(5,10,15,20,25,30)){
  hist(x,breaks=seq(min(x),max(x),length=n+1),main=paste(n,"分割"))
}
Created by Pretty R at inside-R.org


最後の方のinside-R.orgへのリンクが気に入らない人は、HTMLの末尾のアンカータグの部分↓

<p><a href="http://www.inside-r.org/pretty-r" title="Created by Pretty R at inside-R.org">Created by Pretty R at inside-R.org</a></p>

を削除すればいいです。

元のRスクリプトは、カンマの後ろにスペースを入れてないので、モノクロだと読みにくかったりするんですが、色が付くとスペースがなくても見やすくなりました。




2014年5月23日金曜日

Rの複数グラフ表示でヒストグラムの適切な階級幅を見つける(par,mfrow,mfcol)

ヒストグラムの階級幅をどれくらいにするか。小さすぎるとばらつきの影響が出てしまうし、大きすぎると意味が読み取りづらくなってしまいますよね。

パラメタを変えながら何度もhist関数を実行し、表示結果を目で見て判断というのが一般的でしょうか。

どうせなら、いっぺんに描画してまとめて比較しましょう、というやり方です。

# 正規分布に従う乱数を1000個
x <- rnorm(1000)
 
# 表示領域を3行×2列に分割
par(mfrow=c(3,2))
 
# 階級幅の数を変えながら複数のヒストグラムを描画
for (n in c(5,10,15,20,25,30)){
  hist(x,breaks=seq(min(x),max(x),length=n+1),main=paste(n,"分割"))
}
 
上記の表示結果だと、「20分割まで行っちゃうと単峰性が崩れるから、15分割くらいかな」みたいな感じでしょうか。

上記のように mfrow=c(3,2) と指定すれば3行×2列を行方向に並べる、mfcol=c(3,2) と指定すれば3行×2列を列方向に並べる、となります。

length = n+1 となっているのは、5個の階級に分けるためには、6個の境界が必要だからですね。

[蛇足]
関数名やパラメタ名を覚えるにはその語源を知るといいですよね。parはパラメタ(parameter)、rowは行(row)、colは列(column)というのはよさそうですが、mfが何の略だか分からない。「matrix figure」なんじゃないかと個人的には思っていますが、でたらめかも。コメントで教えてくれた人には10ガバスあげます。




2014年5月21日水曜日

Rで繰り返しを含む数列の生成(rep関数、seq関数)

Rで、「123123123...」とか「111222333...」とか「13579...」(奇数)とか「02468...」(偶数)とか、そんな感じの繰り返しを含むような数字の列を作りたい場合、1行でさらっと書けてしまいます。

【繰り返し】
> rep(7, 10)
 [1] 7 7 7 7 7 7 7 7 7 7

「○を○個だけrepeat」って感じで。


【1ずつ増加】
> 1:7
 [1]  1  2  3  4  5  6  7

コロンでつなげば「○から○までインクリメント」


【○ずつ増加】
> seq(1, 7, 2)
[1]  1  3  5  7

sequence関数で、「○から○まで○ずつ増加」


【1ずつ増加する数列を繰り返す】
 > rep(1:7, 3)
 [1] 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7

repeat関数の第1引数に、「1ずつ増加する数列」を指定すればOK。


【○ずつ増加する数列を繰り返す】
> rep(seq(1,7,2), 3)
 [1] 1 3 5 7 1 3 5 7 1 3 5 7

repeat関数の第1引数に、「○ずつ増加する数列」を指定すればOK。


【それぞれ○個】 
> rep(1:7, each=3)
 [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7


【○ずつ増加、それぞれ○個】
> rep(seq(1,7,2), each=3)
 [1] 1 1 1 3 3 3 5 5 5 7 7 7

コロンで指定していたところを、repにするだけなので簡単ですね。


【○ずつ増加、○回繰り返す、それぞれ○個】
> rep(seq(1,7,2), 2, each=3)
 [1] 1 1 1 3 3 3 5 5 5 7 7 7 1 1 1 3 3 3 5 5 5 7 7 7

rep関数の第2引数は繰り返し回数で、「rep(seq(1,7,2), times=2, each=3)」と書いても同じ意味です。こう書いた方が、勘違いは少なそうですね。




2014年5月18日日曜日

Rのsegments関数でグラフを塗りつぶす

やりたいのは↓こういうことだったりします。


51. 低水準作図関数 | R-Source より

やり方も上記リンク先に書いてあるので、ザッツオールです。

でも、領域を複数の台形に分けてpolygon関数で塗りつぶすという上記のやり方って、処理を頭の中で想像するのがけっこう難しい。

もっと単純に、縦棒をびっしり敷き詰めて塗りつぶすというやり方なら、頭の中で想像しやすいんじゃないかと思って、やってみました。



























plot(dnorm, -4, 4)
xvals <- seq(2, 4, length=1000)
dvals <- dnorm(xvals)
segments( xvals,0, xvals,dvals )

最初の3行は、前述のリンク先の流用でして、塗りつぶしのところだけ、線分を描画するsegments関数を使いました。segments(x1, y1,  x2, y2) とやると、(x1,y1) から (x2,y2) まで線分を引いてくれるという関数です。引数がベクトルなら複数の線分を引いてくれます。

length=20 と少なくしてしまうと、この(あまり賢くない)やり方のボロが出てしまいます↓


























なので、塗りつぶされている感じが出る程度まで、lengthを増やす必要があります。たくさんのグラフを描画するような場合は、このやり方は重くなってしまいそうですね。でも、グラフ1枚程度なら、1000本でも10000本でも大したことにはならないです。

塗りつぶしの方法としては理解しやすいでしょ。でも、やり方としてはpolygonで塗りつぶす方が賢いと思います。(実も蓋もない結論ですが・・・)




緒言

統計ツール・言語である「R」を仕事でよく使います。悩んだときに、Webで検索して見つけた皆さんのtipsに助けられることがよくあります。数行のサンプルコードと、その実行結果があるだけで、相当助かるんですよね。

というわけで、私もちょっとしたサンプルコードとその実行結果の痕跡を残していこうと考えました。過去に自分がやったことも忘れちゃうんで・・・