投稿

4月, 2015の投稿を表示しています

Rのhistとggplot(geom_histogram)を比較する

イメージ
ggplot2パッケージって使ったことありませんでした。記法がとっつきにくいというのもあったのですが、あの見た目を素直に「きれい」とは認めたくなかったというか。「ggplotの方がかっこいい」という人への反発みたいな部分も多分にあるのですが。 「白地に線画」みたいなビジュアルが好きなんです。下手に色をつけたらかえって趣味が悪くなったりするじゃないですか。自分自身センスがないことを自覚しているので、色やら塗りつぶしやら飾りっぽいものやらは極力使わないようにしていたんですよね。 でも今回、ggplot2パッケージを使ってみようと思ったのは、ベースとなるレイヤーを用意して、次に重ねたいレイヤーを(場合によっては複数)用意して、最後にバンっとplotする、っていうコンセプトがなんだか使いやすいのではないかという気がしたから。なので、徐々に試していこうかなと思っています。 前置きが長くなりました。 今回はヒストグラムですが、多くの人に馴染みのあるデフォルトのhistでの描画と、ggplot2パッケージを使った場合とで、どんな風に違うのか試してみました。 正規分布に従う乱数を発生させて、普通にヒストグラムを描いてみます。 data1 <- rnorm ( 1000 )   hist ( data1 ) 通常のhist関数で描画 ggplot2パッケージを使うと↓こんな描き方になります。 install.packages ( "ggplot2" ) library ( ggplot2 )   df1 <- data.frame ( data1 ) # データフレームにしておく必要がある   g <- ggplot ( df1 , aes ( x=data1 ) ) # 対象がdata1列であることを指定 g <- g + geom_histogram ( ) # ヒストグラムを描画することを指定 plot ( g ) # 描画 ggplot2パッケージを使ったもの 先ほどの、デフォルトのhistと階級幅を合わせる(=0.5)には、binwidthを指定します↓ g <- g

Rでブートストラップ法

イメージ
「ブートストラップ」という言葉は聞いたことはあるけど、具体的にどんなことなのかはよく知らない、というのが私のレベルです。 が、「統計学入門」(東京大学教養学部統計学教室 編)という教科書↓ 統計学入門 (基礎統計学) の中の演習問題に出てきた(でも、教科書の本編には解説はない)ので、Rでちょっとやってみました。 第3章 練習問題 3.4 <ブートストラップ> 以下の手法をパソコンで行え。 i)  [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。 ii)  i)を11回繰返して実行し、11個のランダムな番号を得たのち、表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、相関係数rを計算せよ。 iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法を ブートストラップ という。 「表3.1、図3.1」というのは、教科書の本編のところに出てきた表と図で、相関関係が認められる例として紹介されていたものです。 統計学入門 (基礎統計学) より 統計学入門 (基礎統計学) より このデータを使う必要があるので、まずは入力しましょう。 > bro <- c(71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66) > sis <- c(69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62) 書籍に載っていた図と同じような感じでプロットしてみました↓ > plot(bro, sis, xlim=c(68,70), ylim=c(59,70), asp=1, pch=16) 書籍に載っていたデータをRでプロットしたもの 書籍に出てくる図とほぼ同じものができました。ただ、点が10個しかなく1個少ないです。これは(70, 65)という全く同じデータが2つあって同じ座標にプロットされているためです。書籍の方では少しずらして、分かるようになっていました。 あと、書籍だと(64,62)あたりに点がありますが(一番左側のもの)、これは間

Rで y = x のグラフを描きたいとき

イメージ
Rでグラフを描きたいときにはcurveを使いますよね。 y = x^2 のグラフなら↓こんな感じ curve(x^2, -1, 1) じゃあ、y = x のグラフのグラフは、↓これでいけるかと思いきや・・・ curve(x, -1, 1)  以下にエラー eval(expr, envir, enclos) :    関数 "x" を見つけることができませんでした xだけじゃ、関数だと思ってくれなかったんでしょうか? とりあえず、↓こうすれば大丈夫ではあります。 f <- function(x){x} # y=x の関数を定義しておく curve(f, -1, 1) でも、なんか面倒臭いですよね。 で、実は括弧で囲めばOKということに最近気付きました。 curve( ( x ) , -1, 1) しかしまあ、「y = x」のグラフを描くことって、滅多にないでしょうけど。

Rの数式描画の中で変数を展開したいときはbquoteを使う

イメージ
ループを使って複数のグラフを描くような状況。 par ( mfrow= c ( 2 , 2 ) ) # 2行2列で作図する   for ( i in 1 : 4 ) { curve ( x^i , - 1 , 1 , asp= 1 , ylab= "y" ) } 複数のグラフを描いてみた(タイトルなし) グラフを区別するためにmainオプションでタイトルを付けましょう↓ par ( mfrow= c ( 2 , 2 ) )   for ( i in 1 : 4 ) { curve ( x^i , - 1 , 1 , asp= 1 , ylab= "y" , main= paste ( "y=x^" , i ) ) } タイトルをつけてみたが格好悪い う~ん、伝わりはするけど、"y=x^2"みたいな見た目が格好悪い。ちゃんと数式(数字の肩に指数が乗っている感じ)で表現したいですよね。 そういうときは、expression関数で、いけるかな・・・? par ( mfrow= c ( 2 , 2 ) )   for ( i in 1 : 4 ) { curve ( x^i , - 1 , 1 , asp= 1 , ylab= "y" , main= expression ( y==x^i ) ) } expressionだと変数iが展開されない ありゃ? 数式の形にはなりましたが、変数iが展開されずそのまま「i」と出ちゃいました。 やっと本題です。 このiを展開(i = 1, 2, 3, 4)したいときには、bquote関数が使えます↓ par ( mfrow= c ( 2 , 2 ) )   for ( i in 1 : 4 ) { curve ( x^i , - 1 , 1 , asp= 1 , ylab= "y" , main= bquote ( y==x^. ( i ) ) ) }

ベンフォードの法則をRでシミュレーションする

イメージ
電気料金の請求書の金額の1桁目をたくさん集めて、その分布を見たらどうなっているでしょう。 数字の出方に法則性があるようにも思えないから、1から9まで一様になっているのでは? という気もしますが、そうではないというのがベンフォードの法則です。 ベンフォードの法則 - Wikipedia ベンフォードの法則は、自然界に出てくる多くの(全てのではない)数値の最初の桁の分布が一様ではない、ある特定のものになっているというものである。 (中略) この直感に反するような結果は、電気料金の請求書、住所の番地、株価、人口の数値、死亡率、川の長さ、物理・数学定数、冪乗則で表現されるような過程(自然界ではとても一般的なものである)など、様々な種類の数値の集合に適用できることがわかっている。 どのように偏って分布しているかというと、↓こんな割合になるとのこと。 ベンフォードの法則 - Wikipedia より 本当にそうなるか、Rでシミュレーションしてみましょう。 数値が満たすべき分布について、ウィキペディアには、 論理的には、数値が対数的に分布しているときは常に最初の桁の数値がこのような分布で出現する。 とあります。いまいちピンとこなかったのですが、あわせて掲載されていた図でなんとなく理解しました。 ふむふむ、対数スケールのグラフ上に、見た目で一様になるように分布すればいいのね。 理解を助けるために、まず、対数の数直線を書いてみましょう。 x <- c ( 1 , 10 , 100 , 1000 ) y <- rep ( 0 , length ( x ) ) plot ( x , y , axes=F , xlab= "" , ylab= "" , log = "x" ) abline ( h= 0 ) text ( x , y , x , pos= 1 ) runif関数で発生させた乱数をそのまま使っちゃだめですよね。それでは数値的な一様になってしまいますから。 数直線に指数表記を追加すると分かりやすいです。 a <- round ( log10 ( x ) , 3 ) # 指数の肩の数字

Rで最上位桁を取り出す

なんでこんなことをしたいのかは、 また別の機会に書く として、やりたいのは数字の最上位桁の数字を取り出したいというものです。123だったら1、31415だったら「3」を取り出すという感じです。 整数商を使えばいけそうです。123から1を取り出すには100で割る、31415から3を取り出すには10000で割る。   123 ÷  100 =  1  余り  23  31415 ÷ 10000 =  3  余り 1415 割る数は、10×(桁数-1)ですね。↓こんな感じでしょうか。 n <- floor(log10(x)) # 桁数-1 msd <- x %/% 10^n    # 整数商 ちなみに、MSDというのはMost Significant Digit(最上位桁、最上位の数字)の略です。 正しく結果が出ているか、乱数で確認してみましょう。 set.seed(1) # 再現性のために x <- sample(1:1000000, 20) n <- floor(log10(x)) # 桁数-1 msd <- x %/% 10^n    # 整数商 data.frame(x, msd) # 並べて比較するためにデータフレームに         x msd 1  265509   2 2  372124   3 3  572853   5 4  908206   9 5  201682   2 6  898386   8 7  944670   9 8  660794   6 9  629110   6 10  61786   6 11 205973   2 12 176555   1 13 687015   6 14 384099   3 15 769831   7 16 497692   4 17 717608   7 18 991890   9 19 380029   3 20 777431   7 うまくいっているようですね。 でも、最上位桁を求めるために、常用対数とfloorと整数商を動員するのか。なんだか仰々しいというか。 要は1文字目を取り出せばいいんでしょ、ってところは後から気づきました。 ↓文字列型に変換して、1文字目を取り出して、それを整数型に変換する、という考え方です。

Rの正規表現で最短マッチ(最短一致)させる

特にRに固有のことというわけではありませんし、正規表現に詳しい人にとっては当たり前のことだとは思いますが、私自身が結構つまづいたので・・・ ↓こんな感じのベクトル(住所リストのイメージ)があったとします。 住所リスト <- c( "佐藤町中村1-2",                  "高橋町新田2-3",                  "佐々木町本町3-4",                  "伊藤町原4-5",                  "山本町町田5-6",                  "長谷川町本郷6-7" ) 住所リスト [1] "佐藤町中村1-2"   "高橋町新田2-3"   "佐々木町本町3-4" [4] "伊藤町原4-5"     "山本町町田5-6"   "長谷川町本郷6-7" で、それぞれから冒頭の「○○町」の部分を取り除きたいとします。 町名の長さも一律ではないので、こういうときは正規表現が便利ですよね。 私が最初に思いついた正規表現が以下です。 sub(" ^.+町 ", "", 住所リスト, useBytes=F) 冒頭( ^ )から、任意の文字( . )が、1つ以上( + )あって、「 町 」で終わるような文字列、というつもりです。 ですが、うまく行きませんでした。↓実行してみると・・・ sub("^.+町", "", 住所リスト, useBytes=F) [1] "中村1-2" "新田2-3" "3-4"      "原4-5"   "田5-6"    "本郷6-7" 町名より下のレベルの住所に「町」という文字を含んでいると、そこまで余分にマッチしてしまったというわけです↓ " 佐々木町本町 3-4" " 山本町町 田5-6" これは正

Rで最頻値(mode)を求める ~連続値データの場合~

イメージ
Rで最頻値を求める方法として、よく見かけるのがtable関数を使った方法です↓ names(which.max(table(x))) 離散値データの場合は上記でうまくいくのですが、連続値データの場合はうまくいきません。 また離散値データであっても、取りうる値の種類に対してサンプルサイズが小さい場合もうまくいきません。 具体例をあげてみます。40人学級でのテストの点数をイメージしました。 # テストの点数っぽいサンプルデータを作る set.seed(4) # (やや恣意的な)再現性のために x <- round(rnorm(n=40, mean=50, sd=15)) x  [1] 53 42 63 59 75 60 31 47 78 77 58 50 56 49 51 53 67 49 48 46 [21] 73 52 70 69 59 46 69 64 36 69 52 66 39 28 63 44 47 64 43 40 上記データに対して、tableを使って最頻値を求めると、 names(which.max(table(x))) [1] "69" 想定している、"50"とかなりずれた値が出てしまいました。実現値がスパースなので、個別に最も現れた値を出すことに意味がないんですね。 ↓階級を1点刻みにしたヒストグラムで見るとよく分かります。 # 1点刻みでヒストグラムを描く hist(x, breaks=seq(20,80,1)) 69点を取った人は3人で最も多いですが、たまたま感がありますよね。50点の周辺の方が頻度が高いように思えます。 すなおに階級に幅を持たせると、予想と近い最頻値が得られます。 # 階級幅はデフォルト(hist関数にまかせる) hist(x) 上記の場合、40~50の階級が最も度数が多いので、この中間値を階級値として「最頻値は45である」という結論は直感に反していないと思われます。 階級値をきりのいい数字にしたい場合は、↓breaksを指定してやればOKです。 hist(x, breaks=seq(25, 85, 10)) 最頻値は50の階級であるという、予想通りの結果となりました。 つまり、連続値(や

所得の分布は対数正規分布に従っているのか?

イメージ
日本の世帯ごとの所得をヒストグラムにすると↓こんな感じになります。 所得金額階級別にみた世帯数の相対度数分布(平成21年調査) 2 所得の分布状況|厚生労働省 より 右側に裾が長く延びているので、正規分布ではないようです。 よく、所得は対数正規分布に従う、なんて聞いたりします。所得額の対数をとってそれで度数分布表を作るような感じですね。 実際に試してみましょう。 グラフで示されていた値を、データ化します。 rate <- 0.01 * c(6.6, 12.7, 13.9, 13.3, 10.0, 8.9, 7.1, 6.2, 5.1, 3.9,                  3.0,  2.1,  1.7,  1.3,  0.9, 0.9, 0.5, 0.4, 0.2, 0.2) income <- seq(50, 1950, 100) df <- data.frame(income, rate) df    income  rate 1      50 0.066 2     150 0.127 3     250 0.139 4     350 0.133 5     450 0.100 6     550 0.089 7     650 0.071 8     750 0.062 9     850 0.051 10    950 0.039 11   1050 0.030 12   1150 0.021 13   1250 0.017 14   1350 0.013 15   1450 0.009 16   1550 0.009 17   1650 0.005 18   1750 0.004 19   1850 0.002 20   1950 0.002 「2000万円以上」の階級はオープンエンドになっているので省きました。所得額の中央の値を階級値(単位:万円)として、rateはパーセントでなく割合で表しています。 そのままplotしてみると↓ plot(df$income, df$rate) 所得ごとの度数をプロット 当然、厚労省のページに載ってたような形が出てきます。 データフレームに、対数をとった所得の列を追加してみます