Open4

統計学まわりの小ネタメモ

畳屋民也畳屋民也

相関係数と回帰係数の関係って、どうなっていたっけ?

設定

変数 X_1, X_2 について考える。
それぞれ平均・分散が(\mu_1, \sigma_1^2), (\mu_2, \sigma_2^2) の正規分布に従い、かつ相関係数は \rho である。

\begin{pmatrix} X_1\\ X_2 \end{pmatrix} \sim N \left( \begin{pmatrix} \mu_1\\ \mu_2 \end{pmatrix}, \begin{pmatrix} \sigma_1^2 & \rho \sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma_2^2 \end{pmatrix} \right)

問い

上記の状況で、単回帰モデル

X_1 = c + \alpha X_2 + \varepsilon, \quad \varepsilon \sim N(0, \sigma_{\varepsilon}^2)

を作成すると、回帰係数 \alpha と誤差項分散 \sigma_{\varepsilon}^2 は、相関係数 \rho を用いてどのように表されるか?

結果

\begin{aligned} \alpha &= \rho \frac{\sigma_1}{\sigma_2}\\ \sigma_{\varepsilon}^2 &= (1-\rho^2) \sigma_1^2 \end{aligned}

導出

X_2 による X_1 への回帰式と条件付き期待値・条件付き分散には、以下のような関係がある:

\begin{aligned} E[X_1 \vert X_2] = c + \alpha X_2\\ V[X_1 \vert X_2] = \sigma_{\varepsilon}^2. \end{aligned}

ここで、多変量正規分布の条件付き確率を思い出す。

ベクトル \boldsymbol{X}^\top = (\boldsymbol{X}_a^\top, \boldsymbol{X}_b^\top) が以下のような正規分布に従っていたとする:

\begin{pmatrix} \boldsymbol{X}_a\\ \boldsymbol{X}_b \end{pmatrix} \sim N \left( \begin{pmatrix} \boldsymbol{\mu}_a\\ \boldsymbol{\mu}_b \end{pmatrix}, \begin{pmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix} \right).

このとき、以下が成り立つ:

\begin{aligned} E[\boldsymbol{X}_a \vert \boldsymbol{X}_b] &= \boldsymbol{\mu}_a + \Sigma_{ab} \Sigma_{bb}^{-1} (\boldsymbol{X}_b - \boldsymbol{\mu}_b)\\ V[\boldsymbol{X}_a \vert \boldsymbol{X}_b] &= \Sigma_{aa} - \Sigma_{ab} \Sigma_{bb}^{-1} \Sigma_{ba} \end{aligned}.

これを用いると、

\begin{aligned} E[X_1 \vert X_2] &= \mu_1 + \rho \sigma_1 \sigma_2 \frac{1}{\sigma_2^2}(X_2 - \mu_2)\\ &= \mu_1 + \rho \frac{\sigma_1}{\sigma_2}(X_2 - \mu_2)\\ E[X_1 \vert X_2] &= \sigma_1^2 - \rho \sigma_1 \sigma_2 \frac{1}{\sigma_2^2} \rho \sigma_1 \sigma_2\\ &= (1 - \rho^2) \sigma^2 \end{aligned}

以上から、

\begin{aligned} \alpha &= \rho \frac{\sigma_1}{\sigma_2}\\ \sigma_{\varepsilon}^2 &= (1-\rho^2) \sigma_1^2 \end{aligned}

が成り立つ。
(なお、c = \mu_1 - \alpha \mu_2

考察

相関係数\rhoと回帰係数\alphaの関係

相関係数\rhoが0の時、回帰係数\alphaも0になる。

相関係数\rho\pm 1のとき、回帰係数\alphaの値はX_1, X_2の標準偏差の比で決まる。

相関係数\rhoと誤差\varepsilonの関係

相関係数\rho\pm 1のとき、誤差項分散は\sigma_\varepsilon^2=0となる。
誤差項はX_1X_2に対して射影した際の残成分を表すので、\rho = \pm 1であればX_1X_2への射影成分のみで表現できていることがわかる。

逆に\rho = 0 であれば誤差項分散とX_1の分散が等しくなる \sigma_1^2 = \sigma_\varepsilon^2 が、
このときは前述のように\alpha=0となることからもX_1X_2に射影して表現することができていないと解釈できる。

逆に、単回帰の係数や誤差項分散の大きさから相関係数について何がわかるか?

\sigma_1>0, \sigma_2>0 を仮定する。

この場合、

  • \alpha=0 ... 無相間 \rho=0
  • \sigma_\varepsilon^2=0 ... \rho = \pm 1

が言える。

一方で、回帰係数の大小から直接相関係数の大小を論じることはできない。
仮に相関係数が 0<\rho <<1と小さくても、\sigma_1 >> \sigma_2 であれば回帰係数 \alphaが大きな値をとりうるからである(逆も同様)。

畳屋民也畳屋民也

尖度と歪度

定義

ある分布に従う確率変数 X の期待値を E[X]、分散を V[X] = E[(X-\mu)^2] と置く。
このとき、この分布の歪度(skewness) \sqrt{}\beta_1 と 尖度(kurtosis) \beta_2 は以下のように定義される:

\sqrt{}\beta_1 = \frac{E[(X - E[X])^3]}{(V[x])^{\frac{3}{2}}},\tag{1}
\beta_2 = \frac{E[(X - E[X])^4]}{(V[x])^2}. \tag{2}

なお、後述のように正規分布では上記の定義における尖度の値が \beta_2=3 になるので、正規分布からのずれを評価する際などでは

\beta_2' = \beta_2 - 3 = \frac{E[(X - E[X])^4]}{(V[x])^2} - 3

を尖度として用いることもある。

意味

歪度(skewness)は、分布が期待値を基準にどれだけ左右非対称かを表す。
\sqrt{}\beta_1>0 では分布の山が左に寄り、分布の裾が右方向に長く伸びる。
逆に \sqrt{}\beta_1<0 では分布の山が右に寄り、分布の裾が左方向に長く伸びる。

尖度(kurtosis)は、分布の尖り具合と裾の広さを表す。
\beta_1 が大きいほど、分布は期待値付近では尖ると同時に裾が広くなっていく。
\beta_1 が0に近づくほど、期待値付近でなだらかな分布になる。

正規分布の場合

確率変数 X が正規分布 P(X=x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\} に従うとき、X の分布の歪度・尖度はそれぞれ以下である:

\begin{aligned} \sqrt{}\beta_1 &= 0,\\ \beta_2 &= 3. \end{aligned}
導出

まず、E[X]=\mu, \, V[X]=\sigma^2 である。

正規分布 N(\mu, \sigma^2)X=\mu について対称なので、E[(X-\mu)^3]= 0 である。したがって、式(1)より歪度 \sqrt{}\beta_1=0

次に、下記より E[(X-\mu)^4] = 3\sigma^4 なので、V[X]=\sigma^2 も用いて式(2)より尖度 \beta_2=3

\begin{aligned} E[(X-\mu)^4] &= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^\infty (x-\mu)^4 \exp \left\{- \frac{(x-\mu)^2}{2\sigma^2} \right\} dx\\ &=\frac{1}{\sqrt{2\pi\sigma^2}} 2 \int_0^\infty z^2 \exp \left(- \frac{z}{2\sigma^2} \right) \frac{1}{2\sqrt{z}}dz \qquad ... \,z = (x-\mu)^2\\ &= \frac{1}{\sqrt{2\pi\sigma^2}} (2\sigma^2)^{5/2} \int_0^\infty z^{5/2-1} \exp (-z) dz\\ &= \frac{1}{\sqrt{\pi}} 4\sigma^4 \Gamma\left(\frac{5}{2}\right)\\ &= \frac{1}{\sqrt{\pi}} 4\sigma^4 \frac{3}{4} \sqrt{\pi} \qquad \because \Gamma(5/2) = \frac{3}{2}\cdot \frac{1}{2} \cdot \Gamma(1/2) = \frac{3}{4} \sqrt{\pi}\\ &= 3\sigma^2 \end{aligned}

Bernoulli 分布の場合

確率変数 X \in \{0,1\} が Bernoulli 分布 P(X=x) = p^x (1-p)^{1-x} に従うとき、X の分布の歪度と尖度は以下のようになる:

\begin{aligned} \sqrt{}\beta_1 &= \frac{1 - 2p}{\sqrt{p(1-p)}},\\ \beta_2 &= \frac{1}{p(1-p)} - 3. \end{aligned}
導出

まず、E[X]=p, \, V[X]=p(1-p) である。
そのうえで、

\begin{aligned} E[(X - E[X])^3] &= (1 - p)^3 p + (0 - p)^3 (1-p)\\ &= p(1-p)(1-2p) \end{aligned}

であるので、式(1)より歪度は

\sqrt{}\beta_1 = \frac{p(1-p)(1-2p)}{\{p(1-p)\}^{3/2}} = \frac{1-2p}{\sqrt{p(1-p)}}.

次に、

\begin{aligned} E[(X-E[X])^4] &= (1-p)^4 p + (0-p)^4 (1-p)\\ &= p(1-p)(1 - 3p + 3p^2)\\ &= p(1-p) - 3 p^2 (1-p)^2 \end{aligned}

であるので、式(2)より尖度は

\beta_2 = \frac{p(1-p) - 3 p^2 (1-p)^2}{p^2(1-p)^2} = \frac{1}{p(1-p)} - 3.

特に、p=1/2 の時、歪度 \sqrt{}\beta_1 = 0、 尖度 \beta_2 = 1 となる。

参考文献

  • 中川重和「正規性の検定 (統計学One Point)」(共立出版、2019)

https://amzn.asia/d/fQlxRcv

畳屋民也畳屋民也

Gamma 分布

以下のような確率密度関数を持つ分布を、Gamma 分布と呼ぶ:

{\rm Ga}(x; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} \exp(-\beta x).

なお、\Gamma(\alpha) は、以下のように定義されたガンマ関数である:

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

上記の表記の他に、k=\alpha, \theta=1/\beta として、

{\rm Ga}(x; k, \theta) = \frac{1}{\Gamma(k) \theta^k} x^{k - 1} \exp(- x / \theta)

と表すこともある。
NumPy, SciPy では後者の表記が用いられており、 scale という引数名で \theta = 1/\beta を指定する。

統計量

期待値

E[X] = \frac{\alpha}{\beta}
導出
\begin{aligned} E[X] &= \frac{\beta^\alpha}{\Gamma(\alpha)}\int_0^\infty x x^{\alpha - 1} \exp( - \beta x) dx\\ &= \frac{\beta^\alpha}{\Gamma(\alpha)}\int_0^\infty x^{(\alpha+1) - 1} \exp( - \beta x) dx\\ &= \frac{\beta^\alpha}{\Gamma(\alpha)}\frac{\Gamma(\alpha+1)}{\beta^{\alpha+1}}\\ &=\frac{1}{\beta}\frac{\alpha \Gamma(\alpha)}{\Gamma(\alpha)} \quad (\because \Gamma(\alpha+1) = \alpha \Gamma(\alpha))\\ &= \frac{\alpha}{\beta} \end{aligned}

分散

V[X] = \frac{\alpha}{\beta^2}
導出
\begin{aligned} V[X] &= E[(X - E[X])^2]\\ &= E[X^2] - (E[X])^2 \end{aligned}

ここで、

\begin{aligned} E[X^2] &= \frac{\beta^\alpha}{\Gamma(\alpha)}\int_0^\infty x^2 x^{\alpha - 1} \exp( - \beta x) dx\\ &= \frac{\beta^\alpha}{\Gamma(\alpha)}\frac{\Gamma(\alpha+2)}{\beta^{\alpha+2}}\\ &=\frac{1}{\beta^2}\frac{\alpha (\alpha + 1) \Gamma(\alpha)}{\Gamma(\alpha)}\\ &= \frac{\alpha(\alpha+1)}{\beta^2} \end{aligned}

より、

\begin{aligned} E[X^2] - (E[X])^2 &= \frac{\alpha(\alpha+1)}{\beta^2} - \frac{\alpha^2}{\beta^2}\\ &=\frac{\alpha}{\beta^2} \end{aligned}

よって、V[X]= \alpha / \beta^2 となる。

歪度(skewness)

{\rm Skew}[X] = \frac{1}{\sqrt{\alpha}}
導出
\begin{aligned} {\rm Skew}[X] = \frac{E[(X - E[X])^3]}{V[X]^{3/2}} \end{aligned}

ここで、

\begin{aligned} E[(X - E[X])^3] &= E[X^3] - 3 E[X^2]E[X] + 3 E[X](E [X])^2 - (E[X])^3\\ &=E[X^3] -3 E[X] V[X] - (E[X])^3\\ &= \frac{\alpha (\alpha+1)(\alpha+2)}{\beta^3} - 3 \frac{\alpha}{\beta} \frac{\alpha}{\beta^2} - \frac{\alpha^3}{\beta^3}\\ &=\frac{2\alpha}{\beta^3} \end{aligned}

である。途中、 E[X^3] = \alpha(\alpha+1)(\alpha+2) / \beta^3 を用いた(導出省略)。

従って、

\begin{aligned} \frac{E[(X - E[X])^3]}{V[X]^{3/2}} &= \frac{2\alpha}{\beta^3} \frac{\beta^3}{\alpha^{3/2}}\\ &= \frac{2}{\sqrt{\alpha}} \end{aligned}

より {\rm Skew}[X] = 2/\sqrt{\alpha} が示された。

尖度(kurtosis)

{\rm Kurt}[X] = \frac{6}{\alpha} + 3
導出
{\rm Kurt}[X] = \frac{E[(X - E[X])^4]}{V[X]^2}

ここで、

E[(X - E[X])^4] = \frac{3\alpha^2 + 6\alpha}{\beta^4}

より(導出省略)、

\begin{aligned} \frac{E[(X - E[X])^4]}{V[X]^2} &= \frac{3\alpha^2 + 6\alpha}{\beta^4}\frac{\beta^4}{\alpha^2}\\ &= \frac{6}{\alpha} + 3 \end{aligned}

となり、{\rm Kurt}[X]=6/\alpha + 3 が示された。

モーメント母関数

E[e^{tX}] = (1 - t / \beta)^{-\alpha} \quad {\rm for }\,t< 1/\beta
導出
\begin{aligned} E[e^{tX}] &= \frac{\beta^\alpha}{\Gamma(\alpha)}\int_0^\infty e^{tx} x^{\alpha - 1} \exp( - \beta x) dx\\ &=\frac{\beta^\alpha}{\Gamma(\alpha)}\int_0^\infty x^{\alpha - 1} \exp\{ - (\beta-t) x\} dx\\ &= \frac{\beta^\alpha}{\Gamma(\alpha)}\frac{\Gamma(\alpha)}{(\beta-t)^{\alpha}}\\ &=(1 - t/\beta)^{-\alpha} \end{aligned}

参考

https://en.wikipedia.org/wiki/Gamma_distribution

畳屋民也畳屋民也

「総額の q % を上位 p % のユーザーが占めている」をどう表すか?

よくアプリの売上や資産額などで、「総額の q %を上位 p % のユーザーが占めている」などといった記述を目にする。

これを数式を用いてどのように表すことができるかを考えてみた。
さらに、Gamma 分布では上記を簡単に計算する方法があることに気がついたのでそれについても触れる。

以下、面倒なので p, q は百分率ではなく小数で表す。

数式で表してみる

売上や資産額を x と表し、その確率密度分布を P(x) とおく。
また、x の最小値・最大値をそれぞれ x_{\rm min}, x_{\rm max} とする。

このとき、

p = \int_z^{x_{\rm max}} P(x)dx \tag{1}

と表せる。ただし、 z は以下を満たすものとする:

\frac{\int_z^{x_{\rm max}} x P(x) dx}{\int_{x_{\rm min}}^{x_{\rm max}} x P(x) dx} = q \tag{2}

Gamma 分布の場合

\begin{aligned} P(x) &= {\rm Ga}(x; \alpha, \beta)\\ &= \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} \exp(-\beta x) \end{aligned}

とする。

このとき、式(2) の左辺は以下のように表される:

\begin{aligned} \frac{\int_z^\infty x {\rm Ga}(x; \alpha, \beta) dx}{\int_{-\infty}^\infty x {\rm Ga}(x; \alpha, \beta) dx} &= \frac{\int_z^\infty x^{(\alpha+1) - 1} \exp(-\beta x)dx}{\int_{-\infty}^\infty x^{(\alpha+1) - 1} \exp(-\beta x) dx}\\ &= \frac{\beta^{\alpha+1}}{\Gamma(\alpha+1)} \int_z^\infty x^{(\alpha+1) - 1} \exp(-\beta x)dx\\ &= \int_z^\infty {\rm Ga}(x; \alpha+1, \beta)dx \end{aligned}

これは、 {\rm Ga}(x; \alpha+1, \beta) の Survival 関数に他ならない。

従って、確率密度分布が {\rm Ga}(x; \alpha+1, \beta) と表される Gamma 分布の累積分布関数を Q_G(x; \alpha+1, \beta) とおくと、z はその逆関数 Q_G^{-1}(y; \alpha+1, \beta)q を用いて以下のように表すことができる:

z = Q_G^{-1}(1-q; \alpha+1, \beta)

これらを用いると、式(1) から p は以下のように求まる:

\begin{aligned} p &= \int_z^\infty {\rm Ga}(x; \alpha, \beta)dx\\ &= 1 - Q_G(z; \alpha, \beta) \end{aligned}

計算例

q=0.99 としたときの p\alpha=0.001, \beta=10^{-5} の場合について Python を用いて計算すると以下のようになる:[1]

from scipy.stats import gamma

alpha = 0.001
beta = 0.00001

q=0.99

z = gamma.ppf(1-q, alpha+1, scale=1/beta)
# もしくは、z = gamma.isf(q, alpha+1, scale=1/beta)
p = gamma.sf(z, alpha, scale=1/beta)

print(p)
# > 0.004020670587092209

(参考)
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html

脚注
  1. 実は、 \beta > 0 には何を入れても構わない。 ↩︎