axjack's blog

### axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz ###

二項分布の正規近似から母比率の検定統計量まで

統計 勉強

二項分布の正規近似は闇が深いっぽいので、深追いせず流れだけを書いておくメモです。

X \sim Bin(n, p)
↓
↓ 正規近似 as n→∞
↓ 実はここで色々厄介な近似をしているらしい。
↓ たとえばスターリングの公式など。
↓
↓ X \sim  N(np, np(1-p) )
↓
↓ Xを標準化(平均値引いて標準偏差で割る)
↓
↓ \displaystyle \frac{ X-np } { \sqrt{ np(1-p) }}  \sim N(0, 1)
↓
↓ 確率変数の分母分子をnで割る。確率変数は実質変わらない。
↓
\displaystyle   \frac{\bar{X} - p}{  \sqrt{   \frac{p(1-p)}{n} } }  \sim N(0,1)

切断正規分布

統計 統計学実践ワークブック

統計学実践ワークブックの問6.1〔4〕は、いわゆる切断正規分布の問題である。このキーワードでググると良い。 なお、期待値を求めるために確率変数 Z|Z>0 のモーメント母関数を計算しようとすると怪我をする？ので、素直に \displaystyle \frac{ f(\rm z) }{ P(Z>0) } から定義通り期待値を求めましょう。

2022-02-19 追記

ワークブックに従い、素直に確率密度関数f(z) = (1/2)×ϕ(z) where z > 0 , 0 where z ≦ 0 と考えた方が楽。

中心極限定理をラフに証明する

統計

厳密さを捨てて大略理解できれば良いぐらいの証明です。

準備

正規分布のモーメント母関数

確率変数X が N(\mu, \sigma^2) に従う時、モーメント母関数 M_{X}[t; \mu , {\sigma}^2 ] は、
\displaystyle M_{X}[t; \mu , {\sigma}^2 ] = \rm exp( \mu t + \frac{1}{2} {\sigma}^2 {t}^2 )
とくに、 \mu = 0, \sigma^2 = 1 のとき、 \displaystyle M_{X}[ t; \mu=0 , {\sigma}^2=1 ] = \rm exp( \frac{1}{2} {t}^2 )

指数関数の底eの極限

\displaystyle \lim_{n \rightarrow \infty } (1+\frac{1}{n} )^n = \rm e

\displaystyle \lim_{n \rightarrow \infty } (1+\frac{a}{n} )^n = \lim_{n \rightarrow \infty } ( (1+\frac{a}{n} )^{\frac{n}{a}} )^a =( \lim_{n \rightarrow \infty } (1+\frac{a}{n} )^{\frac{n}{a}} )^a = \rm e^a

テイラー展開

\displaystyle f(x) = \sum_{k=0}^{\infty} \frac{ x^k }{k!} f^{(k)}(0) = 1 + xf'(0) + \frac{1}{2}x^2 f''(0) + \cdots \approx 1 + xf'(0) + \frac{1}{2}x^2 f''(0)
とくに、 \displaystyle \rm e^x = 1 + x + \frac{1}{2}x^2 + \cdots \approx 1 + x + \frac{1}{2}x^2

逆数

\displaystyle \frac{1}{ \frac{1}{\sqrt{n}} } = \Bigl( \frac{1}{\sqrt{n}} \Bigr) ^{-1} = \sqrt{n}

中心極限定理

示したいこと

確率変数 X_1, X_2, \cdots X_n が平均\mu、分散\sigma^2の分布にi.i.d.で従っている。
このとき、標本平均 \displaystyle \bar{X} を標準化した統計量Tは標準正規分布N(0,1)に分布収束する。

標準化した確率変数の確認

以下を示す

確率変数 X_i はi.i.d. で平均\mu、分散\sigma^2で従っている。
X_iを標準化した確率変数を Z_i とするとき、
\displaystyle \bar{X} を標準化した確率変数は \displaystyle \sqrt{n} \bar{Z} となる。

導出

X_i を標準化した変数を Z_i とすると、 \displaystyle Z_i = \frac{X_i - \mu}{\sigma} 、 E[ Z_i ] = 0, V[ Z_i] = 1 である。

つぎに、 \bar{X}について、\displaystyle E[ \bar{X} ] = \mu , V[ \bar{X} ] = \frac{{\sigma}^2}{n} であるから、

標本平均 \displaystyle \bar{X} を標準化した統計量\rm T は \displaystyle \frac{ \bar{X} - \rm E[\bar{X}] }{ \sqrt{ \rm V[ \bar{X} ]}} = \frac{\bar{X} - \mu }{ \frac{\sigma}{\sqrt{n}}} = \sqrt{n} \Bigl( \frac{\bar{X} - \mu }{\sigma} \Bigr)
となる。

一方、 \displaystyle \bar{Z} = \frac{1}{n} \sum Z_i = \frac{1}{n} \sum (\frac{X_i - \mu}{\sigma} ) = \frac{1}{n \sigma} (n \bar{X}-n\mu) = \frac{\bar{X} - \mu}{\sigma}

となるので、 \displaystyle \sqrt{n} \bar{Z} = \sqrt{n} \frac{\bar{X} - \mu}{\sigma} = T となる。

モーメント母関数

統計量Tのモーメント母関数、すなわち \sqrt{n}\bar{Z}のモーメント母関数を求める。 Z_1, \cdots, Z_nが互いに独立であることに留意すると、

\displaystyle M_{ \sqrt{n}\bar{Z} }(t) = \rm E[ \rm exp( t \sqrt{n}\bar{Z} ) ] = \rm E[ exp( t \sqrt{n} \frac{1}{n} \sum Z_i ) ] = \rm E[ exp( t \frac{1}{\sqrt{n} } \sum Z_i ) ] = \rm \Bigl( E[ exp( t \frac{1}{\sqrt{n} } Z_1 ) ] \Bigr)^n

とかける。ここで、

\displaystyle \rm E[ exp( t \frac{1}{\sqrt{n} } Z_1 ) ] \\ = \displaystyle \rm E[ 1 + \frac{t}{\sqrt{n} } Z_1 + \frac{t^2}{2n} {Z_1}^2 + \cdots ] \\ \displaystyle \approx \rm E[ 1 + \frac{t}{\sqrt{n} } Z_1 + \frac{t^2}{2n} {Z_1}^2 ] \\ = \displaystyle 1 + \frac{t}{\sqrt{n} } \rm E[ Z_1 ] + \frac{t^2}{2n} \rm E [ {Z_1}^2 ] \\ = \displaystyle 1 + \frac{t}{\sqrt{n} } \rm 0 + \frac{t^2}{2n} \rm 1 \\ = \displaystyle 1 + \frac{t^2}{2n}

と近似することができるから、

\displaystyle M_{ \sqrt{n}\bar{Z} }(t) = \rm \Bigl( E[ exp( t \frac{1}{\sqrt{n} } Z_1 ) ] \Bigr)^n = \rm \Bigl( 1 + \frac{t^2}{2n} \Bigr)^n = \rm \Bigl( 1 + \frac{ \frac{1}{2}t^2 }{n} \Bigr)^n = \rm exp( \frac{1}{2}t^2 ) \ ( as \ n \rightarrow \infty )

となって、これは標準正規分布のモーメント母関数に一致する。
よって、標本平均 \displaystyle \bar{X} を標準化した統計量Tは標準正規分布N(0,1)に分布収束する。

クーポンコレクターまたはコンプリートガチャ問題

統計 統計学実践ワークブック

統計学実践ワークブック問5.5より。

問題の概要

4種類のカードを等確率無作為復元抽出で引く。

  • [1] 4種類のカードを全て揃えるまでの、回数の期待値
  • [2] あらたに5種類目のカードが追加されたとする。
    • x: はじめの4種類を集めてから、追加の5種類目を揃えるまでの、回数の期待値
    • y: はじめから5種類を全て揃えるまでの、回数の期待値
    • このとき、xとyの差を求めよ

[1]

幾何分布の復習

成功確率pのベルヌーイ試行において、はじめて成功するまでに起こる「失敗の回数」をXとすると、Xは幾何分布に従う。ここでX = kとなる確率は、\displaystyle P(X=k) = (1-p)^{k}p、k = 0,1,...となる。Xの期待値をもとめるため、Xの確率母関数を計算すると、
\displaystyle \rm E[s^X] = \sum_{k=0}^{\infty} s^{k}(1-p)^{k}p = p\sum_{k=0}^{\infty} (s(1-p))^{k} = \frac{p}{1 - (s(1-p))}

である。

期待値E[X]は確率母関数をsで微分してs=1を代入したものとなる*1ので、
\displaystyle \frac{d}{ds}E[s^X] = \frac{d}{ds} \frac{p}{1 - (s(1-p))} \\= \frac{d}{ds} p(1-s(1-p))^{-1} = p(-1)(1-s(1-p))^{-2}(-(1-p)) = p(1-p)(1-s(1-p))^{-2}
より、s = 1を代入し、
\displaystyle p(1-p)(1-s(1-p))^{-2} |_{s=1} = p(1-p)(1-(1-p))^{-2} = p(1-p)p^{-2} = \frac{1-p}{p}

となる。

ここで、Xは『はじめて成功するまでに起こる「失敗の回数」』であるから、W = X+1とすれば「はじめて成功するまでの回数」の確率変数に変換することができる。Wの期待値は \displaystyle E[W] = E[X+1] = E[X] + 1 = \frac{1-p}{p} + 1 = \frac{1-p}{p} + \frac{p}{p} = \frac{1}{p}

となる。つまり、初めて成功するまでの回数を確率変数とする幾何分布の期待値は、成功確率の逆数となる。

4種類のカードを全て揃えるまでの、回数の期待値

4種類のカードを全て揃えるということは、

確率変数 成功確率 失敗確率 幾何分布に従う時の意味
X_1 \frac{4}{4} 1 - \frac{4}{4} 1種類目を引き当てる
X_2 \frac{3}{4} 1 - \frac{3}{4} 1種類揃い済み、2種類目を引き当てる
X_3 \frac{2}{4} 1 - \frac{2}{4} 2種類揃い済み、3種類目を引き当てる
X_4 \frac{1}{4} 1 - \frac{1}{4} 3種類揃い済み、4種類目を引き当てる

という確率変数X_1 , ... , X_4 についての和の期待値E[X_1 + X_2 + X_3 + X_4]を求めることに等しい。
確率変数ごとに確率が異なる（揃うたびに成功確率は小さくなる）ことに注意しつつ、それぞれ幾何分布に従っていることを踏まえると、求める期待値は、

\displaystyle E[ X_1 + X_2 + X_3 + X_4 ] = E[ X_1] + E[ X_2 ] + E[ X_3 ] + E[ X_4 ]
\displaystyle = (\frac{4}{4})^{-1} + (\frac{3}{4})^{-1} + (\frac{2}{4})^{-1} + (\frac{1}{4})^{-1} = ... = \frac{25}{3}
となる。

[2]

  • x: はじめの4種類を集めてから、追加の5種類目を揃えるまでの、回数の期待値
  • y: はじめから5種類を全て揃えるまでの、回数の期待値

について、yはX_1 \sim Geo(p = \frac{5}{5}) 、X_2 \sim Geo(p = \frac{4}{5} ) 、…、X_5 \sim Geo(p = \frac{1}{5} ) について和の期待値を求めれば良いので、\displaystyle y = (\frac{5}{5})^{-1} + ... + (\frac{1}{5})^{-1} = \frac{127}{12} となる。*2

一方、xについては[1]の期待値に X_5 \sim Geo(p = \frac{1}{5} ) の期待値を加えたものとなるので、
x = \frac{25}{3} + E[X_5] = \frac{25}{3} + ( \frac{1}{5})^{-1} = \frac{40}{3} となる。

以上より、 \displaystyle x - y = \frac{40}{3} - \frac{127}{12} = \frac{23}{12} となる。

*1: \frac{d}{ds}E[s^X] = E[Xs^{X-1}]

*2:Geo(p)は成功確率pを持つ幾何分布の意味。

指数分布の和の分布

統計 統計学実践ワークブック

統計学実践ワークブック問4.2より。指数分布の和の分布を求めた時の教訓・感想です。

オチ・教訓・流れ

  • 指数分布は再生性を持っていない
    • 指数分布のモーメント母関数を求めて掛け算しても、元の分布がよく分からない形に。。
    • なので畳み込み積分で和の分布を求める必要がある
  • ところで、モーメント母関数から導いた和の分布のモーメント母関数と、畳み込み積分で求めた和の分布のモーメント母関数は一致するのだろうか？
    • 計算して確かめよう

やること

  • 1) 畳み込み積分にて和の分布を求める
  • 2) モーメント母関数を使って和の分布のモーメント母関数を求める
    • 再生性がないことを確かめる。
  • 3) 1)で求めた和の分布のモーメント母関数を求め、これが2)と一致することを確認する

問題4.2の概要

1)畳み込み積分にて和の分布を求める

Z=X+Y; 畳み込み積分の公式？より、Zの確率密度関数: h(z) は、
h(z) = \int_{0}^{z}f(t)f(z-t)dt = \int_{0}^{z}\lambda\exp(- \lambda t) \lambda\exp(- \lambda (z-t) )dt
= \lambda^2 \int_{0}^{z} \exp(-\lambda t ) \exp( - \lambda (z-t) )dt = \lambda^2 \int_{0}^{z} \exp(- \lambda z )dt = \lambda^2 \exp(-\lambda z) \int_{0}^{z} 1dt
\Large = \rm \lambda^2 z e^{-\lambda z}

2) モーメント母関数を使って和のモーメント母関数を求める

Xのモーメント母関数M_{X}(t)を求め、XとYは独立であることからM_{X+Y}(t) = M_{X}M_{Y}を計算することにより和の分布のモーメント母関数を求める。
なお積分範囲(0,∞)は記載省略。

Xのモーメント母関数は、
M_{X}(t) = E[exp(tX)] = \int \lambda exp(- \lambda x ) exp(tx)dx = \lambda \int exp(-(t-\lambda )x) dx より、
\large = \lambda [ exp(-(t-\lambda)x) \times \frac{1}{-(t-\lambda)} ]_{0}^{\infty} = \lambda ( \frac{1}{\lambda - t } ) = \frac{\lambda}{\lambda - t }

従って、X+Yのモーメント母関数はM_{X+Y}(t) = M_{X}M_{Y} = {M_{X}}^2
\Large = ( \frac{\lambda}{\lambda - t } )^2


ところで、もし再生性があるのであれば和の分布のモーメント母関数は \Large \frac{2\lambda}{2\lambda - t } となるはずだが、明らかに
\Large ( \frac{\lambda}{\lambda - t } )^2 \neq \frac{2\lambda}{2\lambda - t } である。指数分布は再生性を持たないことが確認できた。

3) 1)で求めた和の分布のモーメント母関数を求め、これが2)と一致することを確認する

1)で求めた、Zの確率密度関数 f(z) = \Large \rm \lambda^2 z e^{-\lambda z} のモーメント母関数を求める。
なお積分範囲(0,∞)は記載省略。

Zのモーメント母関数は、
\rm M_{Z}(t) = E[ exp(tZ) ] = \int exp(tz) \lambda^2 exp(-\lambda z)z dz = \lambda^2 \int exp(-(\lambda - t) z )z dz
\rm = \lambda^2 \int (\frac{1}{-(\lambda - t) } exp(-(\lambda - t )z) )' z dz
\rm = \lambda^2 \Bigl( \Bigl[ \frac{1}{-(\lambda - t )}e^{-(\lambda - t)z} z \Bigr]_{0}^{\infty} - \int (\frac{1}{-(\lambda - t) } exp(-(\lambda - t )z)dz \Bigr)
\rm = \frac{\lambda^2}{\lambda - t } \Bigl( \int (\frac{1}{-(\lambda - t) } exp(-(\lambda - t )z)dz \Bigl)
\rm = \frac{\lambda^2}{\lambda - t } \Bigl( \Bigl[ \frac{1}{-(\lambda - t )}e^{-(\lambda - t)z} \Bigr]_{0}^{\infty}
\rm = \frac{\lambda^2}{ - (\lambda - t)^2 } ( 0 - 1 )
\rm \Large = (\frac{\lambda}{ \lambda - t } ) ^2

となる。これは確かに2)で求めたモーメント母関数に一致する。

2019年6月 統計検定準一級 問8

勉強

問題文

リンクが切れていないなら問題文を参照。

解く時の気持ち

Lassoのことをチョット知っているぐらいのお気持ちで、
機械学習ガチ勢ではなく「統計検定準1級目指してるけど機械学習よくわからん」なスタンスです。

問題の大雑把な概要

  • 中性子の月次カウントデータ217点が与えられている( \displaystyle y_i )
  • Fused Lasso を用いて \displaystyle y_i を平滑化した実数列( \displaystyle \beta_i ) を生成する。
    • Fused Lasso: \displaystyle \hat{\beta} = \arg \underset{\beta \in \mathbb{R}^{217} }\min \frac{1}{2} \sum_{i=1}^{217}{(y_i-\beta_i )}^2 + \lambda \sum_{i=1}^{216}{ | \beta_{i+1} - \beta_i | }
  • 小問〔1〕：
    • \lambda = 500 の場合にFused Lassoで平滑化した結果を表す図を選ぶ。
  • 小問〔2〕：
    • 図が与えられているので適切な平滑化手法(Fused Lasso) を選ぶ。

小問〔1〕：

Lassoは

  • 罰則項 \displaystyle \lambda \sum | \cdots | の 各 | \cdots | を、0にしたりしなかったりするいい感じの回帰分析。
  • λを大きくすると0になる各 | \cdots | は多くなり、λが小さかったら0になる各 | \cdots | は少なくなる。

という考えのもと小問〔1〕を見つめる。

ところで今回の \betaは通常のLassoの回帰係数ではなく「 \displaystyle y_i を平滑化した実数列( \displaystyle \beta_i ) 」ということを念頭におくと、得られるβはyにほぼほぼ近いがいい感じに平滑化されているものだと考えることができる。
また、罰則項の中身が \beta_{i+1} - \beta_i ということは、第i+1項と第i項の差を表している。ということで、

  • | \beta_{i+1} - \beta_i | = 0 ならば隣接する項の差が0 ⇔ グラフの横軸に平行

と解釈できる。従って、λ=500よりFused Lassoで得られるβは、0になる | \beta_{i+1} - \beta_i | が多い→横軸に平行な線が多い。

と考えて、答えは４となる。

小問〔2〕：

Lassoの問題なので、

  • 罰則項はきっと \displaystyle \lambda \sum | \cdots | の形をしているであろう
  • 誤差項は２乗誤差だよねきっと
  • 平滑化と言っているし誤差項の絶対値の中身は、差の形になっている

を踏まえて → 4,5 に絞る。

次に、4と5の違いを考えると

  • 4は
    • \beta_{i+2} - 2\beta_{i+1} + \beta_{i} = ( \beta_{i+2} - \beta_{i+1} )- (\beta_{i+1} - \beta_{i})
    • βの３項間の傾きについて述べている。
  • 5は
    • \beta_{i+3} - 3\beta_{i+2} + 3\beta_{i+1} -\beta_{i} = \bigl( (\beta_{i+3} - \beta_{i+2}) - (\beta_{i+2} - \beta_{i+1})\bigr) - \bigl( (\beta_{i+2} - \beta_{i+1}) - (\beta_{i+1} - \beta_{i})\bigr)
    • βの４項間の傾きについて述べている。


となることから、Lassoで得られたβによる傾き( \beta_{i+1} - \beta_{i} )が３点の間で等しくなりがちな選択肢4が図に適していそうである。
5の場合は \displaystyle \sum | \cdots | の中のパラメータの数が選択肢4よりは複雑なので図ほどシンプルな曲線は描けずもう少し滑らかとなるであろう。

二変量正規分布の条件付き期待値・条件付き分散

統計

以下のような二変量正規分布を考える。

\begin{bmatrix} X \\ Y \end{bmatrix} \sim N( \begin{bmatrix} \mu_{X} \\ \mu_{Y} \end{bmatrix} , \begin{bmatrix} \sigma_{X} & \sigma_{XY} \\ \sigma_{XY} & \sigma_{Y} \end{bmatrix} )

ここで、

Z = X + AY, A = - \sigma_{XY} {\sigma_{Y}} ^{-1}

というZについて、 Cov(Z, Y) を計算すると、

cov(Z,Y) \\ = cov(X+AY,Y) \\ = cov(X,Y) + cov(AY,Y) \\ = \sigma_{XY} + Acov(Y,Y) \\ = \sigma_{XY} - \sigma_{XY} {\sigma_{Y}} ^{-1} \sigma_{Y} \\ = \sigma_{XY} - \sigma_{XY} \\ = 0

となる。 ゆえにX,Y,Zは正規分布に従う確率変数であることから、ZとYは共分散が0⇔ZとYは独立となる。

次に、 E[ X | Y ] を計算する。ここで X = Z - AY を用いると、

E[ X | Y ] \\ = E[ Z - AY | Y ] \\ = E[ Z | Y ] - E[ AY | Y ] \\ = E[ Z ] - A E[ Y | Y ] \\ = E[ X + AY ] - AY \\ = E[ X ] + AE[Y] - AY \\ = \mu_{X} + A\mu_{Y} - AY \\ = \mu_{X} - A( Y - \mu_{Y} ) \\ = \mu_{X} + \sigma_{XY} { \sigma_{Y} }^{-1} ( Y - \mu_{Y} )

となる。*1

最後に、 V[ X | Y ] を計算する。

V[ X | Y] \\ = V[ Z - AY | Y ] \\ = V[ Z | Y ] + V[AY | Y] - 2Cov(Z,AY) \\ = V[ Z | Y ] + V[AY | Y] - 2ACov(Z,Y) \\ = V[ Z | Y ] + 0 + 0\\ = V[ Z ]

より、*2

V[ Z ] \\ = V[ X + AY ] \\ = V[X ] + V[ AY ] + 2Cov(X,AY) \\ = V[X ] + A^2V[ Y ] + 2ACov(X,Y) \\ = \sigma_{X} + A^2\sigma_{Y} + 2A\sigma_{XY} \\ = \sigma_{X} + (- \sigma_{XY} {\sigma_{Y}} ^{-1})^2 \sigma_{Y} + 2(- \sigma_{XY} {\sigma_{Y}} ^{-1}) \sigma_{XY} \\ = \sigma_{X} + {\sigma_{XY}}^2 {\sigma_{Y}}^{-1} - 2{\sigma_{XY}}^2 {\sigma_{Y}}^{-1} \\ = \sigma_{X} - {\sigma_{XY}}^2 {\sigma_{Y}}^{-1}

となる。

参考

stats.stackexchange.com

*1:途中、性質１：XとYが独立ならばE[X|Y] = E[X] と性質２：E[X|X] = X を用いた。

*2:条件付き分散の性質がよくわからない・・・。リンク参照。

正規分布と適合度検定

統計 R

適合度検定 :: 株式会社アイスタット｜統計分析研究所より「適合度の検定（正規性）の結果」をRにて計算してみる。

# パラメータ####
# 平均
m1  <- 64.5
# 標準偏差
sd1 <- 13.41
# 度数の総和
n   <- 40

# 関数 ####
# 区間a<x<bにおける標準正規分布に従うXの確率、を返す関数
f <- function(l,h){ pnorm(h,m1,sd1) - pnorm(l,m1,sd1) }

# データ
# 観測値：observed ####
c(2,4,7,13,10,3,1) -> obs

# 期待度数expected ####
n * c(
pnorm(40, m1, sd1)
,f(40,50)
,f(50,60)
,f(60,70)
,f(70,80)
,f(80,90)
,(1 - pnorm(90, m1, sd1) ) 
) -> expected

# 自由度 = 階級の個数 - 1 - 2
df1 <- 7-1-2

# χ²統計量
chi2 <- sum( (obs - expected)^2/expected )
print( chi2 )
#> print( chi2 )
#[1] 1.382639

# p値
1 - pchisq( chi2, df1) 
#> 1 - pchisq( chi2, df1) 
#[1] 0.8472068
axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.