東京大学 | 大学院医学系研究科 | 国際保健学専攻 | 人類生態学教室TOP | 2nd
社会統計学第2回
「確率と確率分布」(2001年9月27日)
トップ | 更新情報 | 研究と教育 | 業績 | 計算機 | 写真 | 枕草子 | 著者 | 目安箱 | 書評 | 社会統計学目次
最終更新:
2002年1月4日 金曜日 19時14分
全部を一括して読む | 前回へ
講義概要
わかったようでわからない「確率」について,徹底的に考えてみよう
- 確率的な現象と統計的事象
- ▼どういう現象が確率的か?
- サイコロを振ったときの目:振ってみるまでは1から6のどれが出るかはわからない。どの目がでる可能性も等しいから。
- 天気予報:「明日の天気予報は晴れ」といっても「必ず晴れる」とは限らない。「曇ったり雨が降ったりする可能性も少しはあるが,晴れる可能性が高い」ことを意味する。
- 喫煙と肺がんの関係:「タバコを吸うと肺がんになる」というのは,タバコを吸った人と吸わなかった人を比べて,肺がんになった人の割合が吸った人の方で高い,という関係を示す。タバコを吸っても肺がんにならない人もいるし,吸わなくても肺がんになる人もいる。
- 実は自然界のほぼすべての現象は確率的。
- ▼こういう「不確かさ」に潜む法則性(長期間繰り返し観察したり,大集団で観察すると見られる)を考える学問を確率論と呼ぶ。この種の法則性をもつ現象を,「統計的事象」と呼び,統計的事象の確かさの度合いを示すのに便利なモノサシが「確率」
- 標本空間
- ▼統計的事象を捉えるには,「どんなことが起こりうるのか」という範囲を定めることが必要。
- ▼現象は一般に多面的で様々な観察方法がある。以下3点によって統計的現象を捉えた,記号化された結果の集合のことを「標本空間」と呼ぶ。
- 観察を行う面を特定する
- 起こりうる結果の範囲を規定する
- その範囲内の各結果に記号を対応させる
- ▼個々の結果の起こりうる可能性を示す数値(これを「確率」という)を考える。一般には「どの結果も同程度に起こる」と考える。各結果に対応付けられた確率は0から1までの数値であり,各確率の値の総和は1にならねばならない。
- ▼サイコロの目では,標本空間は{1,2,3,4,5,6}
- 事象の確率
- ▼問題は,個々の結果の可能性よりも,いくつかの結果が複合された集合(これを「事象」という)の起こる可能性がどのくらいか,ということ。つまり,事象とは,「標本空間の部分集合」である。
- ▼サイコロの例では,「目が偶数(丁)」とか「目が5以上」とか「目が1」とかいうことが事象である。
- ▼ある事象の確率は,その事象に含まれる各結果の生起確率の和である。従って,各結果の生起確率が等しい場合は,その事象に含まれる結果の場合の数をすべての場合の数で割ると,その事象の確率になる。サイコロの例では,「目が5以上」という事象の確率は,2/6=0.333・・・である。
- 余事象・和事象・積事象・排反事象
- ▼起こりうるすべての結果の集合を「全事象」という。つまり,全事象は標本空間に等しい。
- ▼決して起こらない事象を「空事象」といい,空集合φで表す。
- ▼事象Eに対して,Eが起こらないという事象をEの「余事象」という。サイコロの例では,「目が偶数」という事象の余事象は「目が奇数」である。
- ▼事象EとFの少なくとも一方が起こるという事象を,EとFの「和事象」といい,E∪Fで表す。
- ▼事象EとFの両方が起こるという事象を,EとFの「積事象」といい,E∩Fで表す。
- ▼事象Eが起こればFは決して起こらないとき,EとFは「排反事象」であるという。EとFが排反事象なら, E∩F=φである。
- 相互排反性と加法定理
- ▼サイコロで考えると,1回振ったとき「偶数の目が出る」という事象Eが起こる確率Pr(E)は,{2,4,6}の場合の数3を,{1,2,3,4,5,6}の場合の数6で割った値なのでPr(E)=0.5。
- ▼2回振って「少なくとも1回は偶数の目」の確率は?
- 0.5+0.5=1.0ではないのは自明。
- 『1回目に「偶数の目が出る」事象E1と2回目に「偶数の目が出る」事象E2とは排反ではない』ことに注意。
- 集合から考えると, Pr(E1∪E2)= Pr(E1)+ Pr(E2)-Pr(E1∩E2) =0.5+0.5-9/36=0.75であることが直感的にわかる。この式を「加法法則」と呼ぶ。
- 「2回とも奇数」の余事象ゆえ,1-9/36 =0.75と考えてもよい。
- なお,事象Eと事象Fが排反なら,Pr(E∪F)=Pr(E)+Pr(F)という「加法定理」が成立。
- 事象の独立性と乗法定理
- ▼事象Eが起こっているときに事象Fが起こる確率を,Eが起こったときのFの「条件付き確率」といい,Pr(F|E)と書く。
- ▼「事象Eが起こった」うちの「EもFも起こった」場合を考えるので,Pr(F|E)=Pr(F∩E)/Pr(E)である。
- ▼事象Eと事象Fが互いに無関係(独立)なら,Fの条件付き確率Pr(F|E)は,Pr(F)と一致する。逆にいえば,Pr(F)=Pr(F|E)のときに事象Eと事象Fは互いに独立であるという。独立でないとき「従属である」という。
- ▼上記2つの式から,事象Eと事象Fが独立なら, Pr(F∩E)=Pr(F)×Pr(E)という“乗法定理”が成立。
- 確率の4つの定義
- ▼操作的接近=統計的定義:数多く試したときの相対度数の極限。例えば,事象Eが起こる確率Pr(E)は,n回試したときにn1回事象Eが起こるとして,n1/nという相対度数が,nを無限大にしたときに漸近する値である。
- ▼対称的確率:サイコロの場合,6通りの目の出る確率はどれも等しくなければならず,その和は1でなくてはならないので,例えば1の目が出る確率は1/6となる。限定的かつ循環論法。
- ▼公理的客観確率:Pr(ei)>=0かつPr(e1)+Pr(e2)+・・・+Pr(eN)=1かつPr(E)=ΣPr(ei)を公理とする。数学的に厳密。
- ▼主観確率:観念的にも二度と繰り返すことのできない事象についての「見込み」を扱う。決定理論において重要。
- 大数の法則(操作的接近の根拠)
- ▼Rのプログラムでは,
- a <- c(100,1000,10000,100000)
- op <- par(mfrow=c(2,2))
- for (i in 1:4) { y <- runif(a[i],0,6); hist(y,6) }
- par(op)
とすれば,右図のように,試行回数を増やすとサイコロの特定の目が出る割合はある一定値に近づくことがわかる。数式で書くと,Pr(e1)=pとおけば任意の正の小さな数εに対して,lim(n→∞) Pr{|n1/n-p|<ε}=1ということで,これをベルヌーイ(Bernoulli)の大数の法則という。
- ▼同じことをExcelでやってみるには,=INT(RAND()*6)という式をセルの数だけコピーし,分析ツールの「ヒストグラム」を使えば不可能ではないが,大変面倒である。
- 確率変数と期待値
- ▼例えば,スロットマシンにコインを入れると,ごくたまに,投入金額の何十倍ものコインが出てくる。
- ▼マシン利用者全員に返ってくる賞金の合計を利用回数で割った値が,1回に期待される賞金額である。これを賭け金で割った値を「賞金還元率」と呼ぶ。すべての賭け事は胴元が儲かるようになっているのは,賞金還元率が100%未満だということである。宝くじでは40%,競馬では75%と言われる。
- ▼一般に,賞金額がx1, x2, x3, ・・・で,その賞金が得られる確率がp1, p2, p3, ・・・のように設定されたスロットマシンの期待賞金額Mは,M=x1p1+x2p2+x3p3+・・・で与えられる。
- ▼このスロットマシンのようなものを確率変数といい,期待賞金を一般に期待値と呼ぶ。
- 分散
- ▼期待賞金が同じでも,値動きの幅が小さいと一喜一憂の程度が小さく,逆に幅が大きいと滅多に当たらないが当たったときの喜びは大きくなる。つまり,ギャンブル性は,値動きの幅と,チャンスの大きさに依存している。
- ▼各賞金がどれくらい期待賞金から隔たりがあり,それを獲得できる可能性がどれくらいあるのかを見積もれば,ギャンブル性が表せる。
- ▼V=(マシンのギャンブル性)=Σ(期待値からの隔たり)×(可能性)という値が定義できる。このVを「分散」と呼ぶ。なお,各賞金額xと期待値Mの隔たりは,普通,差のを二乗した値D=(x-M)2で表す。
- 確率変数と確率分布
- ▼一般に,とりうる値の集合x=(x1,x2,x3,・・・)と,それぞれの値が実現する確率p=(p1,p2,p3,・・・)が与えられていて,事象としてxのうちどれか1つの値のみ実現するとき,(x,p)という1セットを,「確率変数」と呼んで,Xで表す。
- ▼期待値はE(X)=μ=Σxipi
- ▼分散はV(X)=σ2=Σ(xi-μ)2pi
- ▼分散の平方根σを標準偏差と呼ぶ。
- ▼横軸にxの各々の値を示す位置に, pの各々の可能性を示す高さの棒を立ててみれば,これが確率変数の「確率分布」ということになる。
- ベルヌーイ試行と2項分布
- ▼1回の実験でSかFかのどちらかが起こり,しかもそれらが起こる可能性が,Pr(S)=p,Pr(F)=1-p=qで何回実験しても変わらないとき,これを「ベルヌーイ試行」という。
- ▼ベルヌーイ試行をn回行って,Sがちょうどk回起こる確率は,Pr(X=k)=nCkpkqn-k
- ▼nCkは2項係数と呼ばれる。このような確率変数Xは,「2項分布に従う」といい,B(n,p)と表す。E(X)=np,V(X)=npqである。
- 正規分布
- ▼nが非常に大きい場合には,2項分布B(n,p)の確率Pr(X=np+d) という値が,1/√(2πnpq)・exp(-d2/(2npq))という値で近似できる。
- ▼一般にこの極限である,Pr(X=x)=1/√(2πσ2)・exp(-(x-μ)2/(2σ2))という形をもつ確率分布を正規分布と呼び,N(μ, σ2)と書く。
- ▼z=(x-μ)/σと置けば,Pr(Z=z)=1/√(2π)・exp(-z2/2)となる。これを標準正規分布と呼び,N(0,1)と書く。
- ▼統計学でよく使われる確率分布であるカイ二乗分布とかt分布とかF分布は,正規分布から導かれる。
- 練習問題
- (1)8頭で出走する競馬のレースがあり,「どの馬が勝つチャンスも等しい」と仮定した場合,ある特定の馬が勝つと予想して当たる確率は1/8となるが,2回のレースの少なくともどちらか一方に当たる確率はいくらか?
- (2)p=0.2として,いろいろなnについて2項分布のグラフはどういう形になるか?
フォロー
- 練習問題(1)の答え
- ▼確率1/8の独立な2つの事象の少なくともどちらかが起こる確率は,それぞれが起こる確率の和から両方とも起こる確率を引いたものなので,1/8+1/8-1/64=15/64
- ▼あるいは,「どちらも起こらない」ことはない,という確率として,1-(7/8*7/8)=15/64
- 2項分布についての補足説明
- ▼正二十面体サイコロをn回(n=4, 10, 20, 50)投げて,1から4までの目が出る回数と考える(2001年12月27日訂正:R-jpメーリングリストでご指摘いただいたが,考えてみれば正五面体というのは存在しないのだった)。1回投げたときに1から4までの目が出る確率は0.2であるとして,試行1000セットの平均値を右図に示す。Rのプログラムは下記の通り。
- n <- c(4,10,20,50)
- op <- par(mfrow=c(2,2))
- for (i in 1:4) {
- nx <- n[i]
- jj <- c(1:(nx+1))
- na <- jj
- for (j in jj) {na[j] <- 0}
- for (k in 1:1000) {chk <- 1
- for (j in 1:nx) {
- if (runif(1,0.0,5.0)<1) {chk <- chk +1}}
- na[chk] <- na[chk]+0.001}
- jj <- jj-1
- plot(jj,na)}
- par(op)
- ▼各nについての理論的な確率分布は,Pr(X=k)=nCk0.2k0.8n-kより右図のようになる。Rのプログラムは下記の通り。
- combi <- function(n, k) {
- x <- 1
- for (i in (n-k+1):n) x <- x*i/(n+1-i)
- if (k<=0) x <- 1
- x}
- n <- c(4,10,20,50)
- jj <- c(1:4)
- op <- par(mfrow=c(2,2))
- for (j in jj) {
- xa <- 1:(n[j]+1)
- jjj <- xa-1
- kk <- 0
- while (kk <= n[j]) { kx <- kk+1;
- xa[kx] <- combi(n[j],kk)*(0.2^kk)*(0.8^(n[j]-kk));
- kk <- kk+1 }
- plot(jjj,xa)}
- par(op)
- (2001年10月17日追記)実は,Rには,組み合わせの関数は組み込まれていた。ここではcombi(n,k)という形で定義したが,もとからchoose(n,k)がその機能をもっていた。また,jjjに0からn[j]までの値を代入してplotするよりも,barplot(xa)とする方が美しい度数分布図が書けることがわかった。
- (2001年12月27日追記)R-jpメーリングリストで,もっと素直でエレガントなプログラムのやり方を教えていただいた。中林@富士通総研さん,森@学芸大さん,佐藤@統数研さんありがとうございます。
- (2002年1月4日追記)群馬大学の青木繁伸先生からも,applyを使って美しくまとめる方法をお教えいただき,一様分布のシミュレーションと尤度を最大にする母比率の実験のプログラムも教えていただいた。ここに掲載することをお許しくださったので,以下に示す。
- 一様分布のシミュレーション
op <- par(mfrow=c(2,2))
x <- apply(matrix(10^(2:5), ncol=1), 1, function(a, n) {y <- runif(a, 0, n); hist(y, n)}, 10)
par(op)
●apply の最後の引数が上限を表す。すなわち,[0,n] の一様分布のシミュレーションになります。
- 二項分布のシミュレーション
draw <- function(n, p, trial=1000){ x <- apply(array(n, trial), 1, function(n) {res <- runif(n, 0, 1) < p; length(res[res])}); hist(x, freq=F, br=seq(-0.5,n+0.5,1)) }
op <- par(mfrow=c(2,2))
x <- apply(matrix(c(4,10,20,50),ncol=1), 1, draw, 0.2, trial=1000)
par(op)
●apply の最後から二番目の引数が母比率,trial は試行回数。
- 二項分布の密度関数
draw <- function(n, p) {x <- 0:n; xa <- choose(n,x)*(p^x)*((1-p)^(n-x)); plot(x,xa)}
op <- par(mfrow=c(2,2))
x <- apply(matrix(c(4,10,20,50), ncol=1), 1, draw, 0.2)
par(op)
●apply の最後の引数は母比率。
- 尤度を最大にする母比率の実験(第5回参照)
bis <- function(n, r, p) { choose(n, r)*p^r*(1-p)^(n-r)}
n <- 20
r <- 2
p <- 0:500/500
y <- bis(n, r, p)
plot(p, y, type="l")
y
全部を一括して読む | 次回へ