The One with ...

思いついたことや作業メモなどを公開しています

StanからWAICを計算するためのメモ

StanからWAICを計算するためのメモ

豊田本（2017）を参考に，Stanの結果からWAICを計算する方法のメモ．

テキスト198ページ記載のコードは略記のため，そのまま書き写して実行してもStanが動きません*1． 公開されたコードを参考に確率分布のパラメータを適切なブロックで定義すれば，Stanが動きます．

まずStanコード内でモデルの対数尤度を生成量（generated quantities）として定義します．

次にR上でlooパッケージを使って，Stanのサンプリング結果から対数尤度を取り出し，waic関数にわたします*2．

それぞれのコードはこんな感じです（ポイントだけ省略して書いています．実際にはデータブロックやパラメータブロックを書かなくてはいけません）．

上記のGitHub Gistからの埋め込みは，Gist上で「Embed」をクリップボードにコピーしたものを ブログのソースに張り付けるだけで完了します．

*1:モデルブロック内で宣言したパラメータが存在しないというエラーメッセージがでます．JAGSならこれでも行ける気がしますが

*2:内部的に対数尤度の総和はtargetに自動記録されるけど，平均計算用に対数尤度をベクトルで記録する必要がある，という理解でいいのだろうか？

TeXstudioショートカットキーのカスタマイズ方法

オリジナルのショートカットキーを設定する方法

TeXstudio便利ですね．私は長い間 秀丸+マクロ を使っていましたが，TeX Live導入後は TeXwokrsとTeXstudioを併用しています．

さて，TeXstudioでオリジナルのショートカットキーを作成したいと思うことが時々あります．

時間がたつと，すぐにやり方を忘れてしまうので，未来の自分のためにメモを残しておきます．

例．ショートカットキーを使って，ワンストロークで

 \begin{align*}　
 &   \\ 
 &　
 \end{align*}

といった改行込みの数式用テンプレ（環境）を入力したい，という状況を考えます．

方法の概略

「メニュー」にまず登録して，その後ショートカットキーを設定する．

手順

  • メニューから「オプション」→「TeXstudioの設定」を選ぶ．

  • 「高度なオプション」にチェックを入れる（チェックしないと「メニュー項目」が出てこない）

  • 左の列から「メニュー」を選ぶ．次に「数式」の項目から適当な一つを選び，右クリックする．

f:id:hamada7418:20170613155650p:plain

  • 右クリックメニューから「（前に）新規メニューアイテムを挿入」を選択する．新しい項目が挿入されたら，適当な名前を入力する（今回の例では「複数行」という名前にしたと仮定します）．

  • 「コマンド」の列に，挿入したいコマンドを書く．例えば

\begin{align}%\ & \ %\ & %\ \end{align}

のように書くと，実行したときにソース上で

 \begin{align*}  　
 &   \\   
 &     
 \end{align*}

が挿入される．メニュー設定のコマンドでは「　%\　」が改行を意味しています． 「　%|　」で挿入後のカーソル位置を指定できます（必要がなければ指定しなくてもOKです）．

  • メニュー登録が完了したら，左の列で「キーボードショートカット」を選択．「数式」の項目を展開すると，今登録した「複数行」があるので「現在のショートカット」のセルをダブルクリックする．

f:id:hamada7418:20170613155654p:plain

  • そこで自分の好きなキーを登録（例として「Ctrl + E 」を登録したと仮定します）．登録後，OKボタンを押して設定を終了する．

  • ソース上で「Ctrl + E 」を押す．

ソースに

 \begin{align*}  　
 &   \\   
 &     
 \end{align*}

が挿入されれば成功です．

うまくいったら，同様のやり方でよく使うコマンドを登録して，入力の手間を省きましょう（もともと登録されているものは，ショートカットキーを編集するだけでOKです．なおメニューに登録していないコマンドをショートカットキー登録する方法は知りません）．

tikz-decoration

TikZのdecorationがうまくいかないので検索してみると，次のような解決策がweb上にあった．

マニュアルに記載の

\draw[decorate,decoration={coil,segment length=4pt}] (0,1) -- (3,1);
\draw[decorate,decoration={coil,aspect=0}] (0,.5) -- (3,.5);

を実行するには，

\usetikzlibrary{decorations}

ではダメで，

\usetikzlibrary{decorations.pathmorphing}

まで指定しないといけないらしい．

マニュアル原典が大部すぎるため，これは確かにわかりにくい．

JAGSでWAICを計算するためのメモ

試しに計算してみたのですが，まだ細部を理解できていないので，検証用にメモを残しておきます．

豊田本2017（『実践ベイズモデリング』）にStanからWAICを計算する手順が解説されていたので，それを参考にして JAGSから作ったサンプルでWAICを計算してみます．

例によってJAGSをRから呼び出すのにパピー4匹本のDBDA2E-utilities.Rを使います．

まずJAGS用のmodel{ }内に

 # Generated quantities
    for (i in 1:N) {
        loglikelihood[i] <- log(dnorm(y[i], mu[i], 1/ (s * s)));
    }

などを追加して，推定したパラメータとデータを用いた対数尤度を計算します． この場合はパラメータはmuとsです．muは説明変数に回帰しているので添え字[i]がついています． 確率密度関数は応答関数にあわせて適宜修正しましょう．

coda.samplesの結果を格納します．

mcmcCoda = coda.samples( jagsModel , variable.names=parameters ,
n.iter=nPerChain , thin=thinSteps )

ここから対数尤度を整形して取り出し，looパッケージのwaic関数に渡します．

こちらの記事を参考にさせていただきました．

ito-hi.blog.so-net.ne.jp

library(loo)
loglik1 <- sapply(1:N, function(i)
    unlist(mcmcCoda[, paste("loglikelihood[", i, "]", sep = "")]))
waic(loglik1)

これでStanから計算した結果とほぼ同じWAICが出力されました． looパッケージのextract_log_lik関数を使えばもっと簡潔になりそう．

TikZ

TikZで簡単かつ綺麗な図がかけるので，これがあれば研究の分析・報告・文章化は全てTeX+beamerだけで事足りるでしょう．

使わないとすぐに忘れそうなので，ポイントだけメモしておきます．

簡易マニュアルでも推奨しているように，最初にグリッドを書いておくと作業が大変ラクです．

これは

\draw[step=1,black!10, very thin] (0,0) grid (12, 8 );

などと入力するだけでOKです．デフォルトの座標単位が1cmなので紙の大きさに合わせて，調整すればよいでしょう．

あとはnodeで数式やプレースフォルダを配置して，drawで矢印を足していく感じ．

nodeの座標は始点位置を指定できるので，例えば

\node [right] at (0.5,7.5 )

のように書いて，調整します．基本的に[right]でいいでしょう．ただしオブジェクトを丸で囲む場合は，サークルの中心点とラベル（ノード）の中心点を合わせましょう．

矢印は

 \draw[thick,-latex] (7.5 ,3 ) -- (9.5 ,2.2 );
    \draw[thick,-latex] (7.5 ,1 ) -- (9.5 ,1.8 );

こんな感じで始点と終点を指定します．

ベイズモデリングでよく使うプレート表現も簡単に描けます．

狙った位置にオブジェクトをビシッと配置できるので気持ちがいいです．

f:id:hamada7418:20170304220554p:plain

まあ，パワポとどっちがラクやねんと言われれば，パワポのほうが手っ取り早い気もしますが，数式の美しさと再現可能性を考慮すればbeamer + TikZがオススメです．

Texworks　ショートカットキーのカスタマイズ

Texworksのショートカットキーの設定を編集すると，任意のTeXコマンド入力を登録できます．

例えばTeXソース編集中に

図 + tab

と入力するだけで，

\begin{figure}[h]
 \begin{center} 
 \includegraphics[scale=0.5]{ .pdf}
 \caption{•} 
 \label{american} 
\end{center}
\end{figure}

というコマンドが挿入できます．

つまり自分が選んだ好きな文字の後にtabキーを押すだけで， 任意のコマンドを挿入できるように設定できます．

これは簡単かつ便利な機能なので是非使ってください．

よく使うコマンドの大半はデフォルトの設定で用意されています．

たとえばギリシア文字の場合

xa+tab

で

\alpha

に変換されます．

設定の方法

まず

\texmf-config\texworks\completion

フォルダ内の

tw-latex.txt

を開きます．

例えば bite+tab という入力に対して

\begin{itemize}
 \item 
 \item 
\end{itemize}

というコマンドの入力を対応させるには， tw-latex.txtの中に

bite:=\begin{itemize}#RET# \item #INS# #RET# \item #RET# #RET#\end{itemize}

という一文を追加するだけでOKです．

#RET#は改行で，#INS#は入力後のカーソル位置です．•を入れておけば，tabでその場所に飛びます．環境内で順番に入力するスペースが決まっている場合はあらかじめ •をいれておくと便利です．

私がよく使っているキーはこんな感じです

v:=\verb|#INS#|
ｖ:=\verb|#INS#|
(:=\left( #INS# \right)
{:=\left{ #INS# \right}
（:=\left( #INS# \right)
｛:=\left{ #INS# \right}
図:=\begin{figure}[h]#RET# \begin{center} #RET# \includegraphics[scale=0.5]{#INS#.pdf}#RET# \caption{•} #RET# \label{american} #RET# \end{center}#RET#\end{figure}
表:=\begin{table}[h]#RET# \caption{}\label{}#RET# \begin{center} #RET# \begin{tabular}{ccc}#RET# #INS# &  & \\ #RET#  &  & \\ #RET#  &  & \\ #RET# \end{tabular} #RET# \end{center} #RET# \end{table}
定義:=\begin{definition}[] #INS# #RET#\end{definition}
命題:=\begin{prop} #INS# #RET#\end{prop}
定理:=\begin{} #INS# #RET#\end{}
系:=\begin{} #INS# #RET#\end{}
行列:=\begin{pmatrix}#RET# &  \\#RET# & \\ #RET#\end{pmatrix}
bite:=\begin{itemize}#RET# \item #INS# #RET# \item #RET# #RET#\end{itemize}
ni:=\begin{nitem}#RET# \item #INS# #RET# \item #RET# #RET#\end{nitem}
e:=\begin{align*}#RET# #INS# & = \\ #RET#  & = #RET#\end{align*}
f:=\frac{#INS#}{•}
%nitem:=\begin{nitem}#RET# \item #INS# #RET# \item #RET#\end{nitem}
r:=\sqrt{#INS#}
s:=\sum_{#INS#}^{•}
4:=$$#INS#
４:=$$#INS#

全角文字でも使えるのがポイントです．日本語でTeXを入力している時，入力モードを切り替えるのが面倒ですよね．そんなとき このショートカットを全角文字に割り当てておけば，いちいち切り替えなくも入力できます．

Income Segregation

www.jstor.org

この論文を読んでいて，Income Segregation と Race Segregation　を組み合わせて表現するコードを思いついたので，メモを残しておきます．

以前から，シェリングのsegregation modelにエージェントの所得情報を別のレイヤーで表示させる，というアイデアを思いついていました． ただ，面倒な割にインプリケーションも少ないかなと思って，今日まで実装するに至りませんでした．

しかし，Reardonのような実証論文があるということは，理論的な挙動を把握するためのプロトタイプを作ってみるのも悪くないかと思い，今回試しにやってみた次第です．

例によってコードはMathematicaでちゃちゃっと書きます*1.

所得情報のArrayPlot

まずエージェントの所得情報をArrayPlotするModuleを考えます．最初にエージェント毎に所得ベクトルを生成します．

(* income　の初期情報を与える *)
income[redn_, bluen_, mred_, mblue_] := Module[{sigma = 1},
   {RandomVariate[LogNormalDistribution[mred, sigma], redn],
    RandomVariate[LogNormalDistribution[mblue, sigma], 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: {{赤の位置情報}，{青の位置情報}} の順で出力する*)

これは以前に解説した関数と同じです．

次に所得情報を位置情報に対応させてArrayPlotするための関数を定義します．

(* 所得情報をplotする関数 *)
visualize2[{redposi_, blueposi_}, incomelist_, size_] :=
  Module[{zero1 = zerogrid[size], redincome, blueincome(* 
    sizeは　Manipulate内で指定*)},
   (* Embedding agents into zero1 for visualization 
   引数のredposi_,blueposi_はそれぞれ「赤」と「青」の現在位置　*)
   redincome = incomelist[[1]];
   blueincome = incomelist[[2]];
   Do[zero1[[ Sequence @@ redposi[[i, {1, 2}]]]] = redincome[[i]], {i,
      Length[redposi]}];
   Do[zero1[[ Sequence @@ blueposi[[i, {1, 2}]]]] = 
     blueincome[[i]], {i, Length[blueposi]}];
   ArrayPlot[zero1, ImageSize -> {300, 300}, AspectRatio -> Automatic,
     Frame -> False, Mesh -> All, Axes -> False]
   ];

ここではゼロのならんだグリッドに，エージェントの位置に応じて所得値を割り当てます．例えばもとの3次正方行列のゼログリッドが

0,0,0,

0,0,0,

0,0,0

だった場合に，3人分の所得{5,3,2}を位置情報{{1,3},{2,1},{3,3}}に割り当てると

0,0,5,

3,0,0,

0,0,2

という出力を得ます．後は所得の大小に応じてグレーに色分けしてプロットする必要があるのですが，ArrayPlotのデフォルト設定をそのまま利用します． 実はデフォルトで数字の大小に応じたグレースケールに変換して，Plotしてくれるのです（便利やー）．

この部分，おもしろくない割に最初から作ると面倒なので，助かります *2

さて，ここで注意すべきは，次の点です．

  • 位置情報の動的更新をせずに，所得値の情報だけ正しい位置に表示できるかを確認する．
  • これができたら，動的更新のModuleに組み込む*3

次にパネルを二つ並べるために， GraphicsRowをManipulateのオブジェクトとして指定します．

GraphicsRow[{
  ArrayPlot[{Table[Random[], {5}]}],
  ArrayPlot[{Table[Random[], {5}]}]
  }]

括弧の数に注意しましょう．できあがりはこんな感じです．

f:id:hamada7418:20161108162851j:plain

右のパネルが所得値を表しています．黒いほど高額所得です．

実行すると二つのパネルが同期して動きます．

人種に関する選好だけを持っていたとしても，人種間の所得分布に差があるとき，結果としてincome segregation(金持ちが金持ち同士で固まって住む)が生じる様子が分かります．これは統計的差別がもたらす意図せざる結果の一種と解釈することも可能です．

ただし，ここまでの結論は計算しなくてもほぼ自明なので，今後は

  1. 片方の人種だけに同人種への選好があっても，income segregationが生じるのか？
  2. 人種選好がincome segregationに及ぼす影響の推定
  3. 逆に人種選好はないが，incomeへの選好だけあるとき，人種間segregationはどの程度生じるのか
  4. 地価上昇ダイナミクスを取り入れると，人種間不平等がどのくらい高まるのか

などを考えてみようと思います．

*1:なお，自分で試してみたい学生は，授業で配布したサンプルコードをいきなり修正せず，追加すべきModuleだけプロトタイプを作り，実験してから組み込むようにしてください

*2:こういう内蔵関数が豊富なところが，Mathematicaのいいところだなあ

*3:はじめてシミュレーションに取り組む学生は，1ステップずつ作業を進めてください．複数のクロージャをまとめて更新すると，どこでエラーが出たか分かりません．慣れるまでは，1つずつエラーを直しながら進めてください．面倒なようですが，このほうが結果的に早いのです．さもなければバグがどこにあるか一日中探さなくてはなりません．はじめのうちはバグを直すより，バグを見つける方が難しいのです