脳内ライブラリアン

脳内ライブラリアン

医療、統計、哲学、育児・教育、音楽など、学んだことを深めて還元するために。

MENU

現代数理統計学の基礎 7章 問8(3)

(3)は尤度比検定を求める問題ですね。

 

まず尤度比を求めます。公式の解答は一発の式変形で何が何やらわかりませんが、もう少し詳しく見てみます。

\frac{L(\hat\alpha, \beta_0)}{L(\hat\alpha, \hat\beta)}=\frac{\hat\alpha^n\frac{1}{\Pi x_i^2}}{(\hat\beta\hat\alpha^{\hat\beta})^n\frac{1}{\Pi x_i^{\hat\beta+1}}}

ここで、(2)の結果より

\hat\beta=\frac{n}{\sum log\frac{x_i}{X_{(1)}}}=\frac{n}{T}

となりますので、これと\hat\alpha=X_{(1)}を代入して

 

\frac{\hat\alpha^n\frac{1}{\Pi x_i^2}}{(\hat\beta\hat\alpha^{\hat\beta})^n\frac{1}{\Pi x_i^{\hat\beta+1}}}\\=X_{(1)}^{n(1-\frac{n}{T})}\Pi x_i^{-(1-\frac{n}{T})}\frac{T^n}{n^n}\\=(expT)^{-(1-\frac{n}{T})}\frac{T^n}{n^n}\\=e^{-T+n}\frac{T^n}{n^n}

 

よって、ここに−2logをつければ

-2nlogT+2nlogn-T+n\gt\chi^2_{1,\alpha}

となり、尤度比検定が求まりました。

現代数理統計学の基礎 7章 問8(2)

問8の(2)をやります。最尤推定量の問題です。

 

まず、αの最尤推定\hat\alphaは尤度関数と条件から決まります。

 

(1)で見たように、条件からαはX>αでないといけないので

X_{(1)}\gt\alpha

となります。

 

同時確率密度関数

\beta^n\alpha^{n\beta}\frac{1}{\Pi x^{\beta+1}}\Pi I_{(x_i\gt\alpha)}

ですので、αができるだけ大きければ尤度も大きくなることがわかります。よって

\hat\alpha=X_{(1)}となります。

 

2017年の統計検定1級問2(1)と似たような問題ですね。検定の問題でも気になっていたのですが、不等号が等号を含んでないのに、解答がイコールで繋いだもので良いのかがずっと気になってます。なぜこれが許されるのか誰かわかったら教えてください・・・。

 

続いてβの最尤推定量です。

 

これは対数尤度関数を微分して、\hat\alphaをαに代入します。

 

\frac{\partial}{\partial\beta}logL=\frac{n}{\beta}+nlog\alpha-\sum logx_i=0とすると

\hat\beta(nlog\alpha-\sum logx_i)=-n\\\hat\beta=\frac{n}{\sum logx_i-nlog\alpha}\\=\frac{n}{\sum log\frac{x_i}{x_{(1)}}}

となります。

現代数理統計学の基礎 7章 問8(1)

最近論文やら発表やらが溜まり気味であんまり更新する余裕がありません、、。そして、今年は人事異動で仕事も増える+神経内科専門医試験もあるので、統計検定1級に受かる自信がじわじわ削がれるわけですが。

ひとまず、あんまり細かく考えずに書ける記事を書いていきます(汗

 

7章問8はパレート分布における十分統計量、最尤推定量と尤度比検定の問題です。

 

(1)から。

 

パラメータの十分統計量を求める問題です。

 

今まで十分統計量ってなんとなく問題は解けるけど概念がよくわからなかったのですが、こちらのページの資料を見て理解がちょっと進みました。

最尤法と十分統計量について

 

どなたが書かれたのかと思ったら京都大学の医療統計の先生だったようです。

土居正明のホームページ

 

十分統計量とは要するにざっくり言えば「データXの関数のうち、対数尤度関数にしてパラメータで微分した時に消えないような部分のこと」のようです。

 

そうした視点で(1)を見てみます。

 

まず同時確率密度関数

\beta^n\alpha^{n\beta}\frac{1}{\Pi x^{\beta+1}}\Pi I_{(x_i\gt\alpha)}

となります。最後のやつは条件を満たす場合を1とする定義関数です。

 

定義関数の十分統計量については、本書中のp119の例(ページ上方です)に書いてあるように、十分統計量はどのデータXを通してαと関連するかによって決まりますので、今回は順序統計量を用いると最も小さいデータXを通してαが決まってきます。

 

よって、X_{(1)}がαに対する十分統計量です。

 

続いてβについては対数尤度関数を見てみると

L(\beta)=nlog\beta+n\beta log\alpha-\beta\sum logx_i-\sum logx_i+\sum logI_{(x_i\gt\alpha)}

となります。

 

これをβで微分したものの中で、残るもののうち、xが入っているものが十分統計量となるわけなので、答えは\beta\sum logx_iとなります。

現代数理統計学の基礎 7章 問7(3)-ガンマ分布とベータ分布の関係性-

引き続いて、問7(3)です。

 

(2)で出てきたZがどんな分布になるかという問題。

 

公式の解答では当然のごとく、ガンマ分布に従う確率変数X、Yであれば\frac{X}{X+Y}\sim B(n,m)となっています。

 

せっかくなので証明をやってみます。

 

まずは変数変換を用います。

\frac{X}{X+Y}=Z, X+Y=W

とします。

 

すると

X=WZ, Y=W(1-Z)となります。

 

ヤコビアン

J=-wz-w(1-z)\\=-w

 

よってXとYの同時確率密度関数にそれぞれを代入して

f_{w,z}\\=\frac{\theta^n}{\Gamma(n)}(wz)^{n-1}exp(-\theta wz)・\frac{\theta^m}{\Gamma(m)}\{w(1-z)\}^{m-1}exp\{-\theta w(1-z)\}w\\=\frac{\theta^{n+m}}{\Gamma(n)\Gamma(m)}w^{n+m-1}z^{n-1}(1-z)^{m-1}exp(-\theta w)

となります。かなりベータ分布らしい形になってきました。

 

最後に不要なwを消して、zの周辺確率密度関数を求めれば良いので、wで積分すると

f_z\\=\frac{1}{\Gamma(n)\Gamma(m)}z^{n-1}(1-z)^{m-1}\int_0^\infty(\theta w)^{n+m-1}\theta exp(-\theta w)dw

\theta w=tと置換すると後半部分はガンマ関数になるので

=\frac{\Gamma(n+m)}{\Gamma(n)\Gamma(m)}z^{n-1}(1-z)^{m-1}\\=\frac{1}{B(n,m)}z^{n-1}(1-z)^{m-1}

となります。

 

後の有意水準はベータ分布を使って解答通り範囲を示せばおしまいです。

現代数理統計学の基礎 7章 問7(1)(2)

代表的な確率分布の特性についてまとめおわったところで、問7です。

 

確率分布の特徴と関連するので、記事を参考にしながらご覧ください。

medibook.hatenablog.com

 

まずは(1)

指数分布の確率関数の和の分布を求める問題ですね。

 

指数分布はガンマ分布の特殊版であることに気づけば容易です。ガンマ分布は和の再生性があるため、足し合わせた結果が簡単にわかります。

 

Ex(\theta)\sim Ga(1,\theta)なので

\sum_{i=1}^n X_i\sim Ga(n, \theta)

となります。

 

※解答も本書の方もGa(1,\frac{1}{\theta})と書かれていると思いますが、二つ目のパラメータの書き方が一般的でないので、単純にθと記載します

 

同様にして

\sum_{i=1}^m Y_i\sim Ga(m, \mu)

 

 

続いて(2)。オーソドックスな尤度比検定の問題です。

 

まず尤度関数は互いに独立であるため、二つの尤度関数の積で表されます。

\theta^n\mu^mexp(-\theta\sum X_i-\mu\sum Y_i)

 

続いて、帰無仮説と対立仮説の時のMLE(最尤推定量)を求めます。

 

①まず帰無仮説が成り立つとき、尤度関数は

L=\theta^{n+m}exp\{-\theta(\sum X_i+\sum Y_i)\}

対数尤度関数にすると

logL=(n+m)log\theta-\theta(\sum X_i+\sum Y_i)

さらにパラメータθで微分して

logL'=\frac{n+m}{\theta}-(\sum X_i+\sum Y_i)

 

=0として、整理すると

\hat\theta_{MLE}=\frac{n+m}{\sum X_i+\sum Y_i}

となります。

 

②次に、対立仮説の場合

各々の最尤推定量を求めます。まず対数尤度関数は

logL=nlog\theta-\theta\sum X_i

パラメータで微分して

logL'=\frac{n}{\theta}-\sum X_i

=0として、整理すると

\hat\theta_{MLE}=\frac{n}{\sum X_i}=\frac{1}{\bar X}

\bar X=\frac{1}{n}\sum X_iとしました。

同様にして

\hat\mu_{MLE}=\frac{m}{\sum Y_i}=\frac{1}{\bar Y}

となります。

 

これで材料は出揃ったのでそれぞれの最尤推定量を尤度関数に代入し、尤度比を求めます。

 

\lambda=\frac{(\frac{1}{\bar X})^n(\frac{1}{\bar Y})^mexp\{-(n+m)\}}{(\frac{n+m}{n\bar X+m\bar Y})^{n+m}exp\{-\frac{n+m}{n\bar X+m\bar Y}(n\bar X+m\bar Y)\}}\\=\frac{(n\bar X+m\bar Y)^{n+m}}{(n+m)^{n+m}(\bar X)^n(\bar Y)^m}

 

expは綺麗に消えました。

あとは−2logをつけると、カイ二乗分布に従うことを利用して、尤度比検定は

-2log\lambda\\=-2(n+m)log(n\bar X+m\bar Y)+2(n+m)log(n+m)+2nlog\bar X+2mlog\bar Y\\=2nlog\frac{\bar X}{n\bar X+m\bar Y}+2mlog\frac{\bar Y}{n\bar X+m\bar Y}+2(n+m)log(n+m)\\=2n\frac{Z}{n}+2m\frac{1-Z}{m}+2(n+m)log(n+m)\gt \chi_{1,\alpha}^2

となります。検定が統計量Zに従うことも示せました。

 

※公式の解答はここがZ/nではなくnZになっているのですが、何回か計算しても公式の方のミスな気がしてなりません。わかる方がいたらご教示ください。

 

(3)は次回に回します。

実際の医学論文から統計を学んでみるⅣ-RMST法(Restrictive Mean Survival Time)-

昨年11月NEJM誌に出た「生体弁に対するDOAC vs ワーファリン」のランダム化比較試験を読みました。

Rivaroxaban in Patients with Atrial Fibrillation and a Bioprosthetic Mitral Valve

 

ここで使われていた生存時間解析の方法はいつものCox比例ハザード回帰ではなく、RMST法なる方法でした。

 

数年前から腫瘍や循環器系などの論文でも使われつつある方法なようで、今後も使われるものは増えていきそうです。ちょっと調べてみた内容をまとめてみます。

 

案の定素人なので疑問点や異なっている点があったら教えてください。

 

 

1、Cox比例ハザード性の問題点

生存時間解析といえば

①カプランマイヤー法で生存曲線を書き

②ログランク検定で仮設検定を行い

③Cox比例ハザード回帰でハザード比を出す

というのが王道な方法でした。

 

ところが、以前から指摘されているように、Cox比例ハザードモデルはそのモデルや解釈に色々問題点があります。

 

問題点①比例ハザード性の仮定は成り立つか

 

そもそも比例ハザード性とは「介入群の人のある瞬間のイベント発症率(ハザード)とコントロール群の人のある瞬間のイベント発症率が常に一定の比に保たれている。」という前提に立ちます。

 

このハザードは時間tに依存する関数であるため、刻一刻と変動します。介入群とコントロール群でその比が常に一定だと仮定するわけです。さて、こんなことはあり得るのでしょうか。

 

具体的な数値で考えてみると、このような例になります。例えば、2群における疾患の発症率が1年目と2年目で以下のような関係が成り立つわけです。

 

f:id:medibook:20210109062329j:plain

 

常に瞬間的な発症率が比例関係にある、というわけです。実際の臨床医学ではそうではない例が多く見受けられます。

 

例えば、脳梗塞に対するスタチンでの再発予防を考えてみましょう。おそらくは内服直後にはまだ有効性を発揮しないため、大きな効果が見られない可能性があり(急性期にも良いという説もある気はしますがそこはさておき)、後からじわじわ差が開いてくる可能性があります。

 

もしそうだとすると2群の1年目、2年目の発症率は以下のような違いを見せます。

 

f:id:medibook:20210109062624j:plain

こうなると、比例ハザード性は保たれていないわけです。このような比例ハザードが保たれない例はJAMAの統計コーナーでも取り上げられており、以前一度記事にまとめました。記事の終わりの方に書いてあります。

Cox比例ハザード回帰モデルについて数式ありで、できるだけわかりやすくまとめる - 脳内ライブラリアン

 

具体的には

①No immediate effect

薬がすぐに効かない場合、上のスタチンの例と同じ

②Immediate and Delayed Effects in Opposite Directions

介入が検査などの場合、検査が有効なら先にあらかた疾患が発見されるので、イベントは増えるが、残った人は見つかりにくくなるので、イベントが減る

③Variation in Disease Susceptibility

介入が効かないタイプの人は早期にイベントを起こし、残ったよく効く人たちはイベントをずっと起こさない

 

の三つが挙げられています。

 

問題点②ハザードの解釈がしにくい

 

各群で出てくる「ハザード」の解釈がしづらいのがもう一つの問題です。上で述べたようにハザードは時間に応じて変動するため、仮に比例ハザード性が成り立って、ハザード比が一定であったとしても、相対リスクは分かりますが、絶対リスクの差は分からないため、臨床的な意味が取りづらいです。

 

例えばハザード比が0.5と出ていたとしても、ハザード1%に対しての0.5%とハザード20%に対しての10%では全然意味合いが違ってきます。

 

また、他に解釈できる数値として生存時間中央値を得ることもできますが、これはイベントの発症が50%まで起こらないと得られないため、数値が得られない場合もあります。

 

以上のことがCox比例ハザード回帰で問題となっていたことでした。それに対してのRMST法という視点で見ていきます。

 

2、RMST法(Restrictive Mean Survival Time)とは

RMST法とは簡単にいえば、生存曲線の曲線下面積を比較する方法です。生存曲線の推定はいつものカプランマイヤー法などで行われます。

 

f:id:medibook:20210109063800j:plain

ある観察期間τ(これは試験デザインで恣意的に決める)までの生存曲線を確認し、介入群とコントロール群におけるこの面積の差を検定していくわけです。

 

確かに治療成績が良ければ、面積は大きくなるわけで、直感的にも納得できるなあと思います。しかも比例ハザード性のような無理のある仮定は不要です。

 

そして、解釈の面でも分かりやすいことがあります。この面積は実はある期間内τにおける平均生存期間を表しています。この理由は数式をみて後述しますが、数式みたくない人はスルーしてください。

 

例えば今回紹介する論文では365日の観察期間で、治療群の平均生存期間(平均的な患者がイベントを起こさない日数)は347.5日でした。つまり多分1年は大丈夫かギリギリぐらい、ということになります。

 

群比較がなくても、群ごとの結果の意味も分かりやすいのはCox比例ハザードと異なる点になります。前述のようにハザード比がわかっても、元のコントロール群のハザードが想定できなければ、具体的な治療効果は想定しにくいですし、介入群だけの数値を見てもどの程度大丈夫なのかは分かりません。

 

3、なぜ曲線下の面積がRMST(ある期間における平均生存期間)になるのか

この項目は数式を用いた話になるので、苦手な方は飛ばしてください。

 

曲線下の面積がなぜRMSTと一致するのでしょうか。日本製薬工業協会がRMST法について素晴らしい資料をまとめてくださっているのでそれを参考にします。(*1 p.79)ページ数多すぎて論文をどう読むかを知りたいだけの人には向かなさそうですが、、、。

 

まず「ある決められた観察期間のうちに観察された生存期間の期待値」(つまりRMSTそのものですね)を数式で定義します。イベントが起きた人の生存期間をT、決められた観察期間をτ(タウ)とすると、以下のようになります。

 

E[min(T, \tau)]

 

min(T, τ)というのはどちらか小さい方が採用されるよ、という意味ですね。

 

続いてはこれを積分で表します。min(T, τ)については場合わけをすることでうまく表現できます。

①T<τのとき

E[min(T, \tau)]=\int_0^\tau tf(t)dt

②T>τのとき

E[min(T, \tau)]=\int_\tau^\infty \tau f(t)dt

となります。なお、f(t)は生存時間Tの確率密度関数です。ハザード関数とはまた別物になります。f(t)と生存関数S(t)の関係性は以下のようになります。

\int f(t)dt=F(t)\\1-F(t)=S(t)

 

先程の①、②の場合をまとめると

 

E[min(T, \tau)]=\int_0^\tau tf(t)dt+\int_\tau^\infty \tau f(t)dt\\=[tF(t)]_0^\tau-\int_0^\tau F(t)dt+[\tau F(t)]_\tau^\infty(第1項に部分積分を使った)\\=\tau F(t)-\int_0^\tau(1-S(t))dt+\tau(1-F(t))(先程の関係式を使った)\\=-\int_0^\tau1dt+\int_0^\tau S(t)dt+\tau\\=\int_0^\tau S(t)dt

 

となります。なんとτから先の積分が消失している!というのは俄かに信じ難いですが、以上の式変形によって平均生存時間=ある期間τ内の生存関数の曲線下面積ということが導き出せました。

 

4、RMST法の特徴と注意点

RMST法の特徴としては先に述べたように、群単体の結果でも理解しやすく、差もわかりやすい点です。差を見れば平均的にイベントの発生をどれだけ遅らせることができるのかが分かります。

 

ただここで出てくる数値は決められた観察期間の影響を受けるので、観察期間が臨床的に妥当な期間がどうかが重要となります。例えば脳梗塞の試験で観察期間1ヶ月だとか、あまりに期間が短いと差が当然出ないので、意味がなくなってしまいます。

 

また、同じ観察期間の研究であれば患者層も同じとすると、研究間での数値比較も容易であるため、なおさら観察期間の統一化が大事となります。

 

5、実際の論文をみてみる

さて、ようやく初めに戻って、NEJMに発表された研究を見てみます。

 

生体弁に対するDOACの使用効果を調べたRCTでRIVER trialとも呼ばれているようです。それなりの規模を集めてますが、それでもまだphase2 trialです。

 

リウマチ熱による弁置換の患者が多いという理由でブラジルで試験は行われています。COIには当然バイエル(リバロキサバンの会社)が絡んでます。

 

試験デザインは多施設ランダム化オープンラベル(アウトカムの判定は盲検化)試験です。そして非劣性試験です。ワーファリンを使用すると調整がいるのでどうしてもオープンラベルになりがちですね。

 

PICOを確認すると

P:心房細動+生体弁で血栓塞栓症予防に抗凝固薬を使用していた(もしくは使用する予定の)18歳以上の成人 合計1005人

 

I:リバロキサバン 20mg/day(腎機能に応じて減量)

 

C:INR 2-3のワーファリン、INRは4週ごとに測定

 

O:複合アウトカム<死亡、心血管イベント(脳卒中TIA、valve thrombosis、全身塞栓症、心不全入院)、重大な出血>

 

そしてこの複合アウトカムの評価をRMST法で行っています。観察期間は1年です。

 

結果としては

リバロキサバン群 347.5日vs ワーファリン群 340.1日

差は7.4日(95%CI -1.4~16.3; 非劣性の基準でp<0.001)

でした。こんな感じで日数で差がわかるのは分かりやすくて面白いところですね。1週間と言われると、まあ差があまりないのかなあという印象を持つかと思います。

 

で、ただこの試験で気になるのはやはり複合アウトカムの部分。アウトカムの中に効果と副作用がごちゃごちゃになっている(出血も含んでいる)点がまず気になります。phase 2 trialで参加者も少ないし、そういうものかなとも思わないでもないですが。

 

さらに気になるのはなぜ心不全入院がアウトカムの一つに入っているのか。プロトコールを見ても明確に書いていません。薬剤の直接的な効用として心不全を予防しうるのはvalve thrombosisが絡む場合ぐらいに見えるのですが、、、。しかも結果としては心不全入院のアウトカムが多いんです。全死亡(リバロキサバン17例、ワーファリン26例、以下同順)や大出血(7例、13例)と同等かそれ以上のアウトカム数(22例、19例)なんです。

 

流石に試験開始後に追加されたものではないですが、プロトコールをみると第5版→第6版で新たに追加されています。ううん、、、非劣性を確実にするために結果を水で薄めた(おそらく差がない心不全入院のイベントで数を稼ぐ)ように見えるのですが、、、。本文中にも記載がないのですが、ここに心不全入院が入る理由がわかる循環器内科に詳しい方がいらっしゃったら教えてください、、。

 

あとは観察期間の1年って妥当なのかどうなのかというところですね。心血管イベントに関してならこれぐらいなのかもしれません。脳梗塞の予防という点でいくと過去のリバロキサバンvsワーファリンの試験であるROCKET-Afなんかは中央値700日ほどは観察していたのでもう少し長くても良い気はします。

 

ただあくまで、phase 2なので、今後より詳しくちゃんとした効果はわかってくるのかもしれません。

 

6、まとめ

・RMST法は、「生存曲線下の面積=平均生存期間」を利用した比較的新しい生存時間解析の方法である

 

・Cox比例ハザードモデルより臨床的な解釈がしやすい、また比例ハザード性の仮定を置かなくても良い

 

・観察期間の決め方が臨床的に妥当かどうかが重要である

 

 

神経内科疾患の論文ではまだあまり見かけていない印象ですが、今後各分野に応用が進みそうですね。

 

7、参考文献

*1

生存時間型応答の評価指標-RMST (restricted mean survival time) を理解する-|治験に関する医薬品評価委員会の成果物|日本製薬工業協会

臨床医が論文を読むために必要な知識はp.13, p.43の表に集約されている気がしますので、そこを読めば良さそうです。

 

*2 『生存時間型変数の要約指標と平均値の推定』

聖マリアンナ医科大学の先生がまとめられています

http://igakukai.marianna-u.ac.jp/idaishi/www/462/46-2-05-Inoue.pdf

代表的な確率分布を覚えやすいようにまとめてみる②-連続型・標本分布-【統計検定1級対策】

前回に引き続いて統計検定1級に出てきやすい連続型確率分布・標本分布について、関係性と式を簡単にまとめてみます。

 

前回の記事はこちら

medibook.hatenablog.com

 

 

 

図にしてみる

 

f:id:medibook:20210109054624j:plain

 

 

各分布の式など

モーメント母関数

M_X(t)=E[e^{tX}]\\M'_X(t)=E[X]\\M''_X(t)=E[X^2]

 

①連続一様分布

f(x;a,b)=\frac{1}{b-a} (a\leq x\leq b)

E[X]=\frac{a+b}{2}

V(X)=\frac{(a-b)^2}{12}

 

②コーシー分布

f(x)=\frac{1}{\pi}\frac{1}{x^2+1}

*自由度1のt分布と等しい

*平均と分散、モーメント母関数は存在しない

 

正規分布

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

*超メジャー分布でありながら、ぱっと見が分かりにくい式ですが、expの手前にあるのは規格化(正規化)定数であり、expの部分を積分したときに、調整のため生じうるものと理解すると何となく意味が掴みやすいです。

 

f:id:medibook:20210109054109j:plain

 

M_X(t)=exp(\mu t+\frac{\sigma^2t^2}{2})

E[X]=\mu

V(X)=\sigma^2

 

N(\mu_X, \sigma^2_X)+N(\mu_Y, \sigma^2_Y)=N(\mu_X+\mu_Y, \sigma^2_X+\sigma^2_Y)

分布の再生性がある

 

④対数正規分布

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

Y\sim(\mu, \sigma^2)の時、exp(Y)=Xと変数変換することで生み出される確率分布

*モーメント母関数は存在しない

 

E[X]=exp(\mu+\frac{1}{2}\sigma^2)

V(X)=exp(2\mu+\sigma^2)\{exp(\sigma^2)-1\} 

*変数変換とガウス積分を使って導出できる

 

⑤多変量正規分布

相変わらず行列が苦手すぎるのと文量も増えそうなのでいったん割愛します

 

⑥ガンマ分布

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

(λでなくβをパラメータとして持ってくる式が多いですが、指数分布と比較しやすくなるようにλにしました。なお、愛用の『現代数理統計学の基礎』ではβの逆数をパラメータにしているため混乱します。)

*次に出てくる指数分布を一般化したもの、としてみると複雑な式も覚えやすい。α=1としたものが指数分布と一致する。

*つまり、ガンマ分布はある一定数の時間においてイベントが複数回起きる確率分布として理解できる。(指数分布は「単位時間あたりλ回起きるイベントが期間xのうちに1回起きる可能性を示す分布」であるので、ガンマ分布は「単位時間あたりλ回起きるイベントが期間xのうちにα回起きる可能性を示す分布である」と言える)

 

以下のブログの説明がわかりやすかったです。

ガンマ分布のはなし - 統計学といくつかのよしなしごと

 

M_X(t)=(\frac{\lambda}{\lambda-t})^\alpha

*定義通りに計算して、ガンマ分布の総和が1になることを利用して導出

E[X]=\frac{\alpha}{\lambda}

V(X)=\frac{\alpha}{\lambda^2} 

*モーメント母関数から導出可

 

Ga(\alpha_1, \lambda)+Ga(\alpha_2, \lambda)=Ga(\alpha_1+\alpha_2, \lambda)

分布の再生性がある

*上述したガンマ分布の意味を考えるとαはバラバラでも良いが、λは一緒でないといけない理由がよくわかる

 

<補足>

ガンマ関数とは→「階乗」を一般化したもの

定義は\Gamma(\alpha)=\int_0^\infty x^{\alpha-1}e^{-x}dx

性質として

\Gamma(\alpha+1)=\alpha\Gamma(\alpha)

\Gamma(\frac{1}{2})=\sqrt\pi, \Gamma(1)=1をもつ。

一つめの性質から、αが自然数を想定すれば階乗になることがわかる。

 

⑦指数分布

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

*ガンマ分布におけるα=1の場合となる

M_X(t)=\frac{\lambda}{\lambda-t}

E[X]=\frac{1}{\lambda}

V(X)=\frac{1}{\lambda^2}

 

*指数分布はハザード関数と関連しており、関係性を示す式変形については以下の記事の一部で以前にまとめました

medibook.hatenablog.com

 

⑧ワイブル分布

f(x; a,b)=abx^{b-1}exp(-ax^b) ただし、x\gt0

*指数分布においてイベントの発生確率と言えるλが、時間経過に伴って増えていくという仮定のもとに\lambda(x)=abx^{b-1}として得られる確率密度関数。λがずっと一定というのは確かに現実的ではないことの方が多いと思います。

*導出過程は以下

f(x)をハザード関数とすると

\lambda(x)=\frac{f(x)}{1-F(x)}

(ここまでの流れは上の記事と同様です)

両辺を積分して

\int_0^x\lambda(t)dt=\int_0^x\frac{f(t)}{1-F(t)}dt

これを変形すると

f(x)=\lambda(x)exp\{-\int_0^x\lambda(t)dt\} 

ここに\lambda(x)=abx^{b-1}を代入すると最初に出てきた確率密度関数の式が得られます。

 

平均と分散の計算は別記事にしました↓ 

medibook.hatenablog.com

 

カイ二乗分布

自由度をnとすると

f(x)=\frac{1}{\Gamma(\frac{n}{2})}(\frac{1}{2})^{\frac{n}{2}}x^{\frac{n}{2}-1}exp\{-\frac{x}{2}\}

*ガンマ分布の特殊形です。自由度nのカイ二乗分布

Ga(\frac{n}{2}, \frac{1}{2})と一致します。

 

*n個の標準正規分布の和は自由度nのカイ二乗分布に従います。過去記事で書きました。

medibook.hatenablog.com

 

E[X]=n

V(X)=2n

*ガンマ分布と同様に考えれば簡単に求められます

Ga(\frac{n}{2}, \frac{1}{2})+Ga(\frac{m}{2}, \frac{1}{2})=Ga(\frac{n+m}{2}, \frac{1}{2})

分布の再生性がある

 

⑩t分布

f(x)=\frac{\Gamma(\frac{m+1}{2})}{\Gamma(\frac{m}{2})}\frac{1}{\sqrt{\pi m}}\frac{1}{(1+\frac{x^2}{m})^{\frac{m+1}{2}}}

モーメント母関数はなし

E[X]=0

V(X)=\frac{m}{m-2}

ただし1≦m<2の時は無限大

 

定義としては

Z\sim N(0,1), U\sim\chi^2_mのとき

T=\frac{Z}{\sqrt{\frac{U}{m}}}

 

⑪F分布

自由度(m,n)のF分布の確率密度関数

 f(x)=\frac{\Gamma(\frac{m+n}{2})}{\Gamma(\frac{m}{2})\Gamma(\frac{n}{2})}\frac{(\frac{m}{n})^\frac{m}{2}x^{\frac{m}{2}-1}}{(1+(\frac{m}{n})x)^{\frac{m+n}{2}}}

 

E[X]=\frac{n}{n-2}

V(X)=\frac{2n^2(m+n-2)}{m(n-2)^2(n-4)}

定義としては

S\sim\chi^2_m, T\sim\chi^2_nのとき

X=\frac{\frac{S}{m}}{\frac{T}{n}}

 

*期待値・分散の導出の方法は以下のサイトが参考になります。定義に戻ってカイ二乗分布の期待値をもとに考えると、E[\frac{1}{T}]を考えればなんとかなることが分かります。自由度のm,nが当ブログと逆なのでご注意ください。

F分布の期待値・分散をカイ二乗分布を用いて導出 | AVILEN AI Trend

 

⑫ベータ分布

f(x; a,b)=\frac{1}{B(a,b)}x^{a-1}(1-x)^{b-1}

*ベータ関数の定義は

B(a,b)=\int_0^1x^{a-1}(1-x)^{b-1}dx

またガンマ関数と以下の式のような関連があります。

B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}

ベータ関数とガンマ関数の関係性の証明は以下にまとめました。 

medibook.hatenablog.com

 

E[X]=\frac{a}{a+b}

V(X)=\frac{ab}{(a+b)^2(a+b+1)}

*導出はややトリッキーでベータ関数の以下の性質が利用されます。

B(a+1, b)=\frac{a}{a+b}B(a, b)

上記の性質は以下のサイトが分かりやすいです。部分積分を使って変形していくことで導出します。

ベータ関数とは? ~ 性質と公式 ~ (証明付) - 理数アラカルト -

期待値と分散の詳しい導出方法は以下のサイトが参考になります。ただサイトの中のパラメータα'、βの「ガンマ分布」と書いてある部分は「ベータ分布」の間違いだと思います。

ベータ分布の期待値・分散の導出 | AVILEN AI Trend

 

参考文献: