3章後半 代表的な確率分布(後半:連続型確率変数の確率分布)
この記事はfMRILab Advent Calendar 2021の記事です!
弊ラボにおける「現代数理統計学の基礎」(久保川,2017)の勉強会資料として作成したものを公開します。
連続分布
確率変数が連続型である場合を扱います
適宜こちらのコードを使って遊びます。
3.2.1 一様分布
確率変数
で与えられることを言う。
平均や分散については次に示す通り。
-
導出
\begin{aligned} Var(X) &= E[X^2]- \{E[X]\}^2 \\ &=\frac{a^2+ab+b^2}{3}-\frac{(a+b)^2}{4}\\ &=\frac{(a-b)^2}{12} \end{aligned}
特に
一様分布の使い道
-
逆関数サンプリング法
wikipediaより(https://upload.wikimedia.org/wikipedia/commons/thumb/2/2d/Inversion_method.svg/450px-Inversion_method.svg.png)
3.2.2 正規分布 ←←めちゃつかう
正規分布は、
- 理論的に扱いやすい分布
- 平均中心で対称な釣り鐘型
- 中心極限定理により標本平均の極限分布になっている
などのことから数理統計学において最も重要な分布の一つ。
確率変数
で与えられることを言い、この分布を
標準化(p.18)して(
と書ける。これは平均
これの分布関数は
で表され、この値は本の巻末で数字の表として与えられることが多い。
確率密度関数の定義から、
が成り立つ。(ガウス積分)
また、
命題3.9
-
[証明]
積率母関数の定義より
\begin{aligned} M_Z(t) &= E[e^{tZ}] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{tz-z^2 /2}dz \\ &= e^{t^2 / 2}\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-(z-t)^2/2}dz = e^{t^2/2} \end{aligned} \begin{aligned} M_Z(t) &= E[e^{tZ}] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{tz-z^2 /2}dz \\ &= e^{t^2 / 2}\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-(z-t)^2/2}dz = e^{t^2/2} \end{aligned} となる。
より、M'_Z(t)=t\exp\{t^2/2 \},M''_Z(t)=(1+t^2)\exp\{t^2/2 \} \begin{aligned} E[Z] &= M'_Z(0)=0\\ Var(Z) &= E[Z^2] = M'_Z(0) = 1 \end{aligned} となる。特性関数については
より求まる。\varphi_Z(t)= M_Z(it) なので、命題3.9を用いるとX=\sigma Z+\mu \begin{aligned} E[X] &= \sigma E[X] + \mu = \mu, \\ Var(X) &= E[\left\{(\sigma Z + \mu ) - \mu \right\}^2 ] = \sigma^2 E[Z^2] = \sigma^2,\\ \varphi_Z(t) &= E[e^{it(\sigma Z + \mu)}] = e^{it\mu} E[e^{(it\sigma)Z}]= e^{it\mu- (t\sigma^2)/2} \end{aligned}
命題3.10
正規分布
で与えられる。
と表される。したがって、
命題3.11
-
[証明]
の特性関数は、aX+b より、\varphi_X(t) = \exp\{\mu it - \sigma^2 t^2 /2 \} \begin{align} \varphi_{aX+b}(t) &= E[e^{(aX+b)it}] = e^{bit} E[e^{(ait)X}] \\ &= e^{bit} e^{\mu ait - \sigma^2 a^2 t^2 /2} = e^{(a\mu + b)it -a^2 \sigma^2 t^2 /2} \end{align} となる。
したがって、定理2.16(特性関数から確率を求めるやつ)と命題3.10より、
.aX+b \sim \mathcal{N}(a\mu + b,a^2 \sigma^2)
3.2.3 ガンマ分布とカイ二乗分布
非負の実数直線上の確率分布の中で代表的な分布がガンマ分布。この分布には特別な場合として指数分布やカイ二乗分布という分布を含んでいて、あとで示すようにカイ二乗分布は正規分布と関連する分布である。
確率変数
で与えられるとき、
ここで
で与えられる。
尺度変換
となり、ガンマ関数の定義から
参考:https://mathtrain.jp/gammadist
ガンマ分布は単位時間あたりの回数が
のイベントが 1/\beta 回起こるまでの時間間隔を表す。 \alpha https://takuyaokada.hatenablog.com/entry/20140812/1407837832
[別の表現]
期間
ごとに1回くらい起こるランダムな事象が \beta 回起こるまでの時間の分布 n
https://manabitimes.jp/math/1386
命題3.12
(3.17)に従う確率変数
-
[証明]
まず積率母関数
\begin{aligned} M_Y(t) &= E[e^{tY}] = \frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}y^{\alpha-1}e^{-(1-t)y}dy\\ &= \frac{1}{(1-t)^{\alpha}} \frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}(1-t)((1-t)y)^{\alpha - 1}e^{-(1-t)y}dy = \frac{1}{(1-t)^{\alpha}} \end{aligned} より、M'(t)=\alpha (1-t)^{-\alpha -1},M''(t) = \alpha (\alpha+1)(1-t)^{-\alpha -2}
(p.21参照)となり、E[Y] = M'_Y(0)=\alpha, E[Y^2]=M''_Y(0)=\alpha(\alpha+1) \begin{align} Var(Y) &= E[Y^2] - \left\{ E[Y] \right\}^2 \\ &= \alpha(\alpha+1) - \alpha^2 \\ &= \alpha \end{align} となる。特性関数は同様にp.21より
\varphi_Y(t) = M_Y(it) = \frac{1}{(1-it)^\alpha} 確率変数
に従う形式に戻すと、X より、X=\beta Y E[X] = \alpha \beta,Var(X) = \beta^2 Var(Y) = \alpha \beta^2 となり、特性関数は
\varphi_X(t) = E[e^{it\beta Y}] = (1-\beta it)^{-\alpha} \tag{3.18}
命題3.13
-
[証明]
(1)について
標準正規分布
の確率密度関数\mathcal{N}(0,1) に平方変換(p.25)を行う。\phi(z) \begin{aligned} 1 &= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-z^2/2}dz = \frac{1}{\sqrt{2\pi}}\int_0^{\infty}\frac{1}{\sqrt{y}}e^{-y/2}dy \\ &= \frac{\Gamma(1/2)}{\sqrt{\pi}} \int_0^{\infty}\frac{1}{\Gamma(1/2)}\frac{1}{2}\left( \frac{y}{2}\right)^{1/2-1}e^{-y/2}dy = \frac{\Gamma(1/2)}{\sqrt{\pi}}. \end{aligned} ※ほかの方法
- ベータ関数を用いたもの
- オイラーの反射光式を用いたもの
- ガンマ関数の定義からごり押しで証明したもの
etc...https://proofwiki.org/wiki/Gamma_Function_of_One_Half
また、
(2)\Gamma(1) = \int_{0}^{\infty} e^{-x} dx =\left[-e^{-x} \right]_{0}^{\infty} = 1 部分積分を用いる。
\begin{aligned} \Gamma(\alpha + 1) &= \int_0^{\infty} x^{\alpha} e^{-x} dx = \int_{0}^{\infty} x^{\alpha}(-e^{-x})'dx \\ &=\left[-x^{\alpha} e^{-x} \right] _0^{\infty} - \int_0^{\infty} (x^{\alpha})' (-e^{-x})dx \\ &= \alpha \int_0^{\infty} x^{\alpha -1} e^{-x}dx = \alpha \Gamma(\alpha). \end{aligned}
ガンマ分布の特殊な場合として、指数分布や '自由度nのカイ二乗分布' が含まれ、それぞれ
正規分布のベイズ推測において、ガンマ分布は正規分布の精度パラメータ
定義3.14
で与えられるとき、自由度
で与えられ、また、
となる。
-
[証明]
特性関数とモーメントに関する式から
\begin{aligned} E[\chi_n^2] &= \left. \frac{1}{i}\frac{d}{dt}\varphi_{\chi_n^2}(t) \right|_{t=0} = \frac{1}{i}\varphi '_{\chi_n^2}(0) \\&= \frac{1}{i}\cdot \frac{-n}{2}\cdot(-2i) = n \end{aligned} 同様に
\begin{aligned} E[(\chi_n^2)^2] &= \left. \frac{1}{i^2}\frac{d^2}{dt^2}\varphi_{\chi_n^2}(t) \right|_{t=0} = \frac{1}{i^2}\varphi''_{\chi_n^2}(0) \\ &= \frac{1}{i^2}\cdot(-2i)\cdot (-n/2)\cdot (-2i)\cdot (-n/2-1) \\ &= n(n+2) \end{aligned} なので、
\begin{aligned} Var(\chi_n^2) &= E[(\chi_n^2)] - \left\{E[\chi_n^2] \right\} \\ &=n(n+2) - n^2 = 2n. \end{aligned}
命題3.15(この先結構出てくるので覚えておくこと!)
確率変数
-
[証明]
標準正規分布の確率密度関数
において、\phi(z)=\frac{1}{\sqrt{2\pi}}\exp\left\{-z^2/2 \right\} という変換を施すと、Y=Z^2 の累積分布関数はY \begin{aligned} F_Y(y) &= P(Y<y) \\ &= P(Z^2<y) \\ &= P(-\sqrt{y}< Z < \sqrt{y}) &= \int_{-\sqrt{y}}^{\sqrt{y}}\phi(z)dz \end{aligned} 累積分布関数の微分によって確率密度関数を得られることから
\begin{aligned} f_Y(y) = \frac{dF_Y(y)}{dy} &= \phi(\sqrt{y})\frac{d}{dy}(\sqrt{y}) - \phi(-\sqrt{y})\frac{d}{dy}(-\sqrt{y})\\ &= \phi(\sqrt{y})\frac{1}{2\sqrt{y}} + \phi(-\sqrt{y})\frac{1}{2\sqrt{y}}\\ &= \frac{1}{\sqrt{2\pi}}\exp(-y/2)\frac{1}{2\sqrt{y}} + \frac{1}{\sqrt{2\pi}}\exp(-y/2) \frac{1}{2\sqrt{y}}\\ &= \frac{1}{\sqrt{2\pi}}y^{-1/2}\exp(-y/2) \\ &= \frac{1}{\Gamma(1/2)}\left(\frac{1}{2}\right)^{1/2}y^{1/2-1} \exp(-y/2) \\ \end{aligned} 式(3.21)
において、n=1となっている場合に相当!したがって自由度1のカイ二乗分布に従う。f_X(x) = \frac{1}{\Gamma(n/2)}\left(\frac{1}{2}\right)^{n/2}x^{n/2-1}\exp\{-x/2\}
参考:https://stats.biopapyrus.jp/probability/chi-squared-distribution.html
カイ二乗分布のさらなる性質は次章以降で取り扱う。
3.2.4 指数分布とハザード関数
ある期間に平均して
で与えられる。
具体例として、webサイトの訪問者、料金所に到着する車と行った事象間の時間分布をモデル化することができる。
分布関数は
指数分布はガンマ分布の特別な場合で
となる。
実際にグラフを作ってみるとたしかに指数分布とガンマ分布の関係性がわかる。
また、
指数分布の性質と無記憶性
指数分布は生存時間の分布として用いられることがある。
となり、それまで生存してきた時間には依存しないことが分かる。
換言すると、指数分布では将来の事象発生までの時間がその過去の経過に依存しない。
幾何分布の時と同様に無記憶性が指数分布についても成立している。これは、指数分布において故障(or死亡)がランダムに起こることを意味し、これに関連した考え方がハザード関数。
ハザード関数
となる。このとき、両辺を
となる。これは、
と書いて、ハザード関数(故障関数)という。
指数分布(3.23)については、常に
-
式
\lambda(x) = \frac{f(x)}{1-F(x)} = \frac{\lambda\exp(-\lambda x)}{1- (1-\exp(-\lambda x))} = \lambda
すなわち、
(3.24)の両辺を積分すると
となるので、
が得られる。非負の連続型確率変数の分布はハザード関数によって特徴づけられることが分かる。
例えば
指数分布(あるいはポアソン分布)を用いたシミュレーションにおけるキーポイントは
例えば交通量やネットワークトラフィックなどは同じ1日という区間でも曜日によって異なるであろうことは容易に想像が出来る。それでもポアソン分布や指数分布を用いたシミュレーションが成立するのは、通常、時間も空間も一様な区間に分割でき、その区間の中だけでなら妥当な仮定となるかららしい。^1
では壊れる可能性が経時的に変化していくような場合を考えてみる。例えば、時間の経過とともに故障しやすくなる場合や故障しなくなる場合には
より、得られる確率密度関数は
となる。これは、ワイブル分布(Weibull distribution)といい、生存解析の分野で基本となる分布。
機械などが壊れる、劣化などの事象がお起こる確率を近似する際に用いられる。
https://recruit.gmo.jp/engineer/jisedai/blog/survival_analysis_covid19/
3.2.5 ベータ分布
ベータ分布は区間
で表される。ここで、
である。ベータ関数はガンマ関数を用いて
なる関係式が成り立つ。(次章の(4.17)で示されるらしい)
これを使って
ベータ分布は次章の例4.14 で扱うベータ・二項分布やベイズ推測における事前分布として用いられる。
-
証明
ベータ分布の確率密度関数は次のように簡略化して記述できる。
-
が自然数の場合a,b \begin{aligned} E[x] &= \int_{0}^{1}xp(x)dx \\ &= \frac{\int_{0}^{1}x^{a}(1-x)^{b-1}dx}{\int_{0}^{1}x^{a-1}(1-x)^{b-1}dx} \end{aligned} 分母と分子はそれぞれ非負の整数m,nに関するベータ関数の積分公式
\int_{0}^{1} x^m (1-x)^n dx = \frac{m!n!}{(m+n+1)!} を用いて変形すると
E[x] =\frac{a!(b-1)}{(a+b)!}\cdot\frac{(a+1-1)!}{(a-1)!(b-1)!}=\frac{a}{a+b} 分散については
E[X^2] - (E[X])^2 = \frac{\int_{0}^{1}x^{a+1}(1-x)^{b-1}dx}{\int_{0}^{1}x^{a-1}(1-x)^{b-1}dx}-\frac{a^2}{(a+b)^2} となる。第1項はベータ関数の積分公式を用いてさらに書き直せる。
\begin{aligned} Var(x)&= \frac{(a+1)!(b-1)!}{(a+b+1)!}\cdot \frac{(a+b-1)!}{(a-1)!(b-1)!}-\frac{a^2}{(a+b)^2}\\ &=\frac{a(a+1)}{(a+b)(a+b+1)}-\frac{a^2}{(a+b)^2}\\ &=\frac{a(a+1)(a+b)-a^2(a+b+1)}{(a+b)^2 (a+b+1)}\\ &=\frac{ab}{(a+b)^2 (a+b+1)} \end{aligned} -
が正の実数である場合a,b
先ほどのように階乗は使えなくなるので、次のベータ関数 の性質をうまく使う。B(a,b) B(a+b) = \int_0^1 x^{a-1}(1-x)^{b-1}\\ B(a+1,b) = \frac{a}{a+b} B(a,b) これらを用いて
\begin{aligned} E[x] &= \int_{0}^{1}xp(x)dx \\ &= \frac{\int_{0}^{1}x^{a}(1-x)^{b-1}dx}{B(a,b)}\\ &=\frac{\int_{0}^{1}x^{a'-1}(1-x)^{b-1}dx}{B(a'-1,b)} \end{aligned} ただし、
とした。このとき、a=a'-1 B(a'-1,b) = \frac{a' + b-1}{a'-1} B(a',b) より、上式は
E[x] = \frac{a'-1}{a'+b-1}\frac{\int_{0}^{1}x^{a'-1}(1-x)^{b-1}dx}{B(a',b)}=\frac{a}{a+b} ここで規格化条件より
であることを用いた。\frac{\int_{0}^{1}x^{a'-1}(1-x)^{b-1}dx}{B(a',b)}=1 分散も同様にして
の計算で出てくるベータ関数をうまく変形していって、最後は規格化条件を用いることで導出できる。E[X^2]
長いので割愛。以下参照。
https://ai-trend.jp/basic-study/beta-distribution/b-parameter-derivation-2/
-
3.3 発展的事項
3.3.1 スタインの等式
次式で示されるスタインの等式(Stein identity)は、モーメントの計算などに役立つらしい
命題3.16
-
[証明]
(3.20)式を用いて
E[(X-\mu)g(X)] = \frac{1}{\sqrt{2\pi}\sigma} \int_{-\infty}^{\infty}g(x)(x-\mu)e^{-(x-\mu)^2/(2\sigma^2)}dx とすると、f'(x)=(x-\mu)e^{-(x-\mu)^2/(2\sigma^2)} として部分積分を用いる。f(x)=-\sigma^2e^{-(x-\mu)^2/(2\sigma^2)} より、f(x) = -\sigma^2 \exp({-(x - \mu)^2/(2\sigma^2)}) \begin{aligned} E[(X-\mu)g(X)]=&\frac{1}{\sqrt{2\pi}\sigma}\left[-\sigma^2 e^{-(x-\mu)^2/(2\sigma^2)}\right]_{-\infty}^{\infty}\\ &+ \frac{\sigma^2}{\sqrt{2\pi}\sigma}\int_{-\infty}^\infty g'(x)e^{-(x-\mu)^2/(2\sigma^2)}dx \end{aligned} となる。右辺第1項は0に収束して、第2項は
のg'(x) に関する期待値に\mathcal{N}(\mu,\sigma^2) をかkたものなので\sigma^2 E[(X-\mu)g(X)]=\sigma^2E[g'(X)]
スタインの等式(3.31)は縮小推定の性質(この場合lassoっぽい?理解が浅いので曖昧です)を議論するときに重要な道具となる。
-
詳細
正則化パラメータの選択方法として本書では3つの手法が紹介されている
- 自由度調整済み決定係数
- 交差検証(よく使うやつ)
- モデルの自由度(モデルのパラメータ数)に基づく評価基準
- 赤池情報量基準
- ベイズ情報量基準
- マローズの
基準C_p
スタインの等式はモデルの自由度に基づく評価基準に関係する。lassoによって推定されたモデルなら、モデルの複雑さは係数パラメータ数だけではなく正則化パラメータ
にも依存しているので、係数パラメータ数を自由度の推定量とするのは適当ではないらしい。スタインの等式を用いていろんな式をこねくり回すと、lassoにおける自由度は係数パラメータについて非0の要素数を用いればよいという結果が出てくる。Zou et al., 2007またはスパース推定法による統計モデリング[川野秀一・松井秀俊・廣瀬慧(2018)]参照。\lambda
モーメントを求める場合の利用についてまとめる。
実数
と変形して、最初の項にスタインの等式を用いると、
と書け、
となり、モーメントの次数を落とせる。以下この操作を繰り返すと、例えば
と書ける。
3.3.2 スターリングの公式
スターリングの公式は、ガンマ関数を近似するときに役立つ。
という具合に近似可能。
-
[証明]
において\Gamma(k+1) = \int_{0}^\infty x^{k+a-1}e^{-1}dx なる変数変換をするとx= k + \sqrt{k}z より、dx=\sqrt k dz \Gamma(k+a) = k^{k+1-1/2}e^{-k}\int_{-\sqrt k}^\infty\left(1 + \frac{1}{\sqrt k}\right)^{k+a-1}e^{-\sqrt{k}z}dz と書き換えられる。積分の中身をまとめると
\exp\{(k+a-1)\ln(1+z/\sqrt k )-\sqrt k z\} となって、対数のテーラー展開より
(k+a-1)\ln\left(1 + \frac{z}{\sqrt k}\right) - \sqrt k z = -\frac{z^2}{2} + o(1) と近似できる。よって、
\left(1 + \frac{z}{\sqrt k}\right)^{k+a-1}e^{-\sqrt k z} = e^{-z^2/2}\{1 + o(1)\} となって、これを積分すると
\int_{-\sqrt k}^\infty \left(1 + \frac{z}{\sqrt k}\right)^{k+a-1}e^{-\sqrt k z}dz = \int_{-\sqrt k}^\infty e^{-z^2/2}dz + o(1) = \sqrt{2\pi} + o(1) が導かれる。以上から、
\Gamma(k+a) = k^{k+a-1/2}e^{-k}\{\sqrt{2\pi} + o(1)\} と近似できるのでスターリングの公式が得られる。
他のバージョンもある。
-
より荒っぽい評価
k! \approx \left(\frac{k}{e}\right)^k -
対数型
\ln k! \approx k\ln k - k + \frac{1}{2}\ln(2\pi k) \ln k!\approx k \ln k -k -
より精密な評価
k!\approx \sqrt{2\pi k}\left(\frac{k}{e}\right)^k\left(1 + \frac{1}{12k}\right)
[参照]https://manabitimes.jp/math/763
情報理論・統計力学の文脈で使われる。エントロピーの導出で出てきます。詳しくはprml上巻の1.6節「情報理論」を参照。
参考文献
- 久保川達也.(2017).現代数理統計学の基礎 . 共立出版 .
- 川野秀一・松井秀俊・廣瀬慧.(2018). スパース推定法による統計モデリング. 共立出版.
- Bishop, Christopher M. Pattern recognition and machine learning. springer, 2006.
Discussion