2008-10-21:Rで統計の問題を解こう - BioEoS::環境解析学研究室::

Rのインストール

後期第1回の研究室ゼミは，Rから始めます.
Rとは統計処理とグラフィックスを行うためのフリーソフトです.

ゼミの前にあらかじめR2.7.2を皆さんのパソコンにインストールしていただいていると思います．

統計の問題を解いてみよう

Rを起動させてください.

皆さんの手元に配られたプリントを見てください．
今日は課題2をRを使って解いてみます.

課題2
①n個のサイコロを振ったときに出た目(X)の和(X1+X2+X3+...+Xn)，または平均値(X1+X2+X3+...+Xn/n)のグラフをつくり，nの数が大きいほど正規分布に近づくことを確認する．
②0から9の乱数を(n個)発生させ，n個の平均のグラフをつくる.nが大きくなるにつれて正規分布に近づくことを確認する.

まず，課題2の②からやってみようと思います．
プリントでは0から9までとなっていますが，0から10の乱数を発生させるとします．

ヒストグラムを作ります．

• hist(runif (100))

ここで，hist()はヒストグラムを作るコマンド，runif()は乱数を発生させるためのコマンドで，runif()の括弧内の数字は繰り返す回数を示しています．

このように描画されて出てきます．

しかし，乱数は通常0から1の数字を示したものなので，0から10の乱数を発生させるには値を10倍する必要があります．

• hist(runif(100)*10)

更に，縦軸の数字の向きを見やすいように変えます．

• par(las=1)

もう一度，描画させて見ましょう．

• hist(runif(100)*10)

縦軸の数字の向きが変わったことがわかると思います．

では，グラフのところどころを変えてみます．

• hist(runif(100)*10,xlab="数値",ylab="頻度",main="(0,10)の乱数100個の頻度分布")

ここで，xlabはx軸のラベル，ylabはy軸のラベル，mainはタイトルを表しています．

ここで，今100となっている数値を1000.10000といった具合に変えてみると
さらに細かいグラフになります.

histのヘルプを見たいときは，

• help(hist)

縦軸の頻度を600にまで増やそうと思えば…

• hist(runif(100)*10,xlab="数値",ylab="頻度",main="(0,10)の乱数100個の頻度分布")

にコンマをうって続けて

• ylim=c(0,600)

とします．

ylimとはy軸のリミット(下限，上限)という意味です．

配列

「means」という名前の配列の箱をつくります.名前は他のものを使っていただいてもかまいません．

• means <- numeric(100)

このコマンドで100個の棚を持つ配列の箱ができました．

• means

上のように入力すると，0という値が100個表示されると思います．

Rにおいて平均値を求めるコマンドを「mean」といいます．

• for(i in 1:100)
• {means[i] <- mean(runif(100)*10}

上のコマンドは「0から10までの乱数を100個作り，その平均値をi番目のmeansの値として入力する」という意味になります．
うまくできているか，確認してみましょう．

• means

これで，値が100個できていれば大丈夫です．

では次に，この値を描画してみましょう．

• hist(means,xlab="数値",ylab="頻度",main="(0,10)の乱数の1000個の平均値の頻度分布")

次に，上で作成したヒストグラムが本当に正規分布しているのかを調べるために，正規分布の式を使って線を引いてみましょう．

まず，x軸刻みを0.1にするために，x軸の設定をします．

• xv <- seq(0,10,0.1)

次に，y軸を計算します．
Rには正規分布を作成するためのコマンド「dnorm」があるので，それを用いて計算させます．

• yv <- dnorm(xv,mean=mean(means),sd=sd(means))

ここで，means(means)とはmeansの平均値を，sd(means)はmeansの標準偏差を表しています．

では，これを描画させてみましょう．

• lines(xv,yv)

さて，ではどれほど経験分布(ヒストグラム)と理論分布(yv)が違っているのかを見るために，表示させてみましょう．
その際，ヒストグラムと線の違いがよくわかるように，色を変えてみましょう．
ヒストグラムを黄色，線を赤にしてみます．

• hist(means,xlab="数値",ylab="頻度",main="(0,10)の乱数の1000個の平均値の頻度分布",col="yellow")
• lines(xv,yv,col="red",lwd=3)

ここで，lwdは線の太さを表しています．

更に見やすくする為に凡例をつけましょう．
凡例をつけるコマンドはlegendといいます．

• legend(4,150,legend=c("経験分布","理論分布"),col=c("yellow","red"),lwd=5)

Rを終えたいときは，

• q()

今回はこれで終わりです．お疲れ様でした．