default eye-catch image.

標本調査に必要なサンプル数の下限を与える2次関数

[mathjax] 2項分布に従う母集団の母平均を推測するために有意水準を設定して95%信頼区間を求めてみた。 母平均のあたりがついていない状況だとやりにくい。 [clink url=\"https://ikuty.com/2019/01/11/sampling/\"] (hat{p})がどんな値であっても下限は(hat{p})の関数で抑えられると思ったので、 気になって(hat{p})を変数のまま残すとどうなるかやってみた。 begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 \\ frac{1.96}{0.05}sqrt{hat{p}(1-hat{p})} le sqrt{n} \\ 39.2^2 hat{p}(1-hat{p}) le n end{eqnarray} 左辺を(f(hat{p}))と置くと (f(hat{p}))は下に凸の2次関数であって、 (frac{d}{dhat{p}}f(hat{p})=0)の時に最大となる。というか(hat{p}=0.5)。 (hat{p}=0.5)であるとすると、これはアンケートを取るときのサンプル数を求める式と同じで、 非常に有名な以下の定数が出てくる。 begin{eqnarray} 1537 * 0.5 (1-0.5) le n \\ 384 le n end{eqnarray} (hat{p})がどんな値であっても、サンプル数を400とれば、 有意水準=5%の95%信頼区間を得られる。 だから、アンケートの(n)数はだいたい400で、となる。 さらに、有意水準を10%にとれば、(n)の下限は100で抑えられる。 なるはやのアンケートなら100、ちゃんとやるには400、というやつがこれ。

default eye-catch image.

標本調査に必要なサンプル数を素人が求めてみる。

[mathjax] ちょっと不思議な計算をしてみる。 仮定に仮定を積み重ねた素人の統計。 成功か失敗かを応答する認証装置があったとする。 1回の試行における成功確率(p)は試行によらず一定でありベルヌーイ試行である。 (n)回の独立な試行を繰り返したとき、成功数(k)を確率変数とする離散確率変数に従う。 二項分布の確率密度関数は以下の通り。 begin{eqnarray} P(X=k)= {}_n C_k p^k (1-p)^{n-k} end{eqnarray} 期待値、分散は、 begin{eqnarray} E(X) &=& np \\ V(X) &=& np(1-p) end{eqnarray} (z)得点(偏差値,つまり平均からの誤差が標準偏差何個分か?)は、 begin{eqnarray} z &=& frac{X-E(X)}{sigma} \\ &=& frac{X-E(X)}{sqrt{V(X)}} \\ &=& frac{X-np}{sqrt{np(1-p)}} end{eqnarray} であり、(z)は標準正規分布に従う。 これを標本比率(hat{p}=frac{X}{n})を使うように式変形する。 begin{eqnarray} z &=& frac{frac{1}{n}}{frac{1}{n}} frac{X-np}{sqrt{np(1-p)}} \\ &=& frac{frac{X}{n}-p}{sqrt{frac{p(1-p)}{n}}} \\ &=& frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} end{eqnarray} (n)が十分に大きいとき、(z)は標準正規分布(N(0,1))に従う。 従って、(Z)の95%信頼区間は以下である。 begin{eqnarray} -1.96 le Z le 1.96 end{eqnarray} なので、 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} le 1.96 end{eqnarray} (hat{p})は(p)の一致推定量であるから、(n)が大なるとき(hat{p}=p)とすることができる。 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{hat{p}(1-hat{p})}{n}}} le 1.96 \\ end{eqnarray} (p)について解くと(p)の95%信頼区間が求まる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} le p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} end{eqnarray} 上記のにおいて、標準誤差(1.96sqrt{frac{hat{p}(1-hat{p})}{n}})が小さければ小さいほど、 95%信頼区間の幅が狭くなる。この幅が5%以内であることを言うためには以下である必要がある。 (有意水準=5%) begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 end{eqnarray} 観測された(hat{p})が(0.9)であったとして(n)について解くと、 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.05 \\ frac{1.96}{0.05} sqrt{0.09} le sqrt{n} \\ 11.76 le sqrt{n} \\ 138.2 le n end{eqnarray} 139回試行すれば、100回中95回は(p)は以下の95%信頼区間に収まる。 つまり95%信頼区間は以下となる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} &le& p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} \\ 0.9-1.96 frac{sqrt{0.09}}{sqrt{139}} &le& p le 0.9 + 1.96 frac{sqrt{0.09}}{sqrt{139}} \\ 0.9-1.96 frac{0.3}{11.78} &le& p le 0.9+1.96 frac{0.3}{11.78} \\ 0.85 &le& p le 0.95 end{eqnarray} (n)を下げたい場合は有意水準を下げれば良い。 統計的に有意水準=10%まで許容されることがある。 有意水準が10%であるとすると、(n)は35以上であれば良いことになる。 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.1 \\ frac{1.96}{0.1} sqrt{0.09} le sqrt{n} \\ 5.88 le sqrt{n} \\ 34.6 le n end{eqnarray} 信頼区間と有意水準の式において(p)を標本から取ってきたけど、 アンケートにおいてYes/Noを答える場合、(p)は標本における最大値(つまり0.5)を 設定して(n)を求める。 つまり、(p)として利用するのは標本比率ではないのかな?と。 このあたり、(hat{p})を変数として残すとどういうことがわかった。 [clink url=\"https://ikuty.com/2019/01/13/sampling_with2/\"]

default eye-catch image.

幾何分布と「過去の結果からは何もわからない話」

[mathjax] いつか起こる大地震がもし昨日起きたとしたら明日起こる確率は下がるのか? 飛行機が昨日落ちたとしたらしばらくは飛行機は落ちないのか? うまくいかない人生が今日またうまくいかなかったとして将来うまくいかない確率は下がるのか? 起こるか起こらないかの確率が変わらないのであれば、将来は過去に影響されないらしい。 影響されるかどうか、と聞かれるとされない、と答えるだろうけど、それを説明することができるらしい。 そんな無記憶性の件読んでみたのでまとめてみる。 統計学入門とこちらを参考にさせて頂きました。 幾何分布 幾何分布 二項分布、ポアソン分布は(n)回の試行のうち(x)回事象(A)が発生したときを話題にしていた。 前提を少し変えて、予め試行する回数を決めないで、 (x)回目の試行で初めて事象(A)が発生した、 という話をすることもできる様子。 確率変数(x)が等間隔に並ぶ時刻であるとすることで、 事象(A)が発生するまでの待ち時間に関する確率分布を作れる。 (x)回目の試行で初めて事象(A)が起こった、ということを、 (x-1)回事象(bar{A})が起こり、次に事象(A)が起きたと考える。 事象(A)の生起確率を(p)、事象(bar{A})の生起確率を(q)とすると、 その確率は以下のようになる。 begin{eqnarray} f(x) = p cdot q^{x-1} end{eqnarray} (x)の増加1回に対して(q)を1回かける構造で、公比(q)の等比数列になっている。 等比数列って英語でgeometric seriesって言うもんだから、幾何分布っていう名前がついてる様子。 (f(x))の確率変数(x)は事象(A)が発生するまでの試行回数(時間)である、と読む様子。 では確率変数(X)が幾何分布に従うとき、期待値はどうかというと(frac{1}{p})となる。 確率の以下の読み方から、(x)回の試行に平均して(frac{1}{p})回かかる、というのは妥当に見える。 事象(A)の生起確率(p)の読み方について、その逆数(frac{1}{p})は、事象(A)が起こるまでの試行回数と読む。 (p)回の試行で初めて事象(A)が起こった、というシーンで事象(A)の生起確率を(p)と置いているので..。 期待値 幾何分布の期待値は(frac{1}{p})。ベルヌーイ試行の確率が(p)であるならば、平均(frac{1}{p})回で事象が起こる。 分散は(frac{1-p}{p^2})。期待値の証明は以下の通り。へぇ。 begin{eqnarray} E(X) &=& sum_x x cdot f(x) \\ &=& sum_x x cdot p q^{x-1} end{eqnarray} (E(X)がfrac{1}{p})であることを示したい。 恒等式を使ってやる奴ではなく、愚直にやる奴を書く。 まず、(q=1-p)として(E(X))を変形しておいてスタート。 begin{eqnarray} E(X) &=& sum_x x cdot p (1-p)^{x-1} \\ &=& p sum_x x cdot (1-p)^{x-1} end{eqnarray} 右辺を生み出すために、(frac{1}{x-1})のテイラー展開を持ち出す。 begin{eqnarray} frac{1}{1-x} &=& 1 + x + x^2 + cdots \\ &=& sum_{k=0}^{infty} x^k end{eqnarray} 左辺を(x)で微分すると以下の通り。 begin{eqnarray} left( frac{1}{1-x} right) frac{d}{dx} = frac{1}{(1-x)^2} end{eqnarray} 右辺を(x)で微分すると以下の通り。 begin{eqnarray} sum_{k=0}^{infty} x^k frac{d}{dx} = sum_{k=1}^{infty} k x^{k-1} end{eqnarray} なので、 begin{eqnarray} frac{1}{(1-x)^2} = sum_{k=1}^{infty} k x^{k-1} end{eqnarray} (x=1-p)として式変形すると、 begin{eqnarray} frac{1}{p^2} = sum_{k=1}^{infty} k (1-p)^{k-1} end{eqnarray} (E(X))にこれらを代入すると、 begin{eqnarray} E(X) &=& p sum_x x cdot (1-p)^{x-1} \\ &=& p cdot frac{1}{p^2} \\ &=& frac{1}{p} end{eqnarray} 本当に(frac{1}{p})になった。 両辺を微分したものが等しいって、なんでだっけ? 無記憶性 どうも、世の中には(n-1)回連続して失敗して(n)回目で初めて成功することを言っているものと、 (n)回連続して失敗して、(n+1)回目で初めて成功することを言っているものがある。 期待値も分散も若干違うものになる。 ここからは(n)回の失敗に続いて(n+1)回目で初めて成功するケースに切り替える。 その時の確率を(P(X=n))とする。で、失敗が(n)回以上連続して起こる確率(P(Xgeq n))を考える。 begin{eqnarray} P(X geq n) &=& P(X=n) + P(X=n+1) + P(X=n+2) + cdots \\ &=& p(1-p)^n + p(1-p)^{n+1} + p(1-p)^{n+2} + cdots \\ &=& p(1-p) left( 1 + (1-p)^1 + (1-p)^2 + cdots right) \\ &=& p(1-p) sum_{k=1}^{infty} (1-p)^{k-1} end{eqnarray} 途中の無限級数は(frac{1}{1-x})の級数展開になっていて、 以下みたいになる。 begin{eqnarray} P(X geq n) &=& p(1-p)^n sum_{k=1}^{infty} (1-p)^{k-1} \\ &=& p(1-p)^n frac{1}{1-(1-p)} \\ &=& p(1-p)^n frac{1}{p} \\ &=& (1-p)^n end{eqnarray} (n)回連続して失敗した上で、さらに連続して(k)回の失敗を重ねる確率を考える。 (n)回連続して失敗する確率は(P(Xgeq n))。 この条件の上でさらに(k)回失敗を重ねる確率は条件付き確率として(P(Xgeq n+k | X geq n))。 条件付き確率の定義と乗法定理から式を展開していく。(ここが難しかった...) begin{eqnarray} P(Xgeq n+k | X geq n) &=& frac{P((X geq n+k)cap (X geq n) )}{P(X geq n)} \\ &=& frac{P(X geq n+k)}{P(X geq n)} \\ &=& frac{(1-p)^{n+k}}{(1-p)^n} \\ &=& (1-p)^k \\ &=& P(X=k) end{eqnarray} ということで、以下が成り立つことがわかる。 begin{eqnarray} P(Xgeq n+k | X geq n) &=& P(X=k) end{eqnarray} よーく見てみると、(n)回連続して失敗した後に(k)回連続して失敗する確率と、 (n)回の失敗無しに、最初から(k)回連続して失敗する確率が同じである、と言っている。 凄まじいことに、(n)回連続して失敗することは、次の(k)回の失敗に全く影響を及ぼさない、と言っている。 何回失敗しようと次に失敗する確率はこれまでの失敗に影響されない。 つまり失敗する確率は過去の影響を受けない。 美しすぎる感じがする。

default eye-catch image.

モーメント母関数と確率密度関数

[mathjax] 期待値、分散、歪度、尖度...、確率分布を形成する確率密度関数の特徴を表す値で、 実は、相互に変換できる値なのだという...。読んでいったら若干感動したのでまとめてみる。 40近いオッさんがはじめてテイラー展開のありがたさを味わう瞬間の記録。 モーメント母関数とモーメント モーメント母関数を以下のように定義。 begin{eqnarray} M_x(t) &=& E(e^{tX}) \\ &=& int_{-infty}^{infty} e^{tx} f(x) dx end{eqnarray} 英語で書くとmoment generating function。モーメントを作る関数。 ここで(f(x))というのが確率密度関数。 もともとこの時点で積分が存在しないかもしれない。 確率密度関数によっては、期待値、分散、歪度、尖度を直接求めるのは難しい。 しかし、モーメント母関数の(r)階導関数からモーメント(mu_r)を解析的に求められる性質から、 期待値、分散、歪度、尖度を求めることができる。 さて...、40近いオッさんは思い出す。 指数関数のテイラー展開は、マクローリン級数を使って以下の通り。 begin{eqnarray} e^x = 1+x+x^2/2!+x^3/3!+cdots end{eqnarray} (tX)だと、 begin{eqnarray} e^{tX} = 1+tX+(tX)^2/2! +(tX)^3/3! + cdots end{eqnarray} 両辺の期待値をとる。右辺全体の期待値はそれぞれの項の期待値の和にできるので、 begin{eqnarray} E(e^{tX}) &=& M_x(t)\\ &=& E(1) + E(tX) + E((tX)^2/2!) + E((tX)^3/3!)) + cdots end{eqnarray} (t)に対する定数を出すと begin{eqnarray} M_x(t) &=& E(e^{tX})\\ &=& E(1) + E(X)t + E(X^2/2!)t^2 + E(X^3/3!)t^3 + cdots \\ &=& 1 + mu_1 t + (mu_2/2!)t^2 + (mu_3/3!)t^3 + cdots end{eqnarray} キター。おわかりだろうか...。 (M_x(t))は(t)に関する展開式の係数に各次数のモーメントを含んでいる。 (t)について1回微分すると0次までの項は消える。 2次以降の項の(t)の次数が1減って残る。1次の項だけ(t)が消える。 そのとき(t=0)とすると、階数の係数(M^{(r)}_X(0)=mu_r)だけ残る! つまり、以下のようなとんでもないことになる。 begin{eqnarray} M_X\'(0) = mu_1 \\ M_X\'\'(0) = mu_2 \\ M_X\'\'\'(0) = mu_3 \\ end{eqnarray} 各次数のモーメントである期待値、分散、歪度、尖度を、 モーメント母関数の(r)階導関数から求められるということになる。 指数分布のモーメント 試しに指数分布でやってみる。 begin{eqnarray} M_x(t) &=& int_{0}^{infty} e^{tx} lambda e^{-lambda x} dx \\ &=& lambda int_{0}^{infty} e^{(t-lambda)x} dx end{eqnarray} 指数関数の積分のところでおっ、と思ったけど、以下となる。 begin{eqnarray} M_x(t) = frac{ lambda }{ lambda -t} end{eqnarray} これ、解析的に微分できるのかな...と思うんだけども高校数学で暗記するやつ。 微分と積分を行ったり来たりするとわかる。 begin{eqnarray} M_x^{(1)}(t) &=& frac{ lambda }{ (lambda -t)^2} \\ M_x^{(2)}(t) &=& frac{ 2 cdot lambda }{(lambda -t)^3} \\ M_x^{(3)}(t) &=& frac{ 2 cdot 3 cdot lambda}{(lambda -t)^4} end{eqnarray} (t=0)とおくと、 begin{eqnarray} mu_1 &=& frac{1}{lambda} \\ mu_2 &=& frac{2}{lambda^2} \\ mu_3 &=& frac{2 cdot 3}{ lambda^3} \\ mu_4 &=& frac{2 cdot 3 cdot 4}{ lambda^4} end{eqnarray} これで、微分の数値計算をしなくても解析的に(mu_1)から(mu_4)が求まった。 そして永遠に微分し続けることで指数分布を形作る指標が決まっていく。 すごいなぁ...。