関数解読normpostsim
Bayesian Computation with R: Second Edition (Use R!)
データ分布
y~N(mu,sigma2)
事前分布
p(mu,sigma2)proportional to (sigma2)^(-1)
の場合の(mu,sigma2)の事後分布からのサンプリング。
この場合，解析的な結果を得ることができるので，
サンプリングしなくても良い。
使用法
normpostsim(data,m=1000)
引数
data:観測値ベクトル
m: 発生させるサンプル数（デフォールト1000）
値mu: シミュレートされた正規分布の平均のベクトル（m次元）
sigma2: シミュレートされた正規分布の分散のベクトル（m次元）
参考
Bayesian Data Analysis, Second Edition (Chapman & Hall/CRC Texts in Statistical Science)
#normpostsimの中身 normpostsim = function (data, m = 1000) { S = sum((data - mean(data))^2) #dataからその平均を引いて，2乗和を計算。Sに代入。 xbar = mean(data) # dataの平均をxbarに代入。 n = length(data) # dataの長さ（要素数）をnに代入。 sigma2 = S/rchisq(m, n - 1) # Bayesian Data Analysis(BDA) p75の数式(3.5)と本質的には同じ。 # 参考 BDA p574-575の一覧表。 mu = rnorm(m, mean = xbar, sd = sqrt(sigma2)/sqrt(n)) # BDA p75の数式(3.3) return(list(mu = mu, sigma2 = sigma2)) # muとsigma2をリスト形式で保存して，出力。 }
例
#データ #データの読み込み # 自家受精させた植物と交雑受精させた植物の高さの差のデータ # 出所：Fisher, R. (1960), Statistical Methods for Research Workers data(darwin) #dataにdarwin$differenceを利用 #mは未指定（デフォールトm=1000をそのまま使う） #結果をsに代入 s=normpostsim(darwin$difference) #平均のヒストグラムを描画 hist(s$mu) #分散のヒストグラムを描画 hist(s$sigma2)
関数解読regroup
Bayesian Computation with R: Second Edition (Use R!)
指定した整数値ごとに行列の行を集計して新しい行列を作成する。
使用法
regroup(data,g)
引数
data:行列
g:1からdataの行数nの間の正整数
値
gごとに行列dataの行を集計した新しい行列
行列dataの行数がgで割り切れなかった場合は，
行列dataの最後尾の余った行（g-1）がその前のg行といっしょにまとめられる（はず）。
値
#regroupの中身 regroup = function (data, g) { d = dim(data) n = d[1] #行数 m = d[2] #列数 N = floor(n/g)#行数をgで割り，小数点以下を切り捨てる。 dataG = array(0, c(N, m))#N行m列の配列を作る。 k = 0 #dataGの行添字kの初期値を設定。 for (j in seq(1, (N - 1) * g + 1, g)) { #行列dataの行を(g+1)行おきに取り出す。 # 例 g=3 なら1行目の次は4行目 # seq(1, N * g, g) とすると問題が起こる？ k = k + 1 # 1だけ増加。 for (i in 0:(g - 1)){ dataG[k, ] = dataG[k,] + data[j+i,] # dataGのk行目とdataのj+i行目を改めて # dataGのk行目に代入。 # i=0,..,(g-1)まで繰り返す。 # 例 j=1でg=3ならi=0,1,2 # dataの1〜3行目を合計することになる。 # apply(data[j:j+(g-1),2,],2,sum) # を使えば，内側のforは不要かな。 } } if (n > N * g) {# nがgで割り切れなかった場合 for (i in (N * g + 1):n){ dataG[N, ] = dataG[N, ] + data[i,] #dataのN*g+1行目からn行目までを #dataGのN行目に合計して代入。 # apply(data[N * g + 1:n],2,sum)が使えそう。 } } return(dataG)#dataGを出力。 }
例
#1から20までの整数を4行5列の行列とする。 data=matrix(c(1:20),nrow=4,ncol=5) data #まとめる行の数を設定。 g=2 #dataを2行ごとまとめて出力。 regroup(data,2) #1から49までの整数を5行5列の行列とする。 data=matrix(c(1:49),nrow=7,ncol=7) data #まとめる行の数を設定。 g=2 #dataを2行ごとまとめて出力。 #dataの行数がg=2で割り切れないので #最後の3行（5行目から7行目）がまとめられている。 regroup(data,2) ###################### #ジーター選手の2004年の154ゲームの打撃成績 data(jeter2004) attach(jeter2004) #ヒット数と打席数？をまとめてdataに代入。 data=cbind(H,AB) #1試合ごとの(H,AB)を5試合ごとに集計してdata1に代入。 #data1の最終行（30行）は9試合（146〜154試合）の集計データ data1=regroup(data,5)
関数解読simcontour
Bayesian Computation with R: Second Edition (Use R!)
グリッド上で定義された一般的2変量密度関数より乱数を発生させる
使用法
simcontour(logf,limits,data,m)
引数
logf:適用する対数密度を定義した関数
limits:グラフを描く領域c(xlo,xhi,ylo,yhi)を指定
xlo:x軸の最小値
xhi:x軸の最大値
ylo:x軸の最小値
yhi:x軸の最大値
data: 適用する対数密度のパラメータ
m:発生させる乱数の数
値
#simcontourの中身 simcontour = function (logf, limits, data, m) { #simcontourの中でのみ使う関数LOGFを定義 # mycontourの場合と同じ LOGF = function(theta, data) { if (is.matrix(theta) == TRUE) { #thetaが行列ならば， val = matrix(0, c(dim(theta)[1], 1)) # ゼロ行列valを作成 # 行数：thetaの行数 # 列数：1 for (j in 1:dim(theta)[1]){ val[j] = logf(theta[j,], data) }#thetaのj行とlogfのパラメータdata # を使って，対数密度logfを計算し，valのj行に代入 #for文を使わずに書けないか。 } else val = logf(theta, data) #thetaが行列でないならば， # thetaとlogfのパラメータdata # を使って，対数密度logfを計算し，valに代入 return(val) #valの値を出力する。 } # (x,y)平面の与えた最小から最大の範囲を49等分 ng = 50 x0 = seq(limits[1], limits[2], len = ng) #最小値がxlo=limits[1],最大値がxhi=limits[2] #長さがngの等差数列をx0に代入。 #この数列はxloとxhiの間をng-1等分する。 y0 = seq(limits[3], limits[4], len = ng) #最小値がylo=limits[3],最大値がyhi=limits[4] #長さがngの等差数列をy0に代入。 #この数列はyloとyhiの間をng-1等分する。 X = outer(x0, rep(1, ng)) #x0とrep(1, ng)（全要素が1の長さngのベクトル） # の外積をXに代入。x軸用のグリッドの準備。 Y = outer(rep(1, ng), y0) #rep(1, ng)とx0 # の外積をYに代入。y軸用のグリッドの準備。 n2 = ng^2 Z = LOGF(cbind(X[1:n2], Y[1:n2]), data) #thetaが行列でないなら，点thetaにおけるlogfの値をZに出力。 #thetaが行列なら，thetaの各行で表される点におけるlogfの値を # n2行1列の行列として # Zに出力。 Z = Z - max(Z) #ZからZの要素の最大値を引いて，新たにZに代入。 # 対数密度の頂点（モード）の高さを0に基準化。 Z = matrix(Z, c(ng, ng)) #再び、Zをng×ngの行列に変換。 ###ここまで関数mycontourと全く同じ。 d = cbind(X[1:n2], Y[1:n2], Z[1:n2]) # 列ベクトルX[1:n2]，Y[1:n2]，Z[1:n2]束ねて行列を作成。 # 各行は領域内のグリッドの座標に対応。 dx = diff(x0[1:2]) #x軸のグリッドの間隔（x0の第1要素と第2要素の差）をdxに代入。 dy = diff(y0[1:2]) #x軸のグリッドの間隔（y0の第1要素と第2要素の差）をdyに代入。 prob = d[, 3] #dの3列目（対数密度の高さ）をprobに代入。 prob = exp(prob) #対数を元の数値に戻してprobに代入。 prob = prob/sum(prob)# 高さの相対比率をprobに代入。 # 結果として，sum(prob)=1となる。 i = sample(2500, m, replace = TRUE, prob = prob) #グリッド上の各点（50×50=2500個）に確率が定義されている。 # それぞれの点を識別する添字（1〜2500）を定義された確率に従って # m回復元抽出する。 return(list(x = d[i,1] + runif(m) * dx - dx/2, y = d[i,2] + runif(m) * dy - dy/2 )) #確率（頻度）に応じた添字iをdの行を表す添字として使う。 #グリッドを中心に半径dx（dy）の範囲で乱数をさらに散らすために # 第2項 runif(m) * dx - dx/2やrunif(m) * dy - dy/2 を追加。 }
例
#平均ベクトルの指定 m=array(c(0,0),c(2,1)) #分散共分散行列の指定 v=array(c(1,.6,.6,1),c(2,2)) #上記の2つをリストにまとめる normpar=list(m=m,v=v) #結果を再現するために（疑似）乱数の種を設定。 set.seed(123) #関数simcontourに従って，乱数の組を500個作成してsに代入。 s=simcontour(lbinorm,c(-4,4,-4,4),normpar,500) #sを散布図で表示。 plot(s$x,s$y) #関数mycontourを使って等高線を追加。 mycontour(lbinorm,c(-4,4,-4,4),normpar,add=T)
関数解読mycontour
Bayesian Computation with R: Second Edition (Use R!)
2変量密度関数の等高線の描画
一般的な2変量密度関数用の等高線グラフ
実際には，対数密度関数を使ってグラフを描く。
密度関数の定数部分は除いても描ける。
モード（最頻値）での密度関数の高さを基準にして
その高さの10%，1%，0.1%の点を結んで等高線を描く。
別のより一般的な関数contourを活用する。
使用法
mycontour(logf,limits,data,...)
引数
logf:適用する対数密度を定義した関数
limits:グラフを描く領域c(xlo,xhi,ylo,yhi)を指定
xlo:x軸の最小値
xhi:x軸の最大値
ylo:x軸の最小値
yhi:x軸の最大値
data: 適用する対数密度のパラメータ
... 関数contourに関するさらなる引数
値
密度関数の等高線グラフの描画
#mycontourの中身 mycontour = function (logf, limits, data, ...) { #mycontourの中でのみ使う関数LOGFを定義 LOGF = function(theta, data) { if (is.matrix(theta) == TRUE) { #thetaが行列ならば， val = matrix(0, c(dim(theta)[1], 1)) # ゼロ行列valを作成 # 行数：thetaの行数 # 列数：1 for (j in 1:dim(theta)[1]){ val[j] = logf(theta[j,], data) }#thetaのj行とlogfのパラメータdata # を使って，対数密度logfを計算し，valのj行に代入 #for文を使わずに書けないか。 } else val = logf(theta, data) #thetaが行列でないならば， # thetaとlogfのパラメータdata # を使って，対数密度logfを計算し，valに代入 return(val) #valの値を出力する。 } # (x,y)平面の与えた最小から最大の範囲を49等分 ng = 50 x0 = seq(limits[1], limits[2], len = ng) #最小値がxlo=limits[1],最大値がxhi=limits[2] #長さがngの等差数列をx0に代入。 #この数列はxloとxhiの間をng-1等分する。 y0 = seq(limits[3], limits[4], len = ng) #最小値がylo=limits[3],最大値がyhi=limits[4] #長さがngの等差数列をy0に代入。 #この数列はyloとyhiの間をng-1等分する。 X = outer(x0, rep(1, ng)) #x0とrep(1, ng)（全要素が1の長さngのベクトル） # の外積をXに代入。x軸用のグリッドの準備。 Y = outer(rep(1, ng), y0) #rep(1, ng)とx0 # の外積をYに代入。y軸用のグリッドの準備。 n2 = ng^2 Z = LOGF(cbind(X[1:n2], Y[1:n2]), data) #thetaが行列でないなら，点thetaにおけるlogfの値をZに出力。 #thetaが行列なら，thetaの各行で表される点におけるlogfの値を # n2行1列の行列として # Zに出力。 #Xはベクトルではないが，X[1:n2]はベクトル。 # 不思議。 # X[1:n2]は行列Xの各列をつなぎ合わせたn2次元のベクトル。 # cbind(X[1:n2], Y[1:n2]) # 列ベクトルX[1:n2]とY[1:n2]束ねて行列を作成。 # 各行は領域内のグリッドの座標に対応。 #上で定義した関数LOGFを利用。 Z = Z - max(Z) #ZからZの要素の最大値を引いて，新たにZに代入。 # 対数密度の頂点（モード）の高さを0に基準化。 Z = matrix(Z, c(ng, ng)) #再び、Zをng×ngの行列に変換。 contour(x0, y0, Z, levels = seq(-6.9, 0, by = 2.3),lwd = 2, ...) #関数contour：等高線を描く汎用関数。 # levels Zの値が-6.9,-4.6,-2.3,0の点を等高線で結ぶ。 # exp(0)=1 # exp(-2.3)=約0.1 # exp(-4.6)=約0.01 # exp(-6.9)=約0.001 # lwd:線の幅の指定。デフォールトは1。 ## （ついでに）lty：線のタイプ ## 0=blank, 1=solid (default), 2=dashed, ## 3=dotted, 4=dotdash, 5=longdash, 6=twodash }
例
関数lbinormを使って等高線を描画
#平均ベクトルの指定 m=array(c(0,0),c(2,1)) #分散共分散行列の指定 v=array(c(1,.6,.6,1),c(2,2)) #上記の2つをリストにまとめる normpar=list(m=m,v=v) #等高線の描画 mycontour(lbinorm,c(-4,4,-4,4),normpar) #x軸y軸を追加 abline(v=0,h=0)
#分散共分散行列のみ変更（相関係数=0） v1=array(c(1,0,0,1),c(2,2)) normpar=list(m=m,v=v1) mycontour(lbinorm,c(-4,4,-4,4),normpar) abline(v=0,h=0)
#分散共分散行列のみ変更（相関係数=0.9） v2=array(c(1,.9,.9,1),c(2,2)) normpar=list(m=m,v=v2) mycontour(lbinorm,c(-4,4,-4,4),normpar) abline(v=0,h=0)
関数解読lbinorm
Bayesian Computation with R: Second Edition (Use R!)
2変量正規分布の密度関数の自然対数の計算
定数部分を除くカーネルの自然対数
使用法
lbinorm(xy,par)
引数
xy: 密度関数の対数を評価する点(x,y)
par:平均ベクトルmと分散共分散行列vをリスト形式にまとめたもの
値
対数密度のカーネルの値
#lbinormの中身 lbinorm=function (xy, par) { m = par$m #リストparからmの部分を取り出し，mに代入。 v = par$v #リストparからvの部分を取り出し，vに代入。 x = xy[1] #xyの第1要素をxに代入。 y = xy[2] #xyの第2要素をyに代入。 zx = (x - m[1])/sqrt(v[1, 1]) zy = (y - m[2])/sqrt(v[2, 2]) #xとyを標準化してzxとzyにそれぞれ代入。 r = v[1, 2]/sqrt(v[1, 1] * v[2, 2]) #分散共分散行列vの要素を使って相関係数rを計算。 return(-0.5/(1 - r^2) * (zx^2 - 2 * r * zx * zy + zy^2)) #2変量正規密度の対数（定数項を除く）を出力。 }
例
#平均ベクトル mean=c(0,0) #分散共分散行列 #対角項が1の対角行列 varcov=diag(c(1,1)) #対数密度を評価する点 value=c(1,1) #平均ベクトルと分散共分散行列をリスト形式で保存 param=list(m=mean,v=varcov) lbinorm(value,param)
関数解読rtruncated
Bayesian Computation with R: Second Edition (Use R!)
切断された確率分布からのサンプリング
累積分布関数と（分位数を出力する）逆累積分布関数（分位数）
が利用可能な場合に適用できる。
Rの場合，例えば正規分布なら
累積分布関数はpnorm
分位数を出力する関数はqnorm
(x <- pnorm(1,mean=1,sd=2)) qnorm(x,mean=1,sd=2)
使用法
rtruncated(n,lo,hi,pf,qf,...)
累積分布関数がpfで分位数を出力する関数がqfの確率分布において
lo <=X <= hi
の範囲からn個のサンプルを出力する。
引数
n:発生させる乱数の数
lo:下限の切断点
hi:上限の切断点
pf:（切断前の）累積分布関数
qf:（切断前の）逆累積分布関数
...:pfやqfで使われる引数
値
この関数に適用する切断された確率分布からのn回分のサンプリング結果
#rtruncatedの中身 #関数の引数に関数を使える。 #runif(n)はn>=2の場合，ベクトルで #pf(lo)やpf(hi)はスカラーだから， #ベクトルのrunifにあわせて（各要素に対して）， #スカラーのpf(lo)やpf(hi)が繰り返し計算に使われる。 rtruncated=function (n, lo, hi, pf, qf, ...){ qf(pf(lo, ...) + runif(n) * (pf(hi, ...) - pf(lo, ...)), ...) }
メモ
切断された新しい累積分布関数のX=xにおける値は
(pf(x)-pf(lo))/(pf(hi)-pf(lo))
である。これをF(x)とおいて
一様分布U[0,1]に従う確率変数をuとするとPr(u<=F(x))=F(x)
が成立する。
u=F(x)は
x = pf^(-1)(pf(lo)+u*((pf(hi)-pf(lo))))
と式変形できる。
ここで，pf^(-1)()はpf()の逆関数だからqf()に等しい。結局，
x = qf(pf(lo)+u*((pf(hi)-pf(lo))))
を得る。
例
# 平均2，標準偏差1の正規分布で3以上の範囲から10個サンプリングする。 n=10 lo=3 # 下限には-Infを使える hi=Inf # 上限にはInfを使える。 rtruncated(n,lo,hi,pnorm,qnorm,mean=2,sd=1)
# ベータ分布beta(2, 5)の0.3<=X<=0.8の範囲から20個サンプリングする。 n=20 lo=0.3 hi=0.8 rtruncated(n,lo,hi,pbeta,qbeta,2,5)
関数解読rigamma
Bayesian Computation with R: Second Edition (Use R!)
逆ガンマ分布IG(a,b)からのサンプリング
ガンマ分布G(a,b)に従う確率変数Xの逆数Y=1/Xは
（定義上）逆ガンマ分布に従うという事実を直接用いる。
使用法
rigamma(n, a, b)
引数
- n: 発生させる乱数の数
- a(>0): 逆ガンマ分布のshapeパラメータ（訳がわからない）。
- b(>0): 逆ガンマ分布のrateパラメータ（訳がわからない）。
確認
X〜G(a,b)のとき，
E(X)=a*b
Var(X)=a*(b)^2
（a>0,b>0ならばE(X)とVar(X)は存在する）Y〜IG(a,b)のとき，
E(Y)=1/((a-1)*b)
Var(Y)=1/((a-1)^2*(a-2)*b^2)
（a>0,b>0ならばYをサンプリングすることはできるが
E(X)が存在するのはa>1の場合
Var(X)は存在するのはa>2の場合）http://en.wikipedia.org/wiki/Gamma_distribution
http://en.wikipedia.org/wiki/Inverse-gamma_distribution
（自分への）注意
Wikipediaの逆ガンマ分布の説明でのパラメータbetaはここでの1/bに対応。
値
逆ガンマ分布IG(a,b)からのn回分のサンプリング結果
#rigammaの中身 rigamma=function (n, a, b) { return(1/rgamma(n, shape = a, rate = b)) # G(a,b)からn回サンプリングして， #各値の逆数を計算して出力。 }
例
a=10 b=5 n=20 ig = rigamma(n,a,b) mean(ig) var(ig) #nを増やすと期待値分散の理論値に近づく。 #期待値，分散が存在しないパラメータ設定で # mean(ig)やvar(ig)を計算すると，nを増やすほど # これらの値は大きくなる（発散する）。
