R言語で二変量正規分布を扱う

2010年08月17日


はじめに

R本体だけでは二変量正規分布(二次元正規分布)をはじめとした多変量正規分布(多次元正規分布)を扱うことはできません。このため、自分で定義するか、MASSあるいはmvtnormというライブラリを読み込む必要があります。単に、多次元の正規乱数を作成するだけなら、MASSパッケージを使うのが楽です。確率密度や累積密度を求めたかったら、mvtnormパッケージを使うのがよいでしょう。

以下では、二変量正規分布を中心に扱いますが、三次元以上も同様に処理することができます。

正規分布とは?

二変量正規分布を扱う前に、もっとも基本的な正規分布、すなわち一次元の正規分布をRでどう扱うかを確認しましょう。

一変量正規分布(一次元正規分布)は、μ,σ2の2つのパラメータを持ちます。μは平均、σ2は分散に相当します。正規分布は、N(μ, σ2)と表します。例えば、平均が4で、分散が9の正規分布は、N(4, 9) と表記されます。正規分布の確率密度函数は以下の通りです。

正規分布の確率密度函数

平均が4で、分散が9の正規分布、すなわち、N(4, 9) の確率密度函数に0.5を代入した場合の値は、以下のようにして求めることができます。

> dnorm(0.5, 4, 3)
[1] 0.0673329

Rでは、正規分布のパラメータとして分散の代わりに標準偏差(すなわち分散の平方根)を用います。このため、dnorm(0.5, 4, 9)とならずに、dnorm(0.5, 4, 3)となります。

同様にして、pnorm()を使うことで累積密度を求めることができます。

> pnorm(0.5, 4, 3)
[1] 0.1216725

qnorm() は、pnorm() の逆で、累積密度から確率点を求めます。

rnorm()という函数を使うと、正規分布に従う乱数を生成することができます。平均が4で、分散が9の正規分布、すなわち、N(4, 9) に従う乱数を6個生成するためには、

rnorm(6, 4, 3)

とコマンドを入れます。ここでも、コマンド入力に当たっては、分散の9でなく標準偏差の3を入力します。

なお、次の方法で標準正規分布 N(0,1)、すなわち平均が0で分散が1の正規分布の確率密度函数のグラフを描くことができます。

# 確率密度函数本体を描く
x <- seq(-3,3,0.01)
plot(x, dnorm(x, mean=0, sd=1), type="n", xlab="", ylab="")
curve(dnorm(x, mean=0, sd=1), type="l", add=T)
# 確率密度函数で、
# x= -3, -2.5, -2, ...., 3 となるところから
# x軸に向かって垂線を引く。
a <- seq(-3, 3, 0.5)
b <- dnorm(a, mean=0, sd=1)
c <- numeric(13)
segments(a, b, a, c, lty=2)
# x軸を描く
abline(h=0)
正規分布のグラフ

多変量正規分布とは?

多変量正規分布についておおざっぱな説明をしたいと思います。先ほど紹介した一変量正規分布では、1系列の確率変数しか出てきませんでした。例えば、先ほど、

rnorm(6, 4, 3)

とコマンドを入れたとき、1×6のデータしか出てきませんでした。要するに、1個正規乱数を出すのを6回やったわけです。

しかし、1回やって正規乱数が2個出てくる場合も想定できます。この場合、6回やるとm2×6=12個のデータが出てきます。つまり、1回ごとに2系列の確率変数が出てくるわけです。

多変量になると、出てくる系列の数が増えます。先ほどの、2×6=12個のデータが出てくる例で言うと、2系列の確率変数が出てきます。これを二変量正規分布と言います。三変量正規分布ならば、3系列の確率変数が出てくるわけです。

多変量正規分布のパラメータとしては、平均のベクトルμと、分散−共分散行列Σを考える必要があります。分散−共分散行列のΣは、対称な行列で、対角成分には分散、非対角成分には共分散が入ります。なお、以下に紹介する方法では、共分散を入力するかわりに相関係数を入力する場合があります。ただ、共分散は、2つの確率変数のそれぞれの分散と、この2つの確率変数の間の相関係数から求めることができますので、これは問題ありません。

自分で定義する

R-tipsでは、第60節で、r2norm() という二変量正規乱数を生成する関数の定義が以下のように記されています。

r2norm <- function(n, mu, sigma, rho) {
  tmp <- rnorm(n)
  x <- mu+sigma*tmp
  y <- rho*x + sqrt(1-rho^2)*rnorm(n)
  return(data.frame(x=x,y=y))
}

この函数は、X, Yがともに、N(μ, σ2)に従い、XとYの相関係数がρであるときの正規乱数としてX, Yの組み合わせを作ります。例えば、X, Yがともに、N(4, 9)に従い、XとYの相関係数が 0.4 のとき、これに従う正規乱数を6個生成するには、以下のようなコマンドを入力します。

r2norm (6, 4, 3, 0.4)

この定義のままでは、XとYのそれぞれの平均・分散が同じでなくてはなりません。つまり、X が N(-3, 12) に、Y が N(8, 21) に従うというような場合はこの定義では扱えません。

MASSパッケージを使う

MASSパッケージには、mvrnorm()という函数が用意されており、簡単に多変量正規乱数を生成することができます。MASSパッケージは、通常Rをインストールするときに一緒に導入されます。このため、特別にMASSパッケージをダウンロードしてきたりする必要はありません。ただし、パッケージですので、

library(MASS)

のようにして読み込む必要があります。

6個の二変量正規乱数を生成する場合を考えましょう。なお、パラメータは以下の通りです。

ある二変量正規分布のパラメータ

この場合、以下のようにコマンドを入力します。

mu    <- c(-1, 7)
Sigma <- matrix(c(4, 3, 3, 25), 2, 2)
mvrnorm(6, mu, Sigma)

上のコマンドで、 mu が平均のベクトルμに、Sigma が分散−共分散行列Σ に相当します。

mvtnorm パッケージを使う

mvtnorm というパッケージを使うことで、多変量正規分布を扱うことができます。このパッケージを使うには、CRAN等を通じて、ダウンロードしてインストールする必要があります。準備できたら、

library(mvtnorm)

として、mvtnorm パッケージを読み込みます。

Rにもとから用意されている一次元の正規分布用の函数に対応する形で、mvtnorm パッケージに多変量正規分布用の函数が用意されています。例えば、一次元の正規乱数を生成する rnorm() に対し、mvtnorm パッケージでは、多次元の正規乱数を生成する rmvnorm() という函数が用意されています。その他の函数の対応関係をまとめると以下の通りになります。

  • dnorm() ― dmvnorm():確率密度を求める
  • pnorm() ― pmvnorm():累積密度を求める
  • qnorm() ― qmvnorm():累積密度から確率点を求める
  • rnorm() ― rmvnorm():乱数を生成する

先ほど、MASSパッケージの mvrnorm() という函数を使って、正規乱数を生成しましたが、同じ事をmvtnorm パッケージの rmvnorm() 函数を使ってやってみましょう。コマンドは以下の通りです。

mu    <- c(-1, 7)
Sigma <- matrix(c(4, 3, 3, 25), 2, 2)
rmvnorm(6, mu, Sigma)

mvtnorm パッケージの用法は三中信宏氏の書いた「二変量正規分布の2D/3Dグラフィクス」というページに詳しくまとめられています。なお、mvtnormパッケージでは、多変量正規分布だけでなく、多変量t分布を扱うこともできます。



Posted by green_idea at 21:54│Comments(0)TrackBack(0)統計・計量・数学 | プログラミング |

この記事へのトラックバックURL

http://trackback.blogsys.jp/livedoor/green_idea/1477686

コメントする

名前
URL
 
  絵文字