政治学方法論 II：多変量のベイズ推定

準備

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

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)

同時分布の形状を詳しく調べるために、三次元の図を描いてみよう。 plot3d() で図を描くと、図を回転させて同時分布の形状を確認することができる。

df$density <- dmvnorm(y, mean = Mu, sigma = Sigma)
plot3d(df$y1, df$y2, df$density)

(3次元の図はブラウザ上には表示されないので、各自で確かめること。)

等高線図 (contour plot) を描いてみよう。

df2 <- expand.grid(df$y1, df$y2)
names(df2) <- c("y1", "y2")
df2$density <- dmvnorm(expand.grid(df$y1, df$y2), mean = Mu, sigma = Sigma)
cont <- ggplot(df2, aes(x = y1, y = y2, z = density)) +
  geom_contour()
cont + labs(x = expression(y[1]), y = expression(y[2]),
            title = "Contour plot of joint distribution")

次に、各変量の周辺分布を確認してみよう。

ggplot(df, aes(y1, ..density..)) + geom_histogram(binwidth = 1) +
    xlab(expression(y[1])) +
    ggtitle(expression(paste("Marginal distribution of ", y[1]))) 

ggplot(df, aes(y2, ..density..)) + geom_histogram(binwidth = .5) + 
    xlab(expression(y[2])) +
    ggtitle(expression(paste("Marginal distribution of ", y[2]))) 

どちらの変量についても、周辺分布が正規分布になっていることがわかる。

最後に、散布図と周辺分布を同時に見てみよう。

mgnl_y1 <- ggplot(df, aes(y1, ..density..)) + geom_histogram() + 
    xlim(min(df$y1), max(df$y1)) + xlab("")
mgnl_y2 <- ggplot(df, aes(y2, ..density..)) + geom_histogram() + 
    xlim(min(df$y2), max(df$y2)) + xlab("") + coord_flip()
empty <- ggplot() + geom_point(aes(0, 0), color = "white") +
    theme(plot.background = element_blank(), 
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          panel.border = element_blank(), 
          panel.background = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks = element_blank() )
grid.arrange(mgnl_y1, empty, mvn1, mgnl_y2,
             ncol = 2, nrow = 2, widths = c(4, 1.5), heights = c(1.5, 4))

変量間に相関がある場合の二変量正規分布

\(p=2\)、すなわち \(y=(y_1, y_2)^T\) で、変量間に相関がある場合を考える。 例として、\(\mathrm{cov}(y_1, y_2) = 0.8^2 = 0.64\) の場合を考えよう。 他の母数の値は、\(\mu = (\mu_1, \mu_2)^T = (0, 0)^T\)、\(\sigma_1 = \sigma_2= 1\) とする。

Rで、\(\mu\) と \(\Sigma\) を設定しよう。

Mu <- c(0, 0)
Sigma <- matrix(c(1, .8, .8, 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])
mvn2 <- 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 = seq(-4, 4, by = 2),
                       labels = c("-4.0", "-2.0", "0.00", "2.00", "4.00")) +
    labs(x = expression(y[1]), y = expression(y[2]))
mvn2 + geom_rug(color = "violet", alpha = .1)

同時分布の形状を詳しく調べるために、三次元の図を描いてみよう。

df$density <- dmvnorm(y, mean = Mu, sigma = Sigma)
plot3d(df$y1, df$y2, df$density)

(3次元の図はブラウザ上には表示されないので、各自で確かめること。)

等高線図 (contour plot) を描いてみよう。

df2 <- expand.grid(df$y1, df$y2)
names(df2) <- c("y1", "y2")
df2$density <- dmvnorm(expand.grid(df$y1, df$y2), mean = Mu, sigma = Sigma)
cont <- ggplot(df2, aes(x = y1, y = y2, z = density)) +
  geom_contour()
cont + labs(x = expression(y[1]), y = expression(y[2]),
            title = "Contour plot of joint distribution")

次に、各変量の周辺分布を確認してみよう。

ggplot(df, aes(y1, ..density..)) + geom_histogram(binwidth = .5) +
    xlab(expression(y[1])) +
    ggtitle(expression(paste("Marginal distribution of ", y[1]))) 

ggplot(df, aes(y2, ..density..)) + geom_histogram(binwidth = .5) +
    xlab(expression(y[2])) +
    ggtitle(expression(paste("Marginal distribution of ", y[2]))) 

どちらの変量についても、周辺分布が正規分布になっていることがわかる。

散布図と周辺分布を同時に見てみよう。

mgnl_y1_2 <- ggplot(df, aes(y1, ..density..)) + geom_histogram() + 
    xlim(min(df$y1), max(df$y1)) + xlab("")
mgnl_y2_2 <- ggplot(df, aes(y2, ..density..)) + geom_histogram() + 
    xlim(min(df$y2), max(df$y2)) + xlab("") + coord_flip()
grid.arrange(mgnl_y1_2, empty, mvn2, mgnl_y2_2,
             ncol = 2, nrow = 2, widths = c(4, 1.5), heights = c(1.5, 4))

逆ガンマ分布 (inverse-gamma distribution)

• 正規分布の尺度母数を推定するときの共役事前分布

• 母数

  • \(\alpha > 0\)：形状母数
  • \(\beta > 0\)：尺度母数

• \(x \sim\) Inv-gamma(\(\alpha, \beta\))

• 確率密度関数 \[p(x |\alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{-(\alpha + 1)} \exp\left(-\frac{\beta}{x} \right)\]

curve(dinvgamma(x, shape = 1, scale = 1), from = 0, to = 3,
      ylim = c(0, 3), lwd = 3, col = "tomato", ylab = "density",
      main = expression(paste("PDF of Inv-gamma(", alpha, ", ", beta, ")")))
curve(dinvgamma(x, shape = 2, scale = 1), from = 0, to = 3,
      ylim = c(0, 3), lwd = 3, col = "royalblue",
      main = "", ylab = "", add = TRUE)
curve(dinvgamma(x, shape = 3.5, scale = 1), from = 0, to = 3,
      ylim = c(0, 3), lwd = 3, col = "seagreen",
      main = "", ylab = "", add = TRUE)
curve(dinvgamma(x, shape = 2, scale = 2), from = 0, to = 3,
      ylim = c(0, 3), lwd = 3, lty = "dashed", col = "orchid4",
      main = "", ylab = "", add = TRUE)
curve(dinvgamma(x, shape = 1, scale = 0.5), from = 0, to = 3,
      ylim = c(0, 3), lwd = 3, lty = "dotted", col = "salmon",
      main = "", ylab = "", add = TRUE)
legend("topright", lwd = 2,
       title = expression(paste("(", alpha, ", ", beta, ")")),
       lty = c(rep("solid", 3), "dashed", "dotted"),
       col = c("tomato", "royalblue", "seagreen", "orchid4", "salmon"),
       legend = c("(1.0, 1.0)", "(2.0, 1.0)", "(3.5, 1.0)",
                  "(2.0, 2.0)", "(1.0, 0.5)"))

スケール（縮小または拡大）された逆カイ二乗分布 (scaled inverse-\(\chi^2\) distribution)

• 正規分布の尺度母数を推定するときの共役事前分布として使う

• 母数

  • \(\nu > 0\)：自由度
  • \(s >0\)：尺度母数

• \(x \sim\) Inv-\(\chi^2\)(\(\nu, s^2\))

• 確率密度関数

\[p(x | \nu, s^2) = \frac{(\nu / 2)^{\nu / 2}}{\Gamma(\nu / 2)} s^\nu x^{-(\nu /2 + 1)} \exp\left( -\frac{\nu s^2}{2x}\right)\]

• 逆ガンマ分布との関係：Inv-\(\chi^2(x | \nu, s^2) =\) Inv-gamma\((x | \alpha = \nu / 2, \beta = (\nu / 2)s^2)\)

nu.vec <- c(2, 4, 7, 4, 2)
s.vec <- c(1, 1, 1, .5, 2)
colors <- c("tomato", "royalblue", "seagreen", "orchid4", "salmon")
ltypes = c(rep("solid", 3), "dashed", "dotted")
plot("n", xlim = c(0, 3), ylim = c(0, 3), xlab = "x", ylab = "density",
     main = expression(paste("PDF of Inv-", chi^2, "(", nu, ", ", s^2, ")")))

## Warning in xy.coords(x, y, xlabel, ylabel, log): NAs introduced by
## coercion

for (i in seq_along(nu.vec)) {
    curve(dinvgamma(x, shape = nu.vec[i]/2, scale = (nu.vec[i]/2) * s.vec[i]^2),
          from = 0, to = 5, lwd = 2, col = colors[i], lty = ltypes[i],
          add = TRUE)
}

無情報事前分布を用いた正規母数の推定

シミュレーション用データセットの作成

シミュレーション練習用にデータセットを作る。 n, mu, sigma の値と乱数の種を変えて様々な状況をq試すこと。

n <- 50
mu <- 8
sigma <- 3
set.seed(2015-05-04)
y <- rnorm(n, mean = mu, sd = sigma)

事後分布の抽出

まず、 \(\sigma^2 | y \sim\) Inv-\(\chi^2\)(\(n-1, s^2\)) によって \(\sigma^2\) を抽出し、 \(\mu | \sigma^2, y \sim\) N(\(\bar{y}, \sigma^2 / n\)) で \(\mu\) を抽出する、

n.sim <- 10000
post.sigma2 <- rinvgamma(n.sim, shape = (n - 1) / 2,
                         scale = ((n - 1) / 2) * var(y))
post.mu <- rnorm(n.sim, mean = mean(y), sd = sqrt(post.sigma2 / n))

is_cpost <- function(x, p = .95) {
    ## given the vector,
    ## return a vector of boolean showin if a value is within
    ## central posterior interval
    cutoff <- quantile(x, prob = c((1 - p) / 2, 1 - ((1 - p) / 2)))
    return(x > cutoff[1] & x < cutoff[2])
}

事後分布のヒストグラムを描いてみよう。95%中心事後区間 (central posterior interval) を赤い線で示す。

df <- data.frame(post.mu = post.mu, post.sigma = sqrt(post.sigma2))
hist.post.mu <- ggplot(df, aes(x = post.mu)) + geom_histogram() +
    labs(x = expression(mu), title = expression(paste("Posterior of ", mu)))
hist.post.mu +
    geom_vline(xintercept = quantile(df$post.mu, prob = c(.025, .975)),
               col = "red") 

hist.post.sigma <- ggplot(df, aes(x = post.sigma)) +  geom_histogram() +
    xlab(expression(sigma)) +
    ggtitle(expression(paste("Posterior of ", sigma))) 
hist.post.sigma + 
    geom_vline(xintercept = quantile(df$post.sigma, prob = c(.025, .975)),
               col = "red") 

授業の内容に戻る