🙆‍♀️

3章後半 代表的な確率分布(後半:連続型確率変数の確率分布)

2021/12/21に公開

この記事はfMRILab Advent Calendar 2021の記事です!
弊ラボにおける「現代数理統計学の基礎」(久保川,2017)の勉強会資料として作成したものを公開します。

連続分布

確率変数が連続型である場合を扱います

適宜こちらのコードを使って遊びます。

3.2.1 一様分布

確率変数X[a,b]上の一様分布(uniform distribution)に従うとは、Xの確率密度関数が

f_X(x|a,b) = \begin{cases} &1/(b-a) & (a \leq x \leq bの場合)\\ &0 & (otherwise) \end{cases} \tag{3.11}

で与えられることを言う。

平均や分散については次に示す通り。

E[X] = \frac{a+b}{2}, E[X^2]=\frac{a^2+ab+b^2}{3}
Var(X) = \frac{(b-a)^2}{12}
  • 導出

    \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}

特にa=0,b=1のときは標準一様分布(standard uniform distribution)と呼ぶ。

一様分布の使い道

3.2.2 正規分布 ←←めちゃつかう

正規分布は、

  • 理論的に扱いやすい分布
  • 平均中心で対称な釣り鐘型
  • 中心極限定理により標本平均の極限分布になっている

などのことから数理統計学において最も重要な分布の一つ。

確率変数Xが平均\mu,分散\sigma^2の正規分布に従うとは、Xの確率密度関数が、

f_X(x| \mu,\sigma^2)=\frac{1}{\sqrt{2\pi}\sigma} \exp\left\{-\frac{(x-\mu)^2}{2\sigma^2} \right \},\,\, -\infty < x < \infty \tag{3.12}

で与えられることを言い、この分布を\mathcal{N}(\mu,\sigma^2)であらわす。

標準化(p.18)して(Z=(X-\mu)/\sigmaで変数変換)すると、

\begin{aligned} \phi(z) = f_Z(z) &= \sigma f_X(\sigma z + \mu | \mu,\sigma^2) \\ &= \frac{1}{\sqrt{2\pi}}\exp\{-z^2 / 2\} \end{aligned} \tag{3.13}

と書ける。これは平均0,分散1の標準正規分布(standard normal distribution)という。
これの分布関数は

\Phi(z) = \int_{-\infty}^{z}\phi(t)dt \tag{3.14}

で表され、この値は本の巻末で数字の表として与えられることが多い。

0<\alpha<1に対して\Phi(z_{\alpha})=1-\alphaとなるz_{\alpha}上側100\alpha%点といい、仮説検定や信頼区間を作るときに使われる。 特にz_{0.025}=1.96,z_{0.05}=1.64の値は頻繁に使われる。

確率密度関数の定義から、\int_{-\infty}^{\infty}\phi(z)dz = 1なので

\int_{-\infty}^{\infty}\exp\{-z^2 / 2 \}dz = \sqrt{2\pi} \tag{3.15}

が成り立つ。(ガウス積分)

また、\phi(z)の対称性より

\Phi(0) = 1/2\\\Phi(-z)=1-\Phi(z)

命題3.9

\mathcal{N}(0,1)に従う確率変数Zの平均は0,分散は1であり、積率母関数はM_Z(t)=\exp \{t^2/2\},特性関数は\varphi_Z(t)=\exp\{-t^2/2\}である。

  • [証明]

    積率母関数の定義より

    \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)より求まる。

    X=\sigma Z+\muなので、命題3.9を用いると

    \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

正規分布\mathcal{N}(\mu,\sigma^2)に従う確率変数Xの平均、分散、積率母関数、特性関数はそれぞれ

\begin{aligned} E[X]&=\mu,\\Var(X)&=\sigma^2,\\ M_X(t) &= \exp \{\mu t + \sigma^2 t^2 /2 \},\\\varphi_X(t)&=\exp \{\mu it - \sigma^2 t^2 /2 \} \end{aligned}

で与えられる。

Xの分布関数は

\begin{aligned} F_X(x) &= P_X(X \leq x) = P((X-\mu)/\sigma \leq (x-\mu)/\sigma) \\ &= P(Z \leq (x-\mu)/\sigma) = \Phi((x-\mu)/\sigma) \end{aligned}

と表される。したがって、P(a< X \leq b)なる確率はP(a < X \leq b) = \Phi((b-\mu))/\sigma) - \Phi((a-\mu)/\sigma)と書けるので、標準正規分布表から求められる。

命題3.11

X\mathcal{N}(\mu,\sigma^2)に従う確率変数とする。定数a,bに対してaX+b\mathcal{N}(a\mu + b, a^2 \sigma^2)に従う。

  • [証明]

    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 ガンマ分布とカイ二乗分布

非負の実数直線上の確率分布の中で代表的な分布がガンマ分布。この分布には特別な場合として指数分布やカイ二乗分布という分布を含んでいて、あとで示すようにカイ二乗分布は正規分布と関連する分布である。

確率変数Xの確率密度関数が

f_X(x|\alpha,\beta) = \frac{1}{\Gamma(\alpha)}\frac{1}{\beta}\left(\frac{x}{\beta}\right)^{\alpha-1}e^{-x/\beta},\, x>0 \tag{3.16}

で与えられるとき、Xガンマ分布(gamma distribution)に従うといい、この分布をGa(\alpha,\beta)であらわす。

ここで\alpha形状母数(shape parameter),\beta尺度母数(scale parameter)と呼ばれ、ともに\alpha>0,\beta>0である。\Gamma(\alpha)ガンマ関数で、

\Gamma(\alpha) = \int_{0}^{\infty} y^{\alpha-1}e^{-y}dy

で与えられる。

尺度変換Y=X/\betaを行うと、命題2.20の中の(\mu,\alpha)(0,\beta)に対応するので、(2.8)と(3.16)を見比べると、Yの分布は\beta f_X(\beta y | \alpha,\beta)すなわち、

f_Y(y|\alpha) = \frac{1}{\Gamma(\alpha)}y^{\alpha-1}e^{-y},\,y>0 \tag{3.17}

となり、ガンマ関数の定義からf_Y(y|\alpha)の積分が1になる。

\begin{aligned} \int_{0}^{\infty}f_Y(y|\alpha)dy &= \int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}y^{\alpha-1}e^{-y}dy \\ &= \frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}y^{\alpha-1}e^{-y}dy \\ &= \frac{1}{\Gamma(\alpha)} \Gamma(\alpha)\\ &=1 \end{aligned}

参考:https://mathtrain.jp/gammadist

ガンマ分布は単位時間あたりの回数が1/\betaのイベントが\alpha回起こるまでの時間間隔を表す。

https://takuyaokada.hatenablog.com/entry/20140812/1407837832

[別の表現]

期間\betaごとに1回くらい起こるランダムな事象がn回起こるまでの時間の分布
https://manabitimes.jp/math/1386

命題3.12

(3.17)に従う確率変数Yについて、平均E[Y]=\alpha,分散Var(T) = \alphaとなり、積率母関数はM_Y(t)=(1-t)^{-\alpha},(|t|<1),特性関数は\varphi_Y(t) = (1-it)^{-\alpha}となる。

  • [証明]

    まず積率母関数

    \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}より、
    E[Y] = M'_Y(0)=\alpha, E[Y^2]=M''_Y(0)=\alpha(\alpha+1) (p.21参照)となり、

    \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)\,\,\Gamma(1/2) = \sqrt{\pi}, \Gamma(1) = 1

(2)\,\,\Gamma(\alpha + 1) = \alpha \Gamma(\alpha),\,(\alpha>0),特にnが自然数の時には\Gamma(n)=(n-1)!

  • [証明]

    (1)について

    標準正規分布\mathcal{N}(0,1)の確率密度関数\phi(z)に平方変換(p.25)を行う。

    \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

    また、

    \Gamma(1) = \int_{0}^{\infty} e^{-x} dx =\left[-e^{-x} \right]_{0}^{\infty} = 1(2)

    部分積分を用いる。

    \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のカイ二乗分布' が含まれ、それぞれGa(1,1/\lambda),Ga(n/2,2)に対応する。

正規分布のベイズ推測において、ガンマ分布は正規分布の精度パラメータ(\beta = 1/\sigma^2)の共役事前分布として使われる。詳しくは6章のベイズ法、あるいはprml上巻のp.95-100.

定義3.14

nを自然数とする。確率変数Xの確率密度関数が

f_X(x) = \frac{1}{\Gamma(n/2)}\left(\frac{1}{2}\right)^{n/2}x^{n/2-1}\exp\{-x/2\}, \, \tag{3.21}

で与えられるとき、自由度nのカイ二乗分布(chi-square distribution with n degrees of freedom)といい、\chi_n^2であらわす。

\chi_n^2の特性関数はGa(n/2,2)と(3.18)より

\varphi_{\chi_n^2}(t)= (1-2it)^{-n/2}

で与えられ、また、

E[\chi_n^2]=n,Var(\chi_n^2)=2n,E[(\chi_n^2)^2]=Var(\chi_n^2) + (E[\chi_n^2])^2 = n(n+2)

となる。

  • [証明]

    特性関数とモーメントに関する式から

    \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(この先結構出てくるので覚えておくこと!)

確率変数Z\mathcal{N}(0,1)に従うとき、Z^2\chi_1^2,すなわち自由度1のカイ二乗分布に従う。

  • [証明]

    標準正規分布の確率密度関数\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)f_X(x) = \frac{1}{\Gamma(n/2)}\left(\frac{1}{2}\right)^{n/2}x^{n/2-1}\exp\{-x/2\}において、n=1となっている場合に相当!したがって自由度1のカイ二乗分布に従う。

参考:https://stats.biopapyrus.jp/probability/chi-squared-distribution.html

カイ二乗分布のさらなる性質は次章以降で取り扱う。

3.2.4 指数分布とハザード関数

ある期間に平均して\lambda回起こる現象が次に起こるまでの期間をXとおく。このときXの確率密度関数は指数分布として

E_X(\lambda)=f_X(x|\lambda)=\lambda e^{-\lambda x}, \, x>0 \tag{3.23}

で与えられる。

具体例として、webサイトの訪問者、料金所に到着する車と行った事象間の時間分布をモデル化することができる。

分布関数はF_X(x)=P(X<x)=1-e^{-\lambda x}

指数分布はガンマ分布の特別な場合でGa(1,1/\lambda)と表されるので、(3.18)から

\begin{aligned}M_X(t) &=\lambda/(\lambda -t) (t<\lambda)\\\varphi_X(t) &= \lambda/(\lambda - it) \end{aligned}

となる。
実際にグラフを作ってみるとたしかに指数分布とガンマ分布の関係性がわかる。

また、E[X] = 1/\lambda,Var(X) = 1/\lambda^2となる。

指数分布の性質と無記憶性

指数分布は生存時間の分布として用いられることがある。

P(X>s)は時間sを超えて生存する確率を表現し、P(X>s) = 1-F_X(s)=e^{-\lambda s}となる。s時間生存していて、さらにt時間生存する確率は

\begin{aligned} P(X \geq s+t | X \geq s) &= \frac{P(X \geq s+t, X \geq s)}{P(X \geq s) } = \frac{P(X \geq s+t)}{P(X\geq s)}\\ &=e^{-\lambda(s+t)}/e^{-\lambda s} =e^{-\lambda t} = P(X \geq t) \end{aligned}

となり、それまで生存してきた時間には依存しないことが分かる。

換言すると、指数分布では将来の事象発生までの時間がその過去の経過に依存しない。

幾何分布の時と同様に無記憶性が指数分布についても成立している。これは、指数分布において故障(or死亡)がランダムに起こることを意味し、これに関連した考え方がハザード関数

ハザード関数

Xが非負の連続型確率変数とし、その密度関数をf(x),分布関数をF(x),Xを故障する時間とする。x時間まで動作していて次の時間x+\Deltaまでに故障する条件付確率は

\begin{aligned} P(x < X \leq x + \Delta | X > x)  &= \frac{P(x<X \leq x + \Delta, X > x)}{P(X > x)} \\  &= \frac{P(x<X \leq x + \Delta)}{P(X>x)} = \frac{F(x+\Delta)- F(x)}{1-F(x)} \end{aligned}

となる。このとき、両辺を\Deltaで割って\Deltaを小さくすると、

\lim_{\Delta \to 0}\frac{1}{\Delta} P(x<X\leq x + \Delta| X+x)= \frac{f(x)}{1-F(x)}

となる。これは、xまで動作している条件の下で次の瞬間に故障する確率密度を表している。これを

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

と書いて、ハザード関数(故障関数)という。

指数分布(3.23)については、常に\lambda(x)=\lambdaで時間xには無関係になる。

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

すなわち、x時間動作していて次の瞬間故障する確率密度は一定で\lambda。これは、指数分布において故障がランダムに起こることからも理解できる。

(3.24)の両辺を積分すると

 \int_0^{x} \lambda(t)dt = \int_0^{x} f(t)/\{1-F(t) \}dt = [-\log(1-F(t))]_0^{x}

となるので、

\begin{align} F(x) &= 1- \exp \left\{-\int_0^x \lambda(t)dt \right\}, \\ f(x) &= \lambda(x)\exp\left\{-\int_0^{x} \lambda(t)dt \right\} \tag{3.25} \end{align}

が得られる。非負の連続型確率変数の分布はハザード関数によって特徴づけられることが分かる。

例えば\lambda(x)=\lambda(\lambda:定数)というふうに置くと、式3.25から指数分布が生ずることが分かる。

指数分布(あるいはポアソン分布)を用いたシミュレーションにおけるキーポイントは\lambda(次の瞬間に壊れる確率)が考慮している期間において一定であること。とはいうものの、ハザード関数が時間の経過に関して一定というのは流石にきつい仮定であるように思われる。

例えば交通量やネットワークトラフィックなどは同じ1日という区間でも曜日によって異なるであろうことは容易に想像が出来る。それでもポアソン分布や指数分布を用いたシミュレーションが成立するのは、通常、時間も空間も一様な区間に分割でき、その区間の中だけでなら妥当な仮定となるかららしい。^1

では壊れる可能性が経時的に変化していくような場合を考えてみる。例えば、時間の経過とともに故障しやすくなる場合や故障しなくなる場合には\lambda(x)=abx^{b-1},a>0,b>0なるハザード関数が考えられる。これを(3.25)に代入すると、

\int_0^{x} abt^{b-1}dt = ax^b

より、得られる確率密度関数は

f(x|a,b) = abx^{b-1}\exp \left\{-ax^b \right\}, \, x>0 \tag{3.26}

となる。これは、ワイブル分布(Weibull distribution)といい、生存解析の分野で基本となる分布。

機械などが壊れる、劣化などの事象がお起こる確率を近似する際に用いられる。
https://recruit.gmo.jp/engineer/jisedai/blog/survival_analysis_covid19/

3.2.5 ベータ分布

ベータ分布は区間(0,1)上の確率分布であり、ガンマ分布と関連する分布である。確率密度関数は

Beta(a,b)=f_X(x|a,b) = \frac{1}{B(a,b)}x^{a-1}(1-x)^{b-1},\,\,0<x<1 \tag{3.27}

で表される。ここで、a>0, b>0であり、B(a,b)ベータ関数

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

である。ベータ関数はガンマ関数を用いて

B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b) \tag{3.29}

なる関係式が成り立つ。(次章の(4.17)で示されるらしい)

これを使って

E[X] = \frac{a}{a+b},\, Var(X) = \frac{ab}{(a+b)^2 (a+b+1)} \tag{3.30}

ベータ分布は次章の例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}

      となる。(参考:https://mathtrain.jp/betadistribution)

    • 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

X\mathcal{N}(\mu,\sigma^2)に従い、g(\cdot)が微分可能でE[|g'(x)|]<\inftyのとき、

E[(X-\mu)g(X)]= \sigma^2 E[g'(X)] \tag{3.31}
  • [証明]

    (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)に関する期待値に\sigma^2をかkたものなので

    E[(X-\mu)g(X)]=\sigma^2E[g'(X)]

スタインの等式(3.31)は縮小推定の性質(この場合lassoっぽい?理解が浅いので曖昧です)を議論するときに重要な道具となる。

  • 詳細

    正則化パラメータの選択方法として本書では3つの手法が紹介されている

    • 自由度調整済み決定係数
    • 交差検証(よく使うやつ)
    • モデルの自由度(モデルのパラメータ数)に基づく評価基準
      • 赤池情報量基準
      • ベイズ情報量基準
      • マローズのC_p基準

    スタインの等式はモデルの自由度に基づく評価基準に関係する。lassoによって推定されたモデルなら、モデルの複雑さは係数パラメータ数だけではなく正則化パラメータ\lambdaにも依存しているので、係数パラメータ数を自由度の推定量とするのは適当ではないらしい。スタインの等式を用いていろんな式をこねくり回すと、lassoにおける自由度は係数パラメータについて非0の要素数を用いればよいという結果が出てくる。Zou et al., 2007またはスパース推定法による統計モデリング[川野秀一・松井秀俊・廣瀬慧(2018)]参照。

モーメントを求める場合の利用についてまとめる。

実数mに対してE[X^m]を求めたいとする。(確率変数Xの原点まわりのm次モーメント)

E[X^m] = E[(X-\mu)X^{m-1}] + \mu E[X^{m-1}] 

と変形して、最初の項にスタインの等式を用いると、

E[(X-\mu)X^{m-1}] = (m-1)\sigma^2 E[X^{m-2}] 

と書け、

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

となり、モーメントの次数を落とせる。以下この操作を繰り返すと、例えばm=3,m=4のときには

\begin{aligned} E[X^3] &= 3\mu \sigma^2 + \mu^3 \\ E[X^4] &= 3\sigma^4 + 6\mu^2 \sigma^2 + \mu^4 \end{aligned}

と書ける。

3.3.2 スターリングの公式

スターリングの公式は、ガンマ関数を近似するときに役立つ。kが大きいとき、\Gamma(k+a) \approx \sqrt{2\pi} k^{k+a-1/2}e^{-k}で近似できる。a=1を代入すると、

k! \approx \sqrt{2\pi}k^{k+1/2}e^{-k}

という具合に近似可能。

  • [証明]

    \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