|
「R」の基本操作
|
ここでは，統計データを取り扱うために必要な「R」の基本操作について，
備忘録的に取りまとめました．
「R」のインストールを含む全般的な解説は他の解説サイトを参照してください．
使用するのは64bit版の「R」(ver. 3.6.2)です．
|
■基本的な演算
ここでは，「R」に関する基本的な演算についてまとめました．
通常の四則演算は下記の通りです．円周率 \(\small \pi\) は「pi」．
平方根は「sqrt()」です．「R」での代入には「＜-」(＜は半角)が用いられるが，
「=」でもかまいません．個人的には「＜-」より「=」がよいと思いますが，
「R」の通常の書式で書いておいた方がよいと思います．
ただし，「＜-」を用いるときは，
不等号やマイナス記号との混乱を避けるため前後に空白を入れた方がよいです．
複数のコマンドは「;」で区切って入力できます．
ネイピア数は「exp(1)」を適当な文字に割り当てておくとよいでしょう．
「( )」で括ると，代入と出力が同時になされます．
四捨五入は「round」，切り捨ては「trunc」，
その値を超えない最大の整数は「floor」です．
ベクトルは「c(a,b,c)」のような感じで入力します．
数列は「seq(start,end,step)」，
連続数は「rep(start:end)」であり，いずれもベクトルとして生成されます．
1刻みの場合は「start:end」だけでもかまいません．
ベクトルを結合したり，共通部分をとったりすることもできます．
- 「append(a,c)」により，「a」の後に「c」が追加されます．
- 「append(a,c,after=2)」により，「a」の2番目の要素の後に「c」が追加されます．
- 「union(b,c)」により，「b」と「c」の和集合が生成されます．
- 「intersect(b,c)」により，「b」と「c」の共通部分が生成されます．
- 「setdiff(a,c)」により，「a」から「c」の要素を除いたベクトルが生成されます．
下図は，上記の操作を行ったものです．
なお，要素を昇順に並べ替えるには「sort(v)」，
降順は「sort(v,decreasing=TRUE)」とします(図は略)．
複数のベクトルをまとめて一つのデータとすることもできます．
長さが同じベクトルを「a」「b」とするとき，
「rbind(a,b)」により2つの行(row)をまとめることができ，
「cbind(a,b)」により2つの行を列(column)にしてまとめることができます．
定義した変数の一覧は「ls()」で表示されます．
変数を消去するには，「remove(変数)，または「rm(変数)」とします．
下図は，表示画面の左半分を切り出したものです．
変数「a,b,c,d」を消去して再度「ls()」を実行しました．
「a,b,c,d」は表示されないので削除されたことが分かります．
|
|
■画面の操作
「R」では空白は無視されますが，
コマンドの綴りは一続きで入れる必要があります．
1行に複数のコマンドを書くには「；」で区切ります．
コマンドが終了しないうちに「Enter」を押すと，
続きの入力を促す「＋」が表示されます．
画面を全部消去するには「ctrl」＋「l」(エル)，または，
トップメニューの「編集」から「コンソール画面を消去」を選択します．
アルファベットの大文字と小文字は区別されるので，
入力には注意が必要です．
コンソール画面では1行ずつ入力しなければならず，煩わしいときがあります．
そのようなときは，コマンドの流れをスクリプトに書いて一気に実行させる
ことができます．
新たなスクリプトの作成は，
「ファイル」から「新しいスクリプト」を選択します．
作成済みのスクリプトを追加・修正するのであれば，
「スクリプトを開く」を選択します．
「Rエディター」が開くので，そこにコマンドを記述します．
実行は，実行するスクリプトを選択して青地表示にしてから
「ctrl」＋「r」を押します．
同様にして，スクリプトの一部分だけ実行させることもできます．
下図は，右側部分を選択して「ctrl」＋「r」を押した状態です．
スクリプトを書いたファイルは「R」の作業フォルダー(working directory)
に保存しておくとよいでしょう．
現在どのフォルダーが作業フォルダーになっているかは「getwd()」で確認できます．
変更するときは，表示されるフォルダーの書式にならって，
「setwd("フォルダー名")」で指定します．
スクリプトに日本語を使用することもできます．
後で振り返るとき，コマンドのいろいろな解説を書いておくとよいでしょう．
行頭に「#」をおくと，その行はコメント行になります．
|
|
■関数の扱い
通常の初等関数を使用することができます．
| sqrt | 平方根 |
| abs | 絶対値 |
| exp | 指数関数 |
| log | 自然対数 |
| log10 | 常用対数 |
| sin | 正弦関数 |
| cos | 余弦関数 |
| tan | 正接関数 |
| floor | ガウス記号 |
| ceiling | 天井関数 |
| trunc | 切り捨て |
| round | 四捨五入 |
双曲線関数は「sinh」「cosh」「tanh」，
逆三角関数は「asin」「acos」「atan」，そして
逆双曲線関数は「asinh」「acosh」「atanh」です．
また，ガンマ関数は「gamma」です．
新しい関数を作るには，
「(関数名) ＜- function(引数){処理)}」(＜は半角)という書式になります．
最後の処理結果は「return(結果)」とします．
たとえば，\(\small f(x)=2(x-2)^2+1\) を細かく分けて定義して
「fun」とするには，次のようにします．
定義した関数を忘れてしまったときは，その定義を再表示させることもできます．
単に関数名を打ち込むだけです．
|
|
■条件分岐・繰り返し
条件分岐の「if-else」は，「if (条件) {処理1} else {処理2}」とします．
たとえば，
\[\small f(x)=
\begin{cases} -1 & (x\lt -1)\\ ~~0 & (-1\leq x\lt 1)\\ ~~1 & (1\leq x)\end{cases}\]
を定義するには，次のようにします．
「＜-」は代入記号なので，
その前後には空白をおくように習慣づけた方がよいでしょう．
また，途中のスペルミス等があったりするので，コンソール画面ではなく，
スクリプトに書いて「ctrl」＋「r」により実行した方がよいでしょう．
プログラム内で用いた変数の値は外部では保持されません．
条件記述のときの真偽判定の論理演算子として，
「等号」は「＝＝」，「不一致」は「！＝」，
「かつ」は「＆＆」，「または」は「||」です．
単に結果を返すだけのときは「ifelse」という書式もあります．
「ifelse(条件，真のときの結果，偽のときの結果)」とします．
これを入れ子にすれば，上記と同じ関数は次のように記述することもできます．
繰り返し処理に「for」が利用できます．書式は
「for (変数 in 範囲) {処理}」です．
たとえば，階乗は次のように計算されます．
任意の実数を代入できて，負の場合は0が，正の場合は整数部分までの階乗の
値が返されます．ただし，\(\small 0!=1\) に関する処理は行っていないので，
実用的ではありません．「return( )」の位置に注意が必要です．
条件が真である間で繰り返すには「while (条件) {処理}」とします．
上記の階乗は，次のようにも定義できます．
同じような繰り返しに「repeat」も利用できます．
書式は「repeat {処理}」です．
ある条件を満たして繰り返しから抜けるには「break」を置きます．
repeatを利用すると，階乗は次のようにも定義できます．
プログラムの簡単のため \(\small 0!=1\) に関する処理は行って
いないので実用的なプログラムではありません．
実用に付すには，その場合の特別処理を追加記述すればよいのです．
|
|
■データの型
数は実数の他に複素数を利用することができます．文字は「" "」で囲います．
虚数単位 \(\small i\) は「1i」として利用できるので，
複素数 \(\small a+bi\) は「a+b\(\small *\)1i」とします．
数や文字を横に並べたものをベクトルと呼びます．
コンマで区切った数をベクトルとして扱うには「c(a,b,c)」を利用します．
規則正しい等差数列は「step(start,end,by=step)」とします．
そのようなベクトルに対する四則演算は，対応する要素どおしの演算になります．
なお，下図は b/a の結果の前半だけを切り取ったものです．
(全部を表示するとスマホでの画面が見づらくなるので．)
単純に一定数だけ増えたり減ったりする数列は簡単に作れますが，
一般的な数列はプログラムを組んで作るしかないと思われます．
たとえば，フィボナッチ数列は次のようになります．
いろいろなベクトルを一つにまとめたものはリストと呼ばれます．
そこでは，数値，文字，論理記号などをまとめて一つのもの(オブジェクト)
として扱うことができます．ベクトルを並べて行列(matrix)とすることもできます．
行列は行と列だけの2次元ですが，
配列(array)を利用すると多次元のものを考えることができます．さらには，
各列の要素が異なる型のデータも扱うことができますが，
ここでは統計解析で必要とするCSVデータの取り扱いに限定して考えます．
|
|
■データの入出力
統計データは列ごとにみます．1行目に列の名前を置くこともできます．
コンマ区切りのデータは「csv」です．
たとえば，
下記のデータ(sample.csv)を読み込むことにします．
csvファイルを読み込むには
「read.csv("filename", header=TRUE)」
または「read.table("filename",header=TRUE,
sep=",")」とします．
1行目に列名がないときは「header=FALSE」とします．
その場合は，各列に「V1, V2, V3」などの列名が自動的に割り振られます．
「read.table」を利用すると，
空白やタブなど多様な区切り記号が利用でき，
いろいろな読み取り方をさせることができますが，
ここでは詳細は省略します．
データ名を「data」とするとき，
「head(data)」とすると出だしの6行が，
「tail(data)」とすると最後の6行が，
「data」と入力すると全データが表示されます．
データ数が多いときは表示のさせ方には注意する必要があります．
特に必要がない限り，headでの表示を習慣づけた方がよいかもしれません．
読み込んだデータの列名だけを表示させるには「names( )」，
列数は「ncol( )」，行数は「nrow( )」とします．
行数は，列名があるときは実データの列数になります．
データの列を指定するには「dataname＄列名」(＄は半角)とします．
このデータ(data)の3列目は「Den」なので「data＄Den」とします．
列名が付されていないときは「data＄V3」とします．
「table(dataname＄列名)」とすると，
そのデータが小さい順に並べ替えられて，
同じ値を持つ度数が表示されます．
離散的な値の場合は度数分布表が表示されることになります．
i行目やj列目を表示させるには「data[i,]」や「data[,j]」とします．
i列とj列を表示させるには「data[,c(i,j)]」とします．
ただし，「head」や「tail」で囲って表示させた方がよいでしょう。
なお，下図では(スマホ表示での見やすさのため)画面右側の一部を省略しました．
|
|
|
ファイルをコンマ区切りのファイル(csv)で保存するには，
「write.csv(データ名, "filename")」とします．
下図は，「write.csv(data,"file.txt")」とした場合です．
上の5行分を示しました．
これから分かるように，単純に保存すると，
列名や行名を引用符で囲んで保存されます．
引用符をつけないようにするには「quote=F」，
そして行名を省略するには「row.names=F」を追加指定します．
ただし，列名を省略する「col.names=F」という指定は，
「write.csv」では無視されて列名が強制出力されます．
これは「R」の仕様のようですが，
「write.csv」というコマンド自体を書き換えて「col.names=F」を
有効にすることもできるようです(
参照)．
この列名も省略して単なるデータ列だけを保存するには
「write.table」を利用します．
その場合は，コンマ区切りであることを指定するために「sep=","」も
追加指定する必要があります．
|
|
■列データの処理
|
サンプルデータ
統計ソフトであるので，
読み込んだデータにいろいろな処理を行うことができます．
ここでは，
サンプルデータとして「sample2.csv」を使用します．
このデータは，Maximaを利用して50個の乱数を発生させたもので，
各列の内容は次の通りです．
| 1列目： | 1〜6の整数からなる乱数． |
| 2列目： | 二項分布Bin(10,1/6)による乱数． |
| 3列目： | 正規分布N(50,10)による乱数． |
| 4列目： | パレート分布Pare(1,1)による乱数． |
| 5列目： | 文字列{A,B,C,D}のランダム抽出． |
なお，正規分布は四捨五入して整数としました．
パレート分布は \(\small 1/x^2~(x\geq 1)\) にしたがう乱数で，
小数第1位までの値としました．
下図は，このファイルを読み込んで先頭の6行を出力させたものです．
1行目に列名があるときは，
「header」による真偽の指定をしなくても列名ごと読み込まれます．
|
列データの操作
それぞれのデータは列データなので「data2」は行列形式のデータです．
したがって，\(\small (i, j)\) 成分は data2[i,j] により参照することができます．
第i行は data2[i, ]，第j列は data2[ ,j] です．
行は列名の下にあるデータ行です．
下図では，(2,5)成分，3行目，そして4列目が出力されています．
1列だけの場合は，横並びのベクトルとして出力されます．
特定の列(たとえば，3列目と4列目)を取り出すには，
列の箇所で c(3,4) を指定します．
逆に，特定の列だけ除くには -c(3,4)を指定します．
行についても同様の指定ができます(図は略)．
いずれも，
「head」をつけないと50行全部が表示されます．
この指定で，c(3,4) をたとえば c(4,3) にすると，
4列目が最初に表示されます．下図では，4列目と2列目を取りだしています．
特定の列を1列だけ取り出すとベクトルとして表示されますが，
列として取り出すには「drop=F」を追加指定します．
下図では2列目の data2$binom を取り出しています．
以上の処理を組み合わせると，幾つかの列データを組み合わせて
新しいデータ行列を生成することができます．
それを行う関数は「cbind」です．
下図では，4行目，2行目，そして3行目からなるデータ行列を data3 として
定義したものです．
ただし，これは data2 の列を入れ替えているだけなので，
実際には，data2[,c(4,2,3)] だけで済みます．
列ではなく行ベクトルを結合するには「rbind」を利用します
．
|
|
基本統計量
データの統計的処理は「suumary」で，
すべての列について一通りの値が出力されます．
下図は，3列目の途中までを切り出したものです．
各列について，
最小値，第1四分位数，メディアン，平均，第3四分位数，そして最大値の順に
表示されます．
列名を指定すると，指定した列についての結果が得られます．
数値データの場合，2列以上の総和は colSums，平均は colMeans
により求めることができます．行については rowSums, rowMeans です．
下図では，文字データである5列を除いた総和や平均を求めています．
この書式で分散を求める関数はないので，必要であれば自分で作ることになります．
個別の統計量として，総和は「sum」，平均は「mean」，
不偏分散は「var」，標準偏差は「sd」，
最小値は「min」，最大値は「max」，分位数は「quantile」，
範囲は「range」，そしてメディアンは「median」です．
ただし，colSums, colMeans を除くと，
標準的な利用では複数列の結果を一度に得ることはできません．
そのような処理が必要なときは，
それを出力させる関数を自分で作ることになります．
または，そのような関数が登録されているパッケージ
(たとえば「fBasics」)を
読み込む必要があります．
また，不偏分散は標本分散ではないので注意が必要です．
標本分散が必要なときは不偏分散を自分で \(\small (n-1)/n\) 倍します．
分位数は％点の値が表示されます．下図では，第3列に対して
平均，分散，そして80％点を求めています．
|
|
度数分布・クロス表
度数分布表は，「table」を利用すると簡単に求めることができます．
下図は，離散データである第2列の場合です．
第1列は1〜6の整数をランダムに出力したものです．
このデータでは「6」が偏って出力されていますが，
文字列からなる第5列と第1列とのクロス表も簡単に出力されます．
第3列は連続分布である正規分布N(50,10)のデータです．
扱いやすいように四捨五入して整数値にしていますが，
このデータを階級に分けて度数分布を作成することも容易です．
個々のデータがどの階級に属するかで区分けするための関数として「cut」があります．
「cut(データ, 区分点のベクトル)」とすると，
個々のデータが該当する階級 (a,b] で置きかえられます．
下記では，最初に data2＄norm の出だしの6行を確認し，
次に「seq」を利用して[10,90] の範囲を10刻みに分け，
その区分点をまとめたベクトルを「kubunten」としています．
そして，その区分点で data2＄norm のデータを区分けして「norm_kubun」とします．
個々のデータが，その値が属する階級 (a,b] で置きかえられています．
その階級データに対して「table」を利用すると，
階級ごとの度数分布表が得られることになります．
下図では右半分は省略されています．
|
データの標準化
データを標準化する関数は「scale」です．
データを指定すると，((データ)-(平均))/(標準偏差) により，
平均が0, 標準偏差が1のデータ列に変換されます．
|
|
■データのグラフ化
1変量のデータ
データ x は，「plot(x)」によりグラフ化することができます．
このとき，データのタイプにより次のようなグラフになります．
- x がベクトル型で，要素がすべて実数のときは，
横軸は何番目かを表す自然数，縦軸には要素の値がプロットされます．
- x の要素が複素数のときは，横軸を実部，縦軸を虚部とするプロットになります．
- x が2列の行列型データのときは，1列目を横軸，2列目を縦軸として
プロットされます．
下図は，data2＄norm のデータをプロットしたものです．
単に「plot」すると，データが並んでいる順に点が丸印でプロットされます．
範囲は自動で設定されます．
配置する点の種類は「点・線種・色」を参照のこと．
横軸の範囲は xlim=c(start,end)，縦軸の範囲は ylim=c(start,end) の形で
指定することができます．
丸印を線で繋ぐには「type="o"」を指定します．
いろいろな結び方は「点・線種・色」を参照のこと．
対数軸とするには，「log="x"」「log="y"」「log="xy"」などを追加指定します．
下図は，data2＄pare (＄は半角) で
縦軸を対数軸にとり，さらに横軸までの垂線を引いたものです．
垂線は，「type="h"」を指定します．
|
1変数関数のグラフ
「R」に登録されている標準関数は，
関数の名称を指定するだけでグラフを描画することができます．
たとえば，正弦関数 \(\small \sin(x)\) のグラフを描くには，
関数の名称と範囲を指定するだけです．
ただし，「sin(x)」とするとエラーとなるので注意が必要です．
範囲は「xlim=c(-pi,pi)」として指定することもできます．
\(\small y\)軸の範囲を指定する場合は，「ylim」について同様にします．
\(\small f(x)\) のグラフを描くには，その関数を定義した上で，
「f」を関数の名称として使用します．「x^3-3\(\small *\)x」 のような
関数の式を直接入力するには「curve」を利用します．
たとえば，\(\small f(x)=x^3-3x\) であれば次のようにします．
グラフを重ね描きするには「add=TRUE」(または「add=T」)を指定します．
線種(line type)は「lty」で指定できます．実線は「lty=1」であり，
デフォルトでは省略されます．サインカーブを実線で，
コサインカーブを破線(lty=2)で描画するには次のようにします．
線種の細かい指定は「点・線種・色」を参照のこと．
グラフを重ねると，\(\small y\) 軸のラベルが重なってしまいます．
軸のラベルは，横軸は「xlab」, 縦軸は「ylab」により指定できます．
たとえば，上記のグラフでコサインカーブのラベルを表示しないように
するには，コサインカーブの箇所で「ylab=""」を指定します(図は略)．
「R」では，グラフ描画に関する多数の関数を利用することができます．
「plot」の他に「curve」もあります．書式は「plot」と同様に
「curve(関数，左端，右端)」ですが，関数の式を直接入力することができます．
下図は，2つの関数 \(\small \sin(x), x-\frac{x^3}{6}\) を重ね描きしています．
グラフを重ね描きするには「add=TRUE」とします．
単に「add=T」だけでもかまいません．
曲線の色は「col」で指定します．下図では2番目の曲線を青(blue)としました．
色の細かい指定は「点・線種・色」を参照のこと．
|
|
点・線種・太さ・色等
離散データの場合に配置する点は，デフォルトでは丸印です．
配置する印は「pch」で指定することができ，デフォルトの丸印は「pch=1」です．
主に使用すると思われる記号として下記のようなものがあります．
なお，7〜14に間には0〜6の記号を組み合わせた記号が登録されています．
| pch | 記号 |
pch | 記号 |
pch | 記号 |
| 0 | □ |
4 | × |
16 | ● |
| 1 | ○ |
5 | ◇ |
17 | ▲ |
| 2 | △ |
6 | ▽ |
18 | ◆ |
| 3 | ＋ |
15 | ■ |
デフォルトでは点が配置されるだけです．点どおしを線で結ぶかどうかは
「type」で指定します．デフォルトは点を置くだけなので「type="p"」
を指定したことになります．
線(line)で結ぶには「plot(norm,type="l")」(エル)とします．
印を置いて線で結ぶには「type="o"」(オー)，
垂線を引くには「type="h"」とします．
下図は，この順に表示されるグラフです．
type="l"
type="o"
type="h"
デフォルトでは実線で描画されます．
線種(line type)は「lty」で指定し，デフォルトの実線(solid)は「lty=1」です．
「lty」の値に応じて，破線(dashed,lty=2)，点線(dotted,lty=3)，
鎖線(dotdash,lty=4)，長破線(longdash,lty5)などで描画されます．
番号の他に，点線であれば「lty="dotted"」のように指定することもできます．
線の太さは「lwd」で指定します．デフォルトの値は「lwd=1」です．
下図は，「type="o"」を指定して太線「lwd=3」で描画しました．
曲線の色は「col」で指定します．具体的な色としては，
黒(1, black)，赤(2, red)、緑(3,green)，青(4, blue)，
水色(5, cyan)，紫(6, purple)，黄(7, yellow)，灰(8,gray)が利用できます．
「col="red"」のように色の名前を指定することもできます(下図)．
ただし，実際にはもっと多数の色が登録済みです．
どのような色があるかは「colors()」を実行してみるとよいでしょう．
657色の名前が表示されます．全部の色の表示を望まないときは，
「head(colors())」か「tail(colors())」として最初と最後だけ表示させます．
下図は，出だしと最後の部分の左半分です．
なお，16進のrgbコードを指定して「col=rgb(赤,緑,青)」のような
指定をすることもできます．
点に配置する記号，線種，太さ，そして色の指定は「curve」でも利用できます．
範囲指定では「xlim=c(statr,end)」「ylim=c(start,end)」も利用でき，
軸の目盛りも細かく指定することができます
(参照)．
|
|
ヒストグラム
データの特徴をつかむには，
代表値を求めるだけではなくグラフ化するのが一番です．
「R」ではヒストグラムを簡単に描画することができ，
「hist(データ)」とするだけです．
このデータはベクトル型のデータである必要があります．
以下は，列データである data2＄norm を「norm」として
ベクトル型に置きかえて「hist(norm)」としたときの図です．
デフォルトでは縦軸は度数(frequency)で表示され，
「freq=TRUE」が省略されています．
個々の階級は \(\small (a,b]\) という右半閉区間であり，
「right=TRUE」が省略されています．\(\small [a,b)\) で考えたいときは
「right=FALSE」とします．
階級の数 (k) はスタージェスの公式 \(\small k\approx 1+\log_2(データ総数)\)
により決められ，それにより階級数や階級幅が設定されます．
この階級を区切りは「breaks」により設定します．
デフォルトでは「breaks="Sturges"」が指定されています．
この分割を決めるアルゴリズムとして，
他に「Scott」「FD」「Freedman-Diaconis」というものもあるようです．
階級数を「breaks=20」のように直接することもできます．
分割点を自分で決めるには，
その分割点を定めるベクトルを指定します．
「norm」はN(50,10)での乱数なので，
たとえば区間 [10,90] を5刻みの階級とするには，
区分点のベクトルとして「seq(10,90,5)」を指定します．
下図は，他に「freq=FALSE」(単に「freq=F」でもよい)として縦軸を確率密度とし，
「density=10」により斜線を引いています．
斜線の密度は，この数の増減で調整します．
色の指定は「col」により行います(図は略)．
階級の区分点は「breaks」に与えるベクトルで指定できます．
その間隔は必ずしも均等である必要はありません．
たとえば，[10,40] と [60,90] は5刻みにして
[50,60] は2刻みにするには，区分点を定めるベクトルを次のように定めます．
下図のベクトル「a」「b」「c」は，それぞれ端点の値が重複するので，
和集合「union」を取っています．
「hist」で「freq=FALSE」の指定をしていないのに，
なぜか縦軸は確率密度で表示されました．
階級幅が均一でないときは，縦軸は自動的に確率密度になるようです．
これに確率密度関数を重ね合わせるには「curve」を利用します．
正規分布N(50,10)の確率密度関数は「dnorm(x,50,10)」です．
横軸の範囲を合わせて，重ね描きするために「add=T」を指定します．
|
|
|
以上は，正規分布をもとにしましたが，
次にベキ分布(パレート分布)のデータを用いてヒストグラムを描画してみましょう．
ベキ分布の確率密度関数は \(\small f(x)=ab^{a}/x^{a+1}~(x\geq b)\)
により定義されます．「data2＄pare」は，\(\small a=1, b=1\) の場合の
ベキ分布にしたがう乱数1000個からなる列ベクトルです．
この場合のベキ分布は平均も分散も存在しないので，
乱数を発生させると，時折，とんでもなく大きな値が発生します．
下図では，最大値を調べて，その箇所までを1刻みにしました．
大部分は20程度までの値ですが，
大きくはずれた値が何個か発生していることが分かります．
なお，真偽値の「TRUE」「FALSE」は，単に T, F だけでもかまいません．
20程度までの部分を詳しく見るため，範囲を限定します．
下図では [1,15] の範囲を表示しました．
下図は，このグラフに確率密度関数である \(\small 1/x^2\) を
「add=TRUE」により重ね描きしたものです．
軸のラベルが重ならないように，密度関数でのラベルの出力を抑制しました．
なお，このベキ分布のように，裾が広い分布のヒストグラムを描画するには
「hist」はちょっと使いにくいところがあります．
最大値がべらぼうに大きな値となるので，標準の仕様では最初に示したような
グラフになります(参照1)．
細かく見るには最大値を調べて分割数を指定することになりますが，
全体像はやはり分かりやすい図ではありません(参照2)．
そのようなヒストグラムを扱うときは「truehist」を利用します．
その関数は，パッケージ「MASS」(Modern Applied Statistics with S)
に含まれているので，まずそれを読み込ませます．
「help(truehist)」とすると書式が表示されます．
「hist」と違い，最大値を指定することなく刻み幅(h)を指定するだけです．
また，デフォルトで縦軸は確率密度で表示されます．
「breaks」も使用できますが，「h」で直接刻み幅を指定する方が簡単です．
度数を表示させたいときは「prob=F」とします．
「hist」と同様に，「xlim」「ylim」「col」「xlab」「ylab」も使用できます．
デフォルトでは「col="cyan"」として表示されます．
下図は，0.5刻みにして表示範囲を [0, 15] とした場合です．
★なお，2つのヒストグラムを重ね合わせるには
「乱数」を，
折れ線で描画するには「ヒストグラム(2)」を参照されたい．
|
|
|
「R」の応用操作
|
ここでは，確率分布を扱うために必要な「R」の応用的操作について取りまとめました．
全般的な説明は他の解説サイトを参照のこと．
|
いろいろな確率分布
「R」は統計処理ソフトなので，いろいろな確率分布について，
その確率密度関数(pdf)，累積分布関数(cdf)，分位関数(quantile)，
そして乱数生成(random)に関する関数が用意されています．
確率密度関数を \(\small f(x)\) とするとき
累積分布関数は \(\small F(x)={\rm P}(X\leq x)=\int_{-\infty}^{x}f(x)\,dx\)
です．定め方から \(\small 0\leq F(x)\leq 1\) です．
また，分位関数は累積分布関数の逆関数であり，
\(\small F(x)=q\) とおくと \(\small x=F^{-1}(q)\) が分位関数であり，
\(\small {\rm P}(X\leq x)=q\) が成り立ちます．
「R」で利用できる主な確率分布と「R」での名称は次の通りです．
他にもありますが，個人的に使いそうにないものは省略しました．
| 確率分布 | 名称 |
| 二項分布 | binom |
| ポアソン分布 | pois |
| 一様分布 | unif |
| 正規分布 | norm |
| コーシー分布 | cauchy |
| カイ二乗分布 | chisq |
| t分布 | t |
| F分布 | f |
| 指数分布 | exp |
| 幾何分布 | geom |
| 超幾何分布 | hyper |
| ガンマ分布 | gamma |
具体的な確率分布に対する確率密度関数等の名称は，
正規分布N(m,s)の場合は次のように表されます．別な確率分布の場合は，
「norm」の箇所を上記の名称で置きかえます．
| 関数名 | Rでの名称 |
| 確率密度関数(pdf) | dnorm(x,m,s) |
| 累積分布関数(cdf) | pnorm(x,m,s) |
| 分位関数(quantile) | qnorm(q,m,s) |
| 乱数生成(random) | rnorm(n,m,s) |
個別にパッケージを読み込ませると，さらに多様な確率分布を
扱うことができます．どのような確率分布を利用できるかは
下記を参照してください。
以下では標準正規分布N(0,1)の場合のグラフを確認します．
その場合の確率密度関数は「dnorm(x)」だけでよく、
「plot」によるグラフは下図の通りです．
plot(dnorm, \(\small -\)4,4)
plot(pnorm,\(\small -\)4,4)
下図は，\(\small {\rm P}(X\leq 1.96)=0.9750021\) であること，
「qnorm」を利用すると逆の関係が得られることが示されています．
デフォルトでは「long.tail=TRUE」が暗黙に指定されています．
「long.tail=FALSE」を指定すると，
\(\small {\rm P}(X\geq x)=q\) の値が返されます．
平均 \(\small m\)，標準偏差 \(\small s\) の正規分布の確率密度関数は，
「dnorm(x,m,s)」と引数を増やすことになりますが，
その関数は \(\small x\) を含んでいるので
「plot」で描画することはできません．その場合は「curve」を利用します．
curve(dnorm(x,50,10),0,100)
個々の確率分布での引数の与え方は下記にまとめられています．
|
|
乱数の利用
正規分布の場合は，「rnorm(n,m,s)」により，
平均\(\small m\)，標準偏差\(\small s\) の正規分布にしたがう乱数が
\(\small n\) 個生成されます．
下図は，
N(50,10)にしたがう乱数100個からなるベクトルを \(\small x\) として，
生成順に線で繋いでグラフ化したものです．
乱数を利用すると，一つの確率分布にしたがう確率変数の
平均がどのような分布をするかをシミュレートすることができます．
試みに，
標準正規分布N(0,1)にしたがう乱数(rnorm)5個の平均を1000個作成して，
そのヒストグラムを描いてみましょう．
まず，そのようなベクトルを発生させる関数を「nor_mean」として定義します．
上記では，「numeric(N)」により0だけからなる長さがNのベクトルを生成して
「nm」とします．そして，「rnorm(n)」により
N(0,1)にしたがう n個の乱数を発生させて，その平均を nm[i] としています．
「nor_mean」を利用して，
「m1」をN(0,1)にしたがう乱数1個だけからなる要素数1000のベクトルとし，
「m5」をN(0,1)にしたがう乱数5個の平均からなる要素数1000のベクトルとします．
それぞれのヒストグラムを重ね描きすると，下図のようになります．
「main」により，グラフ画面上部のタイトルを指定することができます.
ここでの色指定はカラーコードを用いました．
色自体は"#ff9900"で指定されます．
その後の"80"は，色を重ねたときの透明度を指定しています．
念のため、確率密度曲線を重ねてみます．
周知のように，標準正規分布 N(0,1) にしたがう
乱数 \(\small n\) 個の平均は，正規分布 N(0,\(\small 1/\sqrt{n}\))
にしたがいます．
同様にして，N(0,1)にしたがう乱数 n個の平方和を考えます．
周知のように，その分布は自由度 n のカイ二乗分布にしたがいます．
下図では，平方和からなるベクトルを生成する関数を定義して，
1個と10個の和からなるベクトルを生成しています．
次に，それらのヒストグラムを描画して，確率密度関数(dchisq)と
重ね合わせます．
|
|
ヒストグラム(2)
3つ以上のヒストグラムを同時に表示するには，
棒グラフよりも折れ線グラフの方が適しています．
折れ線グラフで描画するには「lines」を利用します．書式は「lines(x, y)」です．
「x」「y」は横軸と縦軸の対応するベクトルか列を指定します．
上記は放物線 \(\small y=x^2\) を，区間 \(\small [-4,4]\) を0.1刻みにして
描画しようとしたものです．
その結果は，事前にどのようなグラフを表示させたかに依存します．
「plot」や「curve」は，グラフ画面の外枠を設定してグラフを描画しますが，
「lines」にグラフ画面の外枠を表示する機能はありません．
すでに設定されているグラフ画面に書き入れるだけなので，
それ以前にグラフを表示していない状態で「lines」を実行すると
「エラー」になります．
グラフを表示していても，範囲が合わないときはグラフは表示されません．
そこで，まず，グラフ画面の外枠だけを表示させて，その上で「lines」を
実行します．
下図における「type="n"」の指定は，外枠だけ設定して何も描画しないように
するための設定です．この指定をしないと，
「plot」の最初で指定した点 (0,0) の箇所に白丸が配置されます．
「plot」「curve」は高水準作図関数と呼ばれるのに対して，
「lines」は低水準作図関数とよばれています．
この「lines」を利用してヒストグラムを折れ線で描画するには，
ヒストグラムの階級値などのデータが必要です．
そのデータは「hist」を実行すると得られます．
たとえば，正規分布N(50,10)にしたがう乱数を1000個発生させて「n1」とする.
最小値と最大値を調べることにより，[20,75] の区間を5刻みに分割して
ヒストグラムを描画して「hn1」に代入すると，
ヒストグラムの情報が「hn1」に保存されます．
グラフも表示されますが，ここでは必要としないので無視します．
「names(hn1)」として1行目の列名を表示させると，
「hn1」は6列からなるデータであることが分かり、
分割点(breaks)，度数(counts)，確率密度(density)，階級値(mids)
などが保存されていることが分かります．
元ファイルのファイル名(xname)も保持されており，
行数の異なるデータになっています．
出だしの6行分を表示させようとして「head(hn1)」とすると，
各列の値が列ごとに全部表示されるので，
分割点を多く取っているときは注意する必要があります．
このデータを利用してヒストグラムを折れ線グラフを描画するには，
「x」には階級値(mids)を，「y」には度数(counts)か確率密度(density)を指定します．
下図では，最初に表示された度数によるヒストグラムに合わせて
度数で表示させました．
なお，下図の「＄」は全角ですが，実際には半角で入力します．
lines(hn1＄mids, hn1＄counts, col="red")
以上のことを利用して，
下図ではN(50,10), N(40,7), N(60,14)にしたがう乱数1000のヒストグラムを，
それぞれの刻み幅を 5, 2, 10 として作成しました．
「n1」と「hn1」はすでに作成済みです．
下図は，それらのヒストグラムのデータを用いて，
3つのヒストグラムをれ線で描画したものです．
階級値の箇所に点を配置したいときは「type="o"」を追記します．
|
|
中心極限定理 [script]
多くの場合，
分析しようとする母集団の統計量がどのような分布をしているのかは不明です．
それを調べようとして，標本抽出して標本平均や標本分散を計算し，
その値をもとに母集団の平均 \(\small \mu\) や分散 \(\small \sigma^2\)
を推測しようとします．
それらが存在する場合，中心極限定理は，標本の数を大きくしていけば，
\[\small \lim_{n\to\infty}
{\rm P} \left(\frac{\bar{X_n}-\mu}{\sigma/\sqrt{n}}\leq x\right)
=\frac1{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{x^2}{2}}\,dx\]
が成り立つこと，つまり，標本平均を標準化した値は標準正規分布 N(0,1)
に収束することを述べています．それは，標本数を大きくしていくと，
標本平均 \(\small \bar{X}\) の分布が
正規分布 N\(\small (\mu, \sigma^2/n)\) に近づいていくことを示しています．
以下では，このことをカイ自乗分布で確かめてみましょう．
■カイ二乗分布
自由度 \(\small n\) のカイ二乗分布は，標準正規分布 N(0,1) にしたがう
互いに独立な \(\small n\) 個の確率変数 \(\small X_1, X_2, \ldots,X_n\)
に対して，\(\small X=X_1^2+X_2^2+\cdots+X_n^2\) の分布であり，
平均と分散はそれぞれ \(\small E(X)=n, V(X)=2n\) です．
「R」でのカイ二乗分布の名称は「chisq」であるので，
確率密度関数は「dchisq(x,自由度)」，乱数生成は「rchisq(個数, 自由度)」です．
下図は，自由度が 1(黒), 2(赤), 5(緑), 10(青)の 場合の
カイ二乗分布の確率密度関数のグラフです．
正規分布とは似ても似つかない分布からの抽出とするため，
自由度 (degree of freedom) が df=1 の場合を考えます．最初に，
各要素が \(\small n\) 個の乱数平均からなるベクトルを生成する関数を定義して，
標本数が 2, 5, 10, 30 個の平均からなる要素数 1000個のベクトルを生成します．
列ベクトルにまとめて「summary」を実行して最小値や最大値を確認します．
「summary」は5数要約のことで，
最小値・第1四分位数・中央値・第3四分位数・最大値の5数
に加えて平均が表示されます．
次に，刻み幅を 0.2 としたヒストグラム情報を取得します．
それをもとに，標本数が \(\small n=2, 5, 10\) の場合のヒストグラムを
折れ線で描画すると次のようになります．
平均をとる標本数を増やしていくと左右対称の山形になっていくのが分かります．
次に，生成したベクトル x2, x5, x10 を標準化したベクトルに対して，
同じようなヒストグラムを描画します．
母集団は自由度が df=1 のカイニ乗分布であるので，
E(X)=n=1, V(X)=2n=2 であり，
標本数 \(\small n\) 個の場合の標本平均は \(\small 2/n\) です．
標準化したベクトルのヒストグラム情報を取得します．
その値をもとに，5, 10, 30個の場合のグラフを折れ線で描画します．
標本数が増えるにしたがい，左右対称の形に近づいていくのが分かります．
標本数30個の場合は，平均が 1，分散は \(\small \frac2{30}\) の正規分布で
近似できます．下図では，このことを標準化しないで確認します．
★平均や分散が存在するような確率分布では，他の分布でも同様にして，
標本平均を標準化すると標準正規分布 N(0,1) に近づいていくことを
確認することができます．
★平均や分散が存在しないような確率分布もあります．
そのような場合の中心極限定理に相当するものとして，
一般化中心極限定理があり，正規分布に代わり安定分布が重要な位置を占めています．
詳しくは「こちら」を
参照してください．
|
|
