研究会用メモ
案の定，2週間ほど前にやった計算をすっかり忘れて，同じコードを書くところであった． ブログにメモを残しておいたおかげで，思い出した．
というわけで
という実験の続きを再開する（10/11のエントリの続きです）． 参照したコードは，例によって Kruschke氏のコードである．
準備はこれだけ
source("DBDA2E-utilities.R") fileNameRoot="Jags-lognormal-direct-sim" graphFileType="png"
まずR内で人工データを発生させる
# "t" means "true value" y0t=10; bt=0.1; pt=0.3; nt=25; mu=log(y0t)+nt*log(1-bt)+log((1+bt)/(1-bt))*nt*pt sigma= sqrt(nt*pt*(1-pt)*(log((1+bt)/(1-bt)))^2) gmc = rlnorm( n=500, mu, sigma)
pt, bt, nt が推定すべき真のパラメータである．
mu, sigmaがこのような決定論的関数で表される理由はgmcで証明が与えられる（ここ重要）
次にJAGSにわたすためのデータを定義する
dataList = list( y = gmc , N = length(gmc), nt=nt, y0t=y0t)#JAGSに渡すパラメータはlist内で再定義する # nt, y0tはR内で定義されているがJAGSにはdataを介して # 渡している．
nt,y0tはスカラである．dataListはRオブジェクトだが，modelstringで参照するために，このようなリストを作っておく
JAGSにわたすモデルを複数定義する．今回は
modelstring0
modelstring1
modelstring2
の三つ．
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(mean.n)#nをポアソン分布で近似 mean.n ~ dunif(1,2*nt) y0 ~ dpois(mean.y0)# mean.y0 ~ dunif(1,2*y0t) # y0t=10; bt=0.2; pt=0.65; nt=25; p ~ dbeta(1,1)#p,bの事前分布にβ分布 b ~ dbeta(1,1)}" # close quote for modelstring
modelstring0は推定するパラメタが一番多い．
つぎは，より理論モデルに忠実な設定．
#modelstring1は理論的に固定すべきパラメータn,y0を定数にする modelstring1 = " model{ for( i in 1 : N ) { y[i] ~ dlnorm(m, 1/(s^2)) } y0 <- y0t n <- nt 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) #prior for hyper parameters p ~ dbeta(1,1) b ~ dbeta(1,1) }" # close quote for modelstring
modelstring1はn,y0を定数として固定する．
ここで，先ほど定義したリスト内のスカラを割り当てる
（最初はRのグローバル変数を参照していたのだが,modelstringの中からは参照できないことが分かった．
まあ，よく考えればこっちはJAGSが読むコードなので，当然である ）
最後のモデル．
# dilectly assign mean of poisson dist. of n,y0 modelstring2 = " 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(nt)#nをポアソン分布で近似 y0 ~ dpois(y0t)# p ~ dbeta(1,1)#p,bの事前分布にβ分布 b ~ dbeta(1,1)}"
最後に，ハイパーパラメータだけ，真の値を代入してその他を推定するモデル
モデルの選択はこんな風にswitchで定義する．（文字型しかダメなのだろうか？）
a<-"1" modelstring<-switch(a, "0"=modelstring0, "1"=modelstring1, "2"=modelstring2)
モデルによって推定するパラメータは直に指定する
parameters0 = c("p","b","n","mean.n","mean.y0","y0","m","s") parameters1 = c("p","b","m","s") parameters2 = c("p","b","n","y0","m","s") parameters<-switch(a, "0"=parameters0, "1"=parameters1, "2"=parameters2)
あとは細かい話だが，plotPostに渡すパラメタがモデルによって異なるので， if文でモデル毎に指定しておく．
分析の度に，モデル選択変数aを指定することに注意する．
#全model共通で推定するパラメタ postInfo = plotPost( mcmcChain[,"p"] , xlab="success p" ) postInfo = plotPost( mcmcChain[,"b"] , xlab="interest b" ) #以下モデルに合わせてplotするパラメタを選択 if(a=="0") postInfo = plotPost( mcmcChain[,"n"] , xlab="time n" ) if(a=="0") postInfo = plotPost( mcmcChain[,"mean.n"] , xlab="mean of n" ) if(a=="0") postInfo = plotPost( mcmcChain[,"mean.y0"] , xlab="mean of y0" ) if(a=="1") postInfo = plotPost( mcmcChain[,"m"] , xlab="mu of y") if(a=="1") postInfo = plotPost( mcmcChain[,"s"] , xlab="sigma of y") if(a=="2") postInfo = plotPost( mcmcChain[,"n"] , xlab="time n" ) if(a=="2") postInfo = plotPost( mcmcChain[,"y0"] , xlab="y0" ) if(a=="2") postInfo = plotPost( mcmcChain[,"m"] , xlab="mu of y") if(a=="2") postInfo = plotPost( mcmcChain[,"s"] , xlab="sigma of y") saveGraph(file=fileNameRoot,type=graphFileType)
モデル1の推定結果
やはり理論的に固定すべき変数を真値で与えたデータの場合，未知のp,bの推定が
正確である．
問題はnt,y0tが現実のデータから推測できないことである．さてどうしよう...
■
Refreshに関するメモ
Mathematica内蔵関数Refreshの引数である，UpdateIntervalは更新を制御できたりできなかったりして，いまいち仕組みが分からない．そこでいろいろと実験してみた．
Manipulate[ Refresh[Random[], UpdateInterval -> If[go, 0.105, Infinity]], {{go, False}, {True, False}}]
このようなコードなら，UpdateIntervalによって更新頻度が制御できる．一方で
x = 1; Manipulate[ Refresh[x = x + 0.01; x, UpdateInterval -> If[go, 0.02, Infinity]], {{go, False}, {True, False}}]
これはUpdateIntervalがInfinityにもかかわらず暴走してしまう
そこでRefreshの中にIfを挟んで代入を制御する.
x = 1; Manipulate[ Refresh[ If[go, x = x + 0.01], UpdateInterval -> If[go, 1, Infinity]], {{go, False}, {True, False}}]
これならオンオフだけ制御できる．しかしインターバルの時間は指定通りにならない．
TrackedSymbols -> {} で制御しようとしてもうまくいかない．
結局，即時割り当てをManipulateの中に入れると，Dynamicに更新されるので，その場合はUpdateIntervalで制御できる範疇を超えてしまうらしい．
オブジェクトを再帰的代入すると，refreshしてもうまくいかないらしい
次のような命令なら代入でも問題ないし，UpdateIntervalによる制御もうまくいく
Manipulate[ Refresh[ If[go, x = Random[]], UpdateInterval -> If[go, 0.8, Infinity]], {{go, False}, {True, False}}, Button["reset", go = False; x = 100] ]
しかし再帰的に代入すると，Dynamicオブジェクトが自分自身を連続的に更新するので UpdateIntervalを指定しても，機能しない （のではないかと推測する）
Manipulate[ Refresh[ If[go, x = x + Random[]], UpdateInterval -> If[go, 0.8, Infinity]], {{go, False}, {True, False}}, Button["reset", go = False; x = 100] ]
グラフィックオブジェクトの動的更新
例えばつぎのようなDynamicオブジェクトはRefreshによって動的に更新される
（ただし止まらない）
Dynamic[Refresh[Histogram[RandomReal[10, 100]], UpdateInterval -> 0.1]]
制御したければManipulate内に入れて，チェックボックスでUpdateIntervalの値を指定する方法がある．
Manipulate[ Refresh[Histogram[RandomReal[10, 1000]], UpdateInterval -> If[x, 0, Infinity]], {{x, False}, {True, False}}]
今まで，ピコピコ系アウトプットはこの方法で作ってきた．
以下のようなアホなコード(めっちゃ計算しているようにみえるが,実は無意味なコード)も簡単に書ける．
Manipulate[Refresh[ Histogram[ {RandomVariate[NormalDistribution[0, 1], n], RandomVariate[NormalDistribution[a, 1/2], n], RandomVariate[NormalDistribution[10, 1], n], RandomVariate[LogNormalDistribution[1, 0.3], n], RandomVariate[NormalDistribution[8, 0.5], n]}, {.25}, PlotRange -> {{-3, 15}, {0, 100}}], UpdateInterval -> If[x, 0, Infinity]], {{x, False, "Run MCMC"}, {True, False}}, {{n, 300}, 1, 500, 1}, {{a, 3, "\[Mu]"}, 1, 10}]
いろいろ実験して分かったことは，ListPlotみたいな即時実行命令をRefreshするとうまくいかない，ということだ．
冗長なユーザ定義関数を割り当て（実際にはクロージャを含まないダミーでよい），遅延割り当てにすると，ちゃんとRefreshが機能する．
だからvisualizeみたいなユーザ関数を使うとうまくいくわけだ．
以上をふまえて作ったコード
(* 関数の定義 *) initialdata[n_, m1_, m2_] := {RandomVariate[NormalDistribution[m1, 5], n], RandomVariate[NormalDistribution[m2, 5], n]}; fightonce[group_] := Module[{power1, power2}, power1 = group[[1]]; power2 = group[[2]]; If[Length[power1] > 0 && Length[power2] > 0, power1 = RandomSample[power1]; power2 = RandomSample[power2]; If[(* 戦闘力を比較して勝者を残す*) RandomChoice[ {power1[[1]]/(power1[[1]] + power2[[1]]), power2[[1]]/(power1[[1]] + power2[[1]])} -> {True, False}], (* True -> 1が勝つ*) power2 = Drop[power2, 1], power1 = Drop[power1, 1]];(* If end*) ];(* If end*) {power1, power2} ]; (* update関数のほうに停止条件を書いておく. dataがなくなったらNullを返す*) visualize[data_, hight_] := Histogram[data, {1}, PlotRange -> {{0, 100}, {0, hight}}]; (* histogram を遅延評価させるために,関数化する*) (*n=50;ylim=25;sample=initialdata[n];*) Manipulate[Refresh[ If[x, sample = fightonce[sample]]; (* If文で実行する動作：UpdateInterval\[Rule]0 のときだけ,update関数を再帰的代入する*) visualize[sample, ylim], UpdateInterval -> If[x, 0, Infinity]], {{x, False}, {True, False}}, {{n, 100}, 50, 1000, 1}, {{m1, 40}, 10, 70}, {{m2, 50}, 10, 70}, Button["reset", x = False; sample = initialdata[n, m1, m2]; ylim = n/5]]
少しだけRefreshの理解が深まった（Dynamicは，まだ謎が多い）．
Dynamic Models of Segregation
Dynamic Models of Segregation (Schelling 1971)をリアルタイムに表現するコードを更新したのでそのメモを残しておきます．
今回はNetLogo風に，エージェントの幸福度をリアルタイムに計算してプロットする部分を追加してみました．
Mathematica ノートブックは以下のサイト内で公開しているseg.nbというファイルをDLしてください．
使い方は，ノートブック内にも書いていますが，
- Module セルを読み込む
- 評価 [RightArrow] 動的更新の実行にチェックを入れる
- Manipulate セルを読み込む
- resetボタンを押す
- run simulationのチェックボックスをon
- 終了すると END!! が出て自動的に停止する
- 条件を変えて4.から繰り返し
です．最終的なアウトプットはこんな感じです．
- 赤と青の個体がトーラス上をランダムに移動します
- 各個体は自分の近傍内に同色の仲間が35%以上存在すると，その位置に定着します
- 左のバーは平均的な幸福度（周りの同色エージェント割合）です
- このシミュレーションは，強い排外意識がなくても，分居が自然に生じる様子を表しています
ノートブックを読み込む段階で動的更新をonにしていると，このUIが表示されるはずですが，エラーが出る場合もあるので 最初はオフにしておき，上記の手順で動的更新を許可すれば問題ないはずです．
このnbはWolfram Demonstrations Projectで公開されているコードを参考にして作りました．
特に，GarbageCollectionByAntsというコードを参考にさせていただきました．
使用関数の定義
まずエージェントの初期配置を決める関数．size×sizeの行列上の番地でエージェントの位置を指定する． 赤エージェントの数は，rednで，赤エージェントの数は，bluenで指定する．
(*赤agent 青agentの初期配置*) initialset[size_, redn_, bluen_] := Module[{a, b, n = redn + bluen}, Developer`ToPackedArray[(* Integer型指定で高速化 *) a = RandomSample[Range[size*size], n];(* 総エージェント数だけ1次元配列から非復元抽出 *) b = Flatten[Table[{i, j}, {i, 1, size}, {j, 1, size}], 1];(* 2次元座標 *) (* 抽出した1次元配列番号を座標に変換する． 前半は赤，後半は青用 *) {Table[b[[a[[i]]]], {i, 1, redn}], Table[b[[a[[i]]]], {i, redn + 1, n}]}, Integer] ]; (* 非復元抽出で初期配置を与える． output: {{赤の位置情報}，{青の位置情報}} の順で出力する initialset[size,redn,bluen][[1]] が 赤初期配置 initialset[size,redn,bluen][[2]] が 青初期配置 *)
初期グリッドに0を並べる手続きです．動的表現を使わない場合は，わざわざ定義する必要はありません．
(*初期グリッド*) zerogrid[size_] := Table[0, {size}, {size}]
エージェントの位置情報を色に変換する関数です．この関数の引数である位置情報がDynamicで更新されるのでエージェントが動いているように見える，という仕組みです．
(* agentの位置情報を検索して，可視化する関数 *) visualize[{redposi_, blueposi_}, size_] := Module[{zero1 = zerogrid[size](*Mnipulate内で指定*)}, (* Embedding agents into zero1 for visualization 引数のredposi_,blueposi_はそれぞれ「赤」と「青」の現在位置 *) Do[zero1[[ Sequence @@ redposi[[i, {1, 2}]]]] = 1, {i, Length[redposi]}]; Do[zero1[[ Sequence @@ blueposi[[i, {1, 2}]]]] = 2, {i, Length[blueposi]}]; ArrayPlot[zero1, ColorFunction -> (If[# == 0, White, If[# == 1, Darker[Red], Darker[Blue]]] &), ColorFunctionScaling -> False, ImageSize -> {300, 300}, AspectRatio -> Automatic, Frame -> False, Mesh -> All, Axes -> False] ];
全セル（トーラス空間）のムーア近傍の座標を出力する関数です．エージェントの位置毎に計算するよりも，グローバルに一回計算して，必要があるやつだけ呼び出したほうが早いんじゃないかと思って定義した関数です．計算は一回きりなので，メモリの節約にはなっていると思うのですが，実際に比較はしてません．ただ，コードの構造上はこの方法が理解しやすいと思います．
(* トーラス上に配置されたエゴのムーア近傍を出力する関数 *) moore[t_(* agent position *), size_] := Module[{m = size, p1, p2, p3, p6, p4, p7, p8, p9, network}, p6 = Table[ If[t[[i]][[2]] != m, t[[i]] + {0, 1}, t[[i]] + {0, -(m - 1)}], {i, Length[t]}]; p9 = Table[Which[(* +-している1を一般化すれば近傍の距離を自由に定義できる*) (t[[i]][[1]] != 1) && (t[[i]][[2]] != m), t[[i]] + {-1, 1}, (t[[i]][[1]] == 1) && (t[[i]][[2]] == m), t[[i]] + {m - 1, -(m - 1)}, (t[[i]][[1]] == 1) && (t[[i]][[2]] != m), t[[i]] + {m - 1, 1}, (t[[i]][[1]] != 1) && (t[[i]][[2]] == m), t[[i]] + {-1, -(m - 1)}], {i, Length[t]}]; p7 = Table[Which[ (t[[i]][[1]] != 1) && (t[[i]][[2]] != 1), t[[i]] + {-1, -1}, (t[[i]][[1]] == 1) && (t[[i]][[2]] == 1), t[[i]] + {m - 1, (m - 1)}, (t[[i]][[1]] == 1) && (t[[i]][[2]] != 1), t[[i]] + {m - 1, -1}, (t[[i]][[1]] != 1) && (t[[i]][[2]] == 1), t[[i]] + {-1, (m - 1)}], {i, Length[t]}]; p3 = Table[Which[ (t[[i]][[1]] != m) && (t[[i]][[2]] != m), t[[i]] + {1, 1}, (t[[i]][[1]] == m) && (t[[i]][[2]] == m), t[[i]] + {-(m - 1), -(m - 1)}, (t[[i]][[1]] == m) && (t[[i]][[2]] != m), t[[i]] + {-(m - 1), 1}, (t[[i]][[1]] != m) && (t[[i]][[2]] == m), t[[i]] + {1, -(m - 1)}], {i, Length[t]}]; p1 = Table[Which[ (t[[i]][[1]] != m) && (t[[i]][[2]] != 1), t[[i]] + {1, -1}, (t[[i]][[1]] == m) && (t[[i]][[2]] == 1), t[[i]] + {-(m - 1), (m - 1)}, (t[[i]][[1]] == m) && (t[[i]][[2]] != 1), t[[i]] + {-(m - 1), -1}, (t[[i]][[1]] != m) && (t[[i]][[2]] == 1), t[[i]] + {1, (m - 1)}], {i, Length[t]}]; p4 = Table[ If[t[[i]][[2]] != 1, t[[i]] + {0, -1}, t[[i]] + {0, (m - 1)}], {i, Length[t]}]; p2 = Table[ If[t[[i]][[1]] != m, t[[i]] + {1, 0}, t[[i]] + {-(m - 1), 0}], {i, Length[t]}]; p8 = Table[ If[t[[i]][[1]] != 1, t[[i]] + {-1, 0}, t[[i]] + {(m - 1), 0}], {i, Length[t]}]; network = Developer`ToPackedArray[ Table[{p6[[i]], p4[[i]], p2[[i]], p8[[i]], p9[[i]], p7[[i]], p3[[i]], p1[[i]]}, {i, Length[t]}], Integer] ];
エージェントの幸福度を計算する関数です．最小値0，最大値1です．位置情報と幸福度をセットにして出力します．
(* ムーア近傍情報に基づき幸福度（ムーア近傍内の同色エージェント割合）を計算, 位置情報に幸福度を追加 引数 agentposiは単色エージェント *) happy[agentposi_, size_] := Module[{neloca, happy3}, neloca = moore[agentposi, size]; happy3 = Table[Plus @@ Table[If[MemberQ[agentposi, neloca[[i]][[j]]], 1, 0], {j, 8}], {i, 1, Length[neloca]}]; happy3 = happy3/8;(* 同色エージェント数を *) Table[Append[agentposi[[i]], happy3[[i]]], {i, Length[agentposi]}] ]; (*入力は {{1,1},{2,1},{3,1}} のような単色エージェントの位置情報．出力は {{1,1,0.2},{2,1,0.3},{3,1,0.8}} のような 位置情報+幸福度 *)
最後にエージェントの移動を管理する関数です．設定した同色エージェント割合（閾値）に達していない場合に移動します．
(* 幸福情報に基づき，移動する visualizeの引数を生成する *) move[agents_, size_, redth_, blueth_, movern_] := Module[{preloca, zerolocate, red0, blue0, red1, blue1, k1, k2, AllocationforRed, l1, l2, AllocationforBlue}, red0 = agents[[1]]; blue0 = agents[[2]]; red1 = happy[red0, size]; blue1 = happy[blue0, size]; preloca = Flatten[{red0, blue0}, 1];(* 埋まっているセル *) zerolocate = Complement[Tuples[Range[size], 2], preloca];(*ゼロセルの番地リスト*) (*red agentからmovern人を選び動かす*非復元抽出* *) k1 = RandomSample[Range[Length[red0]], movern]; k2 = Select[k1, red1[[#, 3]] < redth &]; AllocationforRed = RandomSample[zerolocate, Length[k2]]; red[[k2]] = AllocationforRed; (*Blue agentからmovern人を選び動かす*非復元抽出* *) l1 = RandomSample[Range[Length[blue0]], movern]; l2 = Select[l1, blue1[[#, 3]] < blueth &]; AllocationforBlue = RandomSample[Complement[zerolocate, AllocationforRed], Length[l2]]; blue[[l2]] = AllocationforBlue; {red, blue} ];
さて，以上の関数をforで回すのではなく，Manipulateの中でRefreshします．
するとリアルタイムでエージェントがピコピコ動きます．このような動的表現が，MathematicaのDynamic関数によって実行できます．
以下Manipulateの中身です．
(* ver110:幸福度をlistplotでリアルタイムに更新 自動停止機能を追加，happy level モニタリングをDynamic追加 *) Manipulate[ Refresh[ If[moving, {red, blue} = move[{red, blue}, s, redth, blueth, nomove]](* 更新している部分 *); If[(Min[happy[red, s][[All, 3]]] >= redth) && (Min[happy[blue, s][[All, 3]]] >= blueth), moving = False; state = Text[Style["END !", Red, 20]], moving = moving]; state2 = Mean[happy[red, s][[All, 3]]] // N; state3 = Mean[happy[blue, s][[All, 3]]] // N; (* agentのhappy を 監視して，全員が満足したらmoving をFalseに切り替え, Dynamic[state]表示をendに変える． moving の値がFalseになると {red,blue}の更新が終わり， "run simulation"チェックボックスがoffになる*) visualize[{red, blue}, s], UpdateInterval -> If[moving, 0.01, Infinity]](* reflesh end *), {{s, 20, "length of a side"}, 10, 30, 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{redn, Floor[2*s], "number of red"}, 1, Floor[(s*s)/2], 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{bluen, Floor[2*s], "number of blue"}, 1, Floor[(s*s)/2], 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{redth, 0.35, "red threshold"}, 0, 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{blueth, 0.35, "blue threshold"}, 0, 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{nomove, Floor[(s)/5], "number of mover"}, 1, Min[{redn, bluen}], 1, ImageSize -> Tiny, Appearance -> "Labeled"}, {{moving, False, "run simulation"}, {True, False}}, Button["reset", moving = False; state = "going..."; x = initialset[s, redn, bluen]; red = x[[1]]; blue = x[[2]], ImageSize -> Medium], Dynamic[{"state=", state}], Dynamic[ListPlot[{{{1, state2}}, {{2, state3}}}, Filling -> Axis, PlotRange -> {{0, 3}, {0, 1}}, PlotMarkers -> {Red, Blue}]], SaveDefinitions -> True, TrackedSymbols :> {moving, red, blue, redn, bluen, s}, SynchronousUpdating -> True, ControlPlacement -> Left ]
なにゆえ，このManipulateでエージェントが動くのかを，ざっと解説します．
Manipulateの直接の実行対象となっている関数は Refresh[ ]です．Refresh[ ]内で，上で定義した一連の関数をループさせています．
ループさせている関数の中には，visualizeがあるので，Refreshされる度に，グラフィック出力が変化します．このRefreshするタイミングを
UpdateInterval ->で指定します．intervalが短い（上記コードでは0.01）ので，CPUさえ高性能なら，ヌルヌル動くというわけです*1．
実験としてスライダの下部にDynamic objectを2つ埋め込みました．一つはシミュレーションが停止したかどうかを示すstateです．
停止条件を満たすとここがEND！と表示されます．
もう一つが，今回修正したエージェントの幸福度（リアルタイム）です．色別の幸福度の平均値を計算してListPlotする命令をDynamicに放り込むだけでうまくいきました．
Dynamicは，誤って使用するとPCに負担がかかりすぎてカーネルが落ちるため，なかなか使いこなすのが難しい関数ですが，うまくいけば動的な表現の幅が広がるので是非
使ってみてください．
（今みなおすと，まだ無駄な計算が残っており，もっと効率化できそうです．どのあたりが修正できそうか，気づいた学生さんは授業時間内に教えてください）
*1:厳密に言えば，このオブジェクトは最速で更新されるために，インターバルの指定は無視されます．遅延評価する実行文の場合，インターバルは指定通りになります．これはRefreshの仕様です
SSD換装メモ
ノートPC SSD換装メモ
PCがモッサリしてきたけど，買い換えると再設定が面倒．そんなときOSごとクローンしてSSD換装するわけですが， 数年に1回しかやらないので，手順をすぐに忘れてしまう．そこで，次に作業するとき参照するために記録を残しておきます．
換装手順
- BIOSでCDドライブを起動ドライブに設定
- 移動先SSDをSATA/USB変換ケーブルで接続
- acronis true imageをCDドライブに入れて起動．OSごと内蔵SSDの中身を移動先SSDにクローン
- 電源オフ
- 内蔵SSDを引っ張り出して，移動先SSDをSATA接続．このとき移動先SSDに，取り出し用のテープを貼る．
- BIOSの起動設定を戻す→起動
- ウマー(ﾟДﾟ)
使った道具
- 移動先のSSD
- SATA/USB変換ケーブル
- acronis true image(SSD付属品)
- プラスドライバ
気づいたポイント
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に従うという仮定は妥当か？
といったところ．
続く
