> 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)されている。