関数解読normchi2post

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)


データ分布
y~N(mu,sig2)

事前分布
p(mu,sig2)proportional to (sig2)^(-1)

とする。このときの(mu,sig2)の事後分布（カーネル）の対数に関するある(mu_0,sig2_0)における値。


使用法

normchi2post(theta,data)

引数
theta:パラメータ(M,S2)で構成されるベクトル
data: （標本）観測値のベクトル

値
対数事後分布（カーネル）の値

# 関数normchi2postの中身
normchi2post = function (theta, data)
{
    mu = theta[1]
    　 #thetaの第1要素（平均パラメータ）をmuへ代入。
    sig2 = theta[2]
    　 #thetaの第2要素（分散パラメータ）をsig2へ代入。

    logf = function(y, mu, sig2){

              -(y - mu)^2/2/sig2 - log(sig2)/2

     }# 正規分布N(mu,sig2) の密度関数（カーネル）の対数
      #  ベクトルyの各要素の値について評価したベクトル
    
    z = sum(logf(data, mu, sig2))
      #上で定義した関数logfの引数にnormchi2postの引数dataを使う。
      # ベクトルlogfの要素を合計。
    
    z = z - log(sig2)
      # sig2の無情報事前分布に関する部分を考慮。
      #  対数事後分布カーネルの値をzに代入。
    return(z)  # zを出力。
}


例

parameter=c(25,5)
 #mu=25, sig2=5と設定。
data=c(20, 32, 21, 43, 33, 21, 32)
 #dataの保存。
normchi2post(parameter,data)
 #上記のパラメータとデータに基づく対数事後分布（カーネル）の値。

[1] -61.54247

data(marathontimes)
#20代の男性20人のニューヨークマラソンのタイム（分）
# 変数 マラソンのタイム
#出所     www.nycmarathon.org website.

data = marathonstimes$time
　# マラソン時間をdataに保存。

d=mycontour(normchi2post,c(220,330,500,9000),
  data, xlab="mean",ylab="variance")
　# x軸　220~330
  # y軸　500~9000の範囲で
  # マラソンのタイムを使ったnormchi2postの等高線を描く。

set.seed(66)
post.sim = normpostsim(data,m=1000)
 # dataを前提として，事後分布から(mu,sigma2)をm=1000個発生させて
 #  post.simに保存する。（関数normpostsimを利用。）

points(post.sim$mu,post.sim$sigma2)
 # post.simのmuとsigma2を描画した等高線グラフに追加。

関数解読reg.gprior.post

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)


Zellnerのg-prior

beta (given sigma2) ~N(b0,c0*sigma2*(X'X)^(-1))
p(sigma2)proportional to (sigma2)^(-1)

の場合の，回帰モデルの事後分布のカーネルの対数を計算。

注
　この関数の出力が関数laplaceを用いて対数周辺尤度を計算するための部品になる。
　BCwithR p221

使用法

reg.gprior.post(theta, dataprior)

引数
theta:対数事後分布（のカーネル）を評価するパラメータ(beta,logsigma)の値
dataprior:データ(y,X)とg-priorのパラメータ（b0,co）をまとめたリスト

値
対数事後分布（のカーネル）の値

#関数reg.gprior.postの中身
reg.gprior.post = function (theta, dataprior)
{
    y = dataprior$data$y
    X = dataprior$data$X
     #リストdatapriorから被説明変数と説明変数を取り出す。
     #被説明変数（ベクトル）をyに代入。
     #説明変数（行列）をXに代入。

    c0 = dataprior$prior$c0
    beta0 = dataprior$prior$b0
     #リストdatapriorからg-priorのパラメータを取り出す。
     #betaに関する事前情報の信頼性パラメータをc0に代入。
     #betaの位置に関する事前情報をbeta0に代入。
     # BCwithR p218

    beta = theta[-length(theta)]
     # length(theta) ベクトルthetaの長さ
     #  logsigmaの位置を表す数字として使える。
     # theta[-length(theta)]
     #  thetaからlogsigmaを除いたbetaに関する要素を
     #   betaに代入。
     #  このやり方で任意の長さのbetaに対応できる。

    sigma = exp(theta[length(theta)])
     # theta[length(theta)]
     #  thetaの中のlogsigmaの部分の指数を計算して
     #   sigmaに代入。


    loglike =
    sum(dnorm(y, mean = X %*% as.vector(beta),
         sd = sigma,log = TRUE))
     # 対数尤度を計算してloglikeに代入。
     # y_i~N(x_ibeta,sigma2)

    logprior = dmnorm(beta, mean = beta0, varcov = c0 * sigma^2 *
        solve(t(X) %*% X), log = TRUE)
     # betaに関する対数事前分布を計算してlogpriorに代入。
     # beta~N(beta0,c0*sigma^2*(X'X)^(-1))
     # sigmaに関する事前分布1/sigma^2は省略？

    return(loglike + logprior)
     # 対数事後分布のカーネル（loglike+logprior）を返す。
}

例

data(puffin)
#データの読み込み
#Nest
# 作られた巣の度数？    9平米当たりの穴の数
#Distance
#崖の端からの距離（メートル）
#出所
# Peck, R., Devore, J., and Olsen, C. (2005),
# Introduction to Statistics And Data Analysis,
#  Thomson Learning.

#被説明変数の設定。
y=puffin$Nest
#説明変数の設定。
X=cbind(1,puffin$Distance)

#関数reg.gprior.postの引数を設定
data0=list(y=y, X=X)
#g-priorパラメータの設定
prior0=list(b0=c(0,0), c0=10)
theta0 = c(20,-.5,1)

reg.gprior.post(theta0,list(data=data0,prior=prior0))

[1] -147.6524

関数解読bayesresiduals

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)

Bayesian Computation with R: Second Edition (Use R!)

無情報事前分布を用いた線形回帰モデルに関する（予め定めた）あるカットオフ値をベイズ残差が超える事後確率の計算

使用法

bayesresiduals(lmfit,post,k)

引数
lmfit:関数lmの出力
post：リスト形式で保存された事後分布からシミュレートした値。
　beta: シミュレートされた回帰係数（m×k（説明変数の数））
　sigma: シミュレートされた誤差項の標準偏差（m次元）
　関数blinregの値（出力）がそのまま使える。
k:外れ値を定義するためのパラメータ

値
各サンプルが外れ値である事後確率ベクトル


y = Xbeta + e
y_i = x_i + e_i

y_i: yの第i要素
x_i: Xの第i行列

無情報事前分布を前提とすると，
betaの（y,X,sigma^2を所与とした）事後分布は

beta|y,X,sigma^2 ~ N(betahat,sigma^2*(X'X)^(-1))
betahat = (X'X)^(-1)*X'y

となる。
事後的にはy,Xはすでに定まった値として扱えるので，

ベクトル表記

e|y,X,sigma^2 ~ N(ehat,sigma^2*H)
ただし，
ehat = y - Xbetahat
H = X(X'X)^(-1)X' （hat matrix）

スカラー表記

e_i|y,X,sigma^2 ~ N(ehat_i,sigma^2*h_ii)

ただし，
ehat_i = y_i - x_i*betahat
h_ii = x_i(X'X)^(-1)x_i' （hat matrixの(i,i)要素）

が成立する。
従って，

z_i = (e_i - ehat_i)/(sigma*sqrt(h_ii))

を定義すると

z_i |y,X,sigma^2 ~ N(0,1)

が成立する。

注意
無情報事前分布の場合のみ適用可能な関数。


参考 hat matrixについて
工学のためのデータサイエンス入門　p127
浅野中村　計量経済学[第2版] p63 式(RM-5)
AEwithR p97

工学のためのデータサイエンス入門―フリーな統計環境Rを用いたデータ解析 (工学のための数学)

工学のためのデータサイエンス入門―フリーな統計環境Rを用いたデータ解析 (工学のための数学)

計量経済学 (y21)

計量経済学 (y21)

Applied Econometrics with R (Use R!)

Applied Econometrics with R (Use R!)

#関数bayesresidualsの中身
bayesresiduals = function (lmfit, post, k)
{
    ehat = lmfit$residuals
     # lmの出力から残差を取り出してehatに保存
    h = hat(model.matrix(lmfit))
     # lmの出力からhat matrixの対角項を取り出してhに保存。
     #  diag(X %*% solve(t(X)%*%X) %*% t(X))と結果が一致することを確認。
     #  関数hatと関数model.matrixについて要確認。

    prob = 0 * ehat #ehatと長さが同じのゼロベクトルを設定。

    for (i in 1:length(prob)) {
     #それぞれのサンプルについて

        z1 = (k - ehat[i]/post$sigma)/sqrt(h[i])
        # カットオフ値（上限）
        # z1[1]= (k - ehat[i]/post$sigma[1])/sqrt(h[i])
        
        z2 = (-k - ehat[i]/post$sigma)/sqrt(h[i])
        # カットオフ値（下限）
        # z2[i] = (-k - ehat[i]/post$sigma)/sqrt(h[i])

        prob[i] = mean(1 - pnorm(z1) + pnorm(z2))
        # Pr(e_i >= k*sigma or e_i =< -k*sigma|y,X)を
        #  モンテカルロ積分で計算。BCwithR p208
    }

    return(prob)# probを出力。
}

#例
#説明変数Xと被説明変数の設定
chirps=c(20,16.0,19.8,18.4,17.1,15.5,14.7,17.1,15.4,16.2,15,17.2,16,17,14.1)
temp=c(88.6,71.6,93.3,84.3,80.6,75.2,69.7,82,69.4,83.3,78.6,82.6,80.6,83.5,76.3)
X=cbind(1,chirps)
y=temp

#サンプリング回数　10万回に設定
m=100000

#関数blinregを用いて(beta,sigima)の
#事後分布からサンプリングしてpostに保存。
#（事前分布は無情報事前分布）

post = blinreg(y,X,m)


#yをXで線形回帰した結果をlmfitに保存。
lmfit=lm(y~X-1)

# 外れ値判定基準をsigmaの2倍に設定。
k=2

# 外れ値事後確率の計算。
prob.out = bayesresiduals(lmfit,post,k)

# 上記事後確率が0.05を超えるサンプルかどうかの判定。
out = (prob.out > 0.05)

# 説明変数chirpsと外れ知事後確率の散布図の表示。
plot(chirps,prob.out)
#外れ値を識別する添字を表示。
text(chirps[out],prob.out[out],pos=4)

#散布図と（betaの事後平均を用いた）回帰曲線の表示。
plot(y~chirps)
abline(a=lmfit$coef[1],b=lmfit$coef[2])
#外れ値を識別する添字を表示。
text(chirps[out],y[out],pos=4)
# 3つ外れ値がある模様。


関数解読blinregpred

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)

Bayesian Computation with R: Second Edition (Use R!)

無情報事前分布を用いた線形回帰モデルの
（被説明変数に関する）予測分布からのサンプリング


使用法

blinregpred(X1,theta.sample)

引数
X1：行列（行が説明変数の組に対応）
　　推定に用いなかった説明変数の組でもOKのはず。
theta.sample：リスト形式で保存された事後分布からシミュレートした値。
　beta: シミュレートされた回帰係数（m×k（説明変数の数））
　sigma: シミュレートされた誤差項の標準偏差（m次元）
　関数blinregの値（出力）がそのまま使える。


値
各列がある説明変数の組x*に対応する，
（被説明変数の）予測分布からのシミュレートされた値。


参考
BCwithR p207

#blinregpredの中身

blinregpred = function (X1, theta.sample)
{
    d = dim(X1)# 行列X1の次元（c(行数，列数)）をdに代入。
    n1 = d[1] #  dの第1要素（X1の行数）をn1に代入。
    m = length(theta.sample$sigma)
     # リストtheta.sampleの中のsigmaの要素数をmと設定。

    y1 = array(0, c(m, n1))
    # m行n1列の配列を作成。
    # arrayではなくてmatrixにしても良いのか？

    for (j in 1:n1) {
        y1[, j] = t(X1[j, ] %*% t(theta.sample$beta)) +
        rnorm(m) * theta.sample$sigma
    }
    # 第1項　t(X1[j, ] %*% t(theta.sample$beta))
    #  は 関数linregexpectedと同じ。
    # 第2項　rnorm(m) * theta.sample$sigma
    #  は第i要素が正規分布N(0,($sigma[i])^2)から発生させた値　
    return(y1)  #y1を出力。
}


例

#説明変数Xと被説明変数の設定
chirps=c(20,16.0,19.8,18.4,17.1,15.5,14.7,17.1,15.4,16.2,15,17.2,16,17,14.1)
temp=c(88.6,71.6,93.3,84.3,80.6,75.2,69.7,82,69.4,83.3,78.6,82.6,80.6,83.5,76.3)
X=cbind(1,chirps)
y=temp

#サンプリング回数　10万回に設定
m=100000

#関数blinregを用いて(beta,sigima)の
#事後分布からサンプリングしてtheta.sampleに保存。
#（事前分布は無情報事前分布）

theta.sample=blinreg(y,X,m)

#説明変数の組を設定
covset1=c(1,15)
covset2=c(1,20)

X1=rbind(covset1,covset2)

#被説明変数の予測分布からの値をypredに保存。
ypred = blinregpred(X1,theta.sample)

#被説明変数の期待値の事後分布からの値をEyに保存。
Ey=blinregexpected(X1,theta.sample)

#説明変数ごとに
#被説明変数の予測分布と期待値の事後分布のヒストグラムを描画。
# 赤の垂線　予測分布の最小値と最大値
# 緑の垂線　期待値の事後分布の最小値と最大値
# 当然，赤垂線の間の方が広い。

par(mfrow=c(2,2))
hist(ypred[,1],freq=F,breaks="Scott",col="grey",
main="covset1=c(1,15)",xlim=range(ypred))
abline(v=min(ypred[,1]),col="red")
abline(v=max(ypred[,1]),col="red")
abline(v=min(Ey[,1]),col="green")
abline(v=max(Ey[,1]),col="green")


hist(ypred[,2],freq=F,breaks="Scott",col="grey",
main="covset2=c(1,20)",xlim=range(ypred))
abline(v=min(ypred[,2]),col="red")
abline(v=max(ypred[,2]),col="red")
abline(v=min(Ey[,2]),col="green")
abline(v=max(Ey[,2]),col="green")

hist(Ey[,1],freq=F,breaks="Scott",col="grey",
main="covset1=c(1,15)",xlim=range(ypred))
abline(v=min(ypred[,1]),col="red")
abline(v=max(ypred[,1]),col="red")
abline(v=min(Ey[,1]),col="green")
abline(v=max(Ey[,1]),col="green")

hist(Ey[,2],freq=F,breaks="Scott",col="grey",
main="covset2=c(1,20)",xlim=range(ypred))
abline(v=min(ypred[,2]),col="red")
abline(v=max(ypred[,2]),col="red")
abline(v=min(Ey[,2]),col="green")
abline(v=max(Ey[,2]),col="green")

関数解読blinregexpected

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)

無情報事前分布を用いた線形回帰モデルの（被説明変数の）期待値に関する事後分布からのサンプリング

使用法

blinregexpected(X1,theta.sample)

引数
X1：行列（行が説明変数の組に対応）
　　推定に用いなかった説明変数の組でもOKのはず。
theta.sample：リスト形式で保存された事後分布からシミュレートした値。
　beta: シミュレートされた回帰係数（m×k（説明変数の数））
　sigma: シミュレートされた誤差項の標準偏差（m次元）
　関数blinregの値（出力）がそのまま使える。

値
各列がある説明変数の組x*に対応する
被説明変数の期待値からシミュレートされた値x*beta
であるような行列

ある所与のx*に対して，x*betaはbetaの関数である。

参考
BCwithR p207

#blinregexpectedの中身

blinregexpected = function (X1, theta.sample)
{
    d = dim(X1)# 行列X1の次元（c(行数，列数)）をdに代入。
    n1 = d[1] #  dの第1要素（X1の行数）をn1に代入。
    m = length(theta.sample$sigma)
     # リストtheta.sampleの中のsigmaの要素数をmと設定。
     # betaは行列　sigmaはベクトルだから後者を使った模様。
     
    m1 = array(0, c(m, n1))
    # m行n1列の配列を作成。
    # arrayではなくてmatrixにしても良いのか？

    for (j in 1:n1) {
        m1[, j] = t(X1[j, ] %*% t(theta.sample$beta))
    }
    # X1のj行の説明変数の組を用いた場合の（被説明変数）の期待値
    # の事後分布からのサンプルを配列m1のj列に保存。
    # for文を使わずに書けないか。

    return(m1) #m1を出力。
}


blinregexpectedとblinregpredはどう使い分ける？
BCwithR p213-217の例で確認する。

例

#説明変数Xと被説明変数の設定
chirps=c(20,16.0,19.8,18.4,17.1,15.5,14.7,17.1,15.4,16.2,15,17.2,16,17,14.1)
temp=c(88.6,71.6,93.3,84.3,80.6,75.2,69.7,82,69.4,83.3,78.6,82.6,80.6,83.5,76.3)
X=cbind(1,chirps)
y=temp

#サンプリング回数　10万回に設定
m=100000

#関数blinregを用いて(beta,sigima)の
#事後分布からサンプリングしてtheta.sampleに保存。
#（事前分布は無情報事前分布）

theta.sample=blinreg(y,X,m)


#説明変数の組を設定
covset1=c(1,15)
covset2=c(1,20)

X1=rbind(covset1,covset2)

#被説明変数の期待値の事後分布をEyに保存。
Ey=blinregexpected(X1,theta.sample)

#説明変数ごとに被説明変数の期待値のヒストグラムを描画。
par(mfrow=c(2,1))
hist(Ey[,1],freq=F,breaks="Scott",col="grey",
main="covset1=c(1,15)",xlim=range(Ey))
hist(Ey[,2],freq=F,breaks="Scott",col="grey",
main="covset2=c(1,20)",xlim=range(Ey))

関数解読normpostpred

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)


正規分布のパラメータ（mu,sigma2）の無情報事前分布を前提にして
そのパラメータの事後分布からシミュレーションした上で
その事後分布に基づき，関心のある統計量の事後予測分布を求める。

使用法

normpostpred(parameters,sample.size,f=min)


引数


引数
parameters:
関数normpostsimで出力した結果
muとsigma2のサンプル

sample.size:事後予測分布からのサンプル数

f: 統計量を定義する関数（デフォールトはなぜか最小値）

値
定義した統計量の事後予測分布からのサンプル
（normpostsimの結果と同じ長さ）


参考
Bayesian Data Analysis, Second Edition (Chapman & Hall/CRC Texts in Statistical Science)
p77

#normpostpredの中身
normpostpred2 = function (parameters, sample.size, f = min)
{
    normalsample = function(j, parameters, sample.size){
     rnorm(sample.size,
        mean = parameters$mu[j], sd = sqrt(parameters$sigma2[j]))
        }
    #この関数内でのみ使う関数normalsampleを定義
    # リストparametersのmu[j]を平均に
    # リストparametersのsigma2[j]を分散として
    # sample.sizeだけ正規分布から乱数を発生させて出力する。

    m = length(parameters$mu) #事後分布のサンプルの長さ（要素数）をmに代入。
    
    post.pred.samples =
       sapply(1:m, normalsample, parameters,sample.size)
       # sample.size行×m列の行列
       # j列はN(mu[j],sigma2[j])からsample.sizeだけ発生させたサンプル
       
    stat = apply(post.pred.samples, 2, f)
           # 各列の要素に対して関数fを実行。
           # 結果としてm次元のベクトルが出来るので，statに保存。

    return(stat) #statを出力
}


例

# 未来のサンプル15個に関する最小値の事後予測分布

data(darwin)
s=normpostsim(darwin$difference,m=10000)
sample.size=15
sim.stats=normpostpred(s,sample.size,min)

hist(sim.stats,breaks="Scott",col="grey")


[[]]

関数解読blinreg

Rメモ LearnBayes

Bayesian Computation with R: Second Edition (Use R!)


データ分布

y~N(Xbeta,sigma2*I)

無情報事前分布

p(bea,sigma2)proportional to (sigma2)^(-1)

または
Zellnerのg-prior

beta (given sigma2) ~N(b0,c0*sigma2*(X'X)^(-1))
p(sigma2)proportional to (sigma2)^(-1)

の場合の(beta,sigma)の事後分布からのサンプリング


この場合，解析的な結果を得ることができるので，
サンプリングしなくても良い（はずだが）。


使用法

blinreg(y,X,m,prior=NULL)


引数
y: 被説明変数ベクトル
X: 計画行列（説明変数行列）
m: 発生させるサンプル数
prior:指定しなければ無情報事前分布
リスト(c0,b0)を指定することで，事前分布をZellnerのg-priorにできる。


値
beta: シミュレートされた回帰係数（m×k（説明変数の数））
sigma: シミュレートされた誤差項の標準偏差（m次元）

#blinregの中身

blinreg = function (y, X, m, prior = NULL)
{

    if (length(prior) > 0) {
      #引数にpriorがあれば（g-priorを利用），
      #以下の通り保存。
        c0 = prior$c0 #betaに関する事前情報の信頼性
        beta0 = matrix(prior$b0, c(1, length(prior$b)))
         #betaの位置に関する事前情報
    }

    fit = lm(y ~ 0 + X)
     #yをXで回帰して最小二乗法で推定。
     #その係数などをfitに保存。
     #「0+X」とすることで切片項を別途推定しない（Xに含むから）。
     #「-1+X」とか「X-1」でも良かったはず。

    bhat = matrix(fit$coef, c(1, fit$rank))
     #回帰係数を1行fit$rank列の行列に保存。
     #fit$rankは(X'X)^(-1)の次元（のはず）
     # lmのヘルプ　rank: the numeric rank of the fitted linear model.

    s2 = sum(fit$residuals^2)/fit$df.residual
     #回帰係数を1行fit$rank列の行列に保存。
     #（回帰残差二乗和/自由度）をs2に代入。
     # lmのヘルプ　df.residual: the residual degrees of freedom.
     # 自由度　標本数-（説明変数の数）

    if (length(prior) == 0) {
       #無情報事前分布を用いる場合のサンプリングのためのパラメータを設定。
       # BCwithR p206
       # BDA p 356
        shape = fit$df.residual/2
        rate = fit$df.residual/2 * s2
        beta.m = bhat
        vbeta = vcov(fit)/s2  #要注意 これで(X'X)^(-1)が計算できる？
    }
    else {
       #Zellnerのg-Priorを用いる場合のサンプリングのためのパラメータを設定。
       # BCwithR p219
        shape = length(y)/2
        rate = fit$df.residual/2 * s2 + (beta0 - bhat) %*% t(X) %*%
            X %*% t(beta0 - bhat)/2/(c0 + 1)
        beta.m = c0/(c0 + 1) * (beta0/c0 + bhat)
        vbeta = vcov(fit)/s2 * c0/(c0 + 1)
          #要注意 vcov(fit)/s2で(X'X)^(-1)が計算できる？
    }

    sigma = sqrt(1/rgamma(m, shape = shape, rate = rate))
     # 上記のshapeとrateの設定でガンマ分布からm個サンプリング。
     #  （各要素の）逆数をとって，平方根を取る。
     #  sigmaの事後周辺分布からのサンプルとなる。
    
    beta = rmnorm(m, mean = rep(0, fit$rank), varcov = vbeta)
     # LearnBayesの関数rmnormを使って
     #  N(0,vbata)からm個サンプリング。

    beta = array(1, c(m, 1)) %*% beta.m +
                array(sigma, c(m, fit$rank)) * beta
     #  sigma*beta+beta.mを
     #   標準偏差をを1からsigmaに変えて，平均beta.mだけ平行移動して
     #   結果をbetaに代入。

    return(list(beta = beta, sigma = sigma))
     #betaとsigmaをリスト形式で保存して，出力する。
}

例

#説明変数Xと被説明変数の設定
chirps=c(20,16.0,19.8,18.4,17.1,15.5,14.7,17.1,15.4,16.2,15,17.2,16,17,14.1)
temp=c(88.6,71.6,93.3,84.3,80.6,75.2,69.7,82,69.4,83.3,78.6,82.6,80.6,83.5,76.3)
X=cbind(1,chirps)
y=temp

#10万回サンプリング
m=100000
#結果をsに保存（事前分布は無情報事前分布）
s=blinreg(y,X,m)

#最小二乗法で推定
fit=lm(y~0+X)

#betaの周辺事後分布を表示（burn-in 最初の1000回）
#MCMCでないときはburn-inを考慮しなくても良いのかな。
#赤線はそれぞれの平均
par(mfrow=c(2,1))
hist(s$beta[1001:m,1],breaks="Scott",col="grey")
abline(v=fit$coef[1],col="red",lwd=2)
hist(s$beta[1001:m,2],breaks="Scott",col="grey")
abline(v=fit$coef[2],col="red",lwd=2)

#sigmaを2乗してsigma2のサンプルを求める。
sigma2=(s$sigma[1001:m])^2


#sigma2の周辺事後分布を表示（burn-in 最初の1000回）
#MCMCでないときはburn-inを考慮しなくても良いのかな。
#赤線は平均
#緑線はモード
hist(sigma2,breaks="Scott",col="grey")

shape = fit$df.residual/2
rate = fit$df.residual/2 * s2

abline(v=rate/(shape-1),col="red",lwd=2)
abline(v=rate/(shape+1),col="green",lwd=2)