ヒストグラムの書き方

ヒストグラムの基本

たとえば、正規乱数のヒストグラムを書くには、
> hist( rnorm(1000) )
とすればとい。これで、平均 0, 標準偏差 1 の正規乱数 1000個のヒストグラムを表示してくれる。 階級数をいくつにするかは悩ましい。、
> hist( rnorm(1000), n=100 ) # 階級の数を100にしよう
> hist( rnorm(100000), n=100 ) # データの数を増やしてみよう
階級数n の目安としては、 1 + log2(データ数)[Sturges (1926)]があるが、 私には、多くの場合、少なすぎるように思われる。 試行錯誤をして分布の特徴が分かりやすいものを選ぶと良い。 階級を細かく指定したいときには、breaksというオプションで階級の区分点を指定できる。
> hist( rbinom(200,10,0.4), breaks=-0.5:10.5 ) # (-0.5, 0.5, 1.5,...,10.5)を階級の区切りとして描画
離散分布の場合には、上のように整数を階級の中心とすると、誤解が少ないヒストグラムができる。 階級を指定した場合にどのような違いが生ずるか比較してみよう。
> par( mfrow=c(2,1) )      # 描画領域を2行1列用意する。 
> x <- rbinom(200,10,0.4)  # 同じデータを描画するため乱数をxに取っておく。
> hist( x ) 
> hist( x, breaks=-0.5:10.5 ) # (-0.5, 0.5, 1.5,...,10.5)を階級の区切りとして描画
なお、breaksを指定すると、範囲外の値がデータに含まれるとき、エラーになる。 指定するときには、適切に設定しよう。
> hist( rnorm(200,0,10), breaks=-0.5:10.5 ) 
 hist.default(rnorm(200, 0, 10), breaks = -0.5:10.5) でエラー: 
   'x' の一部分が数えられていません。多分 'breaks' が 'x' の範囲全体をカバーしていません 

描画の基本

表題や日付を図に付けたいときには、以下のようにすればよい。
> hist( rbinom(200,10,0.4), breaks=-0.5:10.5, main="200個の二項乱数 b(10,0.4)" ) 
> text(8,40,date())
ここでは、date()が現在日時の文字列を返してくれることを利用している。 text関数は第1引数の横座標、第2引数の縦座礁を中心として(オプションを与えれば変更できる)、第3引数の文字列を書いてくれる。 タイトルの描画を抑制しておいて、後で書き加えることもできる。
> hist( rbinom(200,10,0.4), breaks=-0.5:10.5, main="" ) 
> title( "200個の二項乱数 b(10,0.4)" )
色も使える。
> hist( rbinom(200,10,0.4), breaks=-0.5:10.5, main="",col=2)
colが色をしていするオプション名で、ここには上のように数字も書けるが、色の名も掛ける。
> hist( rbinom(200,10,0.4), breaks=-0.5:10.5, main="",col="red")
どんな色名があるかは、以下で列挙される。
> colors() 
色番号がどのようなものであるかは、以下のように確認できる。
> par( mfrow=c(1,1) )      # 描画領域を1行1列にもどす。
> n<-25; plot(1:n, 1:n, col=1:n, pch=1:n, cex=3)
>        text(1:n, 1:n-1, 1:n)
ここでは、プロットの記号(pch)も同時に変えている。記号の大きさ(cex)を大き目(3)に設定してみている。

経験分布に理論分布を重ね書きする

関数 hist() を用いて 理論分布と比較する時には probability=T なる オプションを付加すると各棒グラフの棒の(高さ)×(幅)が、それぞれの階級に サンプルが落ちた相対頻度(確率)となる。 ( なお、 probability=T のオプションを追加するときなどに、加えて ylimの指定をするときは、指定する値に注意する必要がある。 大きな区間を指定するとグラフが見えないことがある。) 以下に例を示す。

乱数の個数は固定して例示しているが、いろいろな個数でグラフを描画してみよ。 なお、以下では、cut&paste の便のため、Rのソースの先頭にプロンプト">"を付けない。

連続型分布の場合

data <- rexp(1000,10)           # 平均1/10の指数分布の場合
hist( data, probability=T )     # 確率に直した、頻度分布の作図
xx <- seq(min(data),max(data),length=1000)
                                # 密度関数値を調べるポイント
                                #  初項と末項、項数を指定している
lines( xx, dexp(xx,10) )        # 密度関数の直線の図への追加
title("指数乱数のヒストグラムへの理論確率密度関数の重ね描き")
なお、関数を描画する関数がある。以下でも同じ結果となる。
data <- rexp(1000,10)           # 平均1/10の指数分布の場合
hist( data, probability=T )     # 確率に直した、頻度分布の作図
curve(dexp(x,10),add=T)
title("指数乱数のヒストグラムへの理論確率密度関数の重ね描き")
なお、オプションadd=T は既存の図へ追記してもらうためのものです。

離散型分布の場合

data <- rpois(1000,3)                # 平均3のポアソン分布の場合

hist( data, breaks=seq(-0.5,max(data)+0.5), probability=T )
                                     # 確率に直した、頻度分布の作図
xx <- seq(min(data),max(data),by=1)  # 密度関数値を調べるポイント
                                     #  初項と末項、項差を指定している
lines( xx, dpois(xx,3) )             # 密度関数の直線の図への追加
title("ポアソン乱数のヒストグラムへの理論確率密度関数の重ね描き")

分布関数の描画の方法

理論分布関数は、以下などで描画できる。
curve(pnorm,from=-4,to=4)
curve関数のマニュアルとソースコードを確認しておこう。
?curve
curve
第1引数のexprは関数を受け取る引数になっている。次のような書き方もできる。
curve(pnorm,from=-4,to=4)
curve(pnorm(x,mean=2),from=-4,to=4,add=T,col=2,lwd=2)
これで、色が2(赤)で太さが2倍の線が追記(add=T)されている。
関庸一 Home Page へ