HIV検査で陽性判定を受けた時に本当に陽性である確率をベン図で理解する

　Wikipedia：後天性免疫不全症候群．検査．方法によると、HIV検査（スクリーニング検査）では約0.3%に偽陽性（実際は陰性なのに陽性と判定してしまうこと）が発生するらしい。では、もしもあなたがHIV検査で陽性と判定された時に、実際に陽性である確率はどれほどだろうか。Wikipedia：後天性免疫不全症候群．疫学．世界によると、我が国におけるHIV感染者の割合は0.1%以下となっている。

　これは条件付き確率の問題である。ここで、「実際に陽性である集合」をA、「陽性判定を受けた集合」をBとし、これをベン図で表現すると次のようになる。

　求めたいのは、陽性判定をされた時に実際に陽性である条件付き確率なので、$p(A|B)$である。これは $$p(A|B)=\frac{p(A\cap B)}{p(B)}=\frac{\frac{1}{1000}\cdot \frac{997}{1000}}{\frac{1}{1000}\cdot \frac{997}{1000}+\frac{999}{1000}\cdot \frac{3}{1000}}\approx{0.25}$$
となる。つまり、HIV検査で陽性判定を受けても、実際に陽性である確率は約25%なのである。

　高いと感じますか？低いと感じますか？いずれにしても、一度陽性判定を受けたからといって絶望する必要は（あまり）ないわけです。本当は陰性（偽陽性）である確率が約75%なわけですから。だから複数の検査方法で繰り返し検査（確定診断）していくわけです。

　やんちゃなみなさん、一度陽性だったからって落ち込むことないですよ！確定診断を受けましょう！！

Mac OS X 10.8.3にOctave 3.6.xを入れようとして数時間持っていかれた話

Octaveを入れようとしたら数時間持っていかれたので、ブチ切れながらここに書き残しておく。

Octave for MacOS Xによると、MacPortsを使って入れるのがシンプルだそうなので、まずはMacPortsをインストールする。流れとしては以下のとおり。

1. XcodeをApp Storeからインストール
2. MacPortsをインストール
3. Octaveをインストール
4. ハッピー＼(^o^)／

まず、XcodeをApp Storeで検索して入れる。ここまでは普通に終わる。次に、ここからMountain Lion用のパッケージを落としてインストールするが、

sudo port install octave-devel +atlas+docs

としてもport octave-devel not foundとか言われて何も始まらない。そういえば研究所内から外にアクセスするにはプロキシをかませないと行けないんだった(・ω<)ﾃﾍﾍﾟﾛと思い出して、.profileに

export http_proxy=http://<proxy_server>:<proxy_port>　とか
export rsync_proxy=rsync://<proxy_server>:<proxy_port> 

とか加えたりするけどダメ。ここを参考に/opt/local/etc/macports/sources.confを書き換えたりもしたけど変わらない。イライラしてきたので一旦所内ネットから切り離して外部ネットに繋いでやってみた。←諦め

sudo port -d sync
sudo port list

とすんなりと動いた。 ここでいざ、

sudo port install octave-devel +atlas+docs

をやってみるが上手くいかない。"make"が無いだの"build.cmd"が無いだの言われる。調べてみるとXcodeのインストール時に"Unix Dev Tool"をチェックしてインストールする必要があるようだ。App Storeからインストールしたけど、そんなオプションチェックする所あったかな・・・と思ったが案の定無かった。
なのでXcodeのdisk imageを探してインストールしたみたけどバージョンが古くて使えなかった。Xcodeの[Preferences]-[Downloads]からCommand Line Toolsを入れればなんとかなるという情報を得てためしてみるけど、一覧にCommand Line Toolsが出てこない。おかしいなと思って調べてたらいつの間にか出るようになった。読み込みに時間がかかっていたのか、それならばプログレスバーぐらい出しておいてほしい。Command Line Toolsを入れ終わって再度、

sudo port install octave-devel +atlas+docs

としてみるが、今度は古いバージョンのXcodeが邪魔をしているようだ。そして自分がMacにおけるアンインストールについてよく分かっていないことに気付いた。そこでここを参考にアンインストールした。それで再々度、

sudo port install octave-devel +atlas+docs

何度目かのチャレンジだったが、やっぱり上手くいかない、とりあえずfailedしたことは確かだった。octave-devel(64bit)でなくてもいいのでoctave(32bit)を試してみる。

sudo port install octave

でも上手くいかない。ここで僕はMacPortsを諦めた。homebrewを使うことにする。切り替えが早いのが僕のいいところである。まずはhomebrewのインストールから。ここを見ながら進めていく、

ruby <(curl -fsSk https://raw.github.com/mxcl/homebrew/go)
brew doctor

homebrewをインストールしたら、あとはOctave for MacOS Xに沿って進めていく、

brew tap homebrew/science
brew update && brew upgrade
brew install gfortran
brew install octave

インストールにめちゃくちゃ時間(1時間弱)かかったけれども、無事に入ったようだ。
忘れずに環境変数GNUTERMにX11を設定し、

setenv ("GNUTERM", "X11")

これでやっと動いた。あとはgnuplotとか入れてxtermからｶﾀｶﾀやればグラフとかもちゃんと描ける。

brew install gnuplot

というわけで、計算機に翻弄された一日でした。おつかれ。

状態遷移図で分かる「モンティ・ホール問題」 #2

モンティ・ホール問題を状態遷移図で追っていく。まず、選択肢のうちあたりは1つ、はずれは2つであるため、あたる確率は$\frac{1}{3}$、はずれる確率は$\frac{2}{3}$である。

ここで、実際にははずれは2つあり、そのうち1つがモンティ（司会者）によって消される。

よって、結局選択肢は「はずれ」か「あたり」の2つだけになる。この状態では、選択肢を変更するか変更しないかのどちらか1つになる。

ここで、初めの状態遷移図と整理すると、以下のようになる。

こうなればもう一目瞭然である。以下に示すように、選択肢を変更しないことで当たる確率は$\frac{1}{3}$であるが…

選択肢を変更することで当たる確率は$\frac{2}{3}$となる。

選択肢を変更することで「はずれの確率」を「あたりの確率」にすることができるのである。サヴァントは、回答者が「はずれのドアを開けられる」という情報を得たことによって確率が変化する（事後確率）ことに着目し、正解を得た。

そしてここから学ぶべきことは、どのような強敵（高名な数学者）から批判を受けても、自分が正しいと信じるならばそれを曲げないという信念ではないかと、僕は思っている。なかなか出来ることではないけども。

状態遷移図で分かる「モンティ・ホール問題」 #1

モンティ・ホール問題とは確率論の問題で、アメリカのゲームショー番組を発端として一大論争が巻き起こった。その内容とは、

「プレイヤーの前に3つのドアがあって、1つのドアの後ろには景品の新車が、2つのドアの後ろにはヤギ(はずれを意味する)がいる。プレイヤーは新車のドアを当てると新車がもらえる。プレイヤーが1つのドアを選択した後、モンティが残りのドアの内ヤギがいるドアを開けてヤギを見せる。
ここでプレイヤーは最初に選んだドアを、残っている開けられていないドアに変更しても良いと言われる。プレイヤーはドアを変更すべきだろうか？」 -Wikipedia：モンティ・ホール問題より引用-

ドアを変更する場合としない場合、どちらが当たり（景品の新車）を引く確率が高いのか？という問題である。このゲームのルールを箇条書きにすると以下のようになる。-同上より引用-

1. 3つのドア (A, B, C) に（景品、ヤギ、ヤギ）がランダムに入っている。
2. プレイヤーはドアを1つ選ぶ。
3. モンティ（司会者）は残りのドアのうち1つを必ず開ける。
4. モンティの開けるドアは、必ずヤギの入っているドアである。
5. モンティはプレーヤーにドアを選びなおしてよいと必ず言う。 

直感的に考えると、ドアが2つになることによって当選確率は1/2、つまり当たりもハズレも五分五分になると思ってしまうが、果たして正解はどうなのだろうか。この問題に対し、世界で最も高いIQを有しているとされるマリリン・ボス・サヴァントはこう答えた。-同上より引用-

「正解は『ドアを変更する』である。なぜなら、ドアを変更した場合には景品を当てる確率が2倍になるからだ」

つまり、当選確率は五分五分ではないと主張したのだ。この回答に対し、数多くの数学者から「彼女は間違っている」という意見が寄せられた。しかし、彼女は自分の正しさを疑わず、あの高名な数学者ポール・エルデシュ（エルデシュ数でも有名、ちなみに僕のエルデシュ数は5だった）までもが反論する中、自分の正しさを訴え続けた。そしてコンピュータシミュレーションの結果、最終的に彼女が正しいことが明らかになった。つまり、ドアを変更した場合の当選確率は、ドアを変更しなかった場合の当選確率の2倍になるのだ。こうして論争には決着がついた。さすが天才は凡人にはできないことを平気でやってのけるものである。


前置きはここまで。次回からこの問題を状態遷移図を使って直感的に分かりやすく表現しようと思う。

続く

富士総合火力演習を観覧できる確率を考える

毎年8月に静岡県御殿場市の東富士演習場で行われる富士総合火力演習、通称『総火演』。観覧するには往復はがきまたはインターネットで応募し、恐ろしい倍率の抽選をくぐり抜ける必要がある。ちなみに応募は毎年1人1回のみである。陸上自衛隊：富士総合火力演習応募方法について

応募枠には誰でも応募できる「一般」と29歳以下限定の「青少年」の2種類の枠があり、さらにそれぞれ、車で行く人向けの「駐車券付き」と電車で行く人向けの「駐車券なし」の組み合わせ4種類の応募枠がある。各枠の倍率は、

青少年駐車券なし＜一般駐車券なし≦青少年駐車券付き＜一般駐車券付き

という大小関係となっており、「青少年駐車券なし」が最も当選確率が高い。陸上自衛隊：よくある質問によると、昨年度の倍率は平均約10倍であったらしい。つまり、「青少年駐車券なし」の倍率は10倍以下とみなせる。思ったより倍率が低い。

僕の年齢は現在27歳であり、今年・来年・再来年の3回青少年枠で応募することができる。この3回のうち$k$回当選する確率は、$n=3, p=0.1$の二項分布に従う。ただし、3回とも倍率は10倍とし変化はないものとする。3回のうち少なくとも1回当選するということは、求める確率は$1-P(k=0)$、つまり$1-$(全く当選しない確率)である。
$$1-\binom{3}{0}0.1^0(1-0.1)^{3-0}=0.271$$
よって、僕が青少年枠で総火演を観覧できる確率は約$27\%$となる。

正直言って低い。悲しい。

おそらく、30歳になっても一般枠で応募し続けるだろう。では、「一般駐車券なし」の倍率を、少し悲観して20倍とおき、当選できる確率を$90\%$以上にするためには、何年間応募し続ければいいだろうか。これも先ほどと同様にして、$1-P(k=0)\ge 0.9, p=0.05$である二項分布の$n$がいくつになるのかを考えれば良い。
$$1-\binom{n}{0}0.05^0(1-0.05)^{n-0}\ge 0.9$$
これを解くと$n=45$で当選確率が$90\%$を超える。45年間コツコツ応募し続ければ、きっと夢は叶うのである。75歳の僕には挫けずに頑張ってもらいたい。

富士山関数をMaximaで描く

場合分けされた関数をMaximaで描画するには、if then elseで関数を定義すればいいようだ。

f(x):= if 条件式 then 関数1 else 関数2

 練習のために富士山関数を描いてみる。作る関数は以下の2つ。
$$f(x)= \begin{cases} x^4-x^2+6 & (|x|\le 1)\\ \frac{12}{|x|+1} & (|x|>1) \end{cases}$$
$$g(x)=\frac{1}{2}\cos(2\pi x)+\frac{7}{2}  (|x|\le 2)$$
これを描画する命令は以下。

f(x):= if abs(x)>1 then 12/(abs(x)+1) else x^4-x^2+6;　←f(x)を作る
g(x):= if abs(x)<=2 then 1/2*cos(2*%pi*x)+7/2;　←g(x)を作る
plot2d([f(x), g(x)], [x, -8, 8]);　←適当な範囲で同一平面上に描画

これを実行すると／＾o＾＼ﾌｯｼﾞｻｰﾝ

Prediction in the dark #3

1回目：Prediction in the dark #1

2回目：Prediction in the dark #2

Prediction in the darkの3回目。前回示したアルゴリズムによって予測確率がどのように変化するのかを解析的に求めていく。$n$日目までに自分が間違えた回数を$w$とし、最優秀トレーダー$i'$が予測を間違えた回数を$w'$とする。このときの$w$と$w'$の関係がどのようになっていくのかを解析する。

$t$日目の時点での全トレーダーの重みの合計を$S(t)$とする。

$$\begin{eqnarray} S(t)=\sum_i m(i, t) \end{eqnarray}$$
自分の予測は全トレーダーによる“上がる”または“下がる”の多数決で決めるため、正答を出したグループと誤答を出したグループに別れる。
$$\begin{eqnarray} \frac{S(t)}{2}+\frac{\beta S(t)}{2} \end{eqnarray}$$
予測を間違えたグループは$t$日目よりも$t+1$日目の方が$S$の値が小さくなるので、
$$\begin{eqnarray} S(t+1)&\le& \frac{1+\beta}{2}S(t) \end{eqnarray}$$
が成り立つ。$S(1)=M$なので、自分が$w$回間違えた$n$日目は、
$$\begin{eqnarray} S(n)\le \Bigl(\frac{1+\beta}{2}\Bigr)^w M \end{eqnarray}$$

となる。最優秀トレーダー$i'$の重みは、$w'$回間違えているので、
$$\begin{eqnarray} m(i', n)=\beta^{w'} \end{eqnarray}$$
である。$n$日目での自分の重みと最優秀トレーダー$i'$の重みを整理すると、
$$\begin{eqnarray} S(n)=\sum_i m(i, n)\\ m(i', n)=\beta^{w'} \end{eqnarray}$$
式$(6), (7)$より、$m(i, t)$が全て正の数であることから、和である$S(n)$の方が大きい。すなわち、
$$\begin{eqnarray} \sum_i m(i, n) \ge m(i', n) \end{eqnarray}$$
である。ここで、式$(4), (5), (6), (8)$から、
$$\begin{eqnarray} \beta^{w'} \le \Bigl(\frac{1+\beta}{2}\Bigr)^w M \end{eqnarray}$$
という関係が導き出される。式$(9)$を$w$と$w'$について解く。
$$\begin{eqnarray} w'\log \beta &\le& w\log \Bigl(\frac{1+\beta}{2}\Bigr)+\log M\\ w&\le& \frac{w'\log\beta -\log M}{\log (\frac{1+\beta}{2})}\\ &\le& \frac{w'\log\beta}{\log (\frac{1+\beta}{2})}-\frac{\log M}{\log (\frac{1+\beta}{2})} \end{eqnarray}$$
ここで、$\frac{\log M}{\log (\frac{1+\beta}{2})}$は試行回数に依存しないため、$w$や$w'$が十分に大きければ無視出来る程度の項である。次に、$\beta \rightarrow1$である場合を考える。
$$\begin{eqnarray} \lim_{\beta \rightarrow1}\frac{\log\beta}{\log(\frac{1+\beta}{2})} \end{eqnarray}$$
式$(13)$は不定形となっているため、$\log \beta$および$\log (\frac{1+\beta}{2})$を$1$の回りでtaylor展開すると$\log \beta$は、
$$\begin{eqnarray} \log\beta&=&\beta-1-\frac{(\beta-1)^2}{2}+\frac{(\beta-1)^3}{3}\cdots\\ &=&(\beta-1)\Biggl\{1-\frac{\beta-1}{2}+\frac{(\beta-1)^2}{3}\cdots\\ \end{eqnarray}$$

となり、$\log (\frac{1+\beta}{2})$は、
$$\begin{eqnarray} \log\Bigl(\frac{1+\beta}{2}\Bigr)&=&\frac{\beta-1}{2}-\frac{(\beta-1)^2}{8}+\frac{(\beta-1)^3}{24}\cdots\\ &=&(\beta-1)\Biggl\{\frac{1}{2}-\frac{\beta-1}{8}+\frac{(\beta-1)^2}{24}\cdots \end{eqnarray}$$

となる。ここで、$(\beta-1)$が消え、$\beta\rightarrow 1$とすると以下のようになる。
$$\begin{eqnarray} \lim_{\beta \rightarrow1}\frac{\log\beta}{\log(\frac{1+\beta}{2})}=\frac{1}{\frac{1}{2}}=2 \end{eqnarray}$$

よって最終的に、
$$\begin{eqnarray} w\le 2w' \end{eqnarray}$$
が得られる。したがって、自分が予測を外す回数$w$を最優秀トレーダーが予測を外す回数$w'$の$2$倍以下にすることができる。$Q.E.D.$


…いかがでしたでしょうか。試行回数を十分に大きくする必要がありますが、この威力が分かってもらえたと思います。数学の力は凄いですね、こういうことをゼロから考え出せる人は偉大です。


終わり


2014/02/03 3:17 taylor展開が間違っていたので修正しました。

Prediction in the dark #2

1回目：Prediction in the dark #1


Prediction in the darkの2回目、今回はアルゴリズムの組み立てを追っていく。前回同様デイトレードで例えて考えてみる。

・次の日の値が“上がる”のか“下がる”のかを予測したい

モデル化するには事象を分かりやすく数値で表す必要がある。ここでは$t$日目に値が前日より“上がる”場合を$y(t)=0$、$t$日目に値が前日より“下がる”場合を$y(t)=1$とする。この場合、未来を予測するということは$t$日目に$t+1$日目を予測するということであり、すなわち$y(t+1)=?$を得るということである。

・複数のトレーダーの予測を元にする

$M$人のトレーダーがいて、それぞれが独自に次の日の値を予測する場合を考える。トレーダー$i$人目の$t$日目の予測を$p(i, t)$とする。$p(i, t)$は$0$(=上がる)または$1$(=下がる)のどちらかである。

・全トレーダーに重み（信頼度）を設定する

トレーダー$i$人目に対する$t$日目の信頼度を$m(i, t)$とする。信頼度の初期値は$1$とする。つまり、$m(i, 1)$は全ての$i$において$1$である。

・予測が外れたトレーダーの重みは下げる

トレーダー$i$の$t+1$日目の予測が実際の$t+1$日目の結果と違っていたらトレーダー$i$の重みを下げる。つまり、$p(i, t+1)\ne y(t+1)$ならば、$m(i, t+1)=\beta \cdot m(i, t)$とする。ここで$\beta$は重み下げの係数で$0<\beta <1$の値をとる。

・最終的な予測は、重み付きの多数決で決定する

“上がる”と予測したトレーダーの重みの合計と、“下がる”と予測したトレーダーの重みの合計を比較し、大きい方を採用する。

これがアルゴリズムの全てです。極めてシンプルですがその威力は驚異的です。次回はこの予測確率を解析的に求め、本当に最優秀トレーダー誤答率の2倍以内の予測が可能なのかを確かめます。

続く


2回目：Prediction in the dark #3

Prediction in the dark #1

　その昔、僕が金融工学を勉強しようとネットをウロウロしていた時にDO++：天気予報から機械学習、金融工学までというブログエントリを見つけ、Elad Hazanによる"Prediction in the dark: the multi-armed bandit problem"という機械学習手法を知った。

　当時の僕は（今もだが）この分野はあまり詳しくはなかった、しかし、なんとなく凄まじい威力を感じたので、理解しようと頑張って読み込んでいた。その頃の資料が偶然出てきたので、主に自分の理解度を深めるためにここに書き出してみる。よって上記の元記事以上の情報はおそらくないので悪しからず。

・何ができるのか？

複数の予測者の中で、最も高確率で予測できる予測者が誰なのかを知ること無く、その予測者が予測を外す回数の2倍以内で、自分の予測を出すことができる。

・何がすごいのか？

金融工学っぽく(?)、デイトレードで例えてみる。多数のトレーダーがいる中で、その最優秀トレーダーが予測を外す回数の2倍以内で、自分の予測を出すことができる。

例）100日間、$n$人のトレーダーの意見のもとに、自分の予測を出すとする。その中の最優秀トレーダーは90勝10負であった。この場合、最優秀トレーダーが誰なのかを知ることなく、自分は80勝20負以上の結果を出すことができる。ただし、$n$は十分に大きいとする。

なんとなくワクワクしてきませんか？次回から具体的なアルゴリズムと予測確率の解析を行なっていきます。

続く


2回目：Prediction in the dark #2
3回目：Prediction in the dark #3

競争はどこにでもある

　研究所内の競争的資金に応募してみた。僕はまだペーペーの新米だが、研究所に所属してる人間なら誰でも応募資格があるというのだから、応募しない手はない。うちはあまり外部資金（科研費など）の獲得に積極的じゃないのか(?)、だいたいの所属研究者が応募するらしい。個人の裁量で使える予算が得られるというのは、僕みたいな下っ端研究者からすると非常に魅力的である。

　自分が配属されている研究室のテーマとは全く関係無いテーマを応募してもいい。とは言え、どんなテーマでも自由に応募できるわけでもない。反社会的なのはもちろんダメで、防衛関係も書きようによってはアウトなようだ。一応研究所の憲章に沿った内容である必要がある。憲章はざっくりとした概念が書いてあるだけでどのようにも解釈はできるのだけれども。

　ということで、僕は憲章にももちろん沿い、将来的に研究室のテーマに貢献できるような内容を応募した。これまでの人生でも、自分を選ばせる書類の書き方というのはひと通り経験し、それなりの結果を得られたと思っている。僕のこの腕をもってして、この研究所内のツワモノ相手に勝ち残れるかが試されているわけだ。


　…通るといいな、自分の予算ほしい。