今週は、余り何もする気がしないので、オフィスに出かけてこの本を読んで例題を自分でとくことにする。プログラミングだとか、データ分析というのは、コンセプトの理解と、 hands on で実際にやることを同時にすることが効率がよいからである。休日のオフィスは静かで良い。毎日が日曜日ならよいのに。

この本は内容は比較的平易である、practiceに割り切っていて、余り理論的なこととか証明とかは書いていない、これはある意味とてもよい方針である、 social science系のBayesianの本は、中途半端な数学と、中途半端な実践であまり自分としては読む気がしない。

ただ、この本は一切証明などがないので、証明だとか各々の分布関数の分散（例えばベータの分散てどんな式だったっけ？）を知りたい場合は、別の本を読まないといけない。また、練習問題では普通に、例えばサムはベータ(3,12)を事前分布として考えているのだけど、その分散は？　とか聞いてくるけど、この本の中ではベータの分散などは一切触れていないので数理統計の授業をべつに取った人などには良いだろうけど、最初にこれから始めたら解けないかもしれない。

Casella & BergerのStatistical Inferenceとセットにして使うととてもよいと思える（というか自分がその組み合わせである）、この本の裏にある、確率分布の一覧が大変便利である。予定としてはこの本の例題をできるだけたくさん解きたい。

時系列の分析をしていると、時々月によって標本数が少なかったり、月によってデータの質に多少問題があったりする。そういう際に、半ば自動的にシステマチックに、より信頼度の高い推定する方法が無いかと考えて、いたところ、前から考えていたのだけど、ベイジアンのやり方でやったら良いかなと思いやってみる。

結果から先に載せると、ブルーのラインは、月ごとの(cross sectional surveyである、panel ではない)推定値を単に計算したもの。オレンジのラインは、 posterior である。


教科書どおりに、priorを考える、normalを想定して likelihood (観測データ) からの推定値も normal とすると conjugate になり都合がよいのでそうする。ただ、中心極限定理によれば、平均は正規分布するのでそれほどおかしい仮定ではないので問題はなし。



posterior の平均と分散はこうなる、これも基本的には教科書どおりである。ちなみに分散の部分を良く見て、これを普通の分散 sigma^2/n_o で割ってみると分かるが、この比率は常に１以下になる、すなわち posterior の分散は、観測データだけの分散よりも小さい、これは非常に優れた特性である。お陰で 95%信頼区間を実際に計算してプロットしてみると(Bayesian なら credible region と呼ぶらしいけど)青い縦に走っているバーよりもオレンジの方が狭いことが分かる、これは素晴らしい。



同じデータを使っても、prior を考えるだけで estimator の分散が小さくなるのは面白い、ま当たり前といえば当たり前なのだけど。

新しい本を出版、正確にはそのうちの１章を書いただけだけど。なかなかStataてメジャーにならないんだよな。こんな優れたソフトないと思うんだけど。この北大路書房はプロフェッショナルでよい出版社である。僕も日本で授業するんだったらこれ使うよな。ま、日本でどれだけ、調査データを使う授業があるかは不明だけど。

書誌情報などはこちらに

今週のプロジェクトとして、Call centerの電話の履歴をとってきて分析をしている。ただし、今日は大いにはまって、３時間ほど時間をとられてしまった。

それは何かというとCATIシステムは全ての電話の履歴を管理しており、何月何日の何時何分何秒に電話をかけたという履歴がある、また１つの番号に何度も電話をかけるので、データの形としてはcall level dataすなわち、同じ電話番号が当然何度も登場する。このデータをperson level dataにまず変換するのだけどそれが一苦労だった、基本的には同じ番号のレコードの中から時間が一番新しいものだけを残せばよいのである。CATI systemはデータをMicrosoft SQL serverに保存してあり、ODBCというデータベースとwindows applicationの橋渡しをするプログラムの設定をした後、Stataにデータを読み込める。

Stataにおいて時間変数が例えば time と呼ばれていたのであれば

gen date = day(time) とか Year = yeah(time) とかすれば日にちや年が取り出せるのは知っていたのだけれども、どうやって分や秒を取り出せるかは今まで試したことが無かったので分からなかった。Stataのマニュアルのdate & timeの章は非常に分厚くて読む気になれないが仕方ないので読む、でもまったく意味不明である。

時間変数を普通のreal型で取り出すとその整数部分は 17476である、しばらく考えて365でこれを割ってみると47.879である、マニュアルを見てみるとStataの時間は1960年の１月１日を起点にしているとあるので、２００７年はちょうど４７年後であるからこの変数は単なる「経過日」であることがわかる。この、小数点部分をとってきて　２４を掛ければ時間が得られ、その整数部分を捨ててさらに小数部分だけをとって、６０を掛ければ分が手に入り、同様に小数部分だけにまた６０を掛ければ秒が手に入ることが判明。またこうして得られた秒はSQLサーバーで表示されているものと同じである。

ただし、以上の処理は大変面倒くさいので必ず別のやり方があるはずだと思って調べると、Stataは大まかに分けて６通りの時間の数え方のフォーマットを持っており、経過年は %tdであり、これを秒を単位にした時間や、週や月や四半期に変換も出来ることが判明。それぞれに返還するための関数の変換表も data management manualにある。ここまで理解するのに１時間半ぐらいかかった。このマニュアルは非常に分かりにくいのである。

そうして、計算したもののなぜだか秒がずれてしまう。何度やっても秒がずれる、これはどこかでrounding errorが起きているはずだと思い、SQLからODBCさらにデータを読み込む過程のどこでこれが起きているかなどを一通り見るが分からない。

しばらくして、作った変数を眺めていると、データの型がfloatになっていることに気がついた。Floatでも小数点の38桁まで扱えるはずなので、この程度の誤差なら問題ないはずである。1秒は0.01666分でありそれは、0.0000694444日であり小数点の８桁ぐらいあれば十分なはずであるとPOST ITの裏に計算してみる。それでもと念のためにdoubleで変数を作ってみると解決した。ここまで３時間。Stataは内部では、１９６０年の１月１日からの経過の１０００分の１秒単位で計算しているらしく、double型で変数を作らないと精度が保てないらしい。

tempvar tc_time hour min sec

gen double `tc_time' = cofd(asktime)
gen `hour' = hh(`tc_time')
gen `min' = mm(`tc_time')
gen `sec' = ss(`tc_time')

desc

普通そんなこと思いつかんよ、まったく。だいたいどんな状況で１０００分の１秒単位の分析をStataでやらないといけないのか、エンジニアはC++か何かで計算するだろうし、余り細かすぎるのも困ったものである。

大抵の教科書にのっているbeta-binomialの例がゲルマン本にあったで見てみる。証明自体はCasella & Bergerの本で昔解いたのであるが、久しぶりに紐解いてもう一度やってみる。ゲルマンの例では、事前分布の情報の確かさだけを変えた例を載せているのだけど、せっかくなのでdataの標本数を変えてみる、30,400,1000という３つの組み合わせを試す。

	Prior				data		Posterior
ID	a/(a+b)	a+b	a	b	n-size	note
1	0.515	2	1.03	0.97	30	Prior not informative	0.435
2	0.515	2	1.03	0.97	400		0.445
3	0.515	2	1.03	0.97	1000		0.445
4	0.515	400	206	194	30	prior stronger	0.509
5	0.515	400	206	194	400		0.481
6	0.515	400	206	194	1000	sample stronger	0.465
7	0.515	5000	2575	2425	30	Prior rules	0.515
8	0.515	5000	2575	2425	400		0.510
9	0.515	5000	2575	2425	1000		0.503

これで分かるのは、最初の３つの例、事前分布が単なる一様分布でほとんど何の情報もない場合、事後分布から推定するmedianはdataに依拠したものになる。ヒストグラムだと上段の２つである。一方で、事前分布の標本数が増えるとヒストグラムは下段のようになる。事前分布はbeta分布なのだけど、標本数が増えれば(a+b)当然分散も減っているのがヒストグラムを見ても分かる。また、一様分布というのはbetaの特殊な例であるということもこれから良くわかる。

 ああ、今日は良く勉強した、寝るかな。

StataからC++などを呼び出すことが出来るということで、plug inの使い方のページを見てみたところ余り難しくないので、仕事をとりあえず終えてから数時間取り組んでみたところ、Micro targetingの際に使うモデリングプログラムが吐き出すC++のコードをコンパイルしてStataから呼べるようになった。

これは、なかなかすごいことであるC++でCSVを読んだりする部分や、GUIなどのインターフェースを作るのはだるいなあと思っていたところなので、Stataをフロントエンドに出来るのであればそれに越したことは無い。これはすごいことで、ニューラルネットワークなどのプログラムは大抵Stataには実装はされていないけれども、C++なんかのコードとしては手に入る場合があるのでそれを組み込んでしまえるということだ。

１つ問題としてはlaptopのcygwinに入っているGCCコンパイラで作ったDLLは64bit Windowsで動いている64 bit Stataでは動かないことが判明。GCCは64bit OSの為のバイナリも作れるらしいのだけど、どうもうまくいかないCygwinのGCCのヴァージョンが少し古い(3.4.4)なのが問題なのかもしれない。無料で手に入る Microsoft Visual C++ 2005 Express Editionもあるのでそちらでコンパイルしてみると当然32-bit Windowsでは動くが64 bitは動かない、有料版のVC++ならデフォルトで64bit のアプリが作れるらしいが、無料だとそうは行かないらしい（当然か）

しかし、こんなページを発見How to: Configure Visual C++ Projects to Target 64-Bit Platforms　マイクロソフトも太っ腹じゃないか。今度試してみよう。

1. 来月はAAPORの発表があるので準備をしないといけない
2. 日本の知人と書いている本の原稿をしあげないといけない
3. 論文というか報告書の修正もしないといけない
4. 会社に目を向けてみれば、先週忙しくて調査員データの分析が予定より遅れている
5. MNにおけるターゲティングの分析の修正版Power pointを明日ぐらいには仕上げないといけない
6. また明日10:00から別のbrain stormミーティングがありその為の試案も書かないといけない。

ああ、忙しい。それはともかくThe New Political Targetingを購入。これは先日会社を辞めたAmyが行った会社の社長でもある。Micro targetingというのはこの数年共和党に対抗して民主党陣営も始めることにした政治キャンペーンにおけるマーケティング手法の応用である。要するに米国では選挙人名簿はpublic recordであり手に入るので、それとconsumer databaseをつなげて、投票行動を予測しようというものである。こんなものはマーケティングの人間やクレジットカード会社は何年も前からやっていたことだろうけど、政治の世界でもそういう潮流が生まれているということだ。ミネソタ州全員の有権者名簿とそれ以外の購買データなどをくっつけたものを、もらってきたので(6GigaもありPCに入れると大変なのでサーバーに放り込んだ)、このデータが面白いのは、普通サンプリングするとRDD調査なんかでは誰が住んでいるのかわからないのだけど、全住民の名簿（選挙人登録してない住民のデータもどこから集めてきたのかついてきた）。日本でも始めればよいのにと思う、ただこういうデータは日本では作りにくいだろうな。

今やっている仕事はMNで行われた2006年の知事と、上院議員の選挙の際に作られたそのスコアがどれくらい正確にモデル化できていたかということの検証、オリジナルのモデルは調査データを元にCAHIDなんかで変数選択をしてロジステックモデルで予測値を作っている、この程度のモデル作成というのは、調査方法論の世界では欠測ウェイトの作成の文脈でずいぶん前からされていた。

下の図は別の調査データで聞いた対象の政党帰属意識と、別の会社が作った政党帰属の予測値(MNは選挙人登録の際に、政党を聞かないので支持政党がわからないのである)を回帰してみたもの、実際のフィット(R^2)は余りよくない。Stata9になってからグラフィックが強力になったので、いくつかのグラフを重ねたり、自分で計算した信頼区間なんかも表示できてなかなかよい。点推定から95%信頼区間を計算してそれぞれのカテゴリ毎のバンド幅を計算する。こういのは、もはやお手の物である。

今日は、タンゴ教室に行く、我ながらマルチタスクは苦手で、最初はステップがあまりうまくいかない、順番にマスターしてゆくしかあるまい。アレキサンドリアのold townはなかなか小奇麗で閑静な街である、Dupont Circleのあたりなどよりも気分が良い。今日は単にdrop inでいったのだけど、visiting ｎ先生が来ていて、アルゼンチン人である、話を聞いていると端々に--タンゴっていうのにはかの国の文化が反映されているのだと思う。僕がまだへたくそなのは、これは文化が違うからっていうのもあるだろう、rhythmからして違うんだから。先生はタンゴのCDを聞くことを薦めてくれた。

Bayesian Approaches to Clinical Trials and Health-Care Evaluationという本を少し前に買ってしばらく寝かしておいたのだけど、改めて読んでみると非常に面白い、久しぶりに本を追いながら出来るだけ例 題を、鉛筆で解く。最近やっている仕事で、結構Bayesianのideaが利用できそうなので、この際にしっかり学んでおこうかとおもう。しかし、部分積分のやり方を忘れてしまった自分にがっくし。あせって教科書を取り出す。そういえば、こういうときのためにmapleというすばらしいソフトがあったことをを思い出してインストールすることにする。mapleを使って簡単な例題を自分で解いてみれば、思い出せるというものだ。

後日談：　mapleをインストールして、いろいろな組み合わせを解いてみたら、すっかり思い出した。よかったよかった。MLEが求められなかったら大学院で統計学を学んだとはいえないし。。。

帰りに上記のの本を読みながら思ったことは、研究とは、あまりにホットなものを選んでもしょうがない気がする、なぜならそんなものは、誰でも思いつくし自分がやらなくてもほかの人がやってしまうからである。流行りものを追うというよりかは、自分の好きなテーマを決めて、それを追いかけるのが良い。そうといっても、重要でない事をやってもしょうがないので、中期的にテーマを考えるのが良いか。ただし、企業の中でやる場合は重要でないというテーマはあまりありえない。無論一企業や分野にとって重要な問題が、ほかの世界では重要でないということもある、それは仕方の無いことであると同時に、あまり蛸壺にはいらないような目配りもしておかねばならない。

会社で、自分の専属のアシスタントをつけてもらえるように交渉中である、そうすぐにつけてもらえるわけではない、粘り強く交渉して、証拠というか実績を積み上げながら説得するしかないのである。世界はゆっくりとしか動かない。それで助かるときもあるが、非常にfrustrationを感じる場合もある。

かつて片足を突っ込んでそれ以来離れていた領域、それは調査データで得られる自由回答の分析。修論の１部としてやろうかと思ったが当時はskillが不足していて結局やらなかったんだけど最近仕事で使う必要があって再挑戦である。

STATAは実は正規表現でstringを処理できるので、キーワードをマッチしてデータを拾うことは出来るのである。自由回答の分析には何通りかあって１つは、日本語なら茶筌のような形態素解析ソフトに通して言葉を分割して、全てのuniqueな単語からmatrixを作って分析するというものである。分析自体はMDSのようなものから、クラスター分析、数量化なんかが使われる、ただしこれはマトリックスがでかくなり、このようなsparce matrixの分析は難しい。実際のところ調査データなどで拾える語彙というのは結構少なくて日本のデータの場合でも結構少なかった、そちらの方の数字はすでに失念したけれども、昨日分析した米国のデータだと1000人で8つのopen ended questionのなかから約4000語しかなかった。もちろん質問のframeの仕方でこれは変動するのだけれども、やっぱり少ない。とはいえども、行列データとして 4000 x 4000のデータを扱うのは結構だるい、STATA9になって付属した行列言語の処理システムのMATAなら扱えるけれども、これをそのまま扱うのはあまり賢くない。

もう１つ良く考えられる方法で僕の会社でも使う方法はCODE FRAMEを作っておき、ある単語が（単数でも組み合わせでも）引っかかれば、ある特定争点に言及しているとみなしてコーディングする方法である。この場合のCODE FRAMEの作り方は恣意的であるといっても良いし、研究者のsubstantiveな知識に基づくといっても良い。僕はこの部分を機械化するプログラムはすでに書いた。キーワードを例えば下のようなデータセットに別に作っておいて、STATAに読ませて内部のメモリーに辞書を保存し、これらのキーワードが含まれていたら0-1のダミー変数に返すというものである。１つの争点について複数のキーワードを設定できるようにしてありプログラムは争点ごとにどれだけ独立のキーワードがあるか、またおのおののものが特定のキーワードかそれとも複数かを判別できるようにもしてある。STATAのプログラムはMATAが出来てから非常に柔軟性が高くなりこんなものも比較的簡単に書けるようになった。

Issue	keyword1	Keyword2
war in raq	war	iraq
war in raq	WMD
health care	medicare

先日仕事のために書いた、STATAの回帰モデルのプログラムは実はウェイトは無視していたのだけど、それが推定値と各種指標に反映されるように改造。

基本的にはSRSを仮定しているので、単純にGLSもしくはWLSとして扱う。それで点推定は簡単に求められるのだけれども、R^2などがこの場合どうなるかなどは教科書に書いてなくて、ずいぶん考えた。最終的に、友人にStataのマニュアルをスキャンしてもらったところ、Stataのregコマンドでどう計算しているかが判明、

点推定は簡単。個人ごとの抽出確率のsqrtをとって、対角行列にし、それを独立変数の行列であるXと従属変数のベクトルyにあらかじめかけてやればあとは、通常の(X'X)^-1 X'y　でもとまる。

W = diag(sqrt(weight))
B = W*X
z = W*y
b = invsym(B'B)*B'z

問題はR^2で、これにはSSE(residual error)とSST(total error)を計算してやればいいのだけど、SSTの方の計算が出来なくてずいぶん紙に書きまくって考えたところ、最終的には何とか解けた。ああ、大変だった。次はこのロジステック回帰版を作らないといけないのである。

sst = z'z - (y'*weight)^2/sum(weight)
sse = z'z - b'B'z
ssr = sst - sse

R2 = 1 - sse/sst
aR2 = 1 - (sse/(n-p))/(sst/(n-1))