2016年2月29日 (月)
百年の孤独 (Obra de Garc´ia M´arquez)
[a]
ガブリエル ガルシア=マルケス / 新潮社 / 2006-12
再読だが、20歳前後で読んだきりだから、初めて読んだようなものだ。年始からかなり長いこと読み続けていて、結末に到着したときは、ほっとするような、寂しいような気分であった。
緑衣の女
[a]
アーナルデュル・インドリダソン / 東京創元社 / 2013-07-11
声
[a]
アーナルデュル・インドリダソン / 東京創元社 / 2015-07-29
これもしばらく前の休日に読んだもの。いま流行りの北欧ミステリである。著者はアイスランドの人気作家である由。なかなか面白い。あまり殺人など起こらなさそうなお国柄のせいか、捜査がかなりのんびりしているところが良い。
いと高き貧しさ――修道院規則と生の形式
[a]
ジョルジョ・アガンベン / みすず書房 / 2014-10-25
著者の「ホモ・サケル」シリーズ(という連作なのだそうだ)については、全然わからないのだが、中世の修道院規則の話が面白くて最後まで読み終えることができた。それは規範的であると同時に構成的であるというか、そうやって腑分けすることさえ難しいというか、ただ生きられているとしかいいようがないというか、そういうものだったらしい。ふうん。
哲学・思想(2011-) - 読了:「いと高き貧しさ 修道院規則と生の形式」
ナショナリズム: その神話と論理 (ちくま学芸文庫)
[a]
橋川 文三 / 筑摩書房 / 2015-08-06
著者の橋川文三という人、1983年没だそうだ。もっと昔の人かと思っていた。
面白かったという記憶はあるんだけど、しばらく前に読んだ本なので、ええと、どう面白かったんだっけか...そうそう、この本で明治元年の隠岐コミューンのことを知ったのだった。
神道 (ちくま学芸文庫)
[a]
トーマス カスーリス / 筑摩書房 / 2014-10-08
移動中や待ち時間に読み耽っていた本。ここんところのベストワン、非常に面白い本だった。読み終わってから調べてみたら、著者は比較宗教学の偉い人なのだそうだ。
神道についての解説書というより、実存的 vs. 本質主義的スピリチュアリティという軸を武器に、宗教性とはなにかという大きな問題を解きほぐす内容。なるほどねえ...
憲法と平和を問いなおす (ちくま新書)
[a]
長谷部 恭男 / 筑摩書房 / 2004-04-07
増補新版 法とは何か (河出ブックス)
[a]
長谷部 恭男 / 河出書房新社 / 2015-07-11
年始だったか、時間待ちかなにかで入った書店で何の気なしに新書のほうを買ったのだが、これがまた面白い本で、一気に読了。同じ著者の本をもう一冊探して読んだ。
自分で読もうとまでは思わないけれど、ケルゼンという法学者の話が面白かった。この人にいわせると、違憲の法律というのは無効ではないそうだ。なぜなら、無効ならそれは法律ではないから。ということは違憲の法律も憲法の授権を受けているわけで、従って、国会は憲法から違憲の法律を制定する権限を受け取っていることになる... 分析哲学の本を読んでいるときにも感じることだけど、徹底した議論というのは、そこはかとないユーモアを生むなあ。
ユーロ危機とギリシャ反乱 (岩波新書)
[a]
田中 素香 / 岩波書店 / 2016-01-21
菱田春草 (別冊太陽 日本のこころ 222)
[a]
鶴見香織 / 平凡社 / 2014-09-19
ノンフィクション(2011-) - 読了:「神道」「ユーロ危機とギリシャ反乱」「憲法と平和を問い直す」「法とはなにか」「菱田春草」
私を連れて逃げて、お願い。 3 (ビームコミックス)
[a]
松田 洋子 / KADOKAWA/エンターブレイン / 2016-02-25
完結巻。どん底の境遇にある男女の胸の痛くなるような純愛を描く、傑作であった。
ふつつかなヨメですが! 1 (ビッグコミックス)
[a]
ねむ ようこ / 小学館 / 2016-02-12
可愛らしい兄嫁は実は筋金入りのヤンキーで... というコメディ。風俗や価値観がひと世代上なんじゃないか(特攻服に竹刀はないんじゃないか)と思ったのだが、掲載誌は「スピリッツ」とのこと、読者も高齢化しているだろうし、これでいいのかも。
おひとり様物語(6) (ワイドKC Kiss)
[a]
谷川 史子 / 講談社 / 2016-02-12
おんなの窓 5
[a]
伊藤 理佐 / 文藝春秋 / 2016-01-28
さらば、佳き日 (1) (it COMICS)
[a]
茜田千 / KADOKAWA/アスキー・メディアワークス / 2016-01-14
コミックス(2015-) - 読了:「私を連れて逃げて、お願い。」「ふつつかなヨメですが!」「おひとり様物語」「おんなの窓」「さらば、佳き日」
ほんとはこんなことをしているバアイではないのだが、もう朝なので、片づけるついでに... 最近読んだ本。
忘却のサチコ 5 (ビッグコミックス)
[a]
阿部 潤 / 小学館 / 2016-02-12
プリンセスメゾン 2 (ビッグコミックス)
[a]
池辺 葵 / 小学館 / 2016-02-12
俺の姫靴を履いてくれ (2) (MFコミックス フラッパーシリーズ)
[a]
須河 篤志 / KADOKAWA/メディアファクトリー / 2016-02-23
百姓貴族 (4) (ウィングス・コミックス・デラックス)
[a]
荒川 弘 / 新書館 / 2016-02-25
アップルシードα(2)<完> (モーニング KC)
[a]
黒田 硫黄 / 講談社 / 2016-02-23
昔の士郎正宗のSFコミックス(懐かしい...)のリメイクではなく、その映画のコミカライズである由。それにしても、どこからどうみても黒田硫黄さんのマンガになっている。旧作「大日本天狗党絵詞」を思わせる、独特のリズムと躍動感がたまらない。才能ってあるんだなあ、と感心しながら読んだ。
コミックス(2015-) - 読了:「忘却のサチコ」「プリンセスメゾン」「俺の姫靴を履いてくれ」「百姓貴族」「アップルシード」
しばらく前にメモを取った奴。どうも最近は気ばかり急いて、メモだけ取って放り投げているものが多い。
以前ベイジアン・ネットワークのソフトの講習会に行ったとき、休憩時間に開発元の講師の方を捕まえて、目的変数が二値ひとつ、説明変数が複数という状況で、ナイーブベイズ学習器とロジスティック回帰とどう使い分ければいいんですかね? と尋ねた。その若い方は「さあ...」という返事であった。わからないというよりも、ちょっと一言では答えられない、ということだったのかな。
実のところ、この問いはずっと心に引っかかっている問いで、自分なりの経験的な答えのようなものはなくはないのだが、いちど専門の人に訊ねてみたいものだ、と思っていた。
Ng, A.Y., Jordan, M.I. (2001) On Discriminative vs. Generative classifiers: A comparison of logistic regression and naive Bayes. Advances in Neural Information Processing Systems 14 (NIPS 2001).
そのものずばりのタイトルなので、NIPSなどというオドロオドロシイ出典にも関わらず手に取った。Google様によれば被引用件数1314。まじっすか。
著者いわく。
入力$x$とラベル$y$を受け取って同時確率$p(x, y)$をモデル化し、所与の$x$についてベイズルールで$p(y|x)$を予測し、もっともありそうなラベル$y$をピックアップする。これが生成的分類器である。いっぽう、事後確率$p(y|x)$を直接にモデル化するのが識別的分類器である。
Vapnik(1998)いわく、汝分類課題のためにより一般的な課題を解く勿れ。回り道してないで、さっさと$p(y|x)$をモデル化しろ。とこのように、ふつうは識別的分類器のほうが好まれている。また、必要な事例数はモデルの自由パラメータ数と線形の関係にあることが多いという言い伝えもある。こういう信念は正しいのでしょうか。
いま、「入力とラベルの同時尤度を最適化しましょう」という基準に従う分類器を$h_{Gen}$とする。それと同じパラメトリックなファミリーのモデルで、「条件付き尤度$p(y|x)$を最適化しましょう」ないし「$p(y|x)$に閾値をつけて予測したときの訓練誤差を最小化しましょう」という基準に従う分類器を$h_{Dis}$とする。このペアを生成-判別ペアと呼ぼう。たとえば$p(x|y)$を正規分布、$p(y)$を多項分布とすれば、正規判別分析とロジスティック回帰モデルがペア。入力が離散なら、ナイーブ・ベイズとロジスティック回帰がペアだ。
どっちが優れているだろうか? ナイーブベイズとロジスティック回帰を比べる。
入力変数は$n$個とする。出力ラベルは二値とする。$m$個のiidな事例をとってきて訓練セット$S$とする。
赤コーナー、生成モデルの入場です。ナイーブベイズ分類器。$S$を使って$\hat{p}(x_i | y)$と$\hat{p}(y)$を予測する。
入力変数がすべて二値の場合、$S$におけるイベント発生数を$\sharp s\{\cdot\}$と書くとして、予測式は
$\displaystyle \hat{p}(x_i = 1 | y=b) = \frac{\sharp s\{x_i=1, y=b\}+l}{\sharp s\{y=b\}+2l}$
と書ける。$l=0$とすれば経験的推定値。伝統的には$l$を正の整数にすることが多く(1とか)、これは確率のラプラス・スムージングにあたる。さて、$h_{Gen}(x)=T$という予測値を返すのは、下式が正であるとき、その時に限る。[めんどくさいので省略するけど、$\prod$は$i$から$n$まで回る]
$\displaystyle l_{Gen}(x) = \log \frac{(\prod \hat{p}(x_i|y=T)) \hat{p}(y=T)}{(\prod \hat{p}(x_i|y=F)) \hat{p}(y=F)}$
入力変数がすべて連続の場合は...[めんどくさいので省略。要するに正規判別]
青コーナー、識別モデルの登場です。ロジスティック回帰。予測式は
$p (y=T | x; \beta, \theta) = 1/(1+\exp(-\beta^T x -\theta))$
$h_{Dis}(x)=T$という予測値を返すのは、下式が正であるとき、その時に限る。
$l_{Dis}(x) = \sum \beta_i x_i + \theta$
で...
ここから難しい数学用語が乱舞するので、ほぼ3pにわたってまるきり理解できないのだが、とにかくこういうことがいえるらしい。
- データが無限大なら、ナイーブベイズよりロジスティック回帰のほうがエラーが小さい。
- 識別モデルがエラーの漸近値に到達するのに必要な事例数は、せいぜい$n$のオーダー。
- 生成モデルがエラーの漸近値に到達するのは...それよか速い。(「それよか」のところの説明は!もう全然わかんない!一行たりともわかんない!呪文か!呪文なのかこれは!!)
私のような呑み込みの悪い哀れなアホのために(ちがうかな)、事例が示してある。15個のデータセットで事例数を変えながら成績を比較。ほんとだ、最初はナイーブベイズが勝ち、データが増えるとともにロジスティック回帰がキャッチアップするデータセットが多い。
。。。というわけで、肝心のところは全然理解できていないんだけど、結論のみ学ばせていただきました。へー、知らんかった、そうなんですか。
純粋に素人の感想ですけど、これ、入力空間のなかの学習事例の分布によっても違うんじゃないかと思うんですが、どうなんでしょうか。あくまで経験的な印象だけど、ロジスティック回帰で判別するにあたっては、(もちろんもともと線形判別が可能な課題だとして、) 学習事例の数そのものよりも「閾値の近くにある」学習事例の数が勝負になるような気がする。誰に教わったのかはっきりしないが、「人数は少なくてもいいから、合否ぎりぎりの奴を15人は集めろ」という言い伝えが耳に残っている。その点はナイーブ・ベイズでも同じことなのかなあ。
論文:データ解析(2015-) - 読了:Ng & Jordan (2001) 対決!ナイーブベイズ vs. ロジスティック回帰
2016年2月27日 (土)
しばらく前に読んだ奴のメモだけど、整理の都合上記録しておく。
Li, T., Anderson, J.J. (2009) The vitality model: A way to understand population survival and demograaphid heterogeneity. Theoretical Population Biology, 76, 118-131.
ふつうの生存モデルとは全然違う生存モデル、その名も4パラメータ生命力モデルをご提案します。という論文。
なにかのきっかけで目にして(あ、そうだ。CRANのtask viewにパッケージが載っていたのだ)、気分転換用の読み物としてストックしていた奴。生存モデルは仕事で使うことがあるので、目先を変えて「このたびは生態学での最新のモデルを適用しまして...」なあんて言えたら、ちょっとカッコよくね? なんて色気を出したんだと思う。ははは、センスないな俺。全然ちがう話であった。
どういうモデルかというと...生存曲線を下式とする。
$l(t) = l_a(t) l_v(t)$
$l_a(t)$は生命力とは関係ない外発的な死に対応する。これはポワソン過程に従うと考えて
$l_a(t) = \exp(-kt)$
とりあえず$k$は定数としておく(時変させるのは今後の課題とのこと)。
$l_v(t)$のほうが本命。これは「生命力が0になったことによる死」に対応する。えーと、実は各個体は生命力(vitality)$v$を持っている。開始時点での生命力$v_0$は個体を通じて正規分布する(SDを$\tau$とする)。生命力は時々刻々と次のように変化していく:
$dv/dt = -\rho + \sigma \epsilon_t$
$\rho$は生命力の減衰率の平均、$\sigma$は生命力のばらつき(この2つは定数だと考える)。$\epsilon_t$はホワイトノイズ。こういうのをウィナー過程というのだそうだ。で、$v$が0に到達したら死ぬ。というメカニズムから導出される生存関数が$l_v(t)$である。
ここからの式の導出がまっっっったく理解できないので、数式も省略するけど、結局, 生存関数$l(t)$は次の4パラメータを持つ超ややこしい式で表現される: ドリフト率$r=\rho/\bar{v}_0$, スプレッド率$s=\sigma/\bar{v}_0$, 初期生命力分布の分散$u=\tau/\bar{v}_0$, そして外的死亡力$k$。何だか知らんが、これ、MLで解けるんだそうです。
ふつうの(?)生存分析でも、Coxモデルの共変量のパートに個体差を表す切片項(フレイルティ)を入れることがある。初期生命力$u_0$に個体差を与えるのはそれと同じかな? と思ったのだが、著者に言わせると次の点が異なる。フレイルティ・モデルでは個体のフレイルティは時変しない(低い個体が間引かれるので、分布は変わるけど)。また、ベースライン死亡力からは独立である。いっぽう生命力モデルでは個体の生命力が経験によっても変化するし、生物学的基盤を持つモデルになっている、とのこと。
モデルの特徴についていろいろ説明があったけど、省略して...
生命力モデルが活躍する事例を4つ紹介。
その1、medflyパラドクス。medfly (辞書によればチチュウカイミバエ)はオスのほうが平均寿命が長いが、最後まで生き残るのはたいていメス。これはオスとメスの死亡力が途中で逆転するからであることがわかっている(パラドクスっていうほどのことではないんじゃないの...?)。先行研究ではメスのフレイルティの分散が大きいんだと考えているけれど、このモデルにいわせれば、メスはスプレッド率が大きいのであって、初期生命力の分散はむしろ小さい。前者はきっと多産性のせいだ(現に子供のいないメスは平均寿命が長い由)。後者は...と延々と生物学的な考察があるけど、省略。
その2、食事制限。寿命を延ばす効果があるといわれているが、従来の手法では平均生存時間や生存曲線を比べるくらいしかできない。ショウジョウバエのデータに当てはめてみると、ドリフト率が低下することがわかる。とかなんとか。すいません、適当に読み飛ばしました。
その3、その4は疲れたのでパス。
ディスカッション。古典的な死亡力モデルの文脈ではフレイルティを生命力ということがあるけど、俺たちの生命力モデルとは全然違うからな、気を付けて物を云えよ。生命力モデルは状態依存死亡力モデルともちがうからな、気をつけろよ(←状態依存モデルってのがあるのか...なにそれ...)。死を初通過過程として捉えるってところが最大のポイントだかんな。云々。
。。。いやー、こういう研究分野があるのね。さまざまなtime-to-eventデータのモデリングについての抽象的な話ではなくて、なにかの生物種の生存曲線を、生物学的にリーズナブルな形でどうやって説明するかという、かなり実質科学寄りの話なのであった。心理学でいうと分散学習の効果をどうやってモデル化するかとか、マーケティングでいうとBassモデルをどう拡張するかとか、そういうレベルの話である。半分も理解できてないけど、自分なりに勉強になりましたです。
ちゃんと読んでないせいだと思うけど、わかんなかった点をメモ。(1)観察に打ち切りのあるデータにはどう対応するの? いっかいカプラン・マイヤー推定量を出して、やおら当てはめるのかしらん。それともこの分野では打ち切りなんて問題にならないのかなあ。(2)共変量はどうやってモデルに組み込むの? それとも、この分野にはそもそもそういう問題意識がなくて、観察された生存曲線を記述することだけが関心事なのだろうか。よくわからん。
論文:データ解析(2015-) - 読了:Li & Anderson (2009) 生命力モデルによる生存分析
2016年2月22日 (月)
Carsternsen, B. (2007) Age-period-cohort models for the Lexis diagram. Statistics in Medicine, 26, 2048-3045.
age-period-cohortモデル(いわゆるコウホート分析のこと)について、レキシス図(人口統計の文脈でみかけるわけのわからない図)の観点から解説するという、私にとっては悪夢のような論文。こんなの勉強してる場合じゃないんだけど、仕事の都合でコウホート分析が必要になり、観念して手に取った。
本文では自明視されているけれど、頭の整理をしておくと...
この論文が想定しているのは、「西暦xxxx年(period)、xx歳(age)の、西暦xxxx年生まれの人々(cohort)が、 これこれの率で死にました」というような表の分析である。本文では一貫してtabulationと呼んでいるが、コウホート表と訳す。モデリングもコウホート表から出発する状況を想定している。
レキシス図とは、横軸にperiod、縦軸にageをとった図のこと。本来period=cohort+ageだから、 ある人の生命は右上に向かう傾き1の線分として表現される。
コウホート表の個々の値はレキシス図のある領域に対応している。その領域は正方形だったりそうでなかったりする。たとえば「西暦xxxx年(period)、xx歳(age)の、西暦xxxx年生まれの人々(cohort)」という領域は、左辺が垂直(年始)、上辺ないし下辺が水平(誕生日)である直角二等辺三角形になる。この論文は一般性を重んじているので、領域の形は定めていない(おかげでかなり回りくどい)。なお、本文ではこの領域を一貫してsubsetと呼んでいるので、サブセットと訳す。
いわく。age-period-cohortモデル(以下、APCモデル)を、レキシス図からの観察を記述する一般的モデルとして示そう。
話を先取りすると、本稿のポイントは3つ。(1)表をつくるときにに情報を捨てるな。(2)age, period, cohortは連続変数として捉えろ。(3)率の絶対値を示せ。
まずは、視覚化しましょう。古典的な図は次の4つ。
- 横軸にage, 縦軸に率, periodで層別した折れ線
- 横軸にage, 縦軸に率, cohortで層別した折れ線
- 横軸にperiod, 縦軸に率, ageで層別した折れ線
- 横軸にcohort, 縦軸に率, ageで層別した折れ線
注目点は、率が比例的に動くかどうか。すなわち、縦軸を対数にしたとして、図1,3が直線か (ageごとの率はperiodと比例するか), 図2,4が直線か (cohortと比例するか)。
データをどのように表にすべきか。
レキシス図からみると、あるサブセットが1行のレコード、イベント数とリスク時間[←観察人年のことであろう]がアウトカム、age, period, cohortのサブセット内平均が説明変数である[←ここでいう平均とは、サブセットに属する観察対象者たちの平均ではなく、領域の平均だと思う]。行はできるだけ多くすること。イベント数がゼロの行があってもかまわない。
サブセットのリスク時間はこうやって求める。ageとperiodは年刻みだとしよう。age $a$ のperiod $p$ 初頭での人口サイズを$L_{a,p}$とする。period $p$初頭においてageが$a-1$であった人(=カレンダー年 $p-a$に生まれた人)における、period $p$を通じたリスク時間は、次の2式の和になる。
- age $a-1$において、$y_{a-1, p, p-a} = (1/3) L_{a-1, p} + (1/6) L_{a, p+1}$
- age $a$において、$y_{a, p, p-a} = (1/6) L_{a-1, p} + (1/3) L_{a, p+1}$
[いやー、このくだりで時間がかかった。。。
直感的にはだいたい上の式のような感じだろうなと思ったのだが、こんなにきれいな式になる理由が、考えれば考えるほどわからない。だって、「誕生日に射殺される患者」の年始からの平均余命でしょう? 病死のハザードを定数とし(累積ハザードは時間の指数関数になる)、誕生日の分布を一様として(累積ハザードは時間の一次関数になる)、このふたつの競合リスクの下での平均余命を考えるとなると、それってもっと複雑にならないか?
電車での移動中に宙を眺めてあれこれ考えたが(おかげで傘をどこかに忘れた)、どうにも埒が明かないので付録の証明を読んでみたところ、「年始に$a$歳であり、かつ年内に死んだ人の死亡日の分布」を一様分布と仮定している。この仮定の下では、なるほど、ごく簡単な積分で上式が導出できる。ここにトリックがあるわけだ。本文中には「constant mortalityを仮定する」とあるので、そんなら死亡日の分布は指数分布にすべきだろうと思うのだが、ハザードが低い状況で、1年間くらいならまあ一様分布でいいや、ってことなのだろう。私としてはマーケティング・データへの応用を考えていて、消費者行動のハザードは死亡のハザードよりもはるかに高いから、ここは気を付けないといけない...]
あるサブセットにおけるage, period, cohortの平均について。[←くどいけど、そのサブセットに属する人の平均じゃなくて、領域自体の平均の話。いまなんでこんな話をしているかというと、モデリングにおいてコウホート表から出発するからだ]
ダイアグラム上で、age vs. periodの表は長方形、period vs. cohort の図は左右の辺が垂直な平行四辺形、age vs. cohortの図は上下の辺が水平な並行四辺形になる。いずれにせよ、その範囲のage, period, cohortの平均は、単にその辺の平均である。
しかし、たとえば「1980年の年内に61歳になった人のその誕生日から年末まで」は直角二等辺三角形になる。年齢の平均は61+(1/3)である。こういう風に、平均は辺の1/3とか2/3とかになる。
以下では、任意の形状の領域におけるケース数を$D$、リスク時間の合計を$Y$と呼ぶ。その共変量は、a, p, c(=p-a)のその領域における平均である。
さあ、いよいよモデリングだ。
age $a$, period $p$ における率を$\lambda(a, p)$として、乗算APCモデルは
$\log[\lambda(a, p)] = f(a) + g(p) + h(c)$
入力される$a, p, c$は観察単位における平均である[←たとえばageを年齢で切っていたらaは満年齢+0.5、ということだろう]。率はサブセット内で一定と仮定する(そうせざるを得ない)。
あるサブセットで$(D, Y)$だとして、対数尤度に対する寄与は
$l(\lambda|D,Y) = D \log(\lambda) - \lambda Y$
率が定数だと仮定する場合を別にすれば、これは平均$\lambda Y$のポワソン分布に従う確率変数$D$の、ある観察の対数尤度と同じである。だから独立な観察に対するポワソン回帰のプログラムを使って推定できる。ただし、ケース数がポワソン分布に従っているわけではないことに注意。[...このあたり、ちょっと理解が追い付かない記述があるが省略]
モデル右辺の項を減らし、デビアンス統計量で尤度比検定することも多い(これは記述モデルなんだから検定してもしょうがない、という見方もできるけれど)。a,p,cの間隔を細かくすれば有意になりやすいことに注意。
a,p,cの効果は非線形かもしれないので、factorとして扱い水準ごとにパラメータを推定するのが古典的方法。がんの疫学的研究では5年で区切ることが多かった。パラメトリックなスムージングをかける方法もある。スプライン、自然スプライン、分数多項式など。[パラメータ数の選び方についていろいろ書いている。省略]
線形トレンドを抽出するために、factorモデルで推定した各水準のパラメータを実際のageなりperiodなりcohortなりに回帰することがある。報告の時には傾き、切片、残差も載せる。[←このくだりはこのあとの話への伏線で、いろいろ説明が書いてあるけど大幅省略]
より一般的には、
$\log[\lambda(a, p)] = f(a) + g(p) + h(c)$
の関数f, g, hをどうにかうまく決めて推定する。上の古典的方法はf,g,hをピースワイズ一定とみなしているわけだ。
関数を選ぶ際に目指すのは、(1)意味があること、(2)わかりやすいこと、(3)普通のソフトで推定できること、(4)報告された値から率を再現できること。
まずはage-cohortモデルについて考えてみよう:
$\log[\lambda(a, c)] = f(a) + h(c)$
伝統的には、コホートの参照点$c_0$を決めて$h(c_0)=0$と制約する。$f(a)$はコホート$c_0$における年齢のログレート、$h(c)$はコホート$c_0$と比べたログレート比になる。つまり
$\log[\lambda(a, c)] = \tilde f(a) + \tilde h(c) = (f(a)+\mu) + h(c)-\mu)$
と考え、$\mu=h(c_0)$と定義して、$\tilde h(c)=f(a)-h(c_0)$、$\tilde h(c)= h(c)-h(c_0)$と解するわけだ。
これをage-period-cohortモデルに拡張しよう。
$\log[\lambda(a, p)] = \tilde f(a) + \tilde g(p) + \tilde h(c)$
についてどうにかうまく制約したいわけである。どこか2つの水準と、どれかひとつの傾きを制約する必要がある[←その理由が縷々説明されていたのだけど、よくわかんなかったので今回は断念]。
ひとつの方法は: $h(c_0)$を0に制約。$f(a)$はコホート$c_0$における年齢のperiod調整済ログレートとして、$h(c)$は$c_0$に対する対数RRと解釈できる。さらに、$g(p)$を平均0に制約。$g(p)$はage-cohort予測に対する対数RRと解釈できる(ないし、基準日$p_0$において$g(p_0)=0$と制約。その場合は$a_0=p_0-c_0$における年齢効果が$p_0$の率と等しくなり、periodの効果が$p_0$に対する残差の対数RRになる)[←ぐわー、ややこしくてわかんないー]。
あとは、age-cohortでまず推定してそこにperiodを突っ込む話とか、実例とか。疲れたのでパス。
。。。途中で力尽きて、いい加減な読み方になってしまった。ま、背景知識が足りないので、仕方ない。
想像していたのと全然ちがう内容であった。文脈がよくわかっていないんだけど、 要するに、疫学では間隔を広めにとったfactorモデルで良しとすることが多いので、もっと細かくやるべしという主旨の啓蒙論文だったのではないかと思う。
読んでいてわかったのは、手法の細かいところにdomain-specificな事情が意外に効いているという点。たとえば、レキシス図上で三角形になる領域の集計値の話が多かったのは、生年x死亡時暦年x死亡時満年齢の集計表を分析することが多いという事情からきているんじゃないかなあ。こういう集計、マーケティング・データではちょっと想像しにくい(生年と満年齢のどっちかしか使わないだろう)。うーん、やっぱり社会科学系の解説を読むべきだったかも。
最後まで読んでようやく気が付いたんだけど、著者はRのEpiパッケージの開発者。なるほど、それでEpiパッケージのドキュメントはこの論文を引用しまくってるのか... (気が付くのが遅いぞ!!)
論文:データ解析(2015-) - 読了:Carsternsen(2007) レキシス図で学ぶコウホート分析
2016年2月15日 (月)
昨年末に読んでいたんだけど、記録してなかった本。
東京タラレバ娘(4) (KC KISS)
[a]
東村 アキコ / 講談社 / 2015-12-11
都市部のアラサー未婚女性を厳しく批評する人気恋愛コメディ、4巻目。
こんなことをいうと叱られちゃうかもしれないけど、この人たちの悩みなんて、実にどうでもいいと思うのである。いいじゃないですか、結婚しなくたって。四十代になったら女友達と千葉にバスツアーにいけばいいじゃないですか。一体なんの問題があろうか。
しかし、設定された状況に関心が持てないことを横に置いておけば... どんなくだらない葛藤であっても、それを深く深く掘り下げていくと、人間の普遍的な可笑しさや悲しさに辿りつくものなのだなあ、と感心させられる。その意味でとても面白い物語だ。
トロイラスとクレシダ―シェイクスピア全集〈23〉 (ちくま文庫)
[a]
W. シェイクスピア / 筑摩書房 / 2012-08
小田島訳で読んでいたのだけど、出来心で松岡訳も買っちゃったので、読み直した。
これ、やはり難解な戯曲だと思う。クレシダが何を考えているのかがわからない。たまたま昨年舞台で観る機会があって(鵜山仁演出)、その際は苦境のなかで悪女に変貌せざるを得ないクレシダの悲劇がとてもよくわかる演出だったのだが、こうして読みなおしてみると、最初からちょっとよくわからないところがある人だよな、という気もする。
反知性主義: アメリカが生んだ「熱病」の正体 (新潮選書)
[a]
森本 あんり / 新潮社 / 2015-02-20
キリスト教からみたアメリカ近代史、という感じの本であった。面白かった。
数学が歩いてきた道 (PHPサイエンス・ワールド新書)
[a]
志賀 浩二 / PHP研究所 / 2009-09-19
「私」をつくる――近代小説の試み (岩波新書)
[a]
安藤 宏 / 岩波書店 / 2015-11-21
「日本国憲法」まっとうに議論するために[改訂新版]
[a]
樋口 陽一 / みすず書房 / 2015-09-19
駅をデザインする (ちくま新書)
[a]
赤瀬 達三 / 筑摩書房 / 2015-02-04
著者はかつての営団地下鉄、つくばエキスプレスなどのサイン計画に携わった方。本書の白眉はなんといっても第5章、日本の駅の構造とサインを豊富なカラー写真で示しながら、その無節操さを遠慮会釈なしに批判するところ。
著者いわく、マーケティング・デザインとパブリック・デザインは本来視点が全然違う。後者は誰をも排除しない。なるほどねえ、と感心した。
ノンフィクション(2011-) - 読了:「反知性主義」「『私』をつくる」「駅をデザインする」「数学が歩いてきた道」「『日本国憲法』まっとうに議論するために」
デフレーション―“日本の慢性病"の全貌を解明する
[a]
吉川 洋 / 日本経済新聞出版社 / 2013-01-19
よく知らないけど、著者は反リフレ派の大物経済学者として名前が挙がる人だと思う。細かいところはよくわかんないんですが、貨幣数量説批判のところが面白かった。
学部生のヒマなときに、経済学の講義でも受けておけばよかったな。少しは難しい話にも親しみが沸いたかもしれない。
魚の文化史 (講談社学術文庫)
[a]
矢野 憲一 / 講談社 / 2016-01-09
1983年刊。著者は魚についての著作が多い方だが、本業は伊勢神宮の宮司さんだったのだそうで、神道に関する記述がとても詳しい。楽しい本であった。
中世社会のはじまり〈シリーズ日本中世史 1〉 (岩波新書)
[a]
五味 文彦 / 岩波書店 / 2016-01-21
清朝と近代世界――19世紀〈シリーズ 中国近現代史 1〉 (岩波新書)
[a]
吉澤 誠一郎 / 岩波書店 / 2010-06-19
東芝 不正会計 底なしの闇
[a]
今沢 真 / 毎日新聞出版 / 2016-01-30
ううむ... この内容なら、読むのと同じ時間だけwebで無料記事を読み漁っても同じくらいの情報が得られたんじゃなかろうか... まあ、いいけどさ。
ノンフィクション(2011-) - 読了:「魚の文化史」「デフレーション」「中世世界のはじまり」「東芝 不正会計 底なしの闇」「清朝と近代世界」
昭和元禄落語心中(9) (KCx)
[a]
雲田 はるこ / 講談社 / 2016-02-05
おかあさんの扉5 なにそれ! ?五歳児 (オレンジページムック)
[a]
伊藤理佐 / オレンジページ / 2016-02-02
34歳無職さん (7) (MFコミックス フラッパーシリーズ)
[a]
いけだ たかし / KADOKAWA/メディアファクトリー / 2016-01-23
ダーリンは70歳 (コミックス単行本)
[a]
西原 理恵子 / 小学館 / 2016-01-20
コミックス(2015-) - 読了:「昭和元禄落語心中」「おかあさんの扉」「34際無職さん」「ダーリンは70歳」
健康で文化的な最低限度の生活 3 (ビッグ コミックス)
[a]
柏木 ハルコ / 小学館 / 2016-01-29
小学館青年誌によくある情報マンガのフォーマットを忠実に守っている作品だが(なにも無理にストーリーマンガにしなくてもいいんじゃないか、というような奴)、意外に面白い。演出が良いのだろうか。
ワカコ酒 6 (ゼノンコミックス)
[a]
新久千映 / 徳間書店 / 2016-01-20
ヴィンランド・サガ(17) (アフタヌーンKC)
[a]
幸村 誠 / 講談社 / 2016-01-22
少女終末旅行 3 (バンチコミックス)
[a]
つくみず / 新潮社 / 2016-02-09
コミックス(2015-) - 読了:「健康で文化的な最低限度の生活」「ヴィンランド・サガ」「ワカコ酒」「少女終末紀行」
最近読んだ本、まずはコミックから。
美しの首 (ビームコミックス)
[a]
近藤 ようこ / KADOKAWA/エンターブレイン / 2015-12-25
86~87年の作品の再刊。
めしばな刑事タチバナ 20 (トクマコミックス)
[a]
坂戸佐兵衛,旅井とり / 徳間書店 / 2016-01-30
リバースエッジ 大川端探偵社(6) (ニチブンコミックス)
[a]
ひじかた 憂峰 / 日本文芸社 / 2016-01-29
仲能健児作品集 アジア夜話 (ビームコミックス)
[a]
仲能 健児 / KADOKAWA/エンターブレイン / 2015-12-25
シュトヘル 12 (ビッグ コミックス〔スペシャル〕)
[a]
伊藤 悠 / 小学館 / 2016-01-12
コミックス(2015-) - 読了:「美しの首」「めしばな刑事タチバナ」「リバーズ・エッジ 大川端探偵社」「仲能健児作品集 アジア夜話」「シュトヘル」
2016年2月10日 (水)
Dormann, C.F, et al. (2012) Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography, 36, 27-46.
環境研究の分野における、説明変数間の共線性とその対処に関する解説。すごく引用されているようだし、なにかの足しになるかと思って手に取った。ぱらぱらめくっただけだけど、整理の都合上、読了にしておく。
共線性の診断について、こんな早見表がついていたのでメモ。いろいろ批判がありそうだけど...
- 説明変数間の相関の絶対値。0.7以上が危険。
- 相関行列の行列式。0に近いときに危険。
- condition index。相関行列の固有値の最大値を他の固有値で割ったとき、30以上だったら危険。
- condition number。condition indexの一番大きい奴。
- Kappa。condition numberの二乗。
- VD(variance-decomposition proportions)。ある変数の分散がある固有値で表される割合。0.5以上は危険。
- VIF(variance inflation factor)。10以上は危険。
- Tolerance。1/VIF。
後半はモデリング手法の紹介とシミュレーション。聞いたこともないような手法がいっぱいあった。latent root regressionとか(要するにPLS回帰みたいなものらしいがよくわからん)、OSCARとか(要するに正則化らしいがよくわからん)。シミュレーションはめんどくさいので読み飛ばした。
回帰の文脈で、いつもストレスに感じていることがあるんだけど...
説明変数間に共線性があるとどう困るのか、その理由をきちんと整理しないといけないんじゃないかと思うのである。ひとつの理由は「解が不安定になるのが困る」だけど、もうひとつの理由として、説明変数間の関係について十分な知識がないとき(Pearlさん風にいうと「DAGが描けないとき」)、偏回帰係数の実質科学的解釈がとてつもなく難しくなるのが困る、という点が挙げられると思う。
前者と後者では問題のレベルが異なる。前者はいうなればテクニカルな話であって、サンプルサイズや測定誤差と密接な関連があるし、完璧とは言えないにせよ検出方法もある(VIFとか)。いっぽう、後者はもっと原理的な話、大袈裟に言えば「私たちは偏回帰係数を実質的に解釈するだけの能力を持っているのか」という話だ。
この2つがごっちゃにされているせいでいろんな混乱が起きているように思う。たとえば「偏回帰係数の符号が常識的な因果関係と逆だったらその変数は取り除きなさい」なんていうアドバイスがそうだ。このように信じている人はとても多い。私は静かに暮らしたいだけで、他人を説教する気はないので、そういう人がいると、率直に申し上げてかなり困る。
研究者のみなさま、ぜひ分かりやすい説明を...と願っているのだけれど、研究分野によってはそもそも後者のような問題が起きないだろうから、見境なく期待するわけにもいかない。この論文にしても、共線性についての説明がいきなりVIFの説明に移っちゃうわけで、やっぱ環境科学とかだとそうだよね、うぐぐぐ...と呻きながら読んでいた。
混乱の原因のひとつは、それぞれの現象に良い名前がないという点だと思う。上の2つの問題のうち後者の問題は、伝統的な教科書では抑制(suppression)という見出しの下で解説されていると思うが、なにをもって抑制と呼ぶかも実は相当あいまいらしい。小島隆也さんという先生の一般向け解説書では、前者を「マルチコ」、後者を「マルチコもどき」と呼んでいて、我が意を得たりと思ったのだけれど、考えてみるとマルチコ(multicollinearity; collinearity)という用語自体も、「説明変数間の高い相関のせいで解が不安定になる」現象を指して用いられる場合もあれば、(この論文のように)独立変数間の相関の存在そのものを指して用いられる場合もあるわけで、ここにも混乱の種があると思う。
論文:データ解析(2015-) - 読了:Dormann, et al. (2012) 共線性とその対処
数日前から、朝から晩までモデル・アベレージングについて考えている中年男、それがわたくしである。単にぼんやり考えているだけ。思えば不思議な人生だ。
劉慶豊(2009) モデル平均理論の新展開. 京都大学経済論叢, 183(2).
前半は解説、後半は事例紹介(ARCHモデルとかEGARCHモデルとかOGARCHモデルとかが乱れ飛ぶ。悪夢だ)。いずれも私の理解が及ぶ内容ではなかったが、ま、こんな日もある。
モデル平均の手法の紹介についてメモ:
- ベイジアンモデル平均。データを$y$として、モデル$M_k$に与えるウェイトは$P(M_k | y)$[モデルの事後確率]。Cf. Draper(1995 JRSS), Hoeting et al.(1999 Stat.Sci.), Clyde & George (2004 Stat.Sci.).
- smoothed-BIC。$BIC \approx -2 \log (\lambda_k)$なので[$\lambda_k$とはたぶんモデルの周辺尤度のことだと思う]、ウェイトを$\exp(-BIC_k/2) / \sum_s \exp(-BIC_s/2)$とする。Cf. Buckland et al.(1997 Biometrics).
- smoothed-AIC。上のBICをAICに入れ替える。これもCf. Buckland et al.(1997).
- Mallowの$C_p$を使う方法。丁寧に説明してくださっているのだろうとは思うが、残念ながらさっぱり理解できなかった。Cf. Hansen(2007 Econometrica)。
論文:データ解析(2015-) - 読了:劉(2009) モデル平均理論レビュー
Clyde, M. & George, E.I. (2004) Model Uncertainty. Statistical Science, 19(1), 81-94.
ベイジアン・モデル・アベレージングについてのレビュー。
重回帰などのいくつかのモデルにおける変数選択問題を中心に先行研究を概観。モデル空間の探索手法と事前分布の決め方について概観。木モデルとグラフィカル・モデルでの例について簡単に紹介。周辺尤度が解析的に求められないときの手法を紹介。最後に決定理論的観点からの諸問題を紹介。という構成であった。
多くの話題について、差し迫った関心がないので適当に読み飛ばしてしまったが(←「難しくてわかんなかった」の婉曲表現)、勉強になった箇所をいくつかメモしておく。
- 重回帰の変数選択問題の場合、事前分布の王道は、パラメータについてはゼルナーのg事前分布、モデル空間については独立ベルヌ―イ分布(つまり、各変数がモデルに含まれるかどうかを独立ベルヌーイ試行だと捉える)。
- 独立ベルヌーイ分布のほかに、変数の数についてトランケートつきのポワソン分布を想定するという手もあって、これはスパースなモデルを表現するのに向いている。
- ベルヌーイ事前分布の弱点は、共線性がきついとき、一連のよく似たモデル群に高すぎる事前確率が割り振られちゃうという点。[←そうですよね! この点、Hoetingたちの論文のコメントにも書いてあったんだけどそのときは意味がわからず、実データで試してみてやっと腑に落ちた。モデル・アベレージングって意外に共線性に弱いのだ]
- ベルヌーイとg事前分布の黄金ペアでやるとして、じゃあ無情報だからgはなるべく大きくしよう、という作戦は失敗する。なぜなら、ベイズ・ファクターで見るとナルモデルとすごく極端なパラメータ推定値を持つモデルが有利になっちゃうからだ[←そういうもんか...]。これをLindley-Bartlettパラドクスという。そういうわけで、とりあえずのおすすめは、ベルヌーイ分布のパラメータwは1/2 [←無情報ってことであろうか]、gは10000以下にすること。
- $\sigma^2$が固定されている限り、gとwをどう定めても、事後確率が高くなるモデルはAICやBICも良いことがわかっている。[←えええ? よくわからんけど、なんだかやる気をなくすねえ。George & Foster (2000, Biometrika)をみよとのこと]
論文:データ解析(2015-) - 読了:Clyde & George (2004) ベイジアン・モデル・アベレージング・レビュー
2016年2月 8日 (月)
仕事していて、ふと思うことがあったのである。これって、PLS回帰モデルを協調学習するとうまくいくんじゃないかしらん。ははは、まさかそんな手法はないよね... とためしに検索してみたら、なんと、パッケージがあった... 世の中は広いなあ...
Xiao, N., Cao, D, Xu, Q.S. (2015) enpls: R Package for Ensemble Partial Least Squares Regression.
Rパッケージenplsのvignett。
要するに、訓練データからリサンプリングしてpls回帰モデルをk本つくる。リサンプリングはモンテカルロとブートストラップを用意している。並列処理もできる。
係数の平均をSDで割ったのを重要性として表示してあげるから変数選択に使え。外れ値検出の機能もつけたので使え[詳細略]。 無論予測もできるので使え。k-fold CVもつけたので使え。以上。
そ... そうですか...
論文:データ解析(2015-) - 読了:Xiao, Cao, & Zu (2015) 協調PLS回帰 (Rのenplsパッケージ)
Zeugner, S. & Feldkirchner, M. (2015) Bayesian Model Averaging Employing Fixed and Flexible Priors: The BMS Package for R. Journal of Statistical Software, 68(4)。
Rのベイジアン・モデル・アベレジーング(BMA)のパッケージBMSの紹介。実戦投入前の参詣のような意味合いで目を通した。
RのBMA用パッケージとしては、他に先週読んだ論文の著者HoetingsさんたちによるBMAパッケージとか、その論文にコメントを寄せていたClydeさんによるBASパッケージとかがあり、優劣は良くわかんないんだけど。
BMA概論。ここでは線形回帰のBMAについてだけ考える[BMSはそういうパッケージであるようだ]。
説明変数の候補が$K$個あり、その行列を$X$としよう。そこから$X_\gamma$を選んで、
$y = \alpha_\gamma + X_\gamma \beta_\gamma + \epsilon, \ \ \ \epsilon \sim N(0, \sigma^2 I)$
を推定したい。でも全部入れたモデルを検討するにはサンプルサイズが足りない。
そこでBMAだ。ありうる$2^K$個のモデルについて片っ端から推定し、それぞれにベイズ定理で事後確率を与えて平均する。モデルの事後確率(以下PMPと略記する)は
$p(M_\gamma | y, X) = p(y | M_\gamma, X) p(M_\gamma) / \sum_s p(y|M_s, X) p(M_s)$
分母は分子を全モデルで合計しているだけである。分子のひとつめの$p(y | M_\gamma, X)$はモデルの周辺尤度 [←最大尤度とかじゃないのである。うんうん、ここはHoetingさんたちのおかげで理解できたぞ]。ふたつめの$p(M_\gamma)$はモデルの事前確率。
統計量$\theta$(たとえば回帰係数$\beta$)について、次のようにして事後分布を得ることができる。
$p(\theta | y, X) = \sum_\gamma p(\theta | M_\gamma, y, X) p(M_\gamma | X, y) p(M_\gamma) / \sum_s p(M_s | y, X) p(M_s)$
[←そ、そうなの??? モデルで推定された$\theta$の事後分布$p(\theta | M_\gamma, y, X)$にモデルのPMP$p(M_\gamma | y, X)$を掛けて足しあげればいいんじゃないの? なんでここでまたモデルの事前分布$p(M_\gamma)$が出てくるの?
探してみると、この論文のドラフトであろう BMSパッケージのVignetteでは、素直に
$p(\theta | y, X) = \sum_\gamma p(\theta | M_\gamma, y, X) p(M_\gamma | X, y)$
となっている。これはこの論文の誤記ではなかろうか?]
さて、モデルの周辺尤度$p(y | M_\gamma, X)$、パラメータの事後確率$p(\theta | M_\gamma, y, X)$をどう求めるか。ベイジアン回帰で求めます。[←ここがHoetingさんたちと違うところだ。彼らは個別のモデルの推定は最尤法でいいやと思ってたんじゃないかしら]
パラメータの事前分布は、切片は$p(\alpha_\gamma) \propto 1$、誤差分散は$p(\sigma) \propto \sigma^{-1}$でよしとしよう。問題は$\beta_\gamma$の事前分布。ゼルナーのg事前分布を使うことにします。
[...ここでg事前分布についての定義があるけど、パス。ダレダレの事前分布だなんて、そういう難しい話、できれば一生関わらずに過ごしたい...]
まあとにかく、g事前分布を使うってことは、係数はゼロだよねと思っているということだ。ハイパーパラメータ$g$はその確信の度合いを表す(小さな$g$は「まず間違いなくゼロだよ」に相当)。
g事前分布を使うと、$\beta_\gamma$の事後分布は$t$分布になり、その期待値は、モデル$M_\gamma$のOLS推定量を$\hat\beta_\gamma$として、$\frac{g}{1+g}\hat\beta_\gamma$になる。要するに、$0$と$\hat\beta_\gamma$を$g$で重みづけて足しているわけだ。[←へー!そりゃ面白いなあ。偉いぞゼルナーくん]。
g事前分布を使えば、モデルの周辺尤度もすごく簡単に求められる。[数式省略]
[ここでBMSパッケージの使用例を紹介... 以下、機能の紹介。]
モデルの事前分布について。以下をご用意しております。
- 一様分布。bms(data, mprior="uniform")と指定する。
- 二項モデル事前分布。モデルサイズ(モデルに使う説明変数の数)について考えよう。モデルサイズの事後分布の期待値は、ある説明変数がモデルに含まれる確率(以下PIPと略記)を全変数で合計したものだ。いっぽう、モデルサイズの事前分布は、モデルの事前分布が一様だとすると、問答無用で、$K/2$を中心にした左右対称の山形になっている[←なるほど...]。モデルサイズの事前分布を思い通りに指定できないか? そこで登場するのが「二項モデル事前分布」。モデル$M_\gamma$のモデルサイズを$k_\gamma$として、
$p(M_\gamma) = \theta^{k_\gamma} (1-\theta)^{K-k_\gamma}$
こうすると、モデルサイズの事前分布の期待値が$\bar{m}=K\theta$になる。mprior="fixed", mprior.size=$\bar{m}$とする。 - PIPの直接指定。mprior="pip", mprior.size = c(0.01, ...)という風に指定すると、最初の変数のPIPが0.01になる。
- 二項ベータモデル事前分布。二項モデル事前分布を使うと、モデルサイズの事後分布は$\bar{m}$のあたりでぐわっと高くことが知られている。そうじゃなくて、$\theta$をベータ分布からドローしようという提案がある[ややこしい...]。モデルサイズの事後分布がもっとなだらかになる。mprior="random", mprior.size=$\bar{m}$とする。
- ほかにもいろいろあるよ。
変数の数$K$が大きくて、$2^K$本のモデルを数え上げるだけで寿命が尽きる。そんなあなたのためにMCMCサンプラーをご用意しました。
- birth-deathサンプラー。まずランダムにひとつ変数を選ぶ。もしそれが現在のモデルに含まれていたら落とし、もし含まれていなかったら含める。mcmc="bd"と指定する。
- reversible-jumpサンプラー。50%の確率でbirth-deathサンプリング。残りの50%の確率で、いま含まれている変数と含まれていない変数をひとつづつランダム選びスワップ。mcmc="rev.jump"と指定する。
- MCMCじゃなくて全部数え上げる。mcmc="enumerate"と指定する。
$K$が14までならenumerate、それ以上ならbdがデフォルト。なお、MCMCのときには初期モデルをうまいこと選ぶ機能もつけてあるよ。ほかにもこんな工夫をしているよ[めんどくさいのでパス]。
パラメータの事前分布について。[このパッケージでは、どうやらg事前分布を使うというのはマストらしい]
- unit information prior。ハイパーパラメータ$g$をサンプルサイズ$N$にする。g="UIP"と指定。これがデフォルト。
- Bayesian Risk Information Criterion。$g=max(N, K^2)$とする。g="BRIC"と指定する。
- お好きなgを指定する。たとえばg=5とか。
- モデルごとに$g$を変える。経験ベイズ法(まずOLSで推定しちゃって、そのF値でgを決める)と、$g/(1+g)$がベータ分布に従うと考えてハイパーパラメータを指定する方法を搭載。[気が狂いそうなので大幅に省略。先生方、ちょっと凝りすぎじゃないですか...]
そのほか、予測分布の出し方とか、グラフィック機能の話とか。
... 読みながらパッケージをいじってたんだけど、ちょっと意外なことに気がついた。
ためしに、変数がすごく多くてサンプルサイズがすごく小さい訓練データ(N=15, K=64)でBMAを試してみたら、過学習がすさまじい。それはしょうがないんだけど、意外にも、モデルサイズの事前分布をいじって、事後分布の平均を2くらいまで下げても、やっぱり過学習するのである。つまりモデルサイズにペナルティをかけたってダメだということだ。回帰係数のハイパーパラメータgを下げれば、さすがに過学習は起きなくなるけど... 結局、素直に小さなPLS回帰モデルをひとつつくったほうが、学習性能も般化性能もよかった。
うーん、モデル・アベレージングは過学習に強いと聞いていたのだが、なぜうまくいかないのだろうか。共線性が強すぎるからだろうか、それともMCMCのやりかたがまずいのか...
論文:データ解析(2015-) - 読了:Zeugner & Feldkirchner (2015) RのBMAパッケージ
2016年2月 6日 (土)
Hoeting, J.A., Madigan, D., Raftery, A.E., Volinsky, C.T. (1999) Bayesian Model Averating: A Tutorial. Statistical Science, 14(4), 382-417.
ベイジアン・モデル・アベレージング(BMA)についての解説。仕事の都合で必要に迫られ急遽目を通した。オフィスの本棚には解説書もあるはずなのに、そっちは手に取らず、題名に"A Tutorial"なんて書いてある奴をふらふらと読んでしまう、という...
BMAとはなにか。いきなりフォーマルに書いちゃうと気持ちが萎えますが、いまデータ$D$から関心ある量$\Delta$を予測するモデル$M_1, M_2, \ldots, M_K$があるとして、
$P(\Delta | D) = \sum_k P(\Delta|M_k, D) P(M_k|D)$
を求めましょうよ、という話である。
モデルの事後確率$P(M_k | D)$は、事前確率$P(M_k)$と尤度$P(D|M_k)$の積を、全モデルを通じて合計1になるようにしたものである。
尤度$P(D|M_k)$というのは積分尤度である。つまり、モデルのパラメータを$\theta$として、
$P(D|M_k) = \int P(D | \theta_k, M_k) P(\theta_k | M_k) d\theta_k$
である。[←ここですごく混乱したんだけど、一晩寝て読みなおしたら腑に落ちた。ここでいうモデルの事後確率$P(M_k | D)$とは、「推定されたモデル」が正しい確率ではなくて、なんていえばいいんだろう、モデルの式そのものが正しい確率というか、SASのPROC GLMのMODEL文が正しい確率というか、Rのlm()に与えるformulaが正しい確率というか、そういうものを指しているのだ、きっと。だから、尤度は推定されたモデルの持つ尤度$P(D|\theta_k, M_k)$ではなく、周辺尤度$P(D|M_k)$でなければならないのだ]
BMAには次の問題がある。
- A. モデルの数が多すぎて足し上げきれない。
- B. 積分の計算が難しい。MCMCが登場したいまでも大変。
- C. 事前確率$P(M_k)$をどうするか。
- D. どんなモデルを足し上げるのか。
問題Aについて。これには大きく2つのアプローチがある。
- データによって支持されるモデルだけを足しあげる。これをオッカムの窓という。それぞれのモデルについて$P(M_k|D)$を求め、その最大値と比較する。なにか閾値を決めといて、比が閾値を下回ったら、足し上げせずに捨ててしまう。さらに、「それをもっと単純にしたモデル」が生き残っていて、しかもそいつよりも確率が低い、という図々しいモデルも捨てる。こういうことをやってると、たいていモデルの数はひと桁にまで落ちる。[ああ、なるほど... こないだ読んだ機械学習ユーザ向け箴言集みたいな論文で、BMAなんてのは結局モデルを選んでるだけだ、という悪口が出てきたけど、そういうことか...]
- MCMCモデル合成(MC^3)。えーと、モデル空間$M$を状態空間とするマルコフ連鎖$\{M(t)\}$と、均衡分布$P(M_i|D)$を考えまして... [このへんで頭に霞がかかったようになってきたので中略]。柔軟なやり方だけど、収束するとは限らない。ちょっと似ている方法に確率探索変数選択(SSVS)というのもある。importance samplingを応用するのもある... 云々、云々。
問題Bについて。
- 尤度$P(D|M_k)$がほんとに積分で出せちゃう場合もある(線形回帰の場合とか)。
- 尤度$P(D|M_k)$ を積分なしで近似する方法もいくつかある。BICで近似しちゃうとか。[←ああ、そうだったのか... 各モデルで推定したなにかをBICで重み付けするのってどういう意味だろうかと思っていたのだが、あれは積分尤度の代用品なのか]
- [問題Bからはちょっと話が逸れているんだと思うけど] 事後分布$P(\delta|M_k, D)$のほうを、パラメータの最尤推定$\hat\theta$を使って$P(\Delta | M_k, \hat\theta, D)$で近似しちゃう、というのもある。これをMLE近似という。
問題Cに進む前に、ここで具体的なモデル・クラスを4つ挙げる。
- 例その1、線形回帰。予測子を選択せず、すべての可能な予測子セットを通じて平均する。尤度は閉形式で書けるぞ。俺たちの論文を読め(Raftery, et al., 1997 JASA)。反応変数をBox-Cox変換するとき、そのパラメータを決めずにBMAでどうにかする、という論文も書いたから読め。外れ値をどうこうするっても書いたから読め。[この辺、大幅中略。著者らはどうやらRのBMAパッケージの中の人らしいので、すいません、いずれ勉強いたします]
- 例その2、一般化線形モデル。リンク関数と予測子を変えてBMAする。モデルを$M_k$、ナルモデルを$M_0$として、ベイズファクターは周辺尤度の比
$B_{10} = P(D|M_k) / P(D|M_0)$
事前のオッズを$\alpha_k=P(M_k)/P(M_0)$として、モデルの事後確率は
$P(M_k|D) = \alpha_k B_{k0} / \sum_r \alpha_r B_{r0}$
じゃあ周辺尤度$P(D|M_k)$はどうやって出すのかというと... [なんか超めんどくさい話が続いているので、ええい、省略だ。そういう難しいお話は先生方にお任せしますよ!] - 例その3、生存分析。CoxモデルのBMAはどうやってやるかっているとねえ... [省略!! 落語の「寝床」の丁稚みたいな気分になってきた]
- 例その4、グラフィカル・モデル。[だめだ、視線が紙の上をつつつーっと滑っていく... 省略]
気を取り直して、問題C、モデルの事前確率をどうやって決めるか。
- まず思いつくのは一様分布。
- 予測子セットを動かすBMAで、各予測子について「モデルにこれが含まれる」事前確率を決めておき、積をとってモデルの事前確率にするという方法がある。George & McCullock (1993 Statist. Sinica)を見よ。[←あー、これは面白いわ! 予測子の有用性なら関係者の合議で決められそうだ。応用場面が目に浮かぶなあ...]
- グラフィカルモデルで、リンク有無に事前確率を振る。
- 専門家の合議で情報事前分布を決める。Madigan, Gavrin, & Raftery (1995 Comm.Statist.TheoryMethods), Ibrahim & Laud (1994 JASA)をみよ。
予測成績をどうやって測るか。[このくだりは次の事例1のための説明らしい。事例2では使ってない]
データを学習データ$D^B$と検証データ$D^T$に折半します。成績はこうやって測ります。検証データのすべての実現値 $d$ について
$- \sum_{d \in D^T} \log P(d | M, D^B)$
BMAそのものの予測成績は、logの右側を、それぞれのモデルにおける確率をモデルの事後確率で重み付け合計した奴に変えればよい。
[←これ、対数スコアリング・ルールじゃん。びっくり。この手法はプロパー・スコアリング・ルールになっているんだよ、なんて書いてある。ゲーム理論の文脈で出てくる考え方なんだけど、こういう場面でも使うのかー]
事例紹介をふたつ。
ひとつめ、Cox回帰の共変量を選ばずにBMAでやったという話。[難しくはなさそうだけど、面倒なので読み飛ばした]
ふたつめ、体脂肪率を回帰で予測するときにBMAでやったという話。データは$N=251$、説明変数は年齢、身長、体重など13個。全部使った重回帰で$R^2=0.75$。さて、データを折半し、学習データ側で変数選択ないしBMAをやってみた。BMAは、モデル$2^{13}=8192$個、モデルの事前確率は一様、モデルの事後確率の求め方は俺たちの論文を読め。変数選択手法としては、F値によるステップワイズ、Mallowの$C_p$最小化、調整済$R^2$最小化を試す。[←論文には「良く知られている」なんて書いてあるけど、ごめんなさい、後ろの2つがよく理解できない。$C_p$なり調整済$R^2$なりを基準に最良サブセットを悉皆検索したってこと? 8192本のモデルから? それって良いやり方なの?]
結果:BMAでは、モデルの事後確率は上位3位がそれぞれ12~14%、上位10位で計57%[←ほんとだ、BMAって事実上のモデル選択だ...]。変数選択のほうは、どの手法でも同一の8変数モデルが選ばれた。予測精度はBMAの勝ち。[←いやあ、悪評高きステップワイズ法なんかに勝っても、ねぇ...]
考察。
- 問題Dについて。モデルのクラス(回帰とか)をいろいろとりまぜたBMAというのもあるよ。Draper(1995)を読め。
- 頻度主義的なモデル・アベレージングもあるよ。AICで重みづけするとかね。ほかに機械学習系の、スタッキングとかバギングとかブースティングとか計算学習理論とか[←なにそれ]、いろいろあるけど、そういうのはたいてい点予測に焦点を当てていて、いっぽうBMAでは予測分布に焦点を当てているよ。
- BMAは「モデルのクラスが事前にわかっている」視点でのアプローチだと思われがちだけど、そんなことないよ。オッカムの窓をみてよ、「モデルのクラスが事前にわかってない」という発想に立ち、 いいモデルがみつかったらしょぼいのは捨てているよ。
- BMAの結果は複雑すぎて理解されないっていわれることあるけど、だったら事後確率だけみせればいいじゃん。p値より全然わかりやすいよ。なんならBMAはパラメータについての一種の感度分析だと割り切って、最良のモデルだけ見せればいいんじゃない?[←ああ、そういう思想なんだ...協調学習とは全然ちがうんだな]
- 今後の課題: 事前分布の選択、真のモデルが含まれていなかった時のパフォーマンス、オッカムの窓のような手法のチューニング、いろんなモデルクラスへの適用、など。
この論文には3人の研究者によるコメントと返事がついているので、そっちも一応目を通してみたのだけど、難しい話が多いわ、眠いわ、疲れたわ、途中でうっかり芋焼酎を飲み始めてしまったわで、ほとんど頭にはいらなかった。また次の機会にということで...
論文:データ解析(2015-) - 読了:Hoeting, et al. (1999) ベイジアン・モデル・アベレジーングへの招待
2016年2月 4日 (木)
Dahan, E., Mendelson, H. (2001) A Extreme-value model of concept testing. Management Science, 47(1), 102-116.
仕事の足しになるかと思ってざっと目を通した。
コンセプト・テストについての論文ってのは最近ではちょっと珍しい。えーと、ここでいうコンセプト・テストというのは、ある製品なりサービスなりを開発している途中で、その(実物じゃなくて)コンセプトについて潜在顧客がどう思うかとか、商品化した暁には買ってくれるだろうかとか、そういうのを調べることを指している。こういう市場調査の定番的課題についての議論はとっくに枯れちゃってて、情報を仕入れようにも、かえってろくなものがない。
第一著者は予測市場のマーケティング応用なんかをやっているMITのEly Dahanさん。タイトル通り、コンセプト・テストに極値分布を使ったモデルを使いましょうという話である。先生ってば、またそういう変なことを...
いわく。
新製品開発(NPD)の初期段階(ファジー・フロント・エンド)での手法としては、VOCとかリードユーザ分析とかコンジョイント分析とか狩野法とかPughのコンセプト選択とかいろいろある。コンセプト・テストは最良のデザインとか価格とかを探索する手法であると考えられている。きちんとやれば利益の期待値は上がるがコストもかかる。
コンセプト・テストの先行研究は次の3つに分けられる。
- 探索プロセスのモデル化。テストあたり費用と不確実性によって、同時にテストするコンセプトの最適な数が決まるとか(Nelson,1961;Abernathy & Rosenbloom,1968)、実験計画法を導入すればこんなにいいことがあるよとか(Thomke 1998MgmtSci, 1998Res.Policy)、同時並行的プロトタイピングは一発勝負より儲かるとか(Srinivasan et al., 1997JMR)。
- コンセプトをリアル・オプションとしてみるアプローチ。Hauser(1996 Working Paper), Hauser & Zettelmeyer(1996 Res.Tech.Mgmt.), Baldwin & Clark(2000 "Design Rules")を読め。[←実物資産を金融オプションみたいに評価することをリアル・オプションというそうだ。知らんがな、そんなん。最後のClarkって、藤本&クラークのクラークじゃないかしらん...]
- 開発手法を提案したり調べたりする研究。開発スピードを速める並列的テストとか(Smith & Reinerten, 1995書籍)、同じ失敗するのでも前向きな失敗があるとか(Leonard-Barton 1995書籍)、社内でデザインチームを競争させることの良し悪しとか(Wheelwright & Clark 1992書籍)、技術的問題への同時並行的テストが新技術によって促進されたとか(Thomke et al. 1997Res.Policy)。Dahan & Srinivasan のWeb調査手法とか、[この段階ではまだ論文になってないようだが] Dahanたちの市場メカニズム応用例 STOC なんかも、この路線である。
提案モデル。
次のように想定しよう。それぞれのコンセプトの利益は独立。それぞれのコンセプトは最良の下位コンセプトからなる[←もうrefineされ尽くしているという意味であろう]。利益の分布パラメータは企業と製品カテゴリに依存する。テストにかかるコストはコンセプトあたりで一定。テストしたいすべてのコンセプトを同時並行でテストできるけど、コンセプト数は事前に決めないといけない。利益の期待値が一番でかいやつを上市する。
コンセプト・テストを特徴づけるパラメータってのは、結局次の3つだ。コンセプトあたりのテストのコスト $c$、利益の潜在的不確実性 $b$、そして利益の分布の右の裾の太さ $\alpha$。いま、$n$個のコンセプトを同時にテストするとしよう。コンセプト $i$ の利益を$x_i$とする。
コンセプト数$n$を大きくすれば、利益の最大値$\pi_n$の期待値$E[\pi_n]$が大きくなるけど、コスト$n \cdot c$が増大する。最適なコンセプト数を$n^*$とする。話を簡単にするために、コンセプト数は連続的に動くことにする(実際には整数だけど)。
話を簡単にするために、あるコンセプトについてテストすれば、そのコンセプトから得られる利益$x$が直接にわかっちゃうってことにしよう。
$x$は確率分布$F(x)$に従う確率変数だとする。$n$個の独立したテストの結果の最大値の累積分布は$[F(x)]^n$だ。[←しれっとこう書いてあるけど、ちょ、ちょっと待って... 確率変数$X$の累積確率分布を$G(x)=P(X \leq x)$として、 独立な実現値$x_1, x_2, \ldots, x_n$の最大値が$x$以下である確率は$P(x_1 \leq x) \cdot P(x_2 \leq x) \cdots P(x_n \leq x) = [G(x)]^n$じゃないですか? ってことは、この論文では$F(x)$を確率分布と呼んでいるけど、実は累積確率分布なの?]
その確率密度関数は$n \cdot f(x) \cdot [F(x)]^{n-1}$だ。よって期待される利益は
$E[\pi_n] = n \int^\inf_{-\inf} x \cdot [F(x)]^{n-1} \cdot f(x) dx - c \cdot n$
これを解いて...[中略]... コンセプト数$n$を1増大させることによる限界利益を算出できる。
もし「どれも儲かりそうになかったらなにも上市しない」という選択肢を許すと... [略]
[さあ、ここが本題だ。だんだん読む気が失せてきたけど]
分布$F(x)$は当該コンセプトの利益の不確実性を表している。これは、過去データとかから推定される、そのカテゴリの利益の分布$H(x)$とはちがう。なぜなら、$H(x)$は一般的な利益の分布であり、そこから取り出したありうるコンセプトの巨大な下位集合の、そのまた最大値が、テストされるコンセプトだからだ。つまり、$F(x)$は$H(x)$からランダム・ドローした標本の最大値だと考えられる。$H(x)$の性質に応じて、$F(x)$はフレシェ分布、ガンベル分布、ワイブル分布になるか、あるいはそのどれにもならない。なお、「当たればでかい」製品カテゴリはフレシェ分布、当たっても上限がある製品カテゴリはワイブル分布、だいたい利益が決まっているような製品カテゴリはガンベル分布で表現される。
事例。Inhaleという会社についての実例である。[現在はNektarという社名らしい]
この会社はインシュリンみたいな薬を吸入して肺の奥深くに届けるというシステムを開発している。VCからすごいお金が流れ込んでいる。3つの開発プロジェクトが進んでいる。(A)薬を乾燥した粉にする方法の開発。(B)粉を安く作って梱包する方法の開発。(C)吸入装置の開発。それぞれ他の会社に売ることを考えている。
過去データとかシミュレーションとかで、それぞれのコンセプトの利益の確率分布関数を推定した(その方法はこの論文の本題ではない)。Aはフレシェ分布[←ワイブル分布やガンベル分布と同じく、極値分布の一種なのだそうだ]、Bはワイブル分布、Cはガンベル分布となった。ここから各カテゴリにおける最適コンセプト・テスト数が出せました。
考察。企業のみなさん、新製品の評価をする際には、利益の期待値と分散だけでなく、分布の右裾の太さを考えなさい。テストするコンセプト数をカテゴリごとに最適化しなさい。云々。[他も書いてあるけど省略。「効率の良いテスト手法を考えなさい」なんて、この論文で説教されても困る]
。。。理解できているのか怪しいものだが、要するにこういう論文だったのではないかと思う。
- (1)コンセプト・テストってのは上市後の利益を測るものだよね。
- (2)仮に適当に思いついた製品を上市したら、その製品の利益がどうなるかはわかんないけど、仮にもコンセプト・テストにかけようという製品であればきちんと考え抜かれたものであろうから、その利益の分布はきっと極値分布だよね。それがガンベル型かフレシェ型かワイブル型かは知らねども。
- (3)そのつもりで、ひとつ、コンセプト・テストをやる前に、各製品の利益の分布を極値分布として表現してみなさいな。
- (4)そしたら、何個のコンセプトについてテストすればいいかとかがわかるよ。
- (5)このモデルのおかげでね。
正直、一読して、ぽかーんとしてしまう内容であった。ぽかーん。
多くの研究がそうであるように、この論文もさまざまな想定に基づいているので、どれが非現実的か、またどれがクリティカルか、きちんと考えないといけない。
- まず、(1)コンセプト・テストってのは上市後の利益を測るものだよね、という想定。実際の製品開発では、コンセプト・テストは初期段階に行われるのが普通で、実際の製品・サービスがどうなるかはまだまだわかんないから、「上市後の利益を測っているんだ」というのは相当な違和感があると思う。でもまあ、それはコンセプト・テストになにを期待するかという話であって、上市後の利益がわかるならそれに越したことはないし、利益$x_i$がテストで直接に観察できないとしても、$E[x_i]$が測定されるのだという方向にモデルを拡張することはできる(この点は著者らも述べている)。相当に非現実的だけど、クリティカルな想定ではないと思う。
- (2)テスト対象となるコンセプトの上市後の利益は極値分布に従う、という想定はどうだろうか。うーん、コンセプト開発を「問題空間の下位空間の最大値探索」と捉えることはできるのかしらん? 実際はそうじゃないんじゃないの、と思うんだけど、これは実証的な議論が必要な話だと思う。クリティカルな想定ではあるけれど、非現実的とはいえない。
- (3)コンセプト・テストをやる前に、各コンセプトの利益の分布を表現できる、という想定。多くの人にとってはここが一番びっくりする想定ではないかと思う。「コンセプト・テストをやりたいのね、じゃあその前に、このコンセプトが上市されたときの利益がどうなるか、確率分布で表現して提出しなさい」と命じられたマーケターは、きっと目を白黒させると思うぞ。それがわかんないからテストしたいのに。でもまあ、そういう発想じたいはわからんでもないし(ベイズ推論風に言えば、背景知識によって事前分布を決めろ、といっているわけだ)、データ・リッチなこのご時世、実際にもあり得ない話じゃないと思う。
- (4)テストの最適コストやコンセプト・テスト数が知りたい、という想定。これも、ちょっとレアであるような気はする。ビジネスの規模にもよるけれど、(この論文で想定しているように)仮にコンセプト・テストによって上市後の売上が正確にわかるというならば、テストのコストなど微々たるものだ、といえる場合が多いだろう。でもまあ、非現実的とまではいえない。 ビジネスの規模が小さいとか、テストに莫大な金がかかるとか、そういう事情があるかもしれないし。
こうして整理すると、一番あっけにとられるのは、(5)このモデルでコンセプト数を決める、という点である。
なぜ? なぜここまで苦労して、コンセプト数の最適値を解析的に決めないといけないの? 解を閉形式で書くってのはそんなにすごいことなの?
だって、コンセプト・テストで上市後の利益がわかるんでしょう? その確率分布も事前にわかっているんでしょう? テストのコストを所与として、全体の利益を最大化するためには何個テストすればいいかがわかればいいんでしょう? コンセプトの数は絶対に整数で、いくら多くてもせいぜい一桁でしょう?
。。。いろんなコンセプト数についてモンテカルロ・シミュレーションすればいいじゃん!!!