axjack's blog

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

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

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

統計学実践ワークブック問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

ポアソン分布と適合度検定

統計 R

準1級 例題／解説 の問2より。

ある地域における１日の死亡者数の集計結果表

\begin{array}{c} 死亡者数(人) & 0 & 1 & 2 & 3 & 4 & 5 & 6以上 &計\\ \hline 件数(日数) & 55 & 144 & 140 & 95 & 45 & 15 & 6 &500\\ \end{array}

１日の死亡者数Xがパラメータλのポアソン分布に従うと仮定する。ある日の死亡者数が3人である確率は？

\displaystyle Poisson(x; λ) = \frac{ \lambda^x e^{- \lambda} } { x! }

に於いて、λ=3を代入すればよい。

同分布を仮定した時、E[X2] とλの関係は？

E[X] = \lambda, V[X] = \lambda = E[X^2] - E^2[X]

より、

E[X^2] = \lambda + {\lambda}^2

パラメータλの推定値は？

\hat{ \lambda } = \sum x_i \hat{p}_i

より求める。ここで \hat{p} は 件数を件数総和の500 で割ったものとする。Rで計算すると、

# 人
x <- 0:6
# 件数
v <- c(55,144,140,95,45,15,6)
# 割合
vprop <- v/sum(v)

# 標本平均
lam_hat <- sum(x * vprop) 
print(lam_hat)
#> print(lam_hat)
#[1] 2

より、 \hat{ \lambda} = 2 となる。

パラメータの推定値から期待度数を求めよ

6以上に注意してRで計算すると、

# ポアソン分布(λ=2)からX = 0,1,2,3,4,5までの確率を求める。
dpois05 <- dpois(0:5, lambda = 2)

# 6以上の確率
dpoisGE6 <- 1 - sum(dpois05)

# 期待度数Expected_i = n * p_i。この問題ではnは500
round( 500 * c(dpois05, dpoisGE6), 1)
#> round( 500 * c(dpois05, dpoisGE6), 1)
#[1]  67.7 135.3 135.3  90.2  45.1  18.0   8.3

期待度数をもとに、適合度検定を実行せよ

で適合度検定を実行する。

\displaystyle {\chi_0}^2 = \sum \frac{ { (O_i - E_i) }^2 }{E_i}

より、カイ二乗値は

expected <- round( 500 * c(dpois05, dpoisGE6), 1)
sum( (expected - v )^2/expected )
#[1] 4.498116

ところで自由度は, 7-1-1 より5なので自由度5のカイ二乗分布の上側5%点(⇔下側 95%点)を求めると、

qchisq(1-0.05, 5)
#[1] 11.0705

となるので、4.498116 < 11.0705 より有意水準5%において棄却域に入っていない。よってポアソン分布に従っていることを棄却することはできない*1

*1:ポアソン分布に従っていないとは言えない。言い回しが難しい。。

メモ：Anacondaを使わずにJuliaとJupyterとJupyterLabをインストールする

プログラミング

Anacondaが重かったりアンインストールが面倒くさかったりなんだりで宗教的にAnacondaは受け付けない・・・という人向けのメモです。

  • 環境
  • python
    • /usr/local/bin/python3.8
  • Juliaのインストール
    • 公式サイトからダウンロードしてインストール
    • https://julialang.org/
    • Julia 1.5.3がインストールされた
  • いちどパスの整理
    • python と Juliaのパスが.zshrcに入っているか確認。なければ追加。ないと何かしらのエラーでこける。
    • /Applications/Julia-1.5.app/Contents/Resources/julia/bin:/usr/local/bin/python3.8 などを追加。
  • pythonを使ってJupyterをインストール
    • python3.8 -m pip install jupyter
  • python を使ってJupyterLabをインストール
    • どうやらJupyterをインストールしてもJupyterLabはインストールされないらしい
    • python3.8 -m pip install jupyterlab
  • ターミナルから起動
    • jupyter notebook でjupyterが起動する
    • jupyter lab でjupyterlabが起動する

標本平均からの偏差二乗和をσ²で割ったものが自由度n-1のカイ二乗分布に従う例のアレの導出

勉強 数学 統計

厳密ではないが、カイ二乗分布の再生性やコクランの定理から言えるのではないでしょうかね(無責任)。

導出

Xᵢ ~ N(μ,σ²)に従うiid確率変数とする。
このとき、Zᵢ = (Xᵢ-μ)/σと変換すると、
Z²=ΣZᵢ² over 1 to nはχ²(n) に従う。


ここで、
Σ(Xᵢ-μ)²
= Σ{(Xᵢ-X̅)+(X̅-μ)}²
= Σ{(Xᵢ-X̅)²+2(Xᵢ-X̅)(X̅-μ)+(X̅-μ)²}
= Σ(Xᵢ-X̅)²+2(X̅-μ)Σ(Xᵢ-X̅)+Σ(X̅-μ)²
= Σ(Xᵢ-X̅)²+Σ(X̅-μ)²
= Σ(Xᵢ-X̅)²+n(X̅-μ)²

と変形できることから、

Z²
= Σ{(Xᵢ-μ)/σ}²
= (1/σ²){Σ(Xᵢ-X̅)²+n(X̅-μ)²}
= Σ{(Xᵢ-X̅)/σ}²+{(X̅-μ)/(σ/√n)}²
と表せる。


ところで、X̅~N(μ,σ²/n)であるから
{(X̅-μ)/(σ/√n)}²~χ²(1)である。

よって
Z²~χ²(n)と{(X̅-μ)/(σ/√n)}²~χ²(1)とカイ二乗分布の性質から、
Σ{(Xᵢ-X̅)/σ}²〜χ²(n-1)に従う。
axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.