「R」に関する備忘録

in [数ナビの部屋]

「R」に関する諸操作を備忘録的に取りまとめました.

統計処理ソフト「R」を活用して統計解析を究めよう!

  • モバイル利用(Android)でのメニュー選択は、 SiteMapを利用するか、 または「長押し」から「新しいタブを開く」を選択してください。
■ 「R」に関する備忘録
[御案内] フリーの本格的な統計処理ソフトに「R」があります. ここでは,「R」のパッケージを取り扱うに至るまでの 「R」の利用法について備忘録的に取りまとめました. 「R Studio」は使用していません.


「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」だけでもよい.

ベクトルを結合したり,共通部分をとったりすることもできる.
  1. 「append(a,c)」により,「a」の後に「c」が追加される.
  2. 「append(a,c,after=2)」により,「a」の2番目の要素の後に「c」が追加される.
  3. 「union(b,c)」により,「b」と「c」の和集合が生成される.
  4. 「intersect(b,c)」により,「b」と「c」の共通部分が生成される.
  5. 「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) に近づいていくことを 確認することができる.
★平均や分散が存在しないような確率分布もある. そのような場合の中心極限定理に相当するものとして, 一般化中心極限定理があり,正規分布に代わり安定分布が重要な位置を占める. 詳しくは「こちら」を参照されたい.
▲戻る(トップメニュー)

リンク集
ここでは,「R」の使い方を解説しているサイトを取りまとめた.
  • 統計言語Rの使い方
    「R」の基本的な使い方についてスライドで解説されている。
  • 統計ソフト「R」の使い方
    名古屋市立大学大学院医学研究科の統計学の講義・セミナーで使用している Web教科書で,インストールの仕方や基本的な使い方について解説されている。
  • プログラミング言語R入門
    中京大学での全15回にわたる講義概要で, 「R」の基礎部分から推定・検定や因子分析にいたるまでが解説されている.
  • R言語入門
    「全人類が分かる統計学」の中で、「R言語」についても解説されている.
  • 超初心者向けRガイド
    「Rで重回帰分析をやってみる」という講習会向け資料の概要版で, 具体的な使い方について参考になる.
  • Rの初歩
    三重大学の奥村晴彦先生により、「R」の使い方に関して解説されている.
  • Rでテクニカルデータプレゼンテーション
    Rの基本操作や、いろいろな統計分析についてまとめられている.
  • 統計分析フリーソフト「R」
    「R」のインストールから、検定、回帰分析、判別分析などまで解説されている.
  • プログラマー向けのR言語入門
    「R」のいろいろなコマンドの使い方について、簡潔にまとめられている.
  • Rのページ
    Rについて、計算処理、グラフ、データ構造、プログラミング、オペレーションリサーチなどについて説明されている.
  • R-Tips
    「R」についてある程度知っているとき、 コマンドの使い方などに関する逆引きをしたいときに便利なサイトである. PDFにまとめたファイル(200頁)も登録されている. ( 参照)。
  • Rを利用した統計解析およびデータの視覚化
    Rを用いた統計解析について、基礎編、発展編、 そしてグラフィックスに分けて解説されている.
  • Rのグラフィックス
    Rのグラフィックスについて詳しく解説されています.
  • R、R言語、R環境
    同志社大学文化情報学部のデータサイエンス研究室のサイトのページで, Rの初歩から推定・検定、さらには分散分析・回帰分析・多変量解析等にいたるまで、 解説されている.
  • Rによる統計処理
    「Rによる統計解析(オーム社)」のサポートページで, 書籍を購入していることが前提と思われるが, 「R」の概要を把握した上で見ると 膨大な統計プログラム例が登録されており具体的な使い方について参考になる.
  • Rでプログラミング
    統計解析の知識を前提とした上で、 「R」の初心者でもデータの一括処理やグラフ化のプログラムを 書けるためのチュートリアルである.
  • R基本統計関数マニュアル
    Rの統計関連の種々の関数について、 カテゴリー別に整理した上でヘルプドキュメントが和訳されている.
  • R言語入門
    「R」の基本的な使い方について、13回に分けて動画で解説されている.
  • Rによる統計解析の基礎
    「統計解析」の教科書でいろいろな例に「R」が利用されている. 絶版となったためPDFファイルとして公開されているもので, 通常の推定・検定を含み、分散分析、一般化線形モデル、主成分分析、因子分析、そしてクラスター分析まで扱われている(194頁)。
  • 社会統計演習
    青山学院大学の先生による「社会統計演習」の情報をまとめたページで, 「R」の使い方について多数の動画が登録されている.
▲戻る(トップメニュー)


copyright