Rプログラミングの小ネタ

データがどんな感じで分布しているか見たいとき、とりあえず視覚化してみますよね。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%点 　ひげの左端　　・・・　最大値 慣れないとピンとこないかもしれません。 一次元散布図は、データをそのまま座標にプロットした感じですね。データ数が多くなってくると、データの粗密はなんとなく分かりますが、数はよく分からなくなってしまいます。 こうやって並べてみると、ヒストグラムでいいじゃんって気もしますが、そうでない場合もあります。複数のデータを比較したいような場合は下の２つの方が使いやすい気がします。 ↓こんな感じで、標準偏差の異なる複数のデータを作ってみます。 > sd_05 <- rnorm(100,50, 5) > sd_10 <- rnorm(100,50,10) > sd_15 <- rnorm(100,50,15) > sd_2

グラフの一部を塗りつぶす方法として、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" ) ３点だけだと、全然曲線っぽくないですよね。 点の数を増やしていきます。（サンプルソースの中のnを変える） x座標を4点（polygon関数による描画） x座標を5点（polygon関数による描画） x座標を10点（polygon関数による描画） x座標を20点（polygon関数による描画） ↑20点も取れば、ほとんど曲線に見えますね。 x座標を100点（polygon関数による描画） ↑上記は100点。これ以上増やしても見栄えは変わらないと思います。 次は、polygonへの頂点指定の仕組みを見ていきましょう。 c関数で連結してあるので分かりにくいで

ブログ記事に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の末尾のアン

ヒストグラムの階級幅をどれくらいにするか。小さすぎるとばらつきの影響が出てしまうし、大きすぎると意味が読み取りづらくなってしまいますよね。 パラメタを変えながら何度も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分割くらいかな」みたいな感じでしょうか。 上記のように mf row =c(3,2) と指定すれば3行×2列を 行方向 に並べる、mf col =c(3,2) と指定すれば3行 ×2列を 列方向 に並べる、となります。 length = n +1 となっているのは、5個の階級に分けるためには、6個の境界が必要だからですね。 [蛇足] 関数名やパラメタ名を覚えるにはその語源を知るといいですよね。parはパラメタ（parameter）、rowは行（row）、colは列（column）というのはよさそうですが、mfが何の略だか分からない。「matrix figure」なんじゃないかと個人的には思っていますが、でたらめかも。コメントで教えてくれた人には10ガバスあげます。

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) 」と書いても同じ意味です。こう書いた方が、勘違いは少なそうですね。

やりたいのは↓こういうことだったりします。 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に助けられることがよくあります。数行のサンプルコードと、その実行結果があるだけで、相当助かるんですよね。 というわけで、私もちょっとしたサンプルコードとその実行結果の痕跡を残していこうと考えました。過去に自分がやったことも忘れちゃうんで・・・