2018-11-07

DAアルゴリズムにかんするメモ

RによるDAアルゴリズムの再現

DAアルゴリズム（GSアルゴリズムとも呼ばれます）のコードをRで書きなおして簡単な実験をやってみました．オリジナルのコードは奥村晴彦先生の『C言語による最新アルゴリズム事典』技術評論社を参照しました．

Cではprefix increment，つまりループの中に ++ v を入れることでコードが簡略化されています．同じ記法はMathematicaでも可能だったので，Mathematicaへの翻訳は簡単でした．しかしRには仕様上prefix incrementがないらしく，その部分で多少の変更を余儀なくされました．

DAアルゴリズムそのものの説明は，ここでは省略します（近々刊行予定の本で説明していますので、よろしければそちらを参照してください（宣伝））．

『その問題、数理モデルが解決します』 https://www.amazon.co.jp/dp/4860645685

###################################################
# GS Algorithm　概略
###################################################
#関数gsは引数として，男性の選好(proposer)と女性の選好(receiver)を
#入力すると，GS Algorithmに基づく男女のマッチングを出力します

#二つの引数proposerとreceiverは男女である必要はありませんが，
#先行研究にならって，以下男女の喩えを使います．

# はじめにproposer, receiverの選好を，集団別にlist型で入力します．
# 例えば男性が求婚する時，proposerとして
# list(c(2, 3, 1), c(1, 2, 3), c(3, 1, 2)) を入力する
# これは
#男性1の選好は　　女2>女3>女1
#男性2の選好は　　女1>女2>女3
#男性3の選好は　　女3>女1>女2
#を意味します.
#同様にreceiverは,女性が持つ男性に対する選好のリストです．

関数内のローカル変数の意味は以下の通りです．カウンタがたくさんあるため，構造は結構複雑です．

########################################
# gs内のローカル変数の意味
########################################
#prefb:男性が持つ,女性に対する選好． prefg:女性が持つ,男性に対する選好．
#nextgirl:プロポーズする順番を記録．各男性が独立に持つ．振られる度に1増加
#tempboy:仮マッチした男性のリスト,最初はn+1(番兵)が入る．
#gid:男性がプロポーズする相手女性のid, 
#b:for用のカウンタ．男性1～nに対応
#s:現在プロポーズ中の男性id
#t:振られた男id．この値をsに代入する事で，
#「今回振られた男t」が自動的に次のsになる．

結果以下のようなコードになりました *1．

##################################
# 関数gsの定義
##################################
gs130<- function(proposer,receiver)
{ n<-length(proposer)
prefb = proposer #男性の好みは引数proposerをそのまま入力
z=receiver  #一旦，女性の選好リストを別のzに代入する
for(j in 1:n){
    for(i in 1:n){  
        z[[j]][i]=which(receiver[[j]]==i)}#for j
}#for i
prefg =z #女性の選好を順位表記にする
print("男性の選好") 
print(paste0(matrix(prefb))) 
print(paste0("女性の選好") )
print(paste0(matrix(receiver)) )
for(i in 1:n){
    prefg[[i]] = append(prefg[[i]], n + 1)
}#女性の選好に番兵n+1を追加
nextgirl = rep(1, n)  #個人別のプロポーズ用カウンタ．1位からスタートして，2位，3位と順番に進む
tempboy = rep(n + 1, n)  # while条件用　番人sentinel  
for(b in 1:n){
    s = b
    while(s != n + 1 ){ #sが1,2,3である間ループ*)
        k <- nextgirl[s]
        gid =prefb[[s]][k]#k=1からスタート
        if(prefg[[gid]][s] < prefg[[gid]][tempboy[gid]])#現在の相手順位＜仮の相手順位
        {   t = tempboy[gid]#仮の相手を捨てる
        print(c("振られた男=",t))#初回はt=n+1，入れ替わるとt=以前に仮マッチした男性ID　
        tempboy[gid] = s # 仮の相手を入れ替え
        s = t  # 次のs(proposer)を指定する．
        } #if 条件が真の場合の実行文終わり
        else{#もしsが拒否されたらnextgirl[s]に1を足す
        nextgirl[s] <- nextgirl[s]+1
        }
    } # whileここで終わり
}#　for　ここで終わり

for(i in 1:n)#結果の表示
{print(paste0("男性",tempboy[i], "　と　女性",i))}
}#gs ここで終わり

さて，ランダムな選好を関数

lg <- function(n){replicate(n,list(sample(1:n)))}

で発生させると，関数gsは次のように実行できます．

gs130(lg(10),lg(10))

これで10対10の，マッチング結果が出てきます．n=3の場合の計算結果を確認してみましょう．

gs130(lg(3),lg(3))
[1] "男性の選好"
[1] "c(3, 1, 2)" "c(3, 1, 2)" "c(1, 3, 2)"
[1] "女性の選好"
[1] "c(1, 3, 2)" "c(1, 3, 2)" "c(3, 1, 2)"
[1] "振られた男" "4"         
[1] "振られた男" "4"         
[1] "振られた男" "2"         
[1] "振られた男" "4"         
[1] "男性3　と　女性1"
[1] "男性2　と　女性2"
[1] "男性1　と　女性3"

"振られた男" "4" は番兵（センチネル）です．whileループを抜け出す条件として使います．

このアルゴリズムは相手集団に対する完全な選好が必要なのですが，現実的には全員に順位をつけるのは難しいでしょう．一方でアルゴリズムの原理から，選好によっては10人中5人くらいまで順位をつけておけば事足りる場合もあることは想像できます．

そこで10対10マッチングを実行した場合に，選好リストの下位何番までの情報が必要になるのかを調べてみました．

f:id:hamada7418:20181107180016p:plain

この図は，最も望ましくない相手とペアになったマッチングがどのくらいの確率で生じるのかを表しています． 500回ランダムな選好を発生させてGSアルゴリズムでマッチングを作り，それを100回計測した分布です．例えば一番左は，proposerから見て一番順位の低い人とのペアが実現したマッチングの割合です．全体の3%程度このようなマッチングが存在します．逆に言えば，97%は最後まで選好を表明しなくても安定マッチングができるという意味です．

図の横軸-2は1～2番目に低い順位とのマッチング，-3は1～3番目に順位の低い人とのマッチングが生じた割合を示しています．

DAアルゴリズムはとても単純なアルゴリズムなのに，結果は強力です（提案側最適な安定マッチングを実現する）．

コードを書くよい練習になるので，みなさんも一度再現に挑戦してみてください．

*1:ちょっと細かい話ですが，関数の内部でreceiverの選好をidではなく順位に変換しています．例えば，ある女性の選好がc(2,3,1) 男2>男3>男1 のとき，この選好を（男1は）3番目, （男2は）1番目,　（男3は）2番目　という順位のベクトルに変換します．つまり　c(2, 3, 1)] というreceiverの選好の入力をいったん　 c(3, 1, 2)　と変換します．これは receiver側で相手を留保するかどうかを順位を参照して決定するためです

2018-10-18

Mathematica+TikZ

タイトル通り，MathematicaでつくったグラフをTikZでお化粧してTeXに取り込む方法のメモです．

まず素材となるグラフをMathematicaでちゃっとつくります．

Plot[PDF[LogNormalDistribution[0, 1], x], {x, 0, 5}, 
 PlotStyle -> {Black}, Frame -> True,
 FrameTicks -> {{{0, 0.2, 0.4, 0.6}, None}, {{0, 1, 2, 3, 4, 5}, 
    None}}]

このグラフにlogn_pdfという名前をつけて保存しておきます．

次にTeXのFigure環境の中にTiKZでグリッドを描きます．\includegraphicsで先ほどつくったlogn_pdfを読み込むのですが，tikzpicture環境の中で読み込むのがポイントです．

\begin{figure}[H]
    \centering
\begin{tikzpicture}
\draw[step=1,black!0] (0,0) grid (12,5);
\node[anchor=south west] (image) at (1.5,-0.5)
    {\includegraphics[width=0.7\linewidth]{logn_pdf}};
\node [right] at (4.432,3.5) 
{$\displaystyle \frac{1}{\sqrt{2 \pi \sigma^2} x}\exp \left\lbrace -\frac{(\log x-\mu)^2}{2\sigma^2}\right\rbrace$};
\node [right] at (0.7,4.8){$f(x)$};
\node [right] at (10,0){$x$}; 
\end{tikzpicture}
{\small 対数正規分布の確率密度関数$f(x)$のグラフ. $\mu=0,\sigma=1$}
\end{figure}

できあがりー．

f:id:hamada7418:20181018230630p:plain

グラフの中と周りにTeXの数式がきちんと挿入されました． R MrakdownやJupyter notebookでも同様のことができます．ぜひお試しを．

2018-09-09

TeXの柱とノンブルの設定

jsbookクラスをつかった場合の柱の設定方法を忘れないためにメモを残しておきます．

最終的な出力目標は

偶数頁のヘッダ左側に章タイトル，奇数頁のヘッダ右側にセクションタイトル

偶数頁のフッタ左側に頁番号，奇数頁のフッタ右側に頁番号

であると仮定します．これをデザインAと呼ぶことにしましょう．

\chapterコマンドはデフォルト設定で開始1頁目に\thispagestyle{plain}を実行するので，章のはじめだけ設定したスタイルが表示されないことがしばしばあります．まずはこれを制御します．

直感的に分かりやすいやり方は，\chapterコマンドのすぐ後に自分で\thispagestyle{empty}や\thispagestyle{fancy} を使って，\thispagestyle{plain}を上書きしてやることです．

しかしこの方法だと全ての\chapterに同じ命令を書き込む必要があるのであまりスマートではありません． \newcommandでコマンドを定義すれば多少はラクになりますが，章が多いとはやはり煩雑です．

一番簡単な方法は pLaTeX2e 新ドキュメントクラスの指示通りにプリアンブルに

\makeatletter
\renewcommand{\chapter}{%
  \if@openright\cleardoublepage\else\clearpage\fi
  \global\@topnum\z@
  \secdef\@chapter\@schapter}
\makeatother

と書く方法です．この\renewcommandによって，章開始1頁目は自分で設定した柱のスタイルに一致します．どうしてこれでうまくいくのか，実はよくわかっていません．おそらく， jsbook.clsの\chapterの定義にある\plainifnotemptyの削除により， \thispagestyle{plain}を実行しないのではないかと推測します．\plainifnotemptyというコマンドは全体のpagestyleがemptyでないならplainにせよ，という命令です．ですからこの命令が有効だと，empty以外のpagestyleを使うと，無視されるという現象が生じるのでしょう．

次に，デザインAを実現するには\thispagestyle{fancy}の中身を設定する必要があるので，これを調整します．例としてはこんな感じです．

\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhead{} % clear all header fields
\renewcommand{\chaptermark}[1]{\markboth{第\ \thechapter\ 章\,\, #1}{}}
\fancyhead[RO]{{\footnotesize \rightmark}}
\fancyhead[LE]{{\footnotesize \leftmark}}
\renewcommand{\headrulewidth}{0pt} %ヘッダの罫線を消す 
\fancyfoot{} % clear all footer fields
\fancyfoot[LE,RO]{{\footnotesize \thepage}}%頁番号の表示を偶数頁と奇数頁に分けて指定．
\cfoot{}

\fancyheadとfancyfootのオプションは

R: 右側 L:　左側 O:　奇数頁 E:　偶数頁

を組み合わせて指定します．例えば \fancyfoot[LE,RO]{{\footnotesize \thepage}} は，「LEが偶数頁のフッターの左側」「ROが奇数頁のフッターの右側」に頁数を示す，という指定です．

もし\frontmatterに\chapterを使っている場合は，

\renewcommand{\chaptermark}[1]{\markboth{第\ \thechapter\ 章\,\, #1}{}}

の部分を

\renewcommand{\chaptermark}[1]{#1}{}}

などに改めて定義してください．こうすることで柱を「第0章　イントロダクション」ではなく「イントロダクション」という見出しに変更できます．一方で\mainmatterは「第1章　先行研究のレビュー」のように，{\markboth{第\ \thechapter\ 章\,\, #1}を使います．

要するにポイントは

\renewcommand{\chapter}で章開始頁を修正
\frontmatter用に\pagestyleを定義
\mainmatter用に\pagestyleを定義

です．\frontmatterと\mainmatterの柱が同じなら定義は1種類でOKです．

2018-08-30

tikzのメモ

描き方を忘れないためのメモです．

tikzで関数のグラフを描きたいとき，

\draw [very thick] (始点座標) to [out=5, in=-90] (終点座標);

のようにオプション[out=5, in=-90]で始まりの角度と終わりの角度を指定すると簡単に描けます．効用関数系のグラフだったら，大体これでいけるでしょう．

f:id:hamada7418:20180830200233p:plain

正規分布のような曲線を描きたい場合は

\draw plot[smooth] (\x, {exp(-0.5* \x * \x)/sqrt(2 * pi)});

こんな感じで関数を記述すればtikzの描画機能だけでグラフを作成してくれます．なおオプションの[smooth]を忘れるとカクカクした線になってしまいます．

f:id:hamada7418:20180830200253p:plain

座標を書く場合は\foreachコマンドで繰り返しを記述すると楽できます．

\draw[thick, ->] (-5,0) -- (5,0) node [below] {$x$};%x軸
\draw[thick, ->] (0,-1) -- (0,5) node [left] {$y$};%y軸
\node at (0,0) [anchor=north east] {O};
\foreach \x in {-4,-3,-2,-1,1,2,3,4}
\draw (\x cm,.2) -- (\x cm,0) node[anchor=north] {$\x$};

最近はMathematicaで図を作成→pdf出力→include→tikzでお化粧という方法を多用していました．tikzだけで完結するなら，そのほうが楽かもしれません．

2017-11-30

Agent-based modelのチュートリアル

MesaというPythonライブラリの中に所得分布の生成モデル（Dragulescu & Yakovenko 2000）を題材にしたチュートリアルが含まれている．なかなかよくできているのだが，学部生向けとしてはちょっとコードが複雑すぎるので授業用資料として簡略化してみた*1．

この例でなんとなくPythonの特徴も分かるし，Agent-based simulationの雛形として演習の課題や卒論に使えると思う．

重要なのはコードじゃなくて，モデルのアイデアのほうだ．Dragulescu & Yakovenko のアイデアはシンプルでありながら現象の本質を表現したよいモデルだと思う．

詳細は下記のjupyter notebookをご覧ください

github.com

*1:クラスオブジェクトを使わなければもっと短く書けるのだが，クラスの例を示すために，あえてエージェントをクラス化してある

2017-08-17

ggplot2による事後分布の視覚化

絶対忘れるので自分用にメモしておく．

目的

Stanで推定した結果を分布のまま直接比較したい．

ggplot2用のデータフレームの作り方

流れはこんな感じ．

Stanで推定した結果をfitに記録（CSVで記録すると結果的に遅かった．sample_file からのサンプル列の取り出しは面倒．たいした量じゃないので今回はRDataを保存してloadする）．
stanfitオブジェクトをas.data.frame()でデータフレームに変換
パラメータの各列を1列にスタック（95年と2015年で年齢幅が違うのでindexを小細工したうえでforを回した）
ggplot用のインデックス列を作成（年齢，性別，調査年）．インデックス列は同じラベルがmcmcサンプルの数だけならぶ．今回は1000×3（chain）個．
このときスタックの順番は男女男女とし，性別のインデックスを行方向に対応させる．
plotしたいパラメータとインデックス列をデータフレームにまとめ，ggplotにわたす．
ggplotのfill=で1つの枠内に男女を並べ，facet_wrap(~ age, ncol = 5)で年齢別に並べる．結果年齢別に男女比較が行列型に並ぶ
調査年比較は，パラメータをスタックする際に，2015，1995の順番で並べる．

コードはこんな感じ． fitに2015年度データの推定結果，fit2に1995年度データの推定結果が記録されている．

結果の一例．

f:id:hamada7418:20170817150851p:plain