Kosugitti's BLOG

アンドロイドは正規分布の夢を見るか

この画面は、簡易表示です

2016 / 12月

2016年の10大ニュース

年末になってまいりましたので,一年を振り返ることをしたいと思います。

今年はどんな一年だったんだろ,とブログを読み返して見たら,やっぱりというか結構というか,年齢に伴う精神的な変化をグズグズいうてたり,Stanのコード書いてたり,というそういう一年だったのね。それはまあそうだったと思う。というか,ブログに書いたのはそういうのが多かったってことかしらね。

ということで,改めて今年の自分的トップテンを考えて見たいと思う。

1.引っ越しました。十年ぶりに。

2.娘が入院したのは辛かった。

3.秋の夜の叱咤激励

4.ドラマに共感しまくったね

5.Stanでベイズで遊んだね

6.ディズニーランドに行ったのはいい思い出

7.本が出版されましたよ

8.地味に「体力の衰え」について

9.本を出版する準備をしていますよ,そしてGithub楽しいよ

  1. 足の上にPCを落としたよ。骨はいってないと思うけど。

 

まあやっぱりでかいのは引っ越しですよ。家を建てた云々は昨年のイベントかな。出来上がった家に入ってね,引っ越ししたり生活のリズム整えたり,春はもうずっとそれだったような気がする。家が広くなって嬉しいのは,喧嘩が少なくなることですね。アホみたいなことだけど,子供と喧嘩して「あっちに行け(顔見たくない)」といったときに,そうできるっていうのがね。しばらくして落ち着いたら,仲直りできるしね。子供の部屋を作ってあげられたのは良かったけど,うちの子達があんなにもお片づけできない奴らだとは知らなかった,というのはショッキングなことだけど,徐々に教えて行きますかね。

2番目は末っ子よー。末っ子の肺炎入院は大変でした。出張先で「入院になった」って聞かされたしね。出張に出る前に体調が悪いのはわかってたけど,帰って即病院に行かなあかんのは心理的にも辛かったね。寝泊まりするにしても寝泊まりするスペースないし,とにかく家族の一大事に家族全員で対応しないとなんともならん,そんな日々でした。一週間弱だけど,なかなか「ためされた」ところだったな。

3番目は全く個人的なこと?なんだけど,確実に人生において大事なことだったね。後輩に叱咤激励(正確には叱咤95激励5)されて,論文を書かないと行けないなあという事実を思い知らされたというかね。遊んでいたわけじゃないんですよ,本とか書いてたんですよ。でも本って査読なし,論文じゃないものじゃないですか。しかも著者じゃなくて編者だしね。仕事してないと思われるんだろうな。編者も大変なんだけど,それもこれも含めて,研究者の評価基準はそこじゃない,って教え子に教えられました。最近はちゃんと書いてます。酔った勢いで年度内に3本書くっていうてしまったからねw 俺約束は守る男なんで。

4番目。これはなんていうかなあ。「真田丸」「シン・ゴジラ」「君の名は」「逃げるは恥だが役に立つ」は本当に楽しんだ。12月中旬は真田丸ロスと逃げ恥ロスが一気に来たので,ほんまに辛かったです。シン・ゴジラは二回目を見に行こうとしてまだ行けてない。君の名は,いいなあ。何歳になっても,恋心とか,正義とか,生き様とか,昔は感じなかった色々を感じながら行きてます。シン・ゴジラ見たら公務員のよさもわかるしね。君の名は。をみたら若い恋心の良さもわかるしね。逃げ恥見たら「呪い」から解き放たれた世界の良さを思うしね。ガッキーと星野源のダブルヒロインもよかったしね。源次郎,生きたね・・・。

5番目は研究的な話なんだけど,今年はもう全くベイジアンでした。豊田先生の本を使った読書会も始めたし,学生にもベイズの良さが広がっているよう。自分の昔のデータもベイズでやりなおそうっていうことにしたし,なんだかもう,恋です,恋。最後の恋。たっぷり大恋愛してやりますよ。

6番目は,家族でディズニーに行った思い出です。子供たちが楽しんでくれたのが何よりでした。実はそのちょっと前に下見に行ったりしたんだけど(!)それもこれも込みで,今年はディズニー,たのしかったなあ。また家族旅行したいです。

7番目。遅いようですが,二年前からかかってた仕事,翻訳書の出版ができました。誤訳とか読みにくいところとか,誤字脱字とか色々あるとは思いますよ,思いますけど,なかなか頑張ったんです,いい本を,学界に貢献するような仕事をしたんじゃないかと自負しています。少し古い本なんだけど,内容は今にも生きるものだしね。ほんと,なるべく多くの人に読んでもらいたいです。出版社の方に出版記念祝賀会を開いてもらったこともいい思い出になりました。そのあと叱咤激励に続いたんですが・・・!

8番目。特に何をというわけではないんだけど,体があちこちガタついて来ました。体重は落ちました。友人と誓った三つの約束,「ラーメンのスープは飲まない」「食事と晩酌を分ける」「夜いくら飲んでもいいけどつまみはダメ」の三つで体重はぐんぐん落ちた。最大83キロだったのが,今77ぐらいです。5キロ減といえば十分じゃないかな。一番下がった時は73ぐらいまで行ったんですよ。戻って来たのは誓いの裏技を見つけ始めたからだねw とはいえ,言いたいのはそれではなくて,そういうことをしてもそこまでしか下がらないってことです。今後は「食べない」じゃなくて,「筋肉をつける」という明確な意思を持たないと,痩せて行かないんじゃないかな。基礎代謝がそもそも落ちてるんですよ。こればっかりは年齢のせいだと思うんだなあ。来年どうしよう。

9番目。別の本のプロジェクトが動いてます。多くの人が関わっているので,なかなか自分ノペースだけでできる仕事ではないです。当初の予定では10月出版になってるけどw流石にそんなにジュ長じゃなかった。でも原稿を書いて管理するツールにGithubを使ってるんですが,これの使い方がわかって使い勝手もいいな,ということがわかって来ました。バージョン管理,大事。今後もバージョン管理を意識したプロジェクトづくりができるんじゃないかと思ってます。

さて最後の10番目ですが,これはクリスマスの今日起こったことです。母親のPCがもう壊れたようなので,新しいのを買ってあげようとおもったんですが,前のPCのHDDからデータがサルベージできるかもしれない。できないとしても取り出して壊したほうがいいよね,とネジを外して取り出して・・・という作業をしていた時,あと少しで終わりというタイミングでPC本体が机からぐらりと落ちまして,左の親指にどーん。痛い。痛かった。すぐに冷やしたので今のところなんとかなっているけど,明日の朝,足の指の色がおかしかったら病院に行こうと思います。最悪骨にヒビぐらいいってるかも。うーん,年末に怪我ってなんかやだなあ。親や妹に「疲れた顔してるから・・・」とか言われちゃった。やっぱり俺は寝不足だと怪我する子なのね(今0時半)。おおごとにならないことを祈ってますが,最後にこれをトップ10にいれときましょう。

今年はあと一週間ぐらいあるけど,これ以上でかいイベント・事件は起きないと思うので,この辺で来年に気持ちを向けて生きたいと思います。

 

みなさん良いお年をお迎えください。

 



StanなMCMCで探索的因子分析やってみた

この記事はStan Advent Calender 2016にエントリーしたものです。

やっていることは,ここ数年関わっている探索的因子分析をベイズ推定するっていう話なのですが,最初に言っておきますと,探索的因子分析ってベイズモデリングに向いてないと思うのですよw

ベイズモデリングが何か,という話は一言では終わらないですが,私は「データ生成メカニズムをデザインするモデリング」と理解しています。なので,認知心理学的なモデルなんかにはもってこいだと。そもそも心理学は個人差をどう考えるか,というような世界でもあるので,階層的なモデリングのためにとてもいいと思うのですね。

一方,探索的因子分析はその名の通り「探索的に」データを捉えるわけですから,データがどういう構造しているかわからない,というのはそもそもちょっと違う。もちろん,n因子構造だと仮定してモデルを書くと,ということになるんですが,CFAならともかくEFAだとまあありとあらゆる係数を推定していくことになるのです。なので,ちょっと大変だなあ,と書き終わってから思うたわけです。

 

さてさて,本題です。
因子分析というのは,$$N$$人の個人$$i$$が$$M$$個ある項目$$j$$に反応したデータセット$$x_{ij}$$(ただし$$N \ll M$$,つまり項目数より個人の数が十分に多い)について,$$k$$個の因子という潜在変数を仮定して,因子と項目との関係を表す因子負荷量$$a_{jk}$$と,因子と個人との関係を表す因子得点$$f_{ik}$$との積和の形で表すモデルです。

$$ x_{ij} = \mu_j + \sum_{k=1} a_{jk}f_{ik} + e_{ij}$$

ここの$$e_{ij}$$は誤差項ね。これをベイズ的に考えるときは,基本的に潜在変数が正規分布に従うことを仮定する,というだけなので,それをそのまま素直にモデリングしてあげれば一丁あがりです。

ただ,探索的にやる場合は因子負荷行列に制約をかけてあげる必要があります。因子負荷行列は一般に,直交行列になるような制約を置いてあげるものですが,SEM的には$$a_{jk}=0$$(ただし j < k の場合)という制約をかけることになります(豊田先生の「共分散構造分析【技術編】」を参考)。これをコード化したら次のようになります。

[crayon-5bcac99fd6259356027843/]

ここのtransformed parametersブロックで面倒なことをしています。整数でindexというのを宣言して置いて,このループの中だけでカウントするインデックスを用意しました。もう少しうまいやり方もあると思うんだけど,とりあえずこれでできます。

 

ただし!このやり方はすごく時間がかかります。時間がかかる理由は,データサイズに比例して推定しなければならないパラメタが多いということ。N人A因子あればN*A分の因子得点を推定しますからね!こりゃめんどうです。

実際,一般的な因子分析は,基本モデルにいくつかの条件(因子得点は標準化得点,誤差の平均はゼロ,誤差因子と共通因子の相関はゼロ)を加えて,

$$ \Sigma = A \Phi A’ + \Delta $$

特に直交因子モデルの場合は$$\Phi$$は単位行列なので

$$ \Sigma = AA’ + \Delta $$

とした分散共分散行列で考えることになります。

Stanで多変量正規分布を推定するのはなぜか大変重たいので,ウィシャート分布から推定することにし,因子負荷行列もこれスキー分解したベクトルで与えてやると,Stanコードは次のように大変短くなります。

[crayon-5bcac99fd626c008482944/]

データから平均値を取り去った(中心化した)データセットとして与え,それの積和行列がtransformed dataブロックでspdとして作られており,それがwishart分布に従うようになっています。

これでStanはちゃんと推定してくれます。が!これでお話は終わりではありません。実はここから得られるMCMCサンプルの推定値を考えても,0ばっかり!みたいなことになるのです。それには三つほど問題が関係しています。

一つはMCMCサンプルごとに出てくる因子の順番が違う問題。これはサンプリング後に因子の順番をソートして整えることで対応できます。

もう一つはMCMCサンプルごとに出てくる因子負荷量の符号が違うこと。これは因子毎の負荷量の総和が正になるように整えるとか,任意のMCMCサンプルをつかってその因子毎の負荷量の総和の符号に合わせる,といった恣意的な基準を作って対応できます。

もう一つは因子の回転についてです。stanの中では回転作業をしていませんので,得られたMCMCサンプルを全部回転することで対応できます。

 

・・・とここまで書いておわかりのように,サンプルを得た後の作業が結構大変なんですね。実際にやってみると,サンプリングよりその後の回転やソートのほうが時間がかかったりします。

その作業をするコードを載せても良いのですが,それはRの話でStanの話じゃなくなってしまいますから,この記事はここまでです。Rの関数として,psychパッケージのようにfa,と書いたら探索的ベイジアン因子分析してくれるものを現在開発中ですので,興味がある方は個別にご連絡を。

 

それにしても,確率モデルをほぼそのまま書くだけでできちゃうStanって,本当に便利ですよね。みなさんもEnjoy Stan!

参考文献;

[amazonjs asin=”4254126646″ locale=”JP” title=”共分散構造分析 技術編―構造方程式モデリング (統計ライブラリー)”]



習作;ベイズファクター

岡田(2014)のベイズファクターの計算をStanでやってみた件。

これについては,次のスライドに詳しい。

ただ,途中のStanコードがちょっとわかりにくいのと,計算をstan内部でgenerated quantitiesを使ってやったらいいんじゃないかと思ったので,そのように書き直して見た習作。

[crayon-5bcac99fdb12c128327424/]

推定結果は論文とほぼ一致した。BF1,BF2,PMP1,PMP2を参照。

[crayon-5bcac99fdb13e205362697/]



Stanで体重の推移をみつめてみた(状態空間モデル)

この記事は Stan Advent Calendar 2016 のエントリーです。

Stanいいですね。いいですね,Stan.

私はもっぱらRから触るrstanユーザーなんですが,いろいろなところから触れる間口の広さがいいですね。

もちろん「尤度関数を書かなくても推定してくれる」とか「細かい設定をしなくても大体ちゃんとやってくれる」といったStan本来のありがたみも忘れちゃいけません。

なにより,ベイジアンモデリングのハードルがすごく下がった気がします。頭の中でイメージしていることが,そのまま書けるっていうのは,とてもいいことだと思うのです。

この「頭の中のイメージがそのまま書ける」ということを一番実感したのは,状態空間モデルのお勉強をしているときでした。なので,その時のことを記事にして,Stanへの感謝とエールに代えたいと思います。

時系列データと状態空間モデル

時系列的なデータの場合,一つ前の時点との関係(自己相関)が強いもんだから,そのままの分析じゃダメで,時系列分析という別の考え方が必要になる。これはまあいいです。

でも,従来の(最小自乗法,最尤推定法的な)時系列分析はなんだかわかりにくい。

時系列分析については,いろいろな本も出ています。あるいはこちらのサイトはすごく丁寧にわかりやすく書いてくれていますので,時間かけて読めば何をやっているかはちゃんとわかります。でも,ちゃんと面倒なことがわかります。

というのも,とある時点での推定値がでても,それが期間全体で成立するかどうかを前後にわたってチマチマ調整していかなければならないからです。私の理解したイメージでは,基本的にt時点からt+1時点への回帰で,それがワンステップ前,ワンステップ後でも推定地として大丈夫かどうか,均してからつかうという感じ。すごい力技な感じがするんですね。実際,モデルを決める,パラメタを求めて,フィルタをかけて,スムージング,という順番でやらないといけないのです。

そこで出ました状態空間モデル。これは

  • 目に見える観測の世界
  • 目に見えない状態の世界

を区別して,状態は連綿と続いていく,観測はそれに誤差が乗っかる,というイメージでモデリングしていくのですが(これについてもこちらのサイトが詳しいです!),Stanはこのイメージをモデルに書くときに,本当にそのまんま!という感じなのです。これを実際の例で見ていきましょう。

 

状態空間モデルによる体重の予測

私はなぜか,毎朝起きたら体重と体脂肪を測定して記録する,というのを日課にするようになっているのですが,この習慣が2012年ごろから続いており,時折記録できない日(出張で自宅にいない=いつもの体重計がない場合)なんかもありつつ,2016年7月末の段階で1077レコードありました。

実は昨年末帰省したときに,京都ラーメンが美味しすぎてたくさん食べていたら,体重が最高83kgまで達しまして,これはまずいんじゃないか(当方身長173cm)と,今年に入ってから少し体重を下げる努力をしております。この辺の効果なんかも考えてみたいところです。

 

さて,状態空間モデルは「観測値」と「状態」を区別すると言いました。ここでは,観測値=体重計の目盛りです。でもうちの体重計は,結構数字を出すのに悩むんですね。77.4?77.5?77.3?いやまて77.7・・・?とか考えてから,今日は77.6kgです,と表示したり。体の傾け方とかで数値がぶれるのかしら,と思うのですが,まあこの目盛りは信じるしかない。

こうした測定誤差とは別に,肉がついてるとか落ちたとかいう実際の状態もあるわけで,そこを別なものとしてモデリングします。下のイメージ図参照。

状態空間モデルのイメージ図

さて,これをStanでかくと,モデルの部分はまさにこのイメージのまま,なのです。

[crayon-5bcac99fdcef1167222050/]

いかがでしょうか。観測データ(W[n])は実際の肉(mu[i])に誤差(sdV)がついて変化します。実際の肉(mu[i])は昨日の肉(mu[i-1])に誤差(sdW)がついた形で変化します。
モデルのイメージ図をそのまま書いただけになっているのです!フィルタとかスムージングとかいりません!

欠測と予測

自宅の体重計で毎朝測るのですが,自宅で目覚めない日もあるわけで,そういう日は測定できません。欠損値になります。

状態空間モデル2欠測

それでも肉がなくなるわけではないので,状態は続いていきます。これをモデリングする場合は,測れなかった日の観測値もパラメータとして推定してやれば良いのです。

Stanは残念ながら欠測値に対応していませんから,「これは欠損だよ」「観測データはちゃんと揃ってるよ」という形で渡してあげなければなりません。私はここで,欠測のところに特殊な数字(9999)を入れて,Stanに渡してあげることにしました。

[crayon-5bcac99fdcf04113580679/]

ついでに。
時系列的な変化は,やはり未来に向けての予測がしたいですよね。予測,というのも「(まだ)観測されていないデータ」という意味で,欠損値ですから,同じアルゴリズムが使えます。
状態空間モデル3予測

予測したい日数を欠損としてデータセットに加えてやれば,この問題は解決です。先ほどのコードにちょいと付け足します。

[crayon-5bcac99fdcf0a656931493/]

これでよし。さて,Stanのコードはどうなっているかというと,

[crayon-5bcac99fdcf10632429704/]

観測方程式のところに,missを埋めていくインデックス,nmissというのを作りました。これは0からカウントしていきます。

もし観測(w[i])が9999じゃなかったら,それまでどおり肉の状態(mu[i])に誤差がついた形であらわれているのですが,もし9999だったら=欠測だったら,肉の状態(mu[i])からでてきた推測すべきパラメータ(Miss_W)になるよ,ということです。

そして欠測していようがしていまいが,肉の状態は1-step前の状態から推移し続けます。

これで推定してみた結果は次の通り。

 

[crayon-5bcac99fdcf17830643328/]

上が観測の誤差成分,下が状態の誤差成分です。観測成分による誤差が状態の誤差の約2倍ぐらいあります。
実際のお肉の状態は次のような感じ。

[crayon-5bcac99fdcf1c409318983/]

ほうほう・・・。

さて,この記事は事前に描いて今日公開するよう予約投稿してあるのですが,今年のお正月から7月末までの前半期をつかって,この12月1日の体重を予測しますと!

[crayon-5bcac99fdcf23986939487/]

実は74.1950kgなお肉の状態になっているはずです!当たってるかい?今日の俺よ。
まあ95%確信区間でいくと69.86-78.03kgの間に入って入ればStanの勝ちということです!
確信区間も含めて推移をグラフにしてみると,測定が終わった次の日からパッと幅が広がっていくのがよく見えます。そりゃそうだよねー。

 時系列予測

ともあれ,いわゆる時系列分析をするよりも,Stanでもベイジアンモデリングしたほうが直感的にわかりやすいし,簡単だと思うので,超オススメです。
もちろんこの「状態」のところに,予測するモデルを入れていって,より複雑なモデリングにすることができます。私は曜日によって違うんじゃないかとか,出張の前後で体重が変わるんじゃないかとか,いろいろ細かな変数を入れて推測させて見ましたが,ほとんど効果がなかったので,ここではそうしたモデリングの部分はご紹介しませんでしたが。

 

 

さあEnjoy Stan!




top