色んなものが正規分布っぽい分布をしますが、ＲのMASSライブラリに入っているレプトクラブサスカニの甲羅の長さ（どうやら正規分布っぽいらしい）のデータを使って正規分布の二つのパラメータ（平均と分散）をベイズ推定するのがこの記事です。最後にRのコードも載せてます。最後のコードは分散を所与としたときのコードも最初の方に入れてて、そこで事前分布の影響の大きさを変えて事後分布事前分布に引っ張られる具合の違いをみるってのもやってるので気になる人は見てみてください。適当にコードを実行してグラフを見るだけでもどうぞ。

先ずはデータを見る

このデータってコラムが８個あるんですけどCL（甲羅の長さ）だけに注目します。

まあ大体正規分布ってことにしても問題なさそうですね。甲羅の長さはカウントデータでもないですし。

平均と分散の同時推定

細かい数学は長すぎるので大事な事実だけを書きます。

ベイズ的にいきたいので事前分布とデータから事後分布を求めるわけですが、事前分布と事後分布が同じ分布族（例えば事前も事後も両方パラメータの値だけが違うベータ分布、みたいな）であるのが望ましいことになっているんで正規分布するデータを挟んだ時に平均と分散が同じ分布族になるような事前分布を探さねばなりません。ちなみにそんな事前分布を共役な事前分布とかいいます。

正規分布の場合、それが

• 分散：逆ガンマ分布（大層に聞こえるが、ただのガンマ分布の逆数。つまり、分散の逆数がガンマ分布するということ。）

• 平均：正規分布

です。つまり、分散の事前分布がが逆ガンマ分布で与えられればデータで情報を更新した後の事後分布も逆ガンマ分布となり、平均についてはどっちも正規分布となる、という寸法です。ただ、事前と事後でパラメータが違うというだけです。しかも、凄いのが事後パラメータは事前パラメータとデータから得られる統計量の加重平均になっていたりと無意味な変化ではないので数式を見る時のいらいら感は比較的少ない気がします。

分散$${\sigma^2}$$の事前分布の分散$${\sigma^2_0}$$が$${\nu_0}$$個の事前観測値から得られた分散、平均$${\theta}$$の事前分布の平均$${\mu_0}$$が$${\kappa_0}$$個の事前観測値から得られた平均だとします。すると、

$$
\frac{1}{\sigma^2} : gamma(\frac{\nu_0}{2}, \frac{\nu_0}{2}\sigma^2_0)  \\
\theta|\sigma^2 : normal(\mu_0, \frac{\sigma^2_0}{\kappa_0})
$$

という事前分布にしておくと、$${Normal(\theta, \sigma^2)}$$に独立に従うサンプルサイズ$${n}$$のデータ$${Y}$$を観測した後の事後分布は、$${\nu_n = \nu_0 + n}$$、$${\kappa_n = \kappa_0 + n}$$として

$$
\frac{1}{\sigma^2}|Y: gamma(\frac{\nu_n}{2}, \frac{\nu_n}{2}\sigma^2_n),   ( \sigma^2_n = \frac{1}{\nu_n}[\nu_0\sigma^2_0 + (n + 1)s^2 + \frac{\kappa_0n}{\kappa_n}(\bar{y}-\mu_n)^2]   ) \\
\theta|\sigma^2, Y : normal(\mu_n, \frac{\sigma^2}{\kappa_n}) ,   (
\mu_n = \frac{\kappa_0}{\kappa_n}\mu_0 + \frac{n}{\kappa_n}\bar{y}  )
$$

平均$${\theta}$$の分布は分散$${\sigma^2}$$の条件付き分布になっていのには注意してください。

事前分布の観察

ここでは、事前の情報として適当に甲羅の長さの平均は$${28}$$で分散は$${14^2}$$としました（甲羅の長さは正なので$${0}$$から平均の間の距離を2標準偏差としてみました）。分散はかなり大きめに評価してあります。

上の逆ガンマ分布する分散$${\sigma^2}$$と正規分布する平均$${\theta}$$の分布をモンテカルロ近似させたらこんな感じです。では、一番上のヒストグラムに表される$${n = 200}$$のデータを観察した後の分布はどうなってるのか見てみましょう。

事後分布の観察（事前分布との比較）

分散はかなり大きめに評価してあったので、データを観察したあとの分散の分布はかなり小さい値辺りに分布するようになっています。適当に平均$${28}$$としていたのも、青い事後分布をみると$${30}$$より少し大きいところに分布しているのが見て取れます。

これは$${\sigma^2}$$と$${\theta}$$の事後分布に$${95\%}$$信用区間を赤い線で入れたグラフです。赤い線で囲まれた部分に95%の確率で母集団の真のパラメータが入っているということです。ここで注意するべきは、ここでは信用区間を使っているのであって信頼区間ではないのでこの頻度論的信頼区間の解釈としては起こられてしまう「$${95%}$$の確率で・・・」の解釈で問題ありません。

点が密集しているところほど値が大きいというような見方をする可視化です。何となく点が最も濃い部分と上の事後分布の中心が一致しているのも分かりますよね。

最後に仮説を検証する。

最後に「カニの甲羅の長さの平均の真の値は$${32mm}$$よりも大きいのか」というクソどうでも良い仮説をモンテカルロシミュレーションを使って検証してみます。

>print((sum(mc_theta_post > 32) / length(mc_theta_post)*100))
[1] 1.4853

どうやらこのカニの甲羅の長さが$${32mm}$$以上である確率は$${1.4853\%}$$らしいのでまあ真の平均値は$${32mm}$$以上ではないということでいいでしょう。

今回もつまらない記事を最後まで読んでくれてありがとうございました。Rのコードを最後に置いておきます。

Rのコード

library(MASS); library(ggplot2); library(formattable)
data <- MASS::crabs
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#分散は既知の時の平均の推測 (k.0 = 10)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#Exploratory Analysis
ggplot(data, aes(x = CL)) +
  geom_histogram(color = 'white', fill = 'blue') +
  xlab('Carapace Length (mm)')+ ylab('Frequency')+
  theme_minimal()
#data
sigma2<- var(data$CL)
ybar <- mean(data$CL)
n <- nrow(data)
#Prior of theta
k.0_1 <- 10
mu.0_1 <- 28 ; tau2.0_1 <- sigma2/k.0_1
#Posterior
mu.n_1 <- (k.0_1*mu.0_1 + n*ybar)/(k.0_1 + n)
tau2.n_1 <- 1/((1/tau2.0_1)+(n/sigma2))
#Monte Carlo Sampling fronm prior(theta)
S <- 1000000
mc_theta_pre_1 <- rnorm(S, mu.0_1, sqrt(tau2.0_1))
ggplot(mapping = aes(x = mc_theta_pre_1))+
  geom_histogram(color ='white', fill = 'green')+
  xlab(expression(paste('1000000 MC Sampling from Prior Distribution of ', theta)))
#Prior Stats
mean_theta_pre_1 <- mean(mc_theta_pre_1)
var_theta_pre_1 <- var(mc_theta_pre_1)
std_theta_pre_1 <- sqrt(var(mc_theta_pre_1))
#Monte Carlo Sampling from Posterior
mc_theta_post_1 <- rnorm(S, mu.n_1, sqrt(tau2.n_1))
ggplot(mapping = aes(x = mc_theta_post_1))+
  geom_histogram(color ='white', fill = 'blue')+
  xlab(expression(paste('1000000 MC Sampling from Posterior Distributio of', theta)))
#Comparison between the prior and the posterior
#histograms
ggplot()+
  geom_histogram(mapping = aes(x = mc_theta_pre_1), 
                 color = 'white', fill = 'green', alpha = 0.2)+
  geom_histogram(mapping = aes(x = mc_theta_post_1),
                 color = 'white', fill = 'blue', alpha = 0.2)+
  xlab(expression(paste('Comparison Between the Prior Dist. and the Posterior Dist.')))+
  labs(fill = 'Distribution')+
  scale_fill_manual(values = c('green', 'blue'),
                    labels = c('Prior Dist.', 'Posterior Dist.'))+
  theme_minimal()
#density plots
ggplot()+
  geom_density(mapping = aes(x = mc_theta_pre_1), 
               color = 'transparent', fill = 'green', alpha = 0.2)+
  geom_density(mapping = aes(x = mc_theta_post_1),
               color = 'transparent', fill = 'blue', alpha = 0.2)+
  xlab(expression(paste('Comparison Between the Prior Dist. and the Posterior Dist.')))+
  labs(fill = 'Distribution')+
  scale_fill_manual(values = c('green', 'blue'),
                    labels = c('Prior Dist.', 'Posterior Dist.'))+
  theme_minimal()
#Posterior Stats
mean_theta_post_1 <- mean(mc_theta_post_1)
var_theta_post_1 <- var(mc_theta_post_1)
std_theta_post_1 <- sqrt(var(mc_theta_post_1))
#Table of Stats
stats_table_1 <- data.frame(
  Statistic_1 = c('Mean', 'Variance', 'Standard Deviation'),
  Prior = c(mean_theta_pre_1, var_theta_pre_1, std_theta_pre_1),
  Posterior = c(mean_theta_post_1, var_theta_post_1, std_theta_post_1)
)
formattable(stats_table_1, list(Prior = color_tile('white', 'lightblue'),
                                Posterior = color_text('black', 'red')))
#95% Credible Interval for theta
percentiles_theta_1 <- quantile(mc_theta_post_1, c(0.025, 0.975))
ys_1 <- seq(30, 34, length = 500)
post_dist_1 <- data.frame(ys_1,
                          post_dens_1 = dnorm(seq(30, 34, length = 500), mu.n_1, sqrt(tau2.n_1))
)
ggplot(mapping = aes(x = post_dist_1$ys_1))+
  geom_line(mapping = aes(y = post_dist_1$post_dens_1),
            color = 'blue', size = 1)+
  geom_vline(xintercept = percentiles_theta_1[1], color = 'red')+
  geom_vline(xintercept = percentiles_theta_1[2], color = 'red')+
  xlab(expression(paste('95% Credible Interval for', theta)))+
  ylab('Density')+
  theme_minimal()



#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#分散は既知の時の平均の推測 (k.0 = 100)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#data
sigma2 <- var(data$CL)
ybar <- mean(data$CL)
n <- nrow(data)
#Prior of theta
k.0_2 <- 100
mu.0_2 <- 28 ; tau2.0_2 <- sigma2/k.0_2
#Posterior
mu.n_2 <- (k.0_2*mu.0_2 + n*ybar)/(k.0_2 + n)
tau2.n_2 <- 1/((1/tau2.0_2)+(n/sigma2))
#Monte Carlo Sampling fronm prior(theta)
S <- 1000000
mc_theta_pre_2 <- rnorm(S, mu.0_2, sqrt(tau2.0_2))
ggplot(mapping = aes(x = mc_theta_pre_2))+
  geom_histogram(color ='white', fill = 'green')+
  xlab(expression(paste('1000000 MC Sampling from Prior Distribution of ', theta)))
#Prior Stats
mean_theta_pre_2 <- mean(mc_theta_pre_2)
var_theta_pre_2 <- var(mc_theta_pre_2)
std_theta_pre_2 <- sqrt(var(mc_theta_pre_2))
#Monte Carlo Sampling from Posterior
mc_theta_post_2 <- rnorm(S, mu.n_2, sqrt(tau2.n_2))
ggplot(mapping = aes(x = mc_theta_post_2))+
  geom_histogram(color ='white', fill = 'blue')+
  xlab(expression(paste('1000000 MC Sampling from Posterior Distributio of', theta)))
#Comparison between the prior and the posterior
#histograms
ggplot()+
  geom_histogram(mapping = aes(x = mc_theta_pre_2), 
                 color = 'white', fill = 'green', alpha = 0.2)+
  geom_histogram(mapping = aes(x = mc_theta_post_2),
                 color = 'white', fill = 'blue', alpha = 0.2)+
  xlab(expression(paste('Comparison Between the Prior Dist. and the Posterior Dist.')))+
  labs(fill = 'Distribution')+
  scale_fill_manual(values = c('green', 'blue'),
                    labels = c('Prior Dist.', 'Posterior Dist.'))+
  theme_minimal()
#density plots
ggplot()+
  geom_density(mapping = aes(x = mc_theta_pre_2), 
               color = 'transparent', fill = 'green', alpha = 0.2)+
  geom_density(mapping = aes(x = mc_theta_post_2),
               color = 'transparent', fill = 'blue', alpha = 0.2)+
  xlab(expression(paste('Comparison Between the Prior Dist. and the Posterior Dist.')))+
  labs(fill = 'Distribution')+
  scale_fill_manual(values = c('green', 'blue'),
                    labels = c('Prior Dist.', 'Posterior Dist.'))+
  theme_minimal()
#Posterior Stats
mean_theta_post_2 <- mean(mc_theta_post_2)
var_theta_post_2 <- var(mc_theta_post_2)
std_theta_post_2 <- sqrt(var(mc_theta_post_2))
#Table of Stats
stats_table_2 <- data.frame(
  Statistic = c('Mean', 'Variance', 'Standard Deviation'),
  Prior = c(mean_theta_pre_2, var_theta_pre_2, std_theta_pre_2),
  Posterior = c(mean_theta_post_2, var_theta_post_2, std_theta_post_2)
)
formattable(stats_table_2, list(Prior = color_tile('white', 'lightblue'),
                                Posterior = color_text('black', 'red')))
#95% Credible Interval for theta
percentiles_theta_2 <- quantile(mc_theta_post_2, c(0.025, 0.975))
ys_2 <- seq(29, 33, length = 500)
post_dist_2 <- data.frame(ys_2,
                          post_dens_2 = dnorm(seq(29, 33, length = 500), mu.n_2, sqrt(tau2.n_2))
)
ggplot(mapping = aes(x = post_dist_2$ys_2))+
  geom_line(mapping = aes(y = post_dist_2$post_dens_2),
            color = 'blue', size = 1)+
  geom_vline(xintercept = percentiles_theta_2[1], color = 'red')+
  geom_vline(xintercept = percentiles_theta_2[2], color = 'red')+
  xlab(expression(paste('95% Credible Interval for ', theta)))+
  ylab('Density')+
  theme_minimal()

#!!!!!!!!!!!!!!!!!!!!!!!!!!!
#k0の値による事後分布の比較
#!!!!!!!!!!!!!!!!!!!!!!!!!!!
ys_3 <- seq(29, 34,length = 600)
ggplot()+
  geom_density(mapping = aes(x = mc_theta_pre_1),
               color = 'black')+
  geom_density(mapping = aes(x = data$CL),
               color = 'gray')+
  geom_line(mapping = aes(x = ys_3,
                          y = dnorm(ys_3, mu.n_1, sqrt(tau2.n_1))),
            color = 'green', size = 1)+
  geom_line(mapping = aes(x = ys_3, 
                          y = dnorm(ys_3, mu.n_2, sqrt(tau2.n_2))),
            color = 'blue', size = 1)+
  xlab(expression(paste('Comparison between different ', kappa, '0')))+
  ylab('Posterior Dist.')+
  theme_classic()



#!!!!!!!!!!!!!!!!!!!!!!
#平均と分散の同時推定
#!!!!!!!!!!!!!!!!!!!!!!
#data
sigma2 <- var(data$CL)
ybar <- mean(data$CL)
n <- nrow(data)
#Prior
#Prior of theta~normal(mu.0, tau.0)
mu.0 <- 28; k.0 <- 50
#Prior of sigma2~ inverse-gamma(nu.0/2, (nu.0/2)*sigmma0)
sigma2.0 <- 14**2 ; nu.0 <- 50
 #MC Sampling from the Prior
S <- 1000000
#Sigma
mc_sigma_pre <- 1/rgamma(S, nu.0/2, nu.0*sigma2.0/2)
ggplot(mapping = aes(x = mc_sigma_pre))+
  geom_histogram(color = 'white', fill = 'green')+
  xlab(expression(paste('1000000 MC Sampling from Prior Dist. of ', sigma)))+
  ylab('Frequency')+
  theme_minimal()
#theta
mc_theta_pre <- rnorm(S, mu.0, mc_sigma_pre/k.0)
ggplot(mapping = aes(x = mc_theta_pre))+
  geom_histogram(color = 'white', fill = 'green')+
  xlab(expression(paste('1000000 MC Sampling from Prior Dist of ', theta)))+
  ylab('Frequency')+
  theme_minimal()
#Posterior Prameters
nu.n <- nu.0 + n ;  k.n <- k.0 + n
#sigma
sigma2.n <- (1/nu.n)*(nu.0*sigma2.0 + (n+1)*sigma2 + ((k.0*n)/k.n)*(ybar - mu.0))
#theta
mu.n <- ((k.0*mu.0)/k.n) + ((n*ybar)/k.n)
 #MC Sampling from the Posterior
#sigma
mc_sigma_post <- 1/rgamma(S, nu.n/2, nu.n*sigma2.n/2)
ggplot()+
  geom_density(mapping = aes(x = mc_sigma_pre), 
               color = 'white', fill = 'green', alpha = 0.3)+
  geom_density(mapping = aes(x = mc_sigma_post),
               color = 'white', fill = 'blue', alpha = 0.3)+
  xlab(expression(paste('Prior and Posterior of ', sigma**2)))+
  xlim(0, 400)+
  theme_minimal()
#theta
mc_theta_post <- rnorm(S, mu.n, mc_sigma_post/k.n)
ggplot()+
  geom_density(mapping = aes(x = mc_theta_pre), 
               color = 'white', fill = 'green', alpha = 0.3)+
  geom_density(mapping = aes(x = mc_theta_post),
               color = 'white', fill = 'blue', alpha = 0.3)+
  xlab(expression(paste('Prior and Posterior of ', theta)))+
  xlim(18, 40)+
  theme_minimal()
#95% Credible Interval
#sigma
percentiles_sigma <- quantile(mc_sigma_post ,c(0.025, 0.975))
ggplot()+
  geom_density(mapping = aes(x = mc_sigma_post),
               color = 'white', fill = 'blue')+
  geom_vline(xintercept = percentiles_sigma[1], color = 'red')+
  geom_vline(xintercept = percentiles_sigma[2], color = 'red')+
  xlab(expression(paste('95% Credible Interval for ', sigma**2)))+
  theme_minimal()
#theta
percentiles_theta <- quantile(mc_theta_post, c(0.025, 0.975))
ggplot()+
  geom_density(mapping = aes(x = mc_theta_post),
               color = 'white', fill = 'blue')+
  geom_vline(xintercept = percentiles_theta[1], color = 'red')+
  geom_vline(xintercept = percentiles_theta[2], color = 'red')+
  xlab(expression(paste('95% Credible Interval for ', theta)))+
  theme_minimal()
#Simultaneous Density
#by scatter plot
plot(NA, NA, xlim = c(29.5, 32.6), ylim = c(64, 110),
     xlab = expression(theta), ylab = expression(sigma^2))
points(mc_theta_post[1:10000], mc_sigma_post[1:10000], pch = '.')
title(expression(paste('Simultaneous Posterior Density of ',theta, 'and ', sigma**2)) )