準備

今回の実習に必要なパッケージを(必要ならまずインストールし、)呼び出す。

if (!require(rgl)) {
    install.packages("rgl", dependencies = TRUE)
    library("rgl")
}
## Loading required package: rgl
if (!require(mvtnorm)) {
    install.packages("mvtnorm", dependencies = TRUE)
    library("mvtnorm")
}
## Loading required package: mvtnorm
if (!require(MCMCpack)) {
    install.packages("MCMCpack", dependencies = TRUE)
    library("MCMCpack")
}
## Loading required package: MCMCpack
## Loading required package: coda
## Loading required package: MASS
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2015 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
if (!require(ggplot2)) {
    install.packages("ggplot2", dependencies = TRUE)
    library("ggplot2")
}
## Loading required package: ggplot2
if (!require(gridExtra)) {
    install.packages("gridExtra", dependencies = TRUE)
    library(gridExtra)
}
## Loading required package: gridExtra
## Loading required package: grid

確率分布

同時分布 (joint distributions)

多変量正規分布 (multivariate normal distribution) を例に、同時分布のイメージを掴もう。

確率変数ベクトル \(y=(y_1, \dots, y_p)^T\) があり、\(y\)が多変量正規分布に従うとき、 \[y \sim \mathrm{N}(\mu, \Sigma)\] と書く。ただし、\(\mu\) は長さ\(p\)の位置母数ベクトル、\(\Sigma\)\(p\times p\) の分散共分散行列である。 確率密度関数は、 \[p(y | \mu, \Sigma) = (2\pi)^{-\frac{p}{2}} |\Sigma |^{-\frac{1}{2}}\exp \left[ -\frac{1}{2}(y - \mu)^T \Sigma^{-1}(y - \mu)\right]\] である。

\(y\) の期待値と分散(共分散行列)は、それぞれ、 \[\mathrm{E}(y) = \mu\], \[\mathrm{var}(y) = \Sigma\] と、なる。

変量間に相関がない場合の二変量正規分布

\(p=2\)、すなわち \(y=(y_1, y_2)^T\) で、変量間に相関がない場合を考える。 変量間に相関がないということは、\(\mathrm{cov}(y_1, y_2) = 0\) ということだから、 \[\Sigma = \left( \begin{array}{cc} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \end{array} \right)\] である。

例として、 \(\mu = (\mu_1, \mu_2)^T = (5, 0)^T\)\(\sigma_1 = 3\)\(\sigma_2 = 1\) の場合を考えよう(その他の場合もいくつか自分で確かめること)。

まず、\(\mu\)\(\Sigma\) を設定しよう。

Mu <- c(5, 0)
Sigma <- matrix(c(3, 0, 0, 1)^2, nrow = 2)

この分布から、\(n\) 個の値を無作為抽出し、散布図にしてみよう。

n <- 1000
set.seed(2015-05-03)
y <- rmvnorm(n, mean = Mu, sigma = Sigma)
plot(y, pch = 20, xlab = expression(y[1]), ylab = expression(y[2]))

この図の中心部に点が集まっているので、中心の密度が相対的に高いことが推測できるが、もう少しわかりやすくしてみよう。

df <- data.frame(y1 = y[,1], y2 = y[, 2])
mvn1 <- ggplot(df, aes(x = y1, y = y2)) + geom_point() +
    xlim(min(df$y1), max(df$y1)) + 
    scale_y_continuous(limits = c(min(df$y2), max(df$y2)),
                       breaks = c(-2, 0, 2), labels = c("-2.0", "0.00", "2.00")) +
    labs(x = expression(y[1]), y = expression(y[2]))
mvn1 + geom_rug(color = "violet", alpha = .1)