JAGSによる対数正規分布のパラメータ推定
Kruschke氏のblogでは,対数正規分布のデータを生成するのに，一度正規乱数作ってから基準化して指数化するというやや煩雑な手順を用いている．
doingbayesiandataanalysis.blogspot.jp
最初にコードを見たとき，BUGSのクセとあいまって，何をやっているのか分からなかった． 単純化のため，対数正規分布を使って，直接データを生成してみる．
set.seed(47405) y = rlnorm( n=500, 5.0, 0.5) # Package the data for shipping to JAGS: dataList = list( y = y , N = length(y) , meanOfX = mean(log(y)) ,#事前分布のパラメータに使う sdOfX = sd(log(y)))#事前分布のパラメータに使う
JAGSに渡すコードは以下の通り
modelstring = " model { for( i in 1 : N ) { y[i] ~ dlnorm( mu, 1/sigma^2)#JAGS codeの定義より } sigma ~ dunif( 0.001*sdOfX , 1000*sdOfX ) mu ~ dunif( 0.001*meanOfX , 1000*meanOfX ) meanY <- exp(mu+(sigma^2)/2) modeOfY <- exp(mu-sigma^2) sigmaOfY <- sqrt(exp(2*mu+sigma^2)*(exp(sigma^2)-1)) } " # close quote for modelstring writeLines(modelstring,con="model.txt")
一応これでもMCMC推定はうまくいく．
（今思いついたが，事前分布のパラメータの範囲を指定するのに正規分布を1回経由したほうが便利だから，ああいうコードになっていたのかもしれん）
理論モデルから生成したデータに基づく推定
次の手順として，対数正規分布のパラメータを
y_0, b, p, nの関数とみなして，ダイレクトに対数正規分布の確率密度関数を使って推定してみる．
要するに理論モデル（GMI）に基づいてデータを人工的に発生させて，理論上のパラメータ（潜在変数）をMCMC推定できるのか？ という実験である．
まずデータは次のように生成する
#直接対数正規分布のパラメータを指定する set.seed(47406) y0=20; b=0.2; p=0.65; n=10; mu=log(y0)+n*log(1-b)+log((1+b)/(1-b))*n*p sigma= sqrt(n*p*(1-p)*(log((1+b)/(1-b)))^2) y = rlnorm( n=100, mu, sigma)
JAGS用モデルは次の通り(これをモデルと呼ぶことに依然として抵抗を感じるが今は無視する)．
#modelstring0はパラメータをn,y0,p,b,を全て推定 modelstring0 = " model{ for( i in 1 : N ) { y[i] ~ dlnorm(m, 1/(s^2)) } m <- log(y0)+n*log(1-b)+log((1+b)/(1-b))*n*p s <- sqrt(n*p*(1-p)*(log((1+b)/(1-b)))^2) n ~ dpois(a)#nをポアソン分布で近似 a ~ dunif(13,18) y0 ~ dpois(c)# c ~ dunif(10,18) p ~ dbeta(1,1)#p,bの事前分布にβ分布 b ~ dbeta(1,1)}" # close quote for modelstring
この場合真値をHDIに含む推定をしてくれる．ただ関数gmiは人工データとしては，あまりうまく機能しないようだ．
（自由度に問題がないのか，若干気になる．直感的に言って自由度が高すぎてMCMCが収束しない場合はp以外のパラメータを決定論的に
指定する必要がある．）
再現可能な図
校正をしてて気づいた点をメモしておきます．
図は完全にreproducibleなコードを残しておくべきである．
例えば，Mathematicaで独自関数をplotする場合， パラメータの数値なども図と一緒に出力するようなコードを書いておく．
Module[{plot1, point = 18, t1, t2, \[Rho] = 0.5}, t1 = Graphics[ Text[Style["G of Income", FontSize -> point], {0.6, 0.3}]]; t2 = Graphics[ Text[Style["G of Human Capital", FontSize -> point], {0.6, 0.5}]]; plot1 = Plot[{ giniincome[0.5, s1, 1, \[Rho]], ginihuman[s1] }, {s1, 0, 1}, PlotStyle -> {{Black, Thick, Dashed}, {Gray, Thick}}]; Show[plot1, t1, t2, PlotRange -> {0, 1}, Frame -> True, FrameLabel -> {TableForm[{{"\!\(\*SubscriptBox[\(\[Sigma]\), \(1\)]\ \) of Human Capital"}, {"\[Alpha]=0.5, \ \!\(\*SubscriptBox[\(\[Sigma]\), \(2\)]\)=1, \[Rho]=0.5"}}], "Inequality (Gini index)"}, ImageSize -> Large, LabelStyle -> {FontSize -> point}, AxesStyle -> Directive[point], GridLines -> Automatic, ImageSize -> Medium] ]
ここでは，凡例をplot内に書き込み，パラメータ値をフレーム外に表示させている．ここでTeXやWordなどで 編集している本文のキャプションではなく，計算用ソフト上のグラフィックオブジェクトとして， テキスト文字を埋め込んでしまうのが，ポイントだ．
同じことはggplot2を使えばRでもできる．
こうすれば図を変更した後にキャプションをいちいち修正しなくて済む．なぜなら図を変更したとき コードも修正されているからだ（本当は関数の引数とキャプションの値が同期しているのが望ましい）．
実行例はこんな感じだ．
複数の線が異なるパラメータに対応しているなら，こんな描き方でも分かりやすいと思う．
Mathematicaには出力した図に，オブジェクトをGUI処理で重ね書きする機能があるが， 多少面倒でもグラフィックオブジェクトをコードで書いて，最後にShowでまとめて表示 したほうがよい．
急がば回れ，というやつだろう．要は「未来の自分」が内容を忘れていても， Enterキー一発で，再現できるようなコードを残しておくと，後でラクなのだ．
■
前回の続き．
そもそも分布を導出する理論モデルのパラメータを階層ベイズMCMCで推定できるのか，という確認をやってみる．
まずデータとして，理論モデルからのサンプルをつくる．
今回はJAGSをRから呼び出すので, コードは全部Rで済むように書き換える．
data生成用コードはざっと以下のような感じ
# Generate a person's income under the given condition: set.seed(47406) gmi <- function(y1, b0, p0, n0){ y = y1 for(i in 1:n0){ if(runif(1) < p0){ y <- y + (y*b0)} else{y <- y - (y*b0)} }#for end y } #generate income data for m people m=2000; y1=10; b0=0.2; p0=0.5; n0=15 dat0 <- c() for(i in 1:m){ dat0[i] <- gmi(y1, b0, p0, n0)}
ここで定義しているパラメータが真の値に相当する
次にJAGSに渡すモデルの記述
modelstring0 = " model{ for( i in 1 : N ) { y[i] ~ dlnorm(m, 1/(s^2)) } m <- log(y0)+n*log(1-b)+log((1+b)/(1-b))*n*p s <- sqrt(n*p*(1-p)*(log((1+b)/(1-b)))^2) n ~ dpois(a)#nをポアソン分布で近似 a ~ dunif(1,20) y0 ~ dpois(c)# c ~ dunif(1,20) p ~ dbeta(1,1)#p,bの事前分布にβ分布 b ~ dbeta(1,1)}"
実行には例によってJohn K. Kruschke.さんのラッパー関数DBDA2E-utilities.Rをお借りする． （MCMCの設定はデフォルト）
さて，実行してみると，chainは収束してるみたいだけど，
なんかうまく推定できていない気がする．
よくよく考えたら，nだのy_0だの理論的に定数であるべきパラメータを（非負だからといって）
ポアソン分布に従うとか雑に仮定したのがよくないのでは？
というわけで修正した統計モデルがこちら
modelstring1 = " model{ for( i in 1 : N ) { y[i] ~ dlnorm(m, 1/(s^2)) } m <- log(y1)+n0*log(1-b)+log((1+b)/(1-b))*n0*p s <- sqrt(n0*p*(1-p)*(log((1+b)/(1-b)))^2) #prior for hyper parameters p ~ dbeta(1,1) b ~ dbeta(1,1)}"
p, b以外は定数扱いとして，事前に定義しておく．推定しているのは実質的には
p, bのみである.具体的にはpが投資成功確率でbが投資する資本割合である．
で，結果は
こんな感じでpの事後分布のモードは真値に近いけど，bの方はあまりよくない．
中心極限定理の観点から考えるとpの分散は0でなくてもいいけど（リアプノフの条件から直感的にそう思う），bのほうが理論モデルにあっていないのかもしれない．つまり，そもそもbが確率分布していると，対数正規分布をGMIから導出できないのかもしれない．これでは元も子もない．
ここまでの疑問を整理すると
- ハイパーパラメータはそもそも何個指定できるのか？
- 実際に所得データからパラメータを推定する際に「足りない情報」をどの変数で補うのか？
- 分布Aのパラメータが階層ベイズの仮定として（Aと共役でない）確率分布Bに従う場合，データがAに従うという仮定は妥当か？
といったところ．
続く
■
ずっと気になっている（しかしマイナな）問題が解決しないので，そのメモを記録しておきます．
所得分布の生成モデルは，投資行動の反復から対数正規分布を導出する数理モデルである．仮定は次のとおりである（Hamada 2003; 2004; 2016）.
- 人々が最初に初期資産
y_0を持っている． - 資産のうち割合
bを投資して，確率pで成功すると，資産が増える．失敗すると投資した資産を失う． - 資産の一定割合を所得として得る．人々が投資行動を
n回繰り返すと，資産と所得の分布は対数正規分布にしたがう．
このモデルの結論の一つは，導出された対数正規分布のパラメータがモデルの外生的パラメータであるy_0, b, p, nによって内生的に決まる，というものである．
具体的には，資産の分布は
となる．証明はHamada (2003; 2004; 2016）をみよ．
これは2項分布が正規分布
で近似できる，という性質を利用した命題である（確率変数を基準化した場合にはド・モアブル=ラプラスの定理となる）．
さて
このモデルから得た命題をシミュレーションによって再現できない，というのが問題だった．
どういうことか説明しよう．
まず次の関数によって，確率pで資産y_0 bを増やし，1-pで資産y_0 bを失う仮定を表現する．コードは例によって一番楽に書けるMathematica．
gmi[y0_, b_, p_, n_] := Module[{y = y0}, Do[(* repeat n time *) If[Random[] < p, y = y + (y*b), y = y - (y*b)], {n}]; y];
一方で，理論的に導出される定まる確率分布のパラメータを使って，対数正規分布に従う疑似乱数をコンピュータに作らせる．
(* generate samples from a pdf with parameters derived from Hamada model *) lognsamples[y0_, b_, p_, n_, pop_] := Module[{m, s, a, B }, a = Log[(1 + b)/(1 - b)]; B = Log[y0] + n Log[1 - b]; m = B + a n p; s = Sqrt[n p (1 - p) a^2]; RandomVariate[LogNormalDistribution[m, s], {pop}]];
モデルが正しく対数正規分布を導出できているのであれば，ミクロレベルでの資産生成プロセスを集約した結果と，分布からサンプリングした結果は，ほぼ等しくなるはずである．
modelfit[y0_, b_, \[Gamma]_, n_, pop_] := Module[{dat0, dat1}, (* generate samples from a model *) dat0 = Table[gmi[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*)], {pop(* population *)}]; dat1 = lognsamples[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*), pop]; Histogram[{dat0, dat1}, Automatic, "Probability", ChartStyle -> {Red, Blue}]];
ところが,この関数を実行すると，「理論値から得たサンプルデータ」と「疑似乱数で作ったサンプルデータ」がズレる．
ズレの解消の仕方として，疑似乱数でサンプルデータを作るのではなく，確率密度関数を理論値で得た離散的所得階級値間のレンジで積分して，確率密度を少し「まとめて」やれば，うまくいくことは分かっている．
だが，いちいち積分させるのもコードが煩わしい．
先日，ふと気づいた．疑似乱数から生成したデータベクトルの中には同じ値がdouble精度で1つもないのに，理論に基づいて生成した方は，double精度で同じ値がいくつも含まれているのだ．ようするに，理論値のベクトルの方が疑似乱数のベクトルよりも分散が小さいのである．
これは一体，どういうことなのか？
理論値は理論値で間違っていない．疑似乱数ベクトルも指定したパラメータを持つ分布からランダムにサンプルを抽出している．
分布を導出しているはずの数理モデルから生成したサンプルが，なぜPCで発生させた確率分布疑似乱数と分布がずれてしまうのか？
そこで理論値は厳密すぎて，疑似乱数のランダムさを失っていると考えてみる．
確かに関数gmi[y0_, b_, p_, n_]の中で一様乱数を使っているので投資結果はランダムに変化するのだが，獲得した資本額がある決まった値しかとらない．だから理論値にノイズを足してやればいいのでは？ というわけだ．
別の言い方をすれば疑似乱数の方は連続確率変数の実現値を表現しているのに，理論モデルは離散的確率変数の実現値しか表現できていない．だから確率密度関数を区間毎に積分する操作が，実現値の離散化に対応していたのである．よって修正はただノイズ
RandomVariate[NormalDistribution[0, 0.1]]
を加えるだけでよいのでは，と考えた.
修正したコードはこうなる
lognsamples[y0_, b_, p_, n_, pop_] := Module[{m, s, a, B}, a = Log[(1 + b)/(1 - b)]; B = Log[y0] + n Log[1 - b]; m = B + a n p; s = Sqrt[n p (1 - p) a^2]; RandomVariate[LogNormalDistribution[m, s], {pop}]]; gmi[y0_, b_, p_, n_] := Module[{y = y0}, Do[(* repeat n time *) If[Random[] < p, y = y + (y*b), y = y - (y*b)], {n}]; y]; modelfit[y0_, b_, p_, n_, pop_] := Module[{dat0, dat1}, dat0 = Table[gmi[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*)] + RandomVariate[NormalDistribution[0, 1/Sqrt[n]]], {pop(* population *)}]; (* add randome error since the model is too rigorous *) dat1 = lognsamples[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*), pop]; Histogram[{dat0, dat1}, ChartStyle -> {Red, Blue}]];
これで，だいぶマシにはなる．だがパラメータの値によってノイズの分散に微妙に調整が必要なので 完璧ではない．
多分，記録しておかないと数日したら忘れちゃうので，書き留めておきます．
順番を決めるコード
新学期がはじまるので，演習報告担当者の順番決め用のあみだくじコード（Mathematica）を再掲しておきます．
(* 動的表現 サイズ固定*) Manipulate[ Dynamic[Refresh[RandomSample[Range[n]], UpdateInterval -> x]], {{x, Infinity, "シャッフル"}, {0 -> "Go", Infinity -> "Stop"}}, {{n, 20, "人数"}, 2, 100, 1}, ContentSize -> {600, 150}]
ダイナミックを使ってるので，動的に更新します． この短さでGUIまで自動的に生成するとは驚きである． やっぱりMathematicaで書くとコードがシンプルでええなあ．
探索的分析の自動化実験
内容
データフレームを突っ込むと，データに含まれる全変数の中から1つを順番に応答変数として選び，残りの変数に回するような分析を自動的にやってくれるコードを書いてみた．
「本格的な分析の前に，変数間の関係を探索的・網羅的に知りたい」場合に使うことを想定している． 「こういう頭を使わない機械的な分析はやってはいけない」という例を示すための実験でもある．
目的
データがn個の変数を含む場合に，
が応答変数で，それ以外の全ての変数が説明変数
が応答変数で，それ以外の全ての変数が説明変数
･･･
が応答変数で，それ以外の全ての変数が説明変数
という分析を自動化したい．
コードの作成
基本的なアイデアは次の2つのステップに区分できる．
一つは，テキストで式を自動生成してから，fomulaオブジェクトに変換してlmに突っ込むという方法である．
これは，Rの型キャストを利用した方法であり，オブジェクトの使い方がやや邪道な気もするが，型を意識した
コードの理解には適している．
xnam <- paste("x", 1:25, sep="") fmla <- as.formula(paste("y ~ ", paste(xnam, collapse= "+")))
もう一つは，応答変数を先頭列にもってきて，データフレームごと，lmに突っ込むというワイルド（？）なやり方である．
d2 <- data.frame(y,x1,x2,x3) lm(d2) Call: lm(formula = d2)
この場合，データフレームの先頭列をrotateするコードが別途必要になる．ちなみに上の二つのコードはこちらからの引用である．
上記のコードを元に，自動データ分析関数のプロトタイプを作ってみた．
analyzer <- function(dataframe){ n <- length(dataframe)#data.frame内の変数の個数を取得 no.col <- 1:n #ローテーション用ベクトル n1 <- n-1 #ループエンド用 for(i in 0:n1){ x <- no.col[1] # print(x) no.col[1] <- no.col[i+1] no.col[i+1] <- x #print(no.col) dataframe0 <- dataframe[,no.col] out = lm(formula = dataframe0) print(summary(out)) }}
計算例
なにかテキトーなデータ行列を使って試してみる．
dat0 = data.frame(swiss)#data行列の定義 head(dat0) #中身確認 analyzer(dat0)
アウトプットはこんな感じになる．
Call: lm(formula = dataframe0) Min 1Q Median 3Q Max -15.2743 -5.2617 0.5032 4.1198 15.3213 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 66.91518 10.70604 6.250 1.91e-07 *** Agriculture -0.17211 0.07030 -2.448 0.01873 * Examination -0.25801 0.25388 -1.016 0.31546 Education -0.87094 0.18303 -4.758 2.43e-05 *** Catholic 0.10412 0.03526 2.953 0.00519 ** Infant.Mortality 1.07705 0.38172 2.822 0.00734 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.165 on 41 degrees of freedom Multiple R-squared: 0.7067, Adjusted R-squared: 0.671 F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10 Call: lm(formula = dataframe0) Min 1Q Median 3Q Max -29.958 -9.107 1.287 10.805 24.504 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 135.1861 22.7578 5.940 5.26e-07 *** Fertility -0.7410 0.3027 -2.448 0.018727 * Examination -0.9627 0.5117 -1.881 0.067050 . Education -1.4908 0.4119 -3.619 0.000804 *** Catholic 0.1674 0.0762 2.197 0.033745 * Infant.Mortality -0.3609 0.8637 -0.418 0.678268 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 14.87 on 41 degrees of freedom Multiple R-squared: 0.618, Adjusted R-squared: 0.5715 F-statistic: 13.27 on 5 and 41 DF, p-value: 1.045e-07 Call: lm(formula = dataframe0) Residuals: Min 1Q Median 3Q Max -9.0429 -2.8491 -0.5308 3.0459 6.8874 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 24.57200 8.23969 2.982 0.00480 ** Fertility -0.09523 0.09371 -1.016 0.31546 Agriculture -0.08254 0.04387 -1.881 0.06705 . Education 0.33037 0.12858 2.569 0.01392 * Catholic -0.06778 0.02108 -3.215 0.00254 ** Infant.Mortality 0.09730 0.25297 0.385 0.70250 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.353 on 41 degrees of freedom Multiple R-squared: 0.7346, Adjusted R-squared: 0.7022 F-statistic: 22.7 on 5 and 41 DF, p-value: 7.624e-11
中略．
このように，関数を実行するとlmの結果がずらずら出てくる．
さて，以上の結果をさらにstep関数に渡してモデル選択まで自動でやろうとすると，エラーが出てしまう．
どうやら
formula = dataframe0
の部分が，step関数のオブジェクトとして不適切らしい．
しかたないので，一旦データフレームの列をrotateしてから，順番の入れ替わった変数名リストを取り出し，
テキスト型でy ~ x1 + x2 + ... + xnの形にしてからformulaにキャストする，という荒技(?)
を使ってみる．ついでにlmの実行の度に応答変数が変わるので，名前が分かるようにprintしておく．
#step関数を使い, AIC基準でモデル選択まで実行 analyzer101 <- function(dataframe){ n <- length(dataframe)#data.frame内の変数の個数を取得 no.col <- 1:n#ローテーション用ベクトル n1 <- n-1#ループエンド用 for ( i in 0:n1 ){ x <- no.col[1] # print(x) no.col[1] <- no.col[i+1] no.col[i+1] <- x #ベクトルのrotate #print(no.col) dataframe0 <- dataframe[,no.col]#データフレーム列のrotate variable = names(dataframe0) xnam = variable[2:n] #説明変数のラベル ynam = variable[1] #oucomeのラベル model = as.formula(paste(paste(ynam,"~") , paste(xnam, collapse= "+"))) out1 = lm(model,data=dataframe0) print(paste("Outcome is", names(dataframe0)[1])) out2 <- step(out1,trace=1)#AIC基準でモデル選択 print(out2) }}
結果を確認してみよう. データフレームはさっきと同じ．
> analyzer101(dat0) [1] "Outcome is Fertility" Start: AIC=190.69 Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC - Examination 1 53.03 2158.1 189.86 <none> 2105.0 190.69 - Agriculture 1 307.72 2412.8 195.10 - Infant.Mortality 1 408.75 2513.8 197.03 - Catholic 1 447.71 2552.8 197.75 - Education 1 1162.56 3267.6 209.36 Step: AIC=189.86 Fertility ~ Agriculture + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC <none> 2158.1 189.86 - Agriculture 1 264.18 2422.2 193.29 - Infant.Mortality 1 409.81 2567.9 196.03 - Catholic 1 956.57 3114.6 205.10 - Education 1 2249.97 4408.0 221.43 Call: lm(formula = Fertility ~ Agriculture + Education + Catholic + Infant.Mortality, data = dataframe0) Coefficients: (Intercept) Agriculture Education Catholic Infant.Mortality 62.1013 -0.1546 -0.9803 0.1247 1.0784 [1] "Outcome is Agriculture" Start: AIC=259.3 Agriculture ~ Fertility + Examination + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC - Infant.Mortality 1 38.59 9101.4 257.50 <none> 9062.9 259.30 - Examination 1 782.31 9845.2 261.20 - Catholic 1 1066.74 10129.6 262.53 - Fertility 1 1324.81 10387.7 263.72 - Education 1 2895.46 11958.3 270.33 Step: AIC=257.5 Agriculture ~ Fertility + Examination + Education + Catholic Df Sum of Sq RSS AIC <none> 9101.4 257.50 - Examination 1 816.46 9917.9 259.54 - Catholic 1 1076.03 10177.5 260.76 - Fertility 1 1947.58 11049.0 264.62 - Education 1 3091.51 12192.9 269.25 Call: lm(formula = Agriculture ~ Fertility + Examination + Education + Catholic, data = dataframe0) Coefficients: (Intercept) Fertility Examination Education Catholic 132.6210 -0.7989 -0.9802 -1.5193 0.1681
以下，略．
いまのところnumericなデータにしか使えないが，気が向いたらfactorが応答変数になった場合は自動的に 多項ロジットを実行する機能などを追加してみる．
繰り返すが，こんなAICの使い方は意味が無いし，「有意性を見つけてから理屈を考える」というアドホックな分析は，研究としてやってはいけない． 上記の関数からは，理論的に無意味なアウトプットもたくさん出てくるので，使用はリコードの確認や 探索的分析にとどめよう．
mathjax実験
はてぶろがマークダウンに対応しているのを知ったとき，へえそう，くらいにしか思ってなかった.
が...
いやいや，とっくにMathjax使えるようになっとるやん．
とか
とか．
まあ，Rmarkdownやそのものより使いづらいけど，少なくとも見た目はRpubsより綺麗やん．
これからは全ての資料はで書いて，WEBにアップするときだけマークダウンする方式でいこうかな・・・．
というか，そういう変換コード，誰か既に作ってるよな．．．
