統計学の年

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 01-02-2013

ところで２０１３年は国際統計学の年 International year of statistics。

Range / Inter quartile range

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 04-19-2012

ま、大したことではないんだけど、range って一番大きい値と一番小さい値の差。Interquartile range (IQR) は upper quartile と lower quartile の差。つまりどっちもひとつの数字。

このデータの range は (5, 40) でした、は厳密に言うとたぶん間違い。一番大きい値が４０で一番小さい値が５だったら range は 35。

日本語だと、range は「範囲」で interquartile range は「四分位範囲」。どちらもひとつの数字だ。

カプランマイヤー生存曲線 を描こう（Rで）

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-17-2011

Kaplan-Meier estimates をグラフにするプログラムを書いた。Rだ。

仕事上、何度となくグラフにしていて、その都度せくせくコードを書いていたのだけど、今回、昔書いたプログラムを大幅に改善した。生存曲線を書いて、下のマージンに numbers at risk を書く。

例。セットアップから。

library(survival)
kma <- survfit( Surv(time, status) ~ rx + adhere, data=colon )

この kmaの中にグラフにするべき数字が隠されている。

僕が書いた function はkmplot()という名前で ここ にあるので、

source('https://biostat.app.vumc.org/wiki/pub/Main/TatsukiRcode/RFunctions1.R')

で読み込める。

例の続き：

シンプルにkmplot(kma)で、それなりのグラフは得られるんだけど、細かにいろいろと設定することが可能。

kmplot(kma, mark='',
   xaxis.at=c(0,.5,1:9)*365, xaxis.lab=c(0,.5,1:9),
   lty.surv=c(1,2), lwd.surv=1, col.surv=c(1,1,2,2,4,4),
   col.ci=0,
   group.names=c('Obs ','Obs tumor adh','Lev','Lev tumor adh',
    'Lev+5FU ','Lev+5FU tumor adh'),
   group.order=c(5,3,1,6,4,2),
   extra.left.margin=6, label.n.at.risk=FALSE, draw.lines=TRUE,
   cex.axis=.8, xlab='Years', ylab='Survival Probability',
   grid=TRUE, lty.grid=1, lwd.grid=1, col.grid=grey(.9),
   legend=TRUE, loc.legend='bottomleft',
   cex.lab=.8, xaxs='r', bty='L', las=1, tcl=-.2
)

これでこの図が描ける。

唯一絶対必要なインプットはsurvfit()で作ったオブジェクト（この例ではkma）だけ。それ以降のインプットは・・・

まずmarkは、打ち切りデータをどう表示するか。特に気にしなくて良い。いらない場合はmark=''。

次にxaxis.at。これはn at riskを計算して表示する x軸の値（つまり時）を指定する。何もしていないと適当な所になる。xaxis.labはxaxis.atで表示するラベル。例えば、時間は日にちで指定されているけれど、年を表示したいという時は、xaxis.at=(0:5)*365, xaxis.lab=0:5とすれば、０日、３６５日、７３０日のところに、０，１，２と表示することが出来る。

lty.surv, lwd.surv, col.survは、それぞれ線の種類、太さ、色だ。グループが４つあれば、４つ指定することもできるけれど、４つ未満の場合は繰り返して使われる。Rではいつものことだ。lty.surv=1:2, lty.lwd=1, lty.col=1:4とか。

lty.ci, lwd.ci, col.ciは confidence interval を表示する時に使う。表示したくなければ、lty.ci=0あるいはcol.ci=0とする。

group.namesは別に指定しなくていいんだけど、その場合、kmaオブジェクトで使われているものをそのまま引用する。それだと、たいていの場合うざい。
group.names=c('trt1','trt2','trt3','trt4')とかに書き換えてしまえば良い。 group.orderを指定しないと、ま、そのままグラフになるのだけど、n at riskの表示と、上のグラフがだいたい同じ順序だと見やすい。そこで、表示の順番を変えられる。例えば、group.order=c(2,4,1,3)とすると２番目のグループを１番上に書いて、その次が４番目のグループ、その下に、１番目のグループ、３番目のグループとなる。ま、例を見てください。

group.namesを短く指定してもそれでも見切れる場合はextra.left.marginを指定する。label.n.at.risk=TRUEとすると"Numbers at risk"と表示する。必要ないかな。

grid=TRUEで座標を表示する。細かな設定はcol.gridなどで。

legend=TRUEにすると、グループの名前と線の色、種類なんかを表示できる。場所はloc.legendで変えられる。

その他のインプットはplot()に行くので、cex.label、btyなども指定できる。

図だけで、下の数字はいらないや、という場合はsimple=TRUEで消せる。

他にも ggplot を使って似たようなことをした人は何人かいるけれど、ggplot なんか使わないで済むならそれにこしたことはない、と思う。

団体の大きさの推定

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-02-2011

例えばナッシュビル動物園に来る人は何人グループで来るんだろうか？という疑問を持ったとする。一人で来る人もいるだろうし、恋人同士もいるだろうし、家族４人とか、家族＋親戚９人とか、学校の行事で１２０人ご一行様、とかもあるだろう。

ま、とりあえずグループの平均サイズを知りたいとしよう。理由はなんであれ、ね。

チケットを自動改札機に通す時に（本当はナッシュビル動物園にそんなハイカラなものはないのだけど）、ランダムにピピッとなるようにして、それがなった人に「あなたのグループは何人ですか？」と訊いてデータを集める。

ランダム（無作為）にサンプルを選んでいるので、このサンプルの平均を求めればきっとバイアス（偏り）の無い推定ができるに違いない。

違いない？

★☆☆★

例：とある小さな動物園に昨日やってきたグループは全てで１０グループでした。小さな動物園だし、昨日は雨だったから１０グループしか来なかったのだ。

グループのサイズは 1, 1, 2, 2, 2, 3, 4, 4, 4, 10。一人で来た写真家、雨の中のデート、雨天決行の強行ツアーグループ１０人など。

平均グループサイズは (1+1+2+2+2+3+4+4+4+10)/10 = 3.3 で良い？

もし、平均グループサイズをこう定義すると、前述の「ランダムにピピッ」というサンプルのとり方では、使えるデータは取れない。というのも、グループの中で誰かひとりでもサンプルに入ればそのグループの人数がデータとして記録されるので、大きなグループほどサンプルに入りやすくなってしまうからだ。１０人のグループはサンプルに入る確率が１人グループの１０倍ある。

この弱小動物園の例で計算してみると、この日の入場者数は 33 人だ。その一人ひとりのグループサイズを考慮すると、グループサイズ１が２人、２が６人、３が３人、４が１２人、１０が１０人ということになる。平均だと( 1 + 1 + 2+2 + 2+2 + 2+2 + 3+3+3 + 4+4+4+4 + 4+4+4+4 + 4+4+4+4 + 10+10+10+10+10+10+10+10+10+10 ) / 33 = 5.2。

平均グループサイズは 3.3? 5.2? どっち？

順番の比較　再々訪

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-23-2010

順位の比較 は一段落ついたとして，順番の比較 についてもう少し．

パパが食べた順番に 1. 中トロ，2. コハダ，3. アジ，4. アナゴ，5. 甘エビ，6. しめサバ，7. スズキ，8. ホタテ，と番号を振ると，

となった．どちらがパパの順番に近いか．

順位の比較ではなくて，順番の比較だ．こういう考え方はどうだろう？

「隣り合わせの二つの数字を入れ替える」を何回繰り返したら，パパの順番にたどり着くか．

これだと，慶悟君の順番がほぼパパの順番と同じだということも，移動する数字は１つだけ，ということで何となく表されているような気がする．逆に泰河君の順番はいろいろなところで入れ替え作業をしなくてはいけない．

慶悟君のは，まず，１と８を入れ替えて，23456718 とする．次に１と７を入れ替えて，23456178 とする．ここまでで２手順．この先，入れ替え作業を繰り返して，全部で７回の入れ替え作業で完結する．

泰河君のは，まず，１と２を入れ替えて，31247568．その後，13247568 ⇒ 12347568 ⇒12345768 ，，，と続いて（順番はどうでもいいんだけど）結局，５手で終わる．

あれ？話の展開から行って，慶悟君の方がパパに近いという結論になる指標が出てくるところだったんだけど，またしても泰河君に軍配があがった．

順番の近さを表す指標として「隣どうしを入れ替える作業の回数」はどうなのだろう？

ちょっと落ち付いて考えれば判明することなんだけど，この順番Aから順番Bまでの「入れ替え作業の回数」は，順位の比較 の時に出てきた，総当たり戦の考え方と全く同じ．お寿司の名前がついた８チーム（１位中トロ，２位コハダ．．．）の総当たり戦だとすると，２８試合のうち，慶悟君は７試合外れで泰河君は５試合外れ．順番を考慮にいれたつもりが，順番を考慮に入れない順位の比較方法の最終結論といっしょになってしまった．

▲▼

そんな訳で，パナソニックさん の「最近統計とかの話がないですね」というひとことで書き出した順番の比較だけど，未解決．実はもっともっといろいろと考えたのだけど，それはテクニカル（日本語は？）すぎでここでは書けない．

誰か，順番の比較についての論文をご存知の方は教えてください．

順位の比較

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-20-2010

長文注意．いつになく長い．

⇒順番の比較

サッカー解説者やライターの順位予想はどのくらい当たっているのか？ というページがあった．Jリーグが終わったので，シーズン前の順位予想の答え合わせをしてみよう，ということだ．

予想がどれくらい当たっていたかを判断するのに，予想の順位（１位から１８位）と実際の順位との差を全部足して１８で割る．つまり，平均でどれくらいずれているか．

この「差の平均」は，予想の正しさの指標として妥当なんだろうか？ま，１８で割って平均にするのはややめんどくさいので「差の合計」を指標として使っているとしよう．同じことだ．

優勝したのは，佐藤俊さんで，差の合計は 50，平均は 2.78．２位の清水秀彦さんと３位の野々村芳和さんが同点で，合計 52，平均 2.89．次いで前園真聖さんと渡辺達也さんが同点で，合計 54，平均 3.00．

この計算方法は，まず，１位を当てるのも１０位を当てるのも１８位を当てるのも，同じ重要性としている．もしかしたら，上位をたくさん当てた人の方が予想としては優れているのかもしれない．Jリーグだと，下位のチームがJ2に落ちるので，下の方のチームを正確に予想するのも大切かもしれない．ま，それは差を合計するときになんらかの重みをつけて計算することによって，微調整が可能だけど，今回はシンプルに考える．１位も１０位も１８位も同じ重要性だとする．

外れたときのペナルティーをどう計算するかで，他の指標も考えられる．ここで使われている差（の絶対値）がシンプルだけど，差の２乗というのも容易に考えられる．差の２乗にすると，大きく外した場合のペナルティーがかなり大きくなる．

ちょっと違う考え方をしてみよう．

例として，実際の順位と佐藤さんの予想を使う．

１位から１８位まで予想するということは，結局全てのチームの組み合わせを見て，どっちのチームの方が順位が良いか，という判断を下す必要がある．例えば，優勝した佐藤さんの予想では，１名古屋，２鹿島，３ガンバ大阪，４川崎，５東京，となっている．これで言っていることは（ガンバに注目すると）ガンバ大阪は名古屋と鹿島よりは下だけど，川崎，東京よりは上．名古屋○×ガンバ，鹿島○×ガンバ，ガンバ○×川崎，ガンバ○×東京，ということだ．

１８チームの総当り戦と考えられるので，全部で 153の直接対決があって，そのうち何組で正しく予測出来たかを指標として使えば良い．あるいは間違った予測の数をペナルティーとして足し算すればよい．

例の佐藤さんは，名古屋を正しく１位と予想したので，名古屋に関しては間違いは０．２位だったガンバを３位と予測して，鹿島を２位と予測した．つまりガンバに関しては，ガンバ×◯鹿島だけが間違いで，ペナルティーは１．実際に３位だったセレッソを佐藤さんは１１位と予想した．鹿島，川崎，清水，広島，横浜，浦和，東京，京都の８チームが佐藤さんの予想ではセレッソより上位，実際にはセレッソより下位だったので，ペナルティーは８．

これを全てのチームに関して計算して和を求める．上の表の，総当たり戦の行だ．

予想と実際の順位があまり違わない場合，単純に計算した「差」の方法とあまり違いがないけれど，この総当たり戦順位の最大の利点は意味がストレートだということだ．「順位の上下関係が間違っていたペアの数（かける２）」．
それ以外の利点として，やや引き分けになりにくいかもしれない，というのもある．

佐藤さんの点は 74 だ．例えば清水対磐田に関して，清水の点を計算した時にも磐田の点を計算した時にも数えているので，間違ったペアの数は 74 を 2で割って 37だ．153の直接対決のうち，37の予想が違った．

例えば，実際は 1 >> 2 >> 3 >> 4 >> 5 だった５チームを 5 >> 4 >> 3 >> 2 >> 1 と予想すると，この総当り戦の考え方だと，１ペアも正しく予想できていない．でも，単純に差を計算したら，3位のチームは正確に予想したことになる．それはどうだろう？

「この５チームの中では他のチームとの相対評価は分からないけれど３位だろう」という考え方だったら正解だけど，こういう絶対評価を認めるならば，全てのチームを３位と予想するとか，２位が３チームと４位が２チームという予想も認めるべきだろう．そうしないで，１６チームを１位から１６位までに予想させるなら，相対評価をしていると思われるので，偶然当たってしまっても順位関係が間違いならば，それなりのペナルティーはあるべきで，総当たり戦の考え方がふさわしい．

Jリーグの順位予想をもう一度見てみた．

そういう訳で，総当たり戦の考え方だと，優勝は前園さんだ．懐かしいな，前園さん．

でも「差の平均」が間違っているという訳ではないよ．どんな指標を使うかは主催者の自由だ．差の平均を推す理由は何だかわからないけれど．

総当り戦で勝ち負けがあっている数は『ケンドールの順位相関係数』に通ずる．どうでもいいんだけど．

歴史の問題：　古い順に並べなさい．

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-19-2010

順番の比較．．．問いの続き

以下の事柄を古い順に並べなさい．

１．関が原の戦い
２．猫山虎二郎の乱
３．大政奉還
４．応仁の乱
５．本能寺の変
６．桶狭間の戦い

◯◯

それぞれの事柄が何年に起きたか訊くのは，ただの暗記力のテストみたいになるので，歴史の流れを理解しているかどうかを見るために「古い順に並べなさい」という設問は適切だろう．

応仁の乱は室町時代で，桶狭間は信長が生きていて，本能寺で殺されて，関ヶ原は江戸幕府が始まるちょっと前で，大政奉還で江戸幕府が終わり． 4 >> 6 >> 5 >> 1 >> 3．これで，2がどこに入るかで，点があまり変わらない方が適切な気がする．

どうやって採点するんだろう？

こういう問題を出す先生たちはどうやって採点しているんだろう？

ボーリングのスコア：フレームを並び替え

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-15-2010

このあいだ（と言っても２年ほど前だけど）元同僚のレイフさんとボーリングのスコアの話をした．

スペアのあとに１本しか倒せなかったりして「あぁ，このフレームとあのフレームを入れ替えたい！」と思うことが多々ある．じゃ，フレームを入れ替えて得点を計算し直したらどうだろう，という話．フレームの入れ替え方はたくさんあるので，可能な得点もいろいろある．全て（あるいはとても多く）の入れ替えを考慮したら，倒したピンの数相応の得点の分布が判るという訳だ．

【注意】がいくつかある．

まず，ひとつひとつのフレームを並び替えることに意味があるかどうか，の問題．ボーリングがある程度うまいというのは，クローズドフレーム（スペアとストライク）をまとめて出せる，という技量も含まれているとすると，並び替えてしまうとその辺りの情報は失われる．ひとつひとつのフレームの出方がある程度ランダムという仮定をしなくてはならない．体力の消耗も無視しているな．ま，いいや．

オープンフレームばかりだったら，並び替えても点は変わらない．

１０フレームは特殊（最大で３投できるし，得点の計算もちょっと違う）なので，並び替えはできない．

【注意】おわり．

で，２年ぶりに「そういえばあのボーリングの話どうなった？」と訊かれたので，さらさらさらとボーリングのスコア計算，並び替えプログラムを書いた．例のごとくプログラムは R だ．

とりあえず例：

この 163というスコアは妥当だろうか？ストライクが続いていないので，フレームを並び替えればもっと高い点になるだろう．

９フレームの並べ替えの仕方は全部で 9! = 362,880 通り．並び替えて何点になったか計算してグラフにしてみた．163 点はピンクの所なので，かなり低い．並べ替えで出来るスコアのうち 163点より低かったのは 4.6%．そんな訳でこの人が「あぁ，ストライクが続かなかったから点が低かった！普段はもっと点が高いんだぜ」と言っていたら，ま，信じてもいいだろう．これだけピンを倒していたら最高で 190点まで行く．それにしても，面白い形の分布になったな．

180点バージョンと 190点バージョンのフレームの例：

--

--

だけど，これはあくまでフレームの並び方を換えたらどんな点になったのだろうか？という話．この１フレームから推測される，この人のボーリングスコアの分布はどんなだろうか？という話に続く．ブートストラップですぜ．

順番の比較

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-08-2010

本題は順位の比較だけど，その前置きとして，順番の比較．

泰河君（仮名）と慶悟君（仮名）とパパと３人でお寿司屋さんに行きました．３人ともにぎり寿司８個セットを食べました．パパが食べた順番は
1. 中トロ，2. コハダ，3. アジ，4. アナゴ，5. 甘エビ，6. しめサバ，7. スズキ，8. ホタテ でした．

泰河君は，アジ，コハダ，中トロ，アナゴ，スズキ，甘エビ，しめサバ，ホタテ
慶悟君は，コハダ，アジ，アナゴ，甘エビ，しめサバ，スズキ，ホタテ，中トロの順番で食べました．

さて問題です．パパが食べた順番に近いのは泰河君，慶悟君，どちらでしょうか？

「順番が近い」というのをどうやって判断しようか．

いちいち寿司ネタを書いているとお寿司を食べたくなってしまうので，数字だけにする．パパが食べた順番に，中トロ＝１，コハダ＝２，... とする．

例えば「ぴったり同じ」が多いほうが近い，と決めれば，泰河君は４（アナゴ）と８（ホタテ）がパパとぴったり同じで，慶悟君はぴったり同じは無いから，泰河君の方が近いと言える．

あるいは，最も差が大きいのに注目する，と決めると，泰河君は２，慶悟君は７（中トロが１番と８番）なので，これも泰河君に軍配があがる．

でも，どちらにしても，何となく全体を評価していないような気がするので，例えば「順位の差を全て足す」とかが良さそうだ．パパと泰河君の差は，2, 0, 2, 0, 2, 1, 1, 0 なので，差の合計は 8．パパと慶悟君の差の合計は 14．これでも泰河君の方がパパに近そうだ．

でも，少し詳しく見れば判ることだけど，実は慶悟君の順番はパパの順番とほぼ同じ．中トロを一番最初に食べたパパと一番最後に食べた慶悟君だけど，それ以外は全く同じだ．ふたりとも，コハダ→アジ→アナゴ→甘エビ→しめサバ→スズキ→ホタテ，の順番で食べている．

「順番が似ている」ことをちゃんと評価できる指標はどんなのだろうか？

おまけの問題：はずれあり．

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 04-09-2010

おまけの問題のシンプルな発展で，ハズレがあることにする．

当たりだったら６種類のおまけが同じ割合で入っているのだけど，ハズレもある．

ちょこちょこと細かいところで変わるけれど，大筋では大した違いはない．

当たりの確率を q とすると，k 個買って x 種類の確率は と書ける．前の２項は q がくっついただけで他は全く変わっていない．３項目は，ハズレだったらおまけの種類は増えないので，その時の確率．

これを計算する上では，p(0,0)=1 とする．他にも買った数がおまけの数より少ないケースとか，おまけが -1個のケースとかは確率は全て0 にして計算する．

q=1 の場合はハズレなしなので，シンプルなケースに戻る．

Rでプログラムを書いて解く． coupon.collectors <- function(m,q=1,tol=1e-5){
# m is the number of coupon types
# 1-q = P[ no coupons inside ]

out <- numeric(m+1) ; out[1] <- 1-q ; out[2] <- q
out <- t(matrix(out))
finish <- F ; i <- 2 ; tol <- max(tol, 1e-5)

while(!finish){
qr <- q * (1-(0:m)/m)
out <- rbind( out, out[(i-1),] * (1- qr) + c(0, out[(i-1),] * qr)[1:(m+1)])
finish <- out[i,m+1] > 1-tol ; i <- i + 1
}

out <- data.frame(round(out,abs(log(tol,10))-1)) ; names(out) <- 0:m
out
}

○○○

僕はスパゲッティの問題のときと同じように帰納的に解く方法を書いたのだけど，そうじゃなくて直接的に求める方法もある．いろいろインターネットをふらふらしていたら，そっちの解き方は Weblog on mebius.tokaichiba.jp に判りやすく書いてあったので，ま，いいかということになった．

○○○

おまけの問題：前置き
おまけの問題：問題
おまけの問題：答１
おまけの問題：答２
おまけの問題：答３
おまけの問題：発展問題

おまけの問題：答３

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 04-07-2010

引き続き coupon collector's problem．

問３：全部そろっている確率を 90%以上にするには幾つ買えばいいですか？

○○

平均で云々とかよりも，本当に求めたいのはたぶんこの確率だ．

解き方は問２といっしょ．「○個までにそろっている確率」が判れば解ける．

拡大

○○

おまけが 6種類の時 23個買えば 90%の確率で全種類そろっている．この先は確率がゆっくりしか上がらないので，95%にしたかったら 27個買わなくてはいけないし，99%にしたかったら 36個も買わなくてはいけない．

おまけの問題：前置き
おまけの問題：問題
おまけの問題：答１
おまけの問題：答２
おまけの問題：答３
おまけの問題：発展問題

おまけの問題：答２

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 04-06-2010

引き続き coupon collector's problem．

問２：１５個買ったら，どれくらいの確率で全部そろうか．

１５個というのは，６種類全部そろえるために買わなくてはいけない数の平均 (14.7) に一番近い数，ということ．

１５個買って，６種類全部そろっている確率．いろんな求め方がありそうだ．

ほんとにいろいろあって，ここで書くより簡単でもっとすっきりしているのもあるけれどね．個人的な好みかな．スパゲティの問題 のときと似た考え方だ．

答２：６種類全部そろっている確率だけでなくて，一般的に「m 種類のおまけがある時，k 個買ったら x 種類そろった」という答をめざそう．例えば「６個買った時にちょうど４種類ある」ためには５個買った時点でちょうど３種類あるか，４種類あるかでないといけない．５個買った時点で２種類だったり，５種類だったりしたら，どうがんばっても６個で４種類にはなり得ない．

６個買った時に４種類ある確率は
（５個買った時に３種類ある確率）×（６個目が新しい種類）と
（５個買った時に４種類ある確率）×（６個目が新しくない確率）の和だ．

だから，一歩戻って５個買った時の確率が判らないと６個買った時の確率はわからない．それで，５個買った時の確率は，４個買った時の確率が判らないと判らない．そうやって，さかのぼって行かなくてはいけないので，結局１個買った時に１種類ある確率 (100%) から始めて順番に解いていく．

k 個買った時に x 種類そろっている確率を p_k,x とする．おまけが m 種類あるとすると， と書ける．p_1,1 = 1 だ．x が 0 だったり，x が k より大きかった場合は全て 0 として計算する．

p の前についているのは，k 個目が新しい確率と，k 個目が新しくない確率．

これで，いろいろ計算出来るのだけど，手計算しようとするとかなりめんどうくさい．１５個買ったらうんぬん，という確率を求める場合でも，１個目，２個目，３個目．．．と順に計算しなくてはいけない．

それを計算した結果が，前のページの グラフだ．

１５個買って６種類全部そろっている確率だから「６個目でそろう」から「１５個目でそろう」までの確率を全部足せば１５個目までにそろう確率がわかる．拡大

１５個目までに全部そろっている確率は 64%．全種類そろえるには，平均で約１５個買えばいいのだけど，１５個買って全部そろっている確率は 50%をかなり超えている．前ページのグラフが右にぐでぇと延びているから起こる現象．左右対称だったら，確率は 50%付近になる．

全部そろっている確率を 50%以上にしたいなぁ，と思ったら１３個でよい．確率は 51%．

問３に続く．

おまけの問題：前置き
おまけの問題：問題
おまけの問題：答１
おまけの問題：答２
おまけの問題：答３
おまけの問題：発展問題

おまけの問題：答１

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 04-01-2010

問１：「一個ずつ買って，その場で開けて」という作戦でそろうまで買い続けると，平均で何個買うことになりますか？

答１：
1種類手に入れるには 1個買えば良い．
次に買うのが新たな種類，つまり 1個目と同じでない確率は 5/6．確率が 5/6 であることが起こるまでの試行の平均回数は（話すと長いんだけど）6/5．
次に買うのが新たな種類，つまり 1個目と 2個目と同じでない確率は 4/6．確率が 4/6 であることが起こるまでの試行の平均回数は（話すと長いんだけど）6/4．次に買うのが新たな種類，つまり 1個目と 2個目と 3個目と同じでない確率は 3/6．確率が 3/6 であることが起こるまでの試行の平均回数は（話すと長いんだけど）6/3．

そんな訳で 6種類そろえるために買う数の平均（期待値）は
1 + 6/5 + 6/4 + 6/3 + 6/2 + 6/1 = 14.7．

シンプルに求められるところが面白いけれど，この平均にはどんな意味があるのか．

何個買えば 6種類そろうか，というのをグラフにしてみるとわかるのだけど，右側にぐでーと延びている．こういう左右対称でないことがらの平均というのは，うさんくさいことが多い． 拡大

上のグラフは x回目で 6種類そろう，というのを表している．一番確率が高いのは 11回目で終わる．8.4% だ．（平均に一番近い）15回目で終わる確率は 6.1%．ぐでーと右側に延びているせいで平均値が右に引っ張られてしまっている．

平均で約１５個買えばいいですよ．でも１５個目でぴったり終わる確率というのは，6%でしかない．しかも，１５個目で終わる，というのが一番確率が高い訳でもない．上のグラフから判るように，９～１４個目で終わる確率のほうが高い．

何個目で６種類そろいますか？と訊かれて，答をひとつあげろと言われたら１１個目と答えるのが最善だろう．一番確率が高い．

おまけの問題：前置き
おまけの問題：問題
おまけの問題：答１
おまけの問題：答２
おまけの問題：答３
おまけの問題：発展問題

おまけの問題：問題

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 03-29-2010

限られたある種の世界ではかなり有名な coupon collector's problem．

問題文を作る上で，現実味のあるシチュエーションを作ろうとしたけれど，できなかった．現実的な疑問は持たずに，不義理な算数の問題としてとらえてください．

○○

とあるキャラメルを買うと，おまけのサイコロが入っています．そのサイコロは６種類あって，どれが入っているかはわかりません．どうやら，どのサイコロも同じ確率で（割合で）入っているらしいです．全部そろえたいんだけど何個買えばいいかな？

問１：「一個ずつ買って，その場で開けて」という作戦でそろうまで買い続けると，平均で何個買うことになりますか？

問２：とりあえず今あるお金をはたいて１５個買おう．全部そろっている確率は？

問３：買って中身を見ずに家に帰ります．全部そろっていなかったらかなりショックだから，全部そろっている確率を 90% 以上にするには幾つ買えばいいですか？

○○

この，おまけの問題 coupon collectors problem について書かれた物をインターネット上で検索してさらっと見ると，たいていが「平均で何個買えばそろいますか？」という問題だ．問１．

でも，後で判明することだけど，この平均購入数はややいかがわしい．

おまけの問題：前置き
おまけの問題：問題
おまけの問題：答１
おまけの問題：答２
おまけの問題：答３
おまけの問題：発展問題

おまけの問題：前置き

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 03-27-2010

この前書いた 鶴亀算の話 は，この話の前置き．この話自体も前置きなんだけどね．本当は クーポン収集者の苦悩： coupon collector's problem について書こうと思っていたのだ．

日本語だと「食玩問題」として知られているようだ．食玩なんて言葉使ったことないな．

ある程度有名な問題なんだけど，日本風にアレンジすると「おまけの問題」か．

とあるお菓子を買うとおまけがついてくるんだけど，そのおまけは全部で６種類あってどのおまけが入っているかはわからない．６種類全部そろえるにはいくつ買えばいいんだろう？

ひとつだけ超レアなおまけがあったりするとより現実的だけど，シンプルな問題の方がいいので，どのおまけも同じ確率で（割合で）入っているとする．（同じ確率で，と，同じ割合で，って同じ意味？ これも長い話だなぁ，きっと）

６個買って６種類を集めようとする人は自分の強運さを高く評価し過ぎだし，とりあえず１００個買おうかという人は，たぶんもう少し節約に気を使った方がいいと思う．あなたならいくつ買いますか？というのが問題だ．

○○

６種類そろえるのに必要な購入数の平均を求めればいいのか？でも，たぶん本当に求めたいのは 「これだけ買ったら６種類全部そろう確率がかなり高い」 という数だと思う．たぶん．でも「かなり高い」って？９０％かな？５０％でいいのかな？

○○

現実的な問題を書こうとすると難しい．とりあえず６個買う．その後は１個ずつ買っておまけを確かめて，そろってなかったら買い足す．という作戦をとれば，最初に何個買えばいいのか，という計算は不要だ．それでも，何個買えば（それなりの確率で）全種類そろうか計算する意味はあるか．

とりあえず１００個買っていらないのは ebay で売ればいいじゃん．とか，そもそも最初からおまけだけ ebay で買えばいいじゃん．

とか，現実味のある問題作りについて考えていたら，現実性が全くないものの代表格の鶴亀算について書きたくなったということだ．

前置き終わり．

おまけの問題：前置き
おまけの問題：問題
おまけの問題：答１
おまけの問題：答２
おまけの問題：答３
おまけの問題：発展問題

さいころいっぱい 再訪

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 02-04-2010

さいころをいっぱい振ったときの最大値の 話．再び．

例えば
１：普通の６面のサイコロを４つ振る．一番大きい目が３の確率は？
２：４面，６面，８面のサイコロを２つずつ振る．一番大きい目が５の確率は？

ちなみに答は，１が 2.7%で，２が 16.0%．

この，サイコロ最大値の問題は，計算はやや面倒くさいけれど，考え方もプログラムも極めてシンプル．

サイコロの数を n，それぞれの面の数を m として，最大値が j になる確率は

非常にシンプルなことをめんどくさく書いてみた・・・

（シンプルな）考え方：４面，６面，８面のサイコロを１つずつ振る場合を例にする．

まず最大値が１になる場合を考える．３つのサイコロとも１を出さなくてはいけないので，その確率は
P₁ = (1/4)(1/6)(1/8)
次に最大値が２になる場合：３つのサイコロとも１か２を出さなくてはいけないので，その確率は (2/4)(2/6)(2/8) ぽいのだけど，３つとも１だった場合 （確率は P₁）を除かなくてはいけないので
P₂ = (2/4)(2/6)(2/8) - P₁
最大値が３になる場合：３つとも１か２か３を出さなくていけないんだけど，P₁+P₂ を引かなくてはいけないから
P₃ = (3/4)(3/6)(3/8) - ( P₁+P₂ )

それで，上の式のような 掛け算 - 足し算 になる．最大値が j 以下になる確率は？だったら，マイナス足し算の部分がいらない．

で，min( j, m_i ) なんだけど，これはただ単に，例えば「最大値が５」の確率の時は４面のサイコロはどうがんばっても貢献できないので，(5/4) となってしまう代わりに (4/4) にする処置．例１のように目の数が同じサイコロしか出てこない場合は無視してよい．

これを計算するコードはＲなら２行で書ける．無理すれば１行で済むけれど．given.j <- function(j,v){ prod( pmin(v,j)/v) }
diff( sapply(as.list(0:max(v)), given.j, v=v) ) v が面の数．例えば例題２だったら， v <- rep(c(4,6,8),each=2) としてから走らせればよい．

歩けないのに走ろうと・・・

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 01-24-2010

この前，state of personalized medicine というタイトルの，生物統計グループのチェア（一番偉い人）の講演を聴いた．

personalized medicine （オーダーメイド医療）とは，大まかに言うと，ＤＮＡなどの患者個々の情報を踏まえてその患者に適した治療をしよう，という考え方だ．

それを達成するためには，例えば「こういう DNA 配列の患者さんにはこの治療が効く」という情報がいる．まだそこまでは行き着いていない．

統計学をやっていると，いろんなことにつきまとう誤差というものに直面する機会が多々ある．それを避けては通れない．実験室で観測した値には誤差がつくし，人間がとある治療法にどう反応するかにも誤差がある．厳密に言うと誤差（エラー）ではないのかもしれない．全てをコントロールしたら全く同じように反応するのかもしれないけれど，どこまでをコントロールをして「全て」と呼んでいいのかわからないので，「説明がつかない，コントロールしようがない」ことは誤差として計算にいれる．

personalized medicine の話を聞いていると，この「誤差」の部分の話が欠落していることが多い．全てのことが決定的 (deterministic) で，点と点を直線で結んでいる感じ．実際には無数のごちゃごちゃした点の固まり２つをまっすぐか，あるいは緩やかなカーブを描いている線でなんとなくつないでみようか，というところだ．

で，チェアの話は，ヴァンダービルトで最近 personalized medicine というフレーズをよく聞くけれど，我々は本当はまだそこまで行っていないんじゃない？という話だった．統計をやっている多くの人は personalized medicine に対して（少なくとも今現在では）懐疑的だと思う．

チェアのメッセージのひとつに「簡単なこともできていないのに，難しいことをやろうとしている」というのがあった．従来の臨床試験をより良くする手法はあるけれど，あまり知られていないし使われてもいない．そういう所に目を向けるべきなのに，注目されるのは personalized medicine だ，というのがいけない．

○○

それを聴いていた時に，ジョンレノンの歌が回った．

Children, don't do what I have done.
I couldn't walk and tried to run

というのを，チェアに言ったら，これからこの講演をする時に使おうか，と言っていた．

personalized medicine の話を聴きに行って，ジョンレノンの歌に思いを馳せた，という話だったのに，長かったなぁ．

パイパイチャート

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 08-19-2009

パイチャート．日本語では円グラフ．使い方によっては人目を引くイラストになるので，コマーシャルなどでは使い道が無い訳ではない，と思う．

でも，学術的な目的でデータ（の分布）を示したい場合は，パイチャートは例外なく不適切だ．

理由はいろいろあるのだけど，ま，とりあえず見にくい．角度の比較は長さの比較に比べて格段に難しい．もう少し根本的なところで，もともと一次元のデータ（0と1の間の数字で現せる）を円とか角度とかいう二次元のコンセプトを使って表示したり比べたりするのが無駄だ．

そんな訳なんだけど．．．

今日，素晴らしいパイチャートをみつけた．

出所はタフティーのページだ．

正しい円グラフのあり方．

ビュフォンの針

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 03-05-2009

昔確率を使って円周率を計算する 話 を書いた．これはかなりの力技だった．一般的に円周率と確率と言えばビュフォンの針だ．ビュフォンについては検索すればいろいろみつかる．

平行線がたくさんひいてある紙の上に針を何度も落とす．針が平行線（の一本）と交わる割合から円周率が求められる．

幅が$r$の平行線がいっぱい引いてあるとする．そこに長さが$r$の針を落とす．上の絵で，青い点が針の先が落ちた場所だとする．青い点から見て下にある一番近い平行線からの距離を$y$としている．平行線と垂直の縦の線（薄い灰色）と針が作る角度（鋭角）が$a$より小さい時に，針と下の線が交わる．あるいは，薄い灰色の縦の線と針が作る角度（鋭角）が$b$より小さいと，針と上の線が交わる．

そういう訳で，青い点が与えられた時に，針と上あるいは下の線が交わる確率は$(\; a+b\; )\; /\; \pi$．$\pi$って$180°$のことね．

サインコサインタンジェントを思い出すと$cos(\; a\; )\; =\; y\; /\; r$，$cos(\; b\; )\; =\; (\; r-y\; )\; /\; r$と表せることがわかる．めんどくさいから$r\; =\; 1$にしよう．$cos(\; a\; )\; =\; y$，$cos(\; b\; )\; =\; 1-y$だ．逆関数を使って$a\; =\; cos-1(\; y\; )$，$b\; =\; cos-1(\; 1-y\; )$と書ける．

つまり，青い点が与えられれば，針と平行線の一本とが交わる確率は $(\; a+b\; )\; /\; \pi \; =\; (\; cos-1(\; y\; )\; +\; cos-1(\; 1-y\; )\; )\; /\; \pi$となる．

「青い点が与えられた」という条件は積分を使って取り外してしまえる．要するにまっさらな状態での 「針と平行線が交わる確率」($p$とする)は
となる．

何回も針を落として，針と平行線が交わった割合を$p$として，そこから$\pi \; =\; 2/p$ と計算できる．

実際にやってみた．でも例のごとく，実際にやったのはコンピュータで僕ではない．

Rでやった．

$1,000,000$回試してみたところ，線と交わったのは$636,672$回．$p\; =\; 0.636672$で$2/p\; =\; 3.141335$．ま，円周率だ．

mm

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 03-04-2009

このあいだ，多施設臨床試験の結果の解析について質問されたのをきっかけに mediated moderation について勉強しようと思ったんだけど，うっかり moderated mediation について勉強していた．