バナナでもわかる話

計量経済学・統計学が専門の大学院生です。統計学・経済学・投資理論・マーケティング等々に関する勉強・解説ブログ。ときどき趣味も。極力数式は使わずイメージで説明出来るよう心掛けていますが、時々暴走します。

統計検定1級対策解説~MCMC編1~

統計検定1級の範囲表によると、いつのまにやらMCMCやベイズが追加されていますね。

しかし、公式教科書は2013年出版の物しか出ていないため、その辺の範囲には対応していません。そこで、一通り何個かの記事を通して解説を行い、関連問題も上げていこうと思います。


目次


スポンサーリンク


MCMCとは何か

MCMCは、マルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo method)の略です。

ベイズ統計では、複雑な事後分布の期待値や分散を考える必要が出てきます。
その際、そうした期待値や分散は簡単に計算することが出来ない場合が多くあります。

そうした場合に対応するために、事後分布を定常分布とするマルコフ連鎖を構成することで、近似的な事後分布からのサンプルを得ることで、期待値や分散を計算しようというわけです。

そういった手法群を総称してMCMCと呼びます。

MCMCを考えるには、大雑把に言って次の3つを考える必要があるかと思います。


・マルコフ連鎖はどのような条件の下で収束するのか(定常分布が存在するのか)

・得たい事後分布を定常分布とするマルコフ連鎖をどうやって見つけるのか

・効率の良いアルゴリズムはどのようなものを考えればよいのか


今回は主に上二つに焦点を当てて説明していきます。

マルコフ連鎖

マルコフ連鎖と定常分布

一つ前の状態によって分布が変化する確率過程をマルコフ連鎖と呼びます。

例えば、次のような例を見てみます。

f:id:bananarian:20190712161745p:plain

この表は、状態A~Cから状態A~Cへの遷移行列です。
読み方は、例えば2行1列目には0.1が入っていますが、これは「状態Bから状態Aへ行く確率が0.1である」ことを表します。


ここで、経験的に初期状態で状態Aになる確率が0.3、状態Bになる確率が0.4,状態Cになる確率が0.3であったとしましょう。

この時、先ほどの遷移行列から考えれば

初期状態がAで、Bに遷移する確率は0.3×0.5=0.15であることがわかりますね。

同様にして初期状態BからBに遷移する確率は0.12

初期状態CからBへ遷移する確率は0.03

つまり、2時点目で状態Bである確率は0.3ということになります。

この計算を、全てのパターンで10回やってみることにします。ちょっと大変なので、Rに計算を投げてやると次の通り。

> XX=matrix(c(0.2,0.1,0.6,0.5,0.3,0.1,0.3,0.6,0.3),3,3)
> initV=c(0.3,0.4,0.3)
> res=initV
> 
> for(i in 1:10){
+ res=res%*%XX
+ print(round(res,digit=3))
+ }
     [,1] [,2] [,3]
[1,] 0.28  0.3 0.42
      [,1]  [,2] [,3]
[1,] 0.338 0.272 0.39
      [,1] [,2]  [,3]
[1,] 0.329 0.29 0.382
      [,1]  [,2]  [,3]
[1,] 0.324 0.289 0.387
      [,1]  [,2]  [,3]
[1,] 0.326 0.287 0.387
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386


結果を見てみるとわかる通り、あるところから、確率が0.326 0.288 0.386で動かなくなっています。

この状態を、マルコフ連鎖が収束したと呼び、この時の定常分布は0.326 0.288 0.386です。

定常分布への収束条件

マルコフ連鎖は必ず定常分布へ収束するわけではありません。
次の3条件を満たしている必要があります。

・既約性
ある状態からある状態へ、有限回の遷移でたどり着けるという条件が既約性です。

・正再帰性
任意の要素が、有限期間内で何度も発生しうることを正再帰性と呼びます。
つまり、正再帰性を満たさないとは、70回目の遷移以降、状態Aは発生しなくなるなどと言った状況です。

・非周期性
遷移に周期が存在しないことを非周期性と呼びます。


先ほど用意した遷移行列は、この3条件をしっかり満たしているため、定常分布へ収束したというわけです。


詳細つり合い条件

定常分布への収束条件はわかりました。
しかし、私たちがやりたいことは、ある特定の定常分布へ収束するようなマルコフ連鎖を構成することです。

つまり、ある特定の定常分布 p(x)へ収束するような遷移行列 Cを考えたいということになります。

そこで登場するのが詳細つり合い条件です。
ある状態iから状態jへ遷移する確率を C(j|i)とおくと、詳細つり合い条件は次のようになります。


 C(i|j)p(x=j)=C(j|i)p(x=i)

この条件が全てのi,jで成り立つとき、この遷移行列を用いた、マルコフ連鎖は定常分布 p(x)へ収束することが知られています。

つまり、私たちはこれで、欲しい事後分布を定常分布とするようなマルコフ連鎖を考えることが出来るようになったというわけです。

これらの条件を踏まえると、次に考えるべきは、


じゃあ、事後分布からのサンプルを取ってこれるような効率的なアルゴリズムって何?ってことですよね。

それについてはまた次回。

統計検定1級対策解説~一様最強力検定編~

今回は、一様最強力検定の話をします。

目次

スポンサーリンク


仮説検定

とりあえず、ザックリと仮説検定の話をしていきます。

仮説検定とは、「ある仮説を否定(棄却)出来るかどうかを検討したい」だとか、「ある仮説に基づいて、何かしらの意思決定を行いたい(モデル選択等)」といった場合に利用する、判定方法です。

基本的な枠組みは次の通り。

帰無仮説と対立仮説

まず、棄却出来るかどうかを考えたい仮説として帰無仮説 H_0を考えます。
例えば、ある正規分布に従っているデータを考えた時に、恐らく期待値は0ではないはずだと考えたとします。しかし、これは予想であって、客観的な根拠がありません。そこで、帰無仮説として期待値が0であるという仮説を設けます。この仮説を、仮説検定を使って棄却出来れば、「このデータの期待値は0ではない」ということが出来ますね。

つまり
帰無仮説 H_0 : 期待値が0
対立仮説 H_1 : 期待値は0ではない

というような2種類の仮説をまず、設定するということですね。

誤差の確率

検定をするにあたっては、帰無仮説が正しいのか、対立仮説が正しいのか分からない状態で意思決定を行います。つまり、帰無仮説が正しい状況と対立仮説が正しい状況の両方を考えた上で、妥当な意思決定を検討しなければなりません。

そこで気になるのは次のような数値でしょう。


・帰無仮説が正しい時に棄却される確率 (第一種の誤差)
・対立仮説が正しい時に受容される確率 (第二種の誤差)

両方の確率が0、つまり確率1で正しい仮説を選べるような検定なら最高ですが、少し考えるとわかりますが、そのような状況はほとんど発生しません。どうしてもミスをする確率は出て来てしまいます。

そこで、とりあえず第一種の誤差は小さい値に抑え込もうという方針で妥協します。

そこで利用するのが有意水準です。

つまり、有意水準5%に決めるとは、(帰無仮説内の任意のパラメータにおいて)第一種の誤差を5%以下に抑え込むことを指します。

数学的に書くと、

まず、考える1変量パラメータをもつ確率分布を P_{\theta}とおく。
パラメータ空間全体を \Theta、帰無仮説におけるパラメータの集合を \omegaと置くと、対立仮説におけるパラメータの集合は \Theta-\omegaとなりますが、

棄却域を Rとすると、第一種の誤差は P_{\theta}(R),ただし\theta \in \omega

よって、有意水準 \alphaに決めるとは

 P_{\theta}(R)≦\alpha for all  \theta \in \omega


当然、仮説検定の性能を考えるのであれば、

対立仮説が正しい場合に、きちんと帰無仮説を棄却する確率も知りたいですよね。これを検出力(Power)と呼びます。

つまり、

 P_{\theta}(R) ただし \theta \in \Theta-\omega

を検出力と呼ぶわけです。

一様最強力検定

ここまでが前段で、ここからが本題です。
で、有意水準を使って第一種の誤差は固定しました。

よって、第一種の誤差はとりあえず置いておいて、第二種の誤差をできうる限り小さくしたい(つまり、検出力を大きくしたい)ということに注力することにしましょう。

そこで真っ先に思い浮かぶのが、「対立仮説が正しい場合のパラメータ空間内で一様に検出力が最大になるような検定問題を考えることが出来ないか?」ということです。

つまり、検定の方法は色々考えられるわけですけど、パラメータがどんな値を取ったとしても(対立仮説が正しいもとで)、検出力が最大になる最強の検定があれば、それを使えばいいですよね?


結論から言うと、普通はそんな検定はありません。統計検定1級公式教科書でも、「一般に、そのような検定は存在しない」と書かれています。


ただし、「一般に」成り立たないだけで、特殊なケースでは一様最強力検定が存在します。今回はそれを示すことにします。

単純仮説

単純仮説とは、パラメータ空間全体 \Thetaが2種類の要素しか持たず、帰無仮説 \omega=\{\theta_1\}、対立仮説 \Theta-\omega=\{\theta_2\}のような仮説を単純仮説と呼びます。

この時実は、棄却域として R^*=\{x : \frac{p_{\theta_2}(x)}{p_{\theta_1}(x)}≧k\}を設定することで、これは最強力検定になります。

証明

ここで、任意の標本空間の領域を R,有意水準を \alphaと置くことにします。

 P_{\theta_2}(R^*)-P_{\theta_2}(R)=\int_{\{R^* \cap R^c\}} p_{\theta_2}(x)dx - \int_{\{R^{*c} \cap R\}} p_{\theta_2}(x) dx

これは、 R^*=(R^* \cap R) + (R^* \cap R^c)なのですぐわかります。

ここで、 R^{*} \cap R^c R^*の部分集合なので、 p_{\theta_2}(x)≧k p_{\theta_1}(x)が成り立ちます。

よって、 \int_{\{R^* \cap R^c\}} p_{\theta_2}(x)dx ≧k \int_{\{R^* \cap R^c\}} p_{\theta_1}(x)dx

同様に考えると

  \int_{\{R^{*c} \cap R\}} p_{\theta_2}(x) dx <k  \int_{\{R^{*c} \cap R\}} p_{\theta_1}(x) dx

よって

 P_{\theta_2}(R^*)-P_{\theta_2}(R)≧k(P_{\theta_1}(R^*)-P_{\theta_1}(R))

ここで、有意水準で \alphaで抑えていることから

 P_{\theta_2}(R^*)-P_{\theta_2}(R)≧k(P_{\theta_1}(R^*)-P_{\theta_1}(R))≧0

よって

 P_{\theta_2}(R^*)≧P_{\theta_2}(R)

以上より最強力検定であることが示せた。


統計検定1級対策解説~情報量基準編~

今回はAICです。
AIC周りの話は、カルバックライブラー情報量とAICにどういう結びつきがあるのか概要を知っておくことが重要になってきます。

スポンサーリンク



AICとは

AICは次のような式で表されます。


 AIC=-2 \sum_{i=1}^n log f(X_i,\hat{\theta})+2 dim(\theta)


まず \hat{\theta}は、n個のサンプルから得られるパラメータの最尤推定量です。そして、 dim(\theta)は分布の自由度となっています。

このAICとはなんぞやというのを理解することがこの章で重要な話になってきます。

何の役に立つんだって話

このAIC、何の役に立つのかというと、オーバーフィットを回避するためのモデル選択に利用できます。

オーバーフィットとは

例えば、こんな感じのデータがあったとします。
f:id:bananarian:20190515123427p:plain

x=seq(0,1,0.1)
train=sin(2*pi*x)+rnorm(11,sd=0.3)

で、どういう関数を当てはめたら良いかよくわからないので、最小二乗法(又は最尤法)を使って、多項式を当てはめてみることにするとします。


多項式っていうのはこういうやつですね。

 f(x)=\sum_{m=1}^M \alpha_{m-1}X^{m-1}+\epsilon

これはM-1次の多項式です。
Mはどうやって決めるのか問題がありますが、最尤法であれば、対数尤度最大、最小二乗法であれば二乗誤差を最小にするように、Mを適当に動かしてやると探すことが出来ますね。ちなみに、機械学習界隈ではよくRMSEが利用されますがどれを使おうとさしたる違いはありません。


で、普通に最小二乗法、又は最尤法を使うと次のような曲線が得られます。
f:id:bananarian:20190515124216p:plain

「おぉ~~全ての点を通ってて良い感じ~!」

と満足してはいけません。この曲線には問題があります。

同じ分布から新しいサンプルをこんな風に取ってきたとします。

xtest=seq(0,1,0.01)
test1=sin(2*pi*xtest)+rnorm(length(xtest),sd=0.3)

今の曲線をこの新しいサンプルにあてはめるとこんな感じ。
f:id:bananarian:20190515124809p:plain

変なところにすっ飛んでいってる所があり、予測が見当違いになってしまいます。


このように、既存のデータに極端にフィットしてしまい、新しいデータに対する予測性能が無くなってしまう現象をオーバーフィッティングと呼びます。

パラメータの次元(今回ならM)も最適化しようとした場合に、最尤法や最小二乗法を用いると、このような状況が発生します。それを回避するために利用される手法がAICです。

AICの式の意味

AICの式をもう一度確認してみます。

 AIC=-2 \sum_{i=1}^n log f(X_i,\hat{\theta})+2 dim(\theta)


「このAICを最小化するようなパラメータを選択するのが良い」というのがAICが利用される場合の主張です。

第1項は最大対数尤度ですね。普通の最尤法ではこのAICの第一項を最小化(マイナスがついているので)するような値を選択しているわけですが、それではパラメータの次元を弄り始めるとさっきのような問題が発生します。

そこでAICでは第二項目としてパラメータの自由度を加えます。

つまり、パラメータの自由度を増やせば増やすほど罰則がかかるようになっているわけですね。

AICを使ってみる

では、自由度をいじりながらAICの値を確認してみます。簡単にプロットしてやると次のようになりました。

f:id:bananarian:20190515130147p:plain

最小点は6ですね。つまり、「10次多項式なんていらない、6次で十分だ」とAICで主張できるわけです。

ちなみに6次でさっきのプロットをしてみるとこんな感じ。
f:id:bananarian:20190515130703p:plain
f:id:bananarian:20190515130715p:plain

見た目的にも最小二乗法を使った場合よりも良いですね。

余談(機械学習でよく使われる手法の話)

ちなみに、機械学習や予測の文脈では次のような手法も使われます。

持っているデータを

・最小二乗法(最尤推定)用のトレーニングデータ
・次元Mを決めるためのテストデータ

にわけることで、先ほどの問題を解決するといった手法です。


こんな感じでテストデータとトレーニングデータがあったとします。

x=seq(0,1,0.1)
train=sin(2*pi*x)+rnorm(11,sd=0.3)
test=sin(2*pi*x)+rnorm(11,sd=0.3)

この2つのデータに対して、トレーニングデータで学習した最尤推定量を元に、パラメータの自由度をいじって、その際のRMSEを見ることを考えます。

すると次のような推移になります。
f:id:bananarian:20190515131537p:plain

テストデータを用いた場合は、途中でRMSEの減少が止まっています。
つまり、テストデータを用いた場合、10次元も6次元もいらない。3次多項式で十分だと主張できることになります。

先ほどのAICは、尤度を利用していることからも分かる通り、サンプルサイズが少ない場合十分に機能しない可能性があります。そのため、データセットをこのように分ける手法が取られることもあるわけです。


カルバックライブラー情報量

カルバックライブラー情報量は、2つの分布間の(ある種の)距離を測る考え方です。
次のように定式化します。

2つの密度関数f,gについて

 KL(f|g)=\int f(x)log(\frac{f(x)}{g(x)})dx

2つの密度関数の対数の差の期待値を取っているのが見えますか?
もし、二つの密度が全く同じであれば0になります。

カルバックライブラー情報量とAICの関係性

ある分布 f(y;\theta)から標本Xが得られたとします。
この標本を使って分布のパラメータ \thetaを推定することを考える場合、その推定量 \hat{\theta}を使って尤度 f(y;\hat{\theta})を考えた場合、真の分布 f(y;\theta)に近くならなければならないはずです。


ここで、未知の標本(架空の物、得られていない) Yを使って真の分布と推定量を突っ込んだ分布の距離をカルバックライブラー情報量で考えてやると次のようになります。


 KL(f(y; \theta)|f(y;\hat{\theta}))=\int f(y;\theta)log(f(y; \theta))dy -\int f(y; \theta) log(f(y; \hat{\theta)}) dy

この値を小さくするためには、第1項目は弄りようが無いので、第2項目を考える必要が出てきますね。
つまり、第二項目の

 \int f(y; \theta) log(f(y; \hat{\theta)}) dy

を最大にするような推定量が汎化性能(予測とかの意味で良い性質を持つこと)の意味で良い推定量であると言えます。
これを平均対数尤度と呼んだりします。


もし、ここで新しいデータとして今の分布と同様の分布から Zが得られたとすると、

平均対数尤度の推定値は \frac{1}{n}\sum log(f(z; \hat{\theta}))となるわけですが、 \hat{\theta}はデータXの下での最尤推定量だったので、


 \frac{1}{n}\sum log(f(X; \hat{\theta)}) ≧ \frac{1}{n}\sum log(f(z; \hat{\theta)})

が成り立ちます。つまり、よりよい推定量を考えるためには平均対数尤度を最大化しなければならないのに、それよりも値が大きいXの対数尤度を最大化して満足してしまっていることがわかります。


しかし、そうはいっても「もしZを取ってこれたら」という話をしているに過ぎないので、実際はデータXでなんとかしなければなりません。

そこで、うまく平均対数尤度に近づくように補正する項を持ってきてやろうといことで漸近的に出てくるのが、先ほどAICで出てきた dim(\theta)です。

こうして、平均対数尤度を考える代わりに、既存のデータでどうにかできるようAICが生まれてきたというわけですね。

統計検定1級対策問題集~フィッシャー情報量編~

統計検定1級では、フィッシャー情報量を求める問題が頻出しています。そこで、フィッシャー情報量についてまとめました。

目次


スポンサーリンク


フィッシャー情報量とは


尤度関数 p(x;\theta)を考えます。
この時、この尤度関数について対数を取った対数尤度をパラメータ(ベクトル)について1階微分し、二乗して期待値を取ったものをフィッシャー情報量(行列) I(\theta)と呼びます。つまり

 I(\theta)= E[\{\frac{\partial}{\partial \theta} log( p(x;\theta) ) \}^2]


フィッシャー情報量(行列)は様々な場面で利用されますが、もっとも有名なのはクラメルラオの下限です。

あるパラメータ(ベクトル)における不偏推定量の分散の下限はフィッシャー情報量(行列)の逆数(逆行列)に等しくなります。これをクラメルラオの下限と呼びます。


不偏推定量は、バイアス0の推定量でしたので、不偏推定量のクラスで推定量を考える場合、MSEを最小にするような推定量を考えるには分散を小さくすることだけを考えれば良いので、その下限がわかるのは非常に強力です。


フィッシャー情報量の第二の導出

フィッシャー情報量は定義通り1階微分の二乗の期待値でも求めることは出来ますが、次の計算も適当な正則条件の下で同値になります。


 I(\theta)=E[-\frac{\partial^2}{\partial \theta^2} log(p(x;\theta))]

場合によっては二乗するより二回微分した方が簡単になることもあるので、この関係を知っておくことは重要です。一応簡単に証明しておきます。


 \frac{\partial^2}{\partial \theta^2} log(p(x;\theta))=\frac{\partial}{\partial \theta} \frac{1}{p}\frac{\partial p}{\partial \theta}

 =-\frac{1}{p^2} (\frac{\partial p}{\partial \theta})^2 + \frac{1}{p} \frac{\partial^2 p}{\partial \theta^2}

 =-(\frac{\partial logp}{\partial \theta})^2 +\frac{1}{p} \frac{\partial^2 p}{\partial \theta^2}


ここで E[ \frac{1}{p} \frac{\partial^2 p}{\partial \theta^2} ]=0…(微分と積分の交換が成り立てば)

よって E[- \frac{\partial^2}{\partial \theta^2} log(p(x;\theta)) ]=I(\theta)であることがわかります。


だんだんこの公式で慣れてくると、


「んんん??二乗したやつにマイナスつけるんだっけ?二回微分した奴にマイナスつけるんだっけ??」と混乱してくるかもしれませんが、その時はクラメルラオの下限を思い出すと良いかと思います。

分散の下限になるってことは、値は正になるはずってのが感覚的に普通ですよね。だから二乗したやつにマイナスがつくことなんてあり得ません。

フィッシャー情報量(行列)の具体例


具体的な計算が無いとよくわからないと思うので、1変量の場合と多変量の場合をそれぞれ具体例で確認してみようと思います。

ベルヌーイ分布

独立同一にベルヌーイ分布に従うサンプルを考えた場合、その同時尤度は次のよう。

 p^{\sum x_i}(1-p)^{n-\sum x_i}

そこで、対数を取ってやると

 log(p) \sum x_i + log(1-p) ( n-\sum x_i)

これを pについて二回微分すると次のようになりますね。

 -\frac{\sum x_i}{p^2}-\frac{n-\sum x_i}{(1-p)^2}

ここで、ベルヌーイ分布の期待値は pなので、期待値を取ってマイナスをとるとフィッシャー情報量 I(p)

 I(p)=\frac{n}{p}+\frac{n}{1-p}=\frac{n}{p(1-p)}


正規分布

パラメータが2つあるので、フィッシャー情報行列になります。


 \theta=(\mu,\sigma^2)とおく。

ここで、同時尤度は次のようになります。

 (2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum (x_i-\mu)^2}{2\sigma^2})

ここで、対数を取ると

 \frac{-n}{2} log(2\pi \sigma^2)-\frac{\sum (x_i-\mu)^2}{2\sigma^2}


さて、ここで \muに関して2回微分してやると

 \frac{-n}{\sigma^2}

これについてマイナスをつけて期待値を取ると

 \frac{n}{\sigma^2}

また \sigma^2に関して2階微分してやると

 \frac{n}{2 \sigma^4}-\frac{\sum (x_i-\mu)^2}{\sigma^6}

これについてマイナスをつけて期待値を取ると

 \frac{n}{2 \sigma^4}


 \mu, \sigma^2で1回ずつ微分してやると

 \frac{-\sum (x_i-\mu)}{\sigma^4}

これについてマイナスをつけて期待値を取ると

0


以上より、フィッシャー情報行列 I(\theta)は次のようになります。


 I(\theta)=\begin{pmatrix}
\frac{n}{\sigma^2} & 0 \\
0 & \frac{n}{2 \sigma^4} \\
\end{pmatrix}



ポアソン分布の場合や指数分布の場合も練習になるのでやってみると良いかと思います。

十分統計量に関する小話

ここ最近、統計検定関連の記事が続いていますね。

統計検定の範囲の中に「十分統計量」という単元があります。

前に2記事くらい練習問題記事を書きましたが、今回は十分統計量についてツラツラ思っていることを書いていきます。

完全に個人の感想です。


スポンサーリンク


そもそも統計検定の範囲になっているが

統計検定に関する勉強をしている受験生の方々は、当然今の時代に、統計学を必要としている方々のはずなので、「この十分統計量って何???」と立ち止まってしまったり、「わざわざこんなこと考える必要ある??」なんて感じる人も少なくないかと思います。


「何やら、十分統計量だなんだと難しい計算をわざわざゴリゴリやっていて、ラオブラックウェルの定理でMSEを小さくできる推定量を探せる!嬉しい!」と言われても(これが、統計検定1級の該当範囲)、


「うーん、そんな難しいことをせずとも、乱数を発生させて想定される推定量を比較してやれば良いんじゃないか」

と思う方も多いと思います。



そもそも、十分統計量の前に、不偏推定量についてもなんの意味があるんだろうと素朴な疑問を覚えた方も多いのではないでしょうか。


「わざわざUnbiassedな推定量を探して、そのクラスで分散小さい選手権を開催してる暇があったら、UnbiassedだろうがBiassedだろうが何でもいいから、想定される推定量を片っ端から集めてきて、乱数を100000個くらい発生させて、empiricalなバイアスと分散を推定して、使われる場面に最もフィットするであろう推定量を選べば良いではないか。」


なんて思いますよね。

また、「今の時代、ビッグデータの時代なのだから一致推定量のクラスで分散の小さいものを探したり、MSEが小さくなるようなものを探したほうがよっぽど建設的なのでは?」なんて疑問も持たれるかもしれません。



私も初めて不偏推定量について学んだとき思いましたし、前者については今もそう思っています。


今あげた疑問はコンピュータの計算技術が発展した現代においては(恐らく)真っ当な批判で、当然今はそういうシミュレーションが個人のPCレベルでも簡単に出来てしまうのだから、そうした環境にいる私たちは、小難しい計算を前に(面倒臭がって)そんな文句を言ってしまう訳です。


十分統計量が役に立っていた頃に思いを馳せる

コンピュータシミュレーションが簡単に行える今の時代においては、例えば大学院で指導教官に「僕は十分統計量の研究がしたい」なんて言っても、「うーん。まあ止やしないし、何も出てこないとまでは言わないけど、そこまで面白い話はないと思うよ〜。」と言われてしまうような分野です。

これは、「もうすでに一昔前に多くの研究者が議論し尽くした分野であり」「今や、十分統計量が役に立つ場面は少なくなってしまった」ためです。


そもそも、何で十分統計量や不偏推定量といった分野が、こうも厳密に議論される必要があったのでしょうか。
もう皆さん御察しの通りかと思いますが、これらの概念は、コンピュータで大規模なシミュレーションや計算を行うことが難しい時代に盛んに研究されていた概念で、要は陽に計算できる範疇で、理論上わかる範囲で如何に推定量に対して数学的な正当性を与えるか、妥当性を与えるかが研究されていた時代の考え方なのです。


だから「機械学習」「計算機統計学」「ニューラルネット」などと、どんどんコンピュータのスペックありきの手法が登場する現代で統計学を学ぼうとすると、どうにも十分統計量や不偏推定量のイメージが湧きにくい。この分野、なんか役に立つのか??なんて思ってしまうわけですね。


当然、役立つ分野は今でもあります。特に私の所属する経済学系の研究では、不偏性を持つことは非常に重要な意味を持ちますし、理論的な正当性を完備十分性の観点から考えることも理論研究では必要かもしれません。


ただ、そうした分野に属していないとどうしても「何のための理論なのかよくわからん?」となってしまいます。


統計検定、又は数理統計学の勉強をし始めているという方が本を開いてまず、詰まるのはここだと思いますので、とりあえずは研究の道に進む方以外は、「シミュレーションや大サンプルで正当化するのが難しかった時代は、こんなに細々と理論を積み上げていったのだなあ」と昔話でも読むように、大鏡や源氏物語を読んで、昔の日本人の生活に思いを馳せるかのように理解すれば十分かなと思います。


当然、俺は研究の道に進むので、あるいは教員の道に進むので、しっかり理解しなければならないんだという方であれば、過去の論文を遡ったり、昔から読まれてきている名著を読んで、じっくり格闘する必要はあるかと思います。


十分統計量についてしっかり理解するには、必然測度論が必要になってくるので、まずはそこからでしょうか。


でも、当然勉強し始めの段階では、そういった事情、背景を知らない方がほとんどのはずですので、

「どうか、ここの分野で挫折しないよう、時間との兼ね合いですが、あまり深入りしすぎないようにした方がいいですよ」

と勉強中の皆さんにはお伝えしておきたいです(昔、何も知らずに深入りしすぎて、ひどく時間を溶かした過去が....)

統計検定1級対策問題集~最尤推定量編~

統計検定1級では、最尤推定量を求める問題が頻出しています。そこで、最尤推定量を求める問題についてまとめました。

目次


スポンサーリンク


一様分布

パラメータ1つの場合

最大値が未知パラメータの次のような一様分布を考えます。

 f(x)=\frac{1}{\theta}

 0<x≦\theta

この時の、パラメータ \thetaの最尤推定量を求めます。


n個のサンプルが得られたとして、その同時尤度 l(\theta)

 l(\theta)=\frac{1}{\theta^n}

この尤度 l(\theta)を最大化する推定量 \hat{\theta}を考えます。
ただし、暗黙の条件に次の条件があることに注意します。

 max\{x_1,\cdots,x_n\}≦\theta


不等式制約があるので、正の変数 \lambdaを用いて、次のようなKKT条件を考える。

 L(\theta,\lambda)=\frac{1}{\theta^n}-\lambda (\theta-x_{\{max\}})

 \frac{\partial L(\theta,\lambda)}{\partial \theta} =-\frac{n}{\theta^{n+1}}-\lambda=0…①

 \frac{\partial L(\theta,\lambda)}{\partial \lambda} =-(\theta-x_{\{max\}})≦0…②

 \lambda (\theta-x_{\{max\}})=0…③

 \lambda=0であるとすると、①の等式が成り立たない。

よって \hat{\theta}=x_{\{max\}}


ちなみに、この問題は実際の試験で出題されています。
また、最大値ではなく最小値が未知パラメータである場合も同様の方法で \hat{\theta}=x_{\{min\}}と得られます。


パラメータが2つの場合

最大値も最小値も未知パラメータであるような一様分布として、次のようなものを考えます。

 f(x)=\frac{1}{\theta_2-\theta_1}

 \theta_1<x<\theta_2

この時、n個サンプルを得た時の同時尤度 l(\theta_1,\theta_2)は次のようになります。

  l(\theta_1,\theta_2)=\frac{1}{(\theta_2-\theta_1)^n}

更に、暗黙的に次の条件があることがわかります。

 \theta_1≦x_{\{min\}}

 x_{\{max\}}≦\theta_2

不等式制約なので、正の変数 \lambda_1,\lambda_2を用いて次のようなKKT条件を考えます。

 L(\theta_1,\theta_2,\lambda_1,\lambda_2)=\frac{1}{(\theta_2-\theta_1)^n}-\lambda_1(x_{\{min\}}-\theta_1)-\lambda_2(\theta_2-x_{\{max\}})

 \frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \theta_1}=\frac{n}{(\theta_2-\theta_1)^{n+1}} +\lambda_1=0…①

 \frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \theta_2}=\frac{-n}{(\theta_2-\theta_1)^{n+1}} -\lambda_2=0…②

 \frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \lambda_1}=(x_{\{min\}}-\theta_1)≦0…③

 \frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \lambda_2}=(\theta_2-x_{\{max\}})≦0…④

 \lambda_1(x_{\{min\}}-\theta_1)=0…⑤

 \lambda_2(\theta_2-x_{\{max\}})=0…⑥

⑤、⑥について、もし \lambda_1=0,\lambda_2=0であったとすると、①、②の等式は成り立たない。

よって、それぞれのパラメータの最尤推定量は

 \hat{\theta_1}=x_{\{min\}}

 \hat{\theta_2}=x_{\{max\}}


ベルヌーイ分布

ベルヌーイ分布の確率質量関数 f(x)は次のよう。

 f(x)=p^x(1-p)^{1-x}

ここで、n個のサンプルを考えた場合の同時尤度 l(p)

 l(p)=p^{\{\sum_{i=1}^n x_i\}} (1-p)^{\{n-\sum_{i=1}^n x_i\}}

ここで、 l(p)の対数を取ると次のよう。

 ln(p)=(\sum_{i=1}^n x_i) log(p) +(n-\sum_{i=1}^n x_i) log(1-p)

ここで pの最尤推定量 \hat{p}は次の等式を満たす pに等しい。

 \frac{\partial ln(p)}{\partial p}=\frac{\sum_{i=1}^n x_i}{p}-\frac{n-\sum_{i=1}^n x_i}{1-p}=0

よって
 \hat{p}=\frac{\sum_{i=1}^n x_i}{n}


ポアソン分布

ポアソン分布の確率質量関数 f(x)は次の通り。

 f(x)=\frac{e^{-\lambda} \lambda^{x}}{x!}

サンプルをn個取り出したとすると同時尤度 l(\lambda)は次のよう。

 l(\lambda)=\frac{e^{-n\lambda} \lambda^{\sum_{i=1}^n x_i}}{\prod_{i=1}^n (x_i !)}

ここで、対数を取ってやると

 ln(\lambda)=-log(\prod_{i=1}^n (x_i !))-n\lambda +(\sum_{i=1}^n x_i)log(\lambda)

よって \lambdaの最尤推定量 \hat{\lambda}は次の等式を満たす \lambdaに等しい。

 \frac{\partial ln(\lambda)}{\partial \lambda}=-n+\frac{\sum_{i=1}^n x_i}{\lambda}=0

以上より

 \hat{\lambda}=\frac{\sum_{i=1}^n x_i}{n}


正規分布

正規分布の確率密度関数は次の通り。

 f(x)=\frac{1}{\sqrt{2\pi \sigma^2}} exp(-\frac{((x-\mu)^2)}{2\sigma^2})

同時尤度は
 l(\mu,\sigma)=(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n (x_i-\mu)^2}{2 \sigma^2})

ここで対数を取ると

 ln(\mu,\sigma)=\frac{-n}{2} log(2 \pi \sigma^2)-\frac{\sum_{i=1}^n (x_i-\mu)^2}{2 \sigma^2}


この時、 \mu,\sigmaの最尤推定量 \hat{\mu},\hat{\sigma}は次の等式から得られる。


 \frac{\partial ln(\mu,\sigma)}{\partial \mu}=\frac{\sum_{i=1}^n (x_i-\mu)}{\sigma^2}=0

 \frac{\partial ln(\mu,\sigma)}{\partial \sigma}=\frac{-n}{\sigma}+\frac{\sum_{i=1}^n (x_i-\mu)^2}{\sigma^3}=0

よって、

 \hat{\mu}=\frac{\sum_{i=1}^n x_i}{n}

 \hat{\sigma}=\sqrt{\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n}}


指数分布

指数分布の確率密度関数は次の通り

 f(x)=\lambda exp(-\lambda x)

この時同時尤度は

 l(\lambda)=\lambda^n exp(-\lambda \sum_{i=1}^n x_i)

ここで、対数を取ると

 ln(\lambda)=n log(\lambda) -\lambda \sum_{i=1}^n x_i

この時、 \lambdaの最尤推定量 \hat{\lambda}は次の等式から得られる。

 \frac{\partial ln(\lambda)}{\partial \lambda}=\frac{n}{\lambda}-\sum_{i=1}^n x_i=0

よって、

 \hat{\lambda}=\frac{n}{\sum_{i=1}^n x_i}

統計検定1級対策問題集~十分統計量編2~

十分統計量に関する問題2記事目です。

目次


スポンサーリンク


ラオブラックウェルの定理

ラオブラックウェルの定理とは

初めに完備十分統計量を考える上で重要になってくる「ラオブラックウェルの定理」の証明についてやっておきます。

実際の数理統計の本であれば、ラオブラックウェルの定理を示した後、完備性について解説し、完備十分統計量について話が移ります。

ただ、統計検定1級公式参考書では完備十分統計量までは触れていません。
ラオブラックウェルの定理までの説明で終わっているので、とりあえずこの定理の証明までは確認しておきます。


次のような定理をラオブラックウェルの定理と呼びます。
 T \thetaの十分統計量とする。ここで、 \thetaのある推定量\delta(X)について、次のような推定量 \delta_1(T)をラオブラックウェル推定量と呼ぶことにする。

 \delta_1(T)=E_{\theta}[\delta(X)|T]

そして、ラオブラックウェル推定量が満たす次のような性質をラオブラックウェルの定理と呼ぶ。

 E_{\theta}[(\delta_1(T)-\theta)^2]≦E_{\theta}[(\delta(X)-\theta)^2]

不等式の両サイドは平均二乗誤差になっています。
つまり、この不等式からわかることは

「ある推定量 \deltaを考えた時に、それよりも平均二乗誤差を小さくする(又は同等)推定量を、十分統計量を条件付けることで考えることが出来る」

ということです。単純ですが強力な定理です。

ラオブラックウェルの定理証明

まず、 E_{\theta}[\delta_1(T)]=E_{\theta}[\delta(X)]であることを示します。

 E_{\theta}[\delta_1(T)]=\int_T \int_X \delta(X) dP(X|T) dP(T)

 =\int_X \delta(X) dP(X)=E_{\theta}[\delta(X)]


また、 E_{\theta}[\delta_1(T)^2]≦E_{\theta}[\delta(X)^2]であることも示します。

 E_{\theta}[\delta_1(T)^2]=E_{T;\theta}[ E_{\theta}[\delta(X)|T]^2]

 E_{\theta}[\delta(X)^2]=E_{T;\theta}[E_{\theta}[\delta(X)^2|T]]

更にイェンゼン不等式を用いて、
 E_{T;\theta}[ E_{\theta}[\delta(X)|T]^2]≦E_{T;\theta}[E_{\theta}[\delta(X)^2|T]]

以上より示せた。


最後にラオブラックウェルの定理を示します。
 E_{\theta}[(\delta_1(T)-\theta)^2]=E_{\theta}[\delta_1(T)^2]-2\theta E_{\theta}[\delta_1(T)]+\theta^2

 E_{\theta}[(\delta(X)-\theta)^2]=E_{\theta}[\delta(X)^2]-2\theta E_{\theta}[\delta(X)]+\theta^2


上二つの性質から
 E_{\theta}[(\delta_1(T)-\theta)^2]≦E_{\theta}[(\delta(X)-\theta)^2]

フィッシャーネイマンの分解定理

負の二項分布

負の二項分布の確率質量関数 f(x;p,r)は次のようになります。

 f(x;p,r)=\begin{eqnarray*}
  && {}_{r+x-1} C _x \\
\end{eqnarray*} p^r (1-p)^x


 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、同時分布は

 P(x_1,\cdots,x_n ;n,p,r) = \prod_{i=1}^n \{ \begin{eqnarray*}
  && {}_{r+x_i-1} C _{x_i} \\
\end{eqnarray*} p^r (1-p)^{x_i} \}

 = \{ \prod_{i=1}^n  \begin{eqnarray*}
  && {}_{r+x_i-1} C _{x_i} \\
\end{eqnarray*} \} p^{nr} (1-p)^{ \sum_{i=1}^n x_i }

この時、 T(X)=\sum_{i=1}^n x_iがパラメータ pの十分統計量であることを示します。

フィッシャーネイマンの分解定理より、

 h(X)=\{\prod_{i=1}^n \begin{eqnarray*}
  && {}_{r+x_i-1} C _{x_i} \\
\end{eqnarray*} \}

 g(T(X),p)=p^{nr}(1-p)^{T(X)}

とみると、  T(X)=\sum_{i=1}^n x_iがパラメータ pの十分統計量であることがわかる。


ガンマ分布

ガンマ分布の確率密度関数 f(x;\alpha,\beta)は次のようになります。

 f(x;\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}

 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、

 T_{\alpha}(X)=\prod_{i=1}^n x_i \alphaの十分統計量
 T_{\beta}(X)=\sum_{i=1}^n x_i \betaの十分統計量です。これを示します。


同時分布は

 P(x_1,\cdots,x_n ;\alpha,\beta)=\frac{\beta^{n\alpha}}{\Gamma(\alpha)^n} \{\prod_{i=1}^n x_i\}^{\alpha-1} e^{-\beta\sum_{i=1}^n x_i}


よって、

 h_{\alpha}(X)=e^{-\beta\sum_{i=1}^n x_i}

 g(T_{\alpha}(X),\alpha)=\frac{\beta^{n\alpha}}{\Gamma(\alpha)^n} \{T_{\alpha}(X)\}^{\alpha-1}

とみると、フィッシャーネイマンの分解定理より、 T_{\alpha}(X)=\prod_{i=1}^n x_i \alphaの十分統計量。


また、 h_{\beta}(X)= \frac{1}{\Gamma(\alpha)^n} \{\prod_{i=1}^n x_i\}^{\alpha-1}

 g(T_{\beta}(X),\beta)=\beta^{n\alpha} e^{-\beta T_{\beta}(X)}

と見ると、フィッシャーネイマンの分解定理より、 T_{\beta}(X)=\sum_{i=1}^n x_i \betaの十分統計量


一様分布

下限が0,上限が未知パラメータ \thetaであるときの一様分布を考えます。この時、密度関数は次のようになります。

 f(x;\theta)=\frac{1}{\theta}   (0≦x≦\theta)

これは、見方を変えると次のように書くことも出来ます。

 f(x;\theta)=\frac{1_{\{0≦x≦\theta\}}}{\theta}

ここで、 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、 \thetaの十分統計量は max\{x_1,\cdots , x_n\}となります。

これを示します。

同時分布は次のようになるので

 P(x_1,\cdots x_n;n,\theta)=\frac{1_{\{max\{x_1,\cdots , x_n\}≦\theta\}}}{\theta^n}

フィッシャーネイマンの分解定理より、 \thetaの十分統計量は max\{x_1,\cdots , x_n\}であることは明らか。

統計検定1級対策問題集~十分統計量編1~

今回は十分統計量に関する問題をまとめていきます。少し量が多いので、2回に分けます。
目次


スポンサーリンク


十分統計量とは

標本 Xとその分布のパラメータ \theta を考えます。この時、次の等式が成り立つ統計量 T(X)を十分統計量と呼びます。

 P(X=x|T(X)=t,\theta)=P(X=x|T(X)=t)

これは、どう解釈すれば良いかというと、

「パラメータ \thetaの情報を T(X)は十分に持っている」

と解釈出来ます。

等式を見ていただけるとわかるように \thetaがあろうと無かろうと、 Xの分布には変化がありません。

フィッシャーネイマンの分解定理

フィッシャーネイマンの分解定理

 T(X) \thetaの十分統計量であるとき、確率密度関数(確率質量関数)は次のように分解できる。

 f(x;\theta)=h(x)g(T(x),\theta)

これをフィッシャーネイマンの分解定理(factorization theorem)と呼びます。

この定理の証明は間違いなく出題されませんので、省略します。
というのも、この分解定理の証明には測度論を利用します。これは、1級範囲を逸脱していますし、実際公式教科書でも、しれっと証明を省略しています。

この定理を利用して、ある T(X)が十分統計量であることを証明します。

ベルヌーイ分布

ベルヌーイ分布の確率質量関数 f(x;p)は次のようでした。

 f(x;n,p)=p^x(1-p)^{1-x}

 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、同時分布は

 P(x_1,\cdots,x_n ;n,p)=p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i}


ですが、この時、 T(X)=\sum_{i=1}^n x_iがパラメータ pの十分統計量であることを示します。

 P(x_1,\cdots,x_n ;n,p)=\frac{p}{1-p}^{\sum_{i=1}^n x_i} (1-p)^n=(\frac{p}{1-p})^{T(X)} (1-p)^n

この時 h(x)=1,g(T(x),p)=(\frac{p}{1-p})^{T(X)} (1-p)^nと見ると、フィッシャーネイマンの分解定理から、 T(X)=\sum_{i=1}^n x_iはパラメータ pの十分統計量です。


ポアソン分布

ポアソン分布の確率質量関数 f(x;p)は次のようでした。

 f(x;\lambda)=\frac{e^{-\lambda} \lambda^x}{x!}

 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、同時分布は

 P(x_1,\cdots,x_n ;n,\lambda)=\frac{e^{-n\lambda} \lambda^{\sum_{i=1}^n x_i}}{\prod_{i=1}^n x_i}

ですが、この時、 T(X)=\sum_{i=1}^n x_iがパラメータ \lambdaの十分統計量であることを示します。

 P(x_1,\cdots,x_n ;n,\lambda)=\frac{1}{\prod_{i=1}^n x_i} e^{-n\lambda} \lambda^{\sum_{i=1}^n x_i}=\frac{1}{\prod_{i=1}^n x_i} e^{-n\lambda} \lambda^{T(X)}


 h(x)=\frac{1}{\prod_{i=1}^n x_i} ,g(T(x),\lambda)=e^{-n\lambda} \lambda^{T(X)}と見ると、フィッシャーネイマンの分解定理から T(X)=\sum_{i=1}^n x_iはパラメータ \lambdaの十分統計量です。


正規分布

正規分布の確率密度関数 f(x;\mu,\sigma)は次のようでした。

 f(x;\mu,\sigma)=\frac{1}{\sqrt{2 \pi \sigma^2}} exp(-\frac{(x-\mu)^2}{2\sigma^2})

 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、同時分布は

 P(x_1,\cdots,x_n ;n,\mu,\sigma) = (2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n (x_i-\mu)^2}{2\sigma^2})


この時、 T_{\mu}(X)=\sum_{i=1}^n x_iがパラメータ \muの十分統計量、 (T_{\mu}(X),T_{\sigma}(X))=(\sum_{i=1}^n x_i,\sum_{i=1}^n x_i^2)がパラメータ (\mu,\sigma)の十分統計量であることを示します。


 P(x_1,\cdots,x_n ;n,\mu,\sigma) = (2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n x_i^2-2\mu \sum_{i=1}^n x_i+n\mu^2}{2\sigma^2})

まず、
 P(x_1,\cdots,x_n ;n,\mu,\sigma) =(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n x_i^2-2\mu T_{\mu}(X)+n\mu^2}{2\sigma^2})
 =(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n x_i^2}{2\sigma^2}) exp(-\frac{-2\mu T_{\mu}(X)+n\mu^2}{2 \sigma^2})

このように見ると、 h(x)=(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n x_i^2}{2\sigma^2})であり、
 g(T_{\mu}(X),\mu)=exp(-\frac{-2\mu T_{\mu}(X)+n\mu^2}{2 \sigma^2})のため、フィッシャーネイマンの分解定理から、 T_{\mu}(X) \muに関する十分統計量です。

また、
 P(x_1,\cdots,x_n ;n,\mu,\sigma) =(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{T_{\sigma}(X)-2\mu T_{\mu}(X)+n\mu^2}{2\sigma^2})

であるので、

 h(x)=1,g(T_{\sigma}(X),T_{\mu}(X),\mu,\sigma)=(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{T_{\sigma}(X)-2\mu T_{\mu}(X)+n\mu^2}{2\sigma^2})と見ると、 (T_{\mu}(X),T_{\sigma}(X))=(\sum_{i=1}^n x_i,\sum_{i=1}^n x_i^2)はパラメータベクトル (\mu,\sigma)の十分統計量ベクトルである。


分解定理を使わない例

当然、十分統計量かどうかは分解定理を使わずとも、定義から示すことも可能です。
ただ、計算がしんどいので、普通は分解定理で示します。

ベルヌーイ分布で考える場合、結構計算が簡単なのでこれで確認してみます。


ベルヌーイ分布の確率質量関数 f(x;p)は次のようでした。

 f(x;n,p)=p^x(1-p)^{1-x}

 x_1,\cdots,x_nのサンプルを独立同一に得たとすると、同時分布は

 P(x_1,\cdots,x_n ;n,p)=p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i}


ですが、この時、 T(X)=\sum_{i=1}^n x_iがパラメータ pの十分統計量であることを示します。


 P(x_1,\cdots,x_n |n,p)=p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i}=p^{T(X)}(1-p)^{n-T(X)}

ここで
 P(T(X)=t |n,p)=p^{t}(1-p)^{n-t} 1_{\{T(x)=t\}}

 P(x_1,\cdots,x_n,T(X)=t |n,p)= \sum_{x:\{T(x)=t\}} P(x_1,\cdots,x_n |n,p)=\begin{eqnarray*}
  && {}_n C _t \\
\end{eqnarray*} p^{T(X)}(1-p)^{n-T(X)}


以上より
 P(x_1,\cdots,x_n |n,p,T(X)=t)=\frac{P(T(X)=t |n,p)}{P(x_1,\cdots,x_n,T(X)=t |n,p)}=\frac{1}{\begin{eqnarray*}
  && {}_n C _t \\
\end{eqnarray*}} 1_{\{T(x)=t\}}=P(x_1,\cdots,x_n |n,T(X)=t)

確かに示せました。

統計検定1級対策問題集~ベータ分布編~

統計検定1級対策のために役立ちそうな計算問題をまとめるやつやっていきます。
統計検定前の最終チェックや、統計検定の勉強何をすれば分からないという場合に活用ください。


今回はベータ分布関連。
ガンマ分布の時と同様、部分積分をループさせる計算がいっぱい出てきます。
ベータ分布は部分積分ゲーなので、手を動かして慣れるのが良いかと思います。


スポンサーリンク


目次

ベータ分布の特徴

 f(x) =\frac{x^{\alpha-1} (1-x)^{\beta-1}}{B(\alpha,\beta)}
ただし B(\alpha,\beta)=\int_{0}^1 x^{\alpha-1} (1-x)^{\beta-1}dx

・連続値の分布
 0<x<1
 \alpha,\betaは正

非常にゴチャついていて、嫌になるかもしれませんが、よく見てください。
 B(\alpha,\beta)は単なる正規化定数(積分したらうまく1になるよう調整するための定数)に過ぎず、分布の本体は x^{\alpha-1} (1-x)^{\beta-1}であることがわかります。そう思うと、とっつきにくさは多少和らぐのではないでしょうか。

正規化定数の計算

ベータ分布の期待値や分散の導出、その他様々な計算で、次の性質を利用します。

 B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}

 \Gamma()はガンマ関数です。ガンマ分布の記事で嫌というほど使いましたね笑

まず、この性質について確認します。

 B(\alpha,\beta)=\int_{0}^1 x^{\alpha-1} (1-x)^{\beta-1}dx

 =\int_{0}^1 (\frac{1}{\alpha})^{'} (1-x)^{\beta-1} dx

 =\int_{0}^1 \frac{\beta-1}{\alpha} x^{\alpha} (1-x)^{\beta-2}dx

=\int_{0}^1 \frac{(\beta-1)(\beta-2)}{\alpha(\alpha+1)} x^{\alpha+1}(1-x)^{\beta-3}dx

 \cdots

 =\int_{0}^1 \frac{ \Gamma(\beta) \Gamma(\alpha) }{ \Gamma(\alpha+\beta-1) }  x^{\alpha+\beta-2} dx

 =\frac{\Gamma(\beta) \Gamma(\alpha)}{\Gamma(\alpha+\beta)}

モーメント周りの計算

積率母関数は、存在するのですが導出しません。
というのも、ベータ分布の積率母関数はウィキか何かで調べてもらえればわかる通り、複雑すぎて役に立ちません。
実際公式テキスト(2015年出版時点)でも、ベータ分布の積率母関数は一切触れられず、スルーされています。
まあ、なので導出する必要もないだろうというわけで省略します。

期待値の導出

定義に従った計算

定義に従って期待値を求めてみます。
 E[x]=\frac{1}{B(\alpha,\beta)}\int_{0}^1 x^{\alpha} (1-x)^{\beta-1}dx

これも、さっき導出した B(\alpha,\beta)の計算と同様の手順をひたすら繰り返すと、

 = \frac{1}{B(\alpha,\beta)} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\alpha}{\alpha+\beta}

 = \frac{\alpha}{\alpha+\beta}


スポンサーリンク


分散の導出

定義に従った計算

 Var[x]=E[x^2] -(E[x])^2

 E[x^2] =\frac{1}{B(\alpha,\beta)}\int_{0}^1 x^{\alpha+1} (1-x)^{\beta-1}dx

これも、次数がズレただけでさっきと同じ部分積分の繰り返しですね。打ち込むのが大変なので省略します笑
一回は出しておいた方が良いと思います。

計算してやると

 Var[x]=\frac{\alpha \beta}{(\alpha+\beta)^2 (\alpha+\beta+1)}

になるはずです。


ベータ分布の導出

実は、ベータ分布は2つの独立なガンマ分布に従う確率変数を用いて導出出来ます。
 x_1 ~ Ga(\alpha_1,\beta)
 x_2 ~ Ga(\alpha_2,\beta)

について、
 X=\frac{x_1}{x_1+x_2}

 Y=x_1+x_2

と置きます。

この時、逆変換した際のヤコビアンは Yなので


 f(X,Y) = Y \frac{\beta^{\alpha_1}}{\Gamma(\alpha_1)} (XY)^{\alpha_1-1}exp(-\beta (XY)) \frac{\beta^{\alpha_2}}{\Gamma(\alpha_2)} (Y(1-X))^{\alpha_2-1} exp(-\beta Y(1-X))

 = \frac{\beta^{\alpha_1+\alpha_2}}{\Gamma(\alpha_1)\Gamma(\alpha_2)} X^{\alpha_1-1}(1-X)^{\alpha_2-1} Y^{\alpha_1+\alpha_2-1} exp(-\beta Y)

 \frac{\beta^{\alpha_1+\alpha_2}}{\Gamma(\alpha_1+\alpha_2)}  Y^{\alpha_1+\alpha_2-1} exp(-\beta Y)   \frac{X^{\alpha_1-1}(1-X)^{\alpha_2-1}}{B(\alpha_1,\alpha_2)}

はい、見事にガンマ分布の密度関数とベータ分布の密度関数の積に分解することが出来ました。あとはXに関して周辺分布を考えてやればよく、ガンマ分布の密度関数は全範囲で積分すると1になるので

Xはベータ分布に従います。

ちなみに、Yがガンマ分布に従うのは、ガンマ分布の再生性からある意味で自明ですね。




ベータ分布と二項分布の関係

ベータ分布の上側確率は二項分布の確率関数の和と解釈出来ます。

 \int_{p}^{1} \frac{x^{k-1} (1-x)^{n-k}}{B(k,n-k+1)}dx =\begin{eqnarray*}
  && {}_n C _{k-1} \\
\end{eqnarray*} p^{k-1} (1-p)^{n-k+1} + \int_{p}^{1} \frac{x^{k-2}(1-x)^{n-k+1}}{B(k-1,n-k+2)} dx

 \cdots

=\sum_{z=0}^{k-1} \begin{eqnarray*}
  && {}_n C _{z} \\  \end{eqnarray*} p^z (1-p)^{n-z}


リンク

統計学を勉強するためのオススメ本

www.bananarian.net

2017年数理1級の解説記事

www.bananarian.net

その他の問題記事

統計検定1級対策問題集~二項分布編~

www.bananarian.net

統計検定1級対策問題集~ポアソン分布編~

www.bananarian.net

統計検定1級対策問題集~負の二項分布編~

www.bananarian.net

統計検定1級対策問題集~正規分布編~

www.bananarian.net

統計検定1級対策問題集~指数分布編~

www.bananarian.net

統計検定1級対策問題集~ガンマ分布編~

www.bananarian.net

統計検定1級対策問題集~ガンマ分布編~

統計検定1級対策のために役立ちそうな計算問題をまとめるやつやっていきます。
統計検定前の最終チェックや、統計検定の勉強何をすれば分からないという場合に活用ください。


今回はガンマ分布関連。
ガンマ関数の処理に慣れるまでは難しいかもしれません。


スポンサーリンク


目次

ガンマ分布の特徴

 f(x) =\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}exp(-\beta x)

・連続値の分布
 x >0
 \alpha,\betaは正
 \Gamma(\alpha)=\int_{0}^{\infty} t^{\alpha-1}exp(-t) dt

ガンマ関数は、階乗を一般化したものです。


モーメント周りの計算

積率母関数の導出

積率母関数の定義は次の通りでした。
 E[exp(tx)]

計算していきます。

 E[exp(tx)]=\int_{0}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} exp(tx-\beta x) dx

 = \int_{0}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha) (\beta-t)^{\alpha-1}} ((\beta-t) x)^{\alpha-1} exp(tx-\beta x) dx

 = \frac{\beta^{\alpha}}{\Gamma(\alpha) (\beta-t)^{\alpha-1}} \int_{0}^{\infty}  ((\beta-t) x)^{\alpha-1} exp(-(\beta-t) x) dx

 =\frac{\beta^{\alpha}}{\Gamma(\alpha) (\beta-t)^{\alpha}}  \int_{0}^{\infty}  ((\beta-t) x)^{\alpha-1} exp(-(\beta-t) x) d((\beta-t)x)

 =\frac{\beta^{\alpha} \Gamma(\alpha)}{\Gamma(\alpha) (\beta-t)^{\alpha}}

 =(\frac{\beta}{\beta-t})^{\alpha}


期待値の導出

定義に従った計算

定義に従って期待値を求めてみます。
 E[x]=\int_{0}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha} exp(-\beta x) dx

 =\frac{1}{\Gamma(\alpha)} \int_{0}^{\infty} (\beta x)^{\alpha} exp(-\beta x) dx

 = \frac{1}{\beta \Gamma(\alpha)} \int_{0}^{\infty} (\beta x)^{\alpha} exp(-\beta x) d(\beta x)

 =\frac{\Gamma(\alpha+1)}{\Gamma(\alpha) \beta}

 =\frac{\alpha}{\beta}


スポンサーリンク


分散の導出

定義に従った計算

 Var[x]=E[x^2] -(E[x])^2

 E[x^2] =\int_{0}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha+1} exp(-\beta x) dx

 = \int_{0}^{\infty} \frac{1}{\beta \Gamma(\alpha)} (\beta x)^{\alpha+1} exp(-\beta x) dx

 = \frac{1}{\beta^2 \Gamma(\alpha)} \int_{0}^{\infty}  (\beta x)^{\alpha+1} exp(-\beta x) d(\beta x)

 = \frac{\Gamma(\alpha+2)}{\beta^2 \Gamma(\alpha)}

 \frac{\alpha(\alpha+1)}{\beta^2}

 Var[x] = \frac{\alpha(\alpha+1)}{\beta^2} - \frac{\alpha^2}{\beta^2}=\frac{\alpha}{\beta^2}


ガンマ分布の再生性について

 x_1 ~ Ga(\alpha_1,\beta)
 x_2 ~ Ga(\alpha_2,\beta)

に関して、ガンマ分布は再生性がある。これは先ほど導出した積率母関数を考えると明らかで、

 (\frac{\beta}{\beta-t})^{\alpha_1} (\frac{\beta}{\beta-t})^{\alpha_2}=(\frac{\beta}{\beta-t})^{\alpha_1+\alpha_2}

積率母関数と分布は1対1対応するため、 x_1+x_2もガンマ分布することがわかる。

ガンマ分布とポアソン分布の関係

ガンマ分布の上側確率は適当な仮定のもとでポアソン確率関数の和と解釈出来ます。
この性質はポアソン過程を考える際に用います。


まず、 \alphaを1以上の正整数とする。そして正の実数 \omegaを用意して、次のようなものを考える。

 \int_{\omega}^{\infty} f(x) dx =\int_{\omega}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}exp(-\beta x) dx

 = -\int_{\omega}^{\infty} \frac{1}{\Gamma(\alpha)} (\beta x)^{\alpha-1} (exp(-\beta x))^{'} dx

 =-\{ [ \frac{(\beta x)^{\alpha-1}}{ exp(\beta x)\Gamma(\alpha)} ]_{\omega}^{\infty} - \int_{\omega}^{\infty} \frac{\alpha-1}{\Gamma(\alpha)}  (\beta x)^{\alpha-2} exp(-\beta x) dx \}

 =\frac{(\beta \omega)^{\alpha-1}  exp(-\beta \omega)}{ (\alpha-1)!} + \int_{\omega}^{\infty} \frac{1}{\Gamma(\alpha-1)}  (\beta x)^{\alpha-2} exp(-\beta x) dx

 =\frac{(\beta \omega)^{\alpha-1}  exp(-\beta \omega)}{ (\alpha-1)!} + \frac{(\beta \omega)^{\alpha-2}  exp(-\beta \omega)}{ (\alpha-2)!} + \cdots +\frac{(\beta \omega)^{\alpha-\alpha}  exp(-\beta \omega)}{ (\alpha-\alpha)!}

 =\sum_{k=0}^{\alpha-1} \frac{(\beta \omega)^k  exp(-\beta \omega)}{ k!}

これはパラメータ \beta \omegaのポアソン確率関数の和。


リンク

統計学を勉強するためのオススメ本

www.bananarian.net

2017年数理1級の解説記事

www.bananarian.net

その他の問題記事

統計検定1級対策問題集~二項分布編~

www.bananarian.net

統計検定1級対策問題集~ポアソン分布編~

www.bananarian.net

統計検定1級対策問題集~負の二項分布編~

www.bananarian.net

統計検定1級対策問題集~正規分布編~

www.bananarian.net

統計検定1級対策問題集~指数分布編~

www.bananarian.net

統計検定1級対策問題集~ガンマ分布編~

www.bananarian.net