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

ところで2013年は国際統計学の年 International year of statistics

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|統計学コメント(0)

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

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

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 木曜日|統計学コメント(0)

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('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/TatsukiRcode/TatsukiRcodeKMplot.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.labxaxis.atで表示するラベル。例えば、時間は日にちで指定されているけれど、年を表示したいという時は、xaxis.at=(0:5)*365, xaxis.lab=0:5とすれば、0日、365日、730日のところに、0,1,2と表示することが出来る。

lty.surv, lwd.surv, col.survは、それぞれ線の種類、太さ、色だ。グループが4つあれば、4つ指定することもできるけれど、4つ未満の場合は繰り返して使われる。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)とすると2番目のグループを1番上に書いて、その次が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.labelbtyなども指定できる。

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 土曜日|統計学コメント(2)

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

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

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

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

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

違いない?

★☆☆★

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

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

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

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

この弱小動物園の例で計算してみると、この日の入場者数は 33 人だ。その一人ひとりのグループサイズを考慮すると、グループサイズ1が2人、2が6人、3が3人、4が12人、10が10人ということになる。平均だと( 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? どっち?

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 金曜日|統計学コメント(2)

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

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

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

泰河君32147568
慶悟君23456781

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

順位の比較ではなくて,順番の比較だ.こういう考え方はどうだろう?

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

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

慶悟君のは,まず,1と8を入れ替えて,23456718 とする.次に1と7を入れ替えて,23456178 とする.ここまでで2手順.この先,入れ替え作業を繰り返して,全部で7回の入れ替え作業で完結する.

泰河君のは,まず,1と2を入れ替えて,31247568.その後,13247568 ⇒ 12347568 ⇒12345768 ,,,と続いて(順番はどうでもいいんだけど)結局,5手で終わる.

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

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

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

▲▼

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 木曜日|統計学コメント(2)

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

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

順番の比較

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

予想がどれくらい当たっていたかを判断するのに,予想の順位(1位から18位)と実際の順位との差を全部足して18で割る.つまり,平均でどれくらいずれているか.

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

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

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

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

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

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

チーム実際の順位佐藤さんの予想総当たり戦
名古屋1100
ガンバ2311
セレッソ31188
鹿島4222
川崎5411
清水6602
広島7924
横浜8713
新潟91677
浦和10824
磐田111766
大宮121204
山形131304
仙台141515
神戸151415
東京1651111
京都171077
湘南181800
-- ------5074

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

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

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

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

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

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

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

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

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

名前差の合計順位総当り順位
佐藤俊 501744
清水秀彦 522744
野々村芳和522662
前園真聖 544621
渡辺達也 544766
猪狩真一 566723
前田秀樹 5877810
早野宏史 5877810
中西哲生 587766
熊崎敬587766
G. エンゲルス60117810
木村浩吉 60118013
加部究 60118414
西部謙司 6011766
福西崇史 62158615
後藤健生 62158816
金田喜稔 66179619
山内雄司 66179417
浅田真樹 66179417
木崎伸也 782010220

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

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 月曜日|統計学コメント(2)

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

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

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

1.関が原の戦い
2.猫山虎二郎の乱
3.大政奉還
4.応仁の乱
5.本能寺の変
6.桶狭間の戦い

◯◯

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

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

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 日曜日|統計学コメント(0)

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

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

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

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

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

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

10フレームは特殊(最大で3投できるし,得点の計算もちょっと違う)なので,並び替えはできない.

【注意】おわり.

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

とりあえず例:

9/3/x809/x71x3/8/x
133351597997105125143163

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

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

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

9/3/xxx809s713/8/x
13336391109117134142160180

--

713/9/xxx9/803/8/x
8274777106126144152170190

--

713/803/9/9/xxx8/x
82634537292122150170190

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|統計学コメント(2)

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

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

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

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

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

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

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

泰河君32147568
慶悟君23456781

例えば「ぴったり同じ」が多いほうが近い,と決めれば,泰河君は4(アナゴ)と8(ホタテ)がパパとぴったり同じで,慶悟君はぴったり同じは無いから,泰河君の方が近いと言える.

あるいは,最も差が大きいのに注目する,と決めると,泰河君は2,慶悟君は7(中トロが1番と8番)なので,これも泰河君に軍配があがる.

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

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|統計学コメント(7)

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

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

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

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

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

これを計算する上では,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 に判りやすく書いてあったので,ま,いいかということになった.

○○○

おまけの問題:前置き
おまけの問題:問題
おまけの問題:答1
おまけの問題:答2
おまけの問題:答3
おまけの問題:発展問題

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 金曜日|統計学コメント(0)

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

引き続き coupon collector's problem.

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

○○

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

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

拡大

○○

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

おまけの問題:前置き
おまけの問題:問題
おまけの問題:答1
おまけの問題:答2
おまけの問題:答3
おまけの問題:発展問題

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|統計学コメント(2)

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

引き続き coupon collector's problem.

問2:15個買ったら,どれくらいの確率で全部そろうか.

15個というのは,6種類全部そろえるために買わなくてはいけない数の平均 (14.7) に一番近い数,ということ.

15個買って,6種類全部そろっている確率.いろんな求め方がありそうだ.

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

答2:6種類全部そろっている確率だけでなくて,一般的に「m 種類のおまけがある時,k 個買ったら x 種類そろった」という答をめざそう.例えば「6個買った時にちょうど4種類ある」ためには5個買った時点でちょうど3種類あるか,4種類あるかでないといけない.5個買った時点で2種類だったり,5種類だったりしたら,どうがんばっても6個で4種類にはなり得ない.

6個買った時に4種類ある確率は
(5個買った時に3種類ある確率)×(6個目が新しい種類)と
(5個買った時に4種類ある確率)×(6個目が新しくない確率)の和だ.

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

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

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

これで,いろいろ計算出来るのだけど,手計算しようとするとかなりめんどうくさい.15個買ったらうんぬん,という確率を求める場合でも,1個目,2個目,3個目...と順に計算しなくてはいけない.

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

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

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

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

問3に続く.

おまけの問題:前置き
おまけの問題:問題
おまけの問題:答1
おまけの問題:答2
おまけの問題:答3
おまけの問題:発展問題

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 火曜日|統計学コメント(0)

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

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

答1:
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%.ぐでーと右側に延びているせいで平均値が右に引っ張られてしまっている.

平均で約15個買えばいいですよ.でも15個目でぴったり終わる確率というのは,6%でしかない.しかも,15個目で終わる,というのが一番確率が高い訳でもない.上のグラフから判るように,9~14個目で終わる確率のほうが高い.

何個目で6種類そろいますか?と訊かれて,答をひとつあげろと言われたら11個目と答えるのが最善だろう.一番確率が高い.

おまけの問題:前置き
おまけの問題:問題
おまけの問題:答1
おまけの問題:答2
おまけの問題:答3
おまけの問題:発展問題

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 木曜日|統計学コメント(0)

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

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

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

○○

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

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

問2:とりあえず今あるお金をはたいて15個買おう.全部そろっている確率は?

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

○○

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

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

おまけの問題:前置き
おまけの問題:問題
おまけの問題:答1
おまけの問題:答2
おまけの問題:答3
おまけの問題:発展問題

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 月曜日|統計学コメント(0)

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

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

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

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

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

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

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

○○

6種類そろえるのに必要な購入数の平均を求めればいいのか?でも,たぶん本当に求めたいのは 「これだけ買ったら6種類全部そろう確率がかなり高い」 という数だと思う.たぶん.でも「かなり高い」って?90%かな?50%でいいのかな?

○○

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

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

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

前置き終わり.

おまけの問題:前置き
おまけの問題:問題
おまけの問題:答1
おまけの問題:答2
おまけの問題:答3
おまけの問題:発展問題

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 土曜日|統計学コメント(0)

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

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

例えば
1:普通の6面のサイコロを4つ振る.一番大きい目が3の確率は?
2:4面,6面,8面のサイコロを2つずつ振る.一番大きい目が5の確率は?

ちなみに答は,1が 2.7%で,2が 16.0%.

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

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

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

(シンプルな)考え方:4面,6面,8面のサイコロを1つずつ振る場合を例にする.

まず最大値が1になる場合を考える.3つのサイコロとも1を出さなくてはいけないので,その確率は
P1 = (1/4)(1/6)(1/8)
次に最大値が2になる場合:3つのサイコロとも1か2を出さなくてはいけないので,その確率は (2/4)(2/6)(2/8) ぽいのだけど,3つとも1だった場合 (確率は P1)を除かなくてはいけないので
P2 = (2/4)(2/6)(2/8) - P1
最大値が3になる場合:3つとも1か2か3を出さなくていけないんだけど,P1+P2 を引かなくてはいけないから
P3 = (3/4)(3/6)(3/8) - ( P1+P2 )

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

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 木曜日|統計学コメント(0)

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

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

personalized medicine (オーダーメイド医療)とは,大まかに言うと,DNAなどの患者個々の情報を踏まえてその患者に適した治療をしよう,という考え方だ.

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

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

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

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

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

○○

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

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

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 日曜日|統計学コメント(2)

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

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

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

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

そんな訳なんだけど...

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

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

正しい円グラフのあり方

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|統計学コメント(0)

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

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

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

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

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

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

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

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

何回も針を落として,針と平行線が交わった割合をpとして,そこからπ = 2/p と計算できる.

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

Rでやった.

1,000,000回試してみたところ,線と交わったのは636,672回.p = 0.6366722/p = 3.141335.ま,円周率だ.

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 木曜日|統計学コメント(0)

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

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

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|統計学コメント(0)

ブログ検索

カレンダー

JanFebMar AprMay JunJulAugSepOctNovDec 2005200620072008200920102011
- - - - - - 1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 - - - - - -

最近のコメント

履歴

2017年04月 (1)

2017年03月 (1)

2017年02月 (1)

2017年01月 (1)

2016年12月 (1)

2016年11月 (1)

2016年09月 (1)

2016年08月 (1)

2016年07月 (1)

2016年05月 (1)

2016年04月 (1)

2016年03月 (1)

2016年02月 (1)

2016年01月 (1)

2015年12月 (1)

2015年07月 (1)

2015年06月 (2)

2015年05月 (2)

2015年04月 (3)

2015年03月 (2)

2015年02月 (1)

2015年01月 (3)

2014年12月 (1)

2014年11月 (1)

2014年10月 (2)

2014年09月 (2)

2014年08月 (2)

2014年07月 (1)

2014年06月 (1)

2014年05月 (2)

2014年04月 (1)

2014年03月 (2)

2014年02月 (1)

2014年01月 (2)

2013年12月 (1)

2013年11月 (3)

2013年10月 (2)

2013年09月 (1)

2013年08月 (4)

2013年07月 (1)

2013年06月 (2)

2013年05月 (2)

2013年04月 (3)

2013年03月 (1)

2013年02月 (1)

2013年01月 (4)

2012年12月 (1)

2012年11月 (3)

2012年10月 (1)

2012年09月 (1)

2012年08月 (3)

2012年07月 (3)

2012年06月 (2)

2012年05月 (6)

2012年04月 (2)

2012年03月 (8)

2012年02月 (2)

2012年01月 (1)

2011年12月 (6)

2011年11月 (5)

2011年10月 (4)

2011年09月 (6)

2011年08月 (9)

2011年07月 (5)

2011年06月 (5)

2011年05月 (5)

2011年04月 (6)

2011年03月 (17)

2011年02月 (6)

2011年01月 (10)

2010年12月 (10)

2010年11月 (4)

2010年10月 (6)

2010年09月 (5)

2010年08月 (11)

2010年07月 (8)

2010年06月 (8)

2010年05月 (3)

2010年04月 (8)

2010年03月 (11)

2010年02月 (4)

2010年01月 (8)

2009年12月 (6)

2009年11月 (6)

2009年10月 (6)

2009年09月 (7)

2009年08月 (6)

2009年07月 (10)

2009年06月 (10)

2009年05月 (10)

2009年04月 (6)

2009年03月 (7)

2009年02月 (9)

2009年01月 (12)

2008年12月 (6)

2008年11月 (10)

2008年10月 (8)

2008年09月 (9)

2008年08月 (12)

2008年07月 (8)

2008年06月 (12)

2008年05月 (12)

2008年04月 (12)

2008年03月 (11)

2008年02月 (10)

2008年01月 (10)

2007年12月 (12)

2007年11月 (14)

2007年10月 (13)

2007年09月 (11)

2007年08月 (16)

2007年07月 (10)

2007年06月 (10)

2007年05月 (6)

2007年04月 (10)

2007年03月 (13)

2007年02月 (10)

2007年01月 (8)

2006年12月 (13)

2006年11月 (15)

2006年10月 (9)

2006年09月 (8)

2006年08月 (18)

2006年07月 (14)

2006年06月 (16)

2006年05月 (23)

2006年04月 (20)

2006年03月 (12)

2006年02月 (14)

2006年01月 (21)

2005年12月 (20)

2005年11月 (17)

2005年10月 (18)

2005年09月 (16)

2005年08月 (10)

2005年07月 (6)