axjack's blog

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

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

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

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


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

1日の死亡者数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)に従う。

忘れては思い出したいポアソン分布

スマホで文字打ち込んでいるのでLaTeXに出来ませんがまぁご愛嬌。

ポアソン分布の確率関数

Pr(X=k) = exp(-θ)θ^k/k!
ただしk は0以上の整数

ポアソン分布の期待値

E[X] = θ

ポアソン分布の期待値の導出

overはΣの添字が走る範囲です。またexp(θ)=Σθ^h/h! over h≧0を用います。

E[X]
= Σkexp(-θ)θ^k/k! over k≧0
= Σkexp(-θ)θ^k/k! over k≧1
= exp(-θ)Σθ^k/(k-1)! over k≧1
= θexp(-θ)Σθ^(k-1)/(k-1)! over k-1≧0
= θexp(-θ)Σθ^h/h! over h≧0
= θexp(-θ)exp(θ)
= θ

Rで階層的クラスタリング

あけましておめでとうございます今年も統計学とRの勉強を継続します。

ということで、いつもの通り(?)Yuyaさんの動画にて階層的クラスタリングを見たので、Rを使ってクラスタリングしてみます。

www.youtube.com

コード

# データの準備
dd <- data.frame(
  ht = c(2,1,2,4,6)
  ,sp = c(1,4,6,1,5)
)
row.names(dd) <- c('A','B','C','D','E')

# 距離行列
dd.dst <- dist(dd)

# 階層的クラスターの作成
hc.comp <- hclust(dd.dst,method = "complete")  # 最遠隣法
hc.sing   <- hclust(dd.dst,method = "single") # 最近隣法
hc.wd     <- hclust(dd.dst,method = "ward.D2") # ウォード法

# グラフの作図
dev.off()
par(mfrow = c(3,1) )
par(family = "HiraKakuProN-W3")
for(i in 1:3){
  plot(hc.comp,main = "最遠隣法")
  plot(hc.sing,main = "最近隣法")
  plot(hc.wd, main = "ウォード法")
}

グラフ

f:id:axjack:20210103213509p:plain

参考

Rによる階層的クラスタリング

データサイエンス100本ノックをRでやってみた

github.com

1日1~2時間×1週間ぐらいで終えました。

感想

  • tidyverseに感謝
  • 日付の計算はstrptimeよりlubridate使った方が簡単
  • inner_joinはbyを指定しないときnatural joinになって便利
  • 正規表現を使う問題はstr_subやstr_detectを用いても代替可能(問題の趣旨に反するが・・・。)
  • slice_*系でデータ取り出し、が失敗することもあるのでそんな時はheadを利用

つまずいた問題たち

供養と復習を兼ねて書いておきます。どの問題も実用的なのでスラスラと書けるようになりたいですね。

R-029

R-029: レシート明細データフレーム(df_receipt)に対し、店舗コード(store_cd)ごとに商品コード(product_cd)の最頻値を求めよ。

最頻値ってどうやって出すのっていうか最頻値出す関数無くない?と一時停止した問題です。inner_joinを駆使してなんとかひねり出した自分の答えに対してAnswerの解法はエレガントでした。勉強も兼ねてrecodeしたものが以下です。

df_receipt %>%
  group_by(store_cd, product_cd) %>%
  summarise(cnt = n(), .groups = 'drop') %>%
  group_by(store_cd) %>%
  filter( cnt == max(cnt) ) # Ans の n == n %>% max() は、 n == max(n) と同じ。

R-058

R-058: 顧客データフレーム(df_customer)の性別コード(gender_cd)をダミー変数化し、顧客ID(customer_id)とともに抽出せよ。結果は10件表示させれば良い。

ダミー変数化って要するにどうなるんだっけ?と思考停止した問題です。変数の取りうる値の範囲が事前に分かっていればdecodeすれば良いのだけど、そうではない時にはAnswerにある通りdummyVars()を使うそうです。

# ダミー変数化がわからなかったのでAnsを参照した。
# すると、gender_cd = {0,1,9} を gender_cd0,gender_cd1,gender_cd9と3つの変数に分離することらしい。
# 横持ち形式ということかなぁ。
# で、dummyVars()というのは{caret}パッケージの関数であるもよう。

# if_else で書けるといえば書けるのだが、
# よく考えると「取りうる値」が事前に分かっていないとコーディングできないね。
# ということで、「取りうる値」が分からない場合はdummyVars()を使うしかないわね。
df_customer %>% 
    mutate(
        gender_cd0 = if_else(gender_cd == 0,1,0)
        ,gender_cd1 = if_else(gender_cd == 1,1,0)
        ,gender_cd9 = if_else(gender_cd == 9,1,0)
        ) %>%
    select(customer_id, starts_with("gender")) %>%
    slice(1:10)

R-087

R-087: 顧客データフレーム(df_customer)では、異なる店舗での申込みなどにより同一顧客が複数登録されている。名前(customer_name)と郵便番号(postal_cd)が同じ顧客は同一顧客とみなし、1顧客1レコードとなるように名寄せした名寄顧客データフレーム(df_customer_u)を作成せよ。ただし、同一顧客に対しては売上金額合計が最も高いものを残すものとし、売上金額合計が同一もしくは売上実績の無い顧客については顧客ID(customer_id)の番号が小さいものを残すこととする。

SQLだったらwindow関数とpartition by でうまくいくのになぁとモヤモヤした問題です。Answerはとてもエレガントでdistinct関数をうまく利用した解法です。自分で出した解答は確かに問題文の題意をそのまま愚直に書き下した感がありエレガントではないでした。一番悔しかった問題です。

df_receipt %>% 
    group_by(customer_id) %>% 
    summarise(sum_amount = sum(amount) , .groups = 'drop' ) -> dfr

df_customer %>% left_join(dfr) %>% 
    mutate(sum_amount = replace_na(sum_amount,0) ) %>%
    group_by(customer_name, postal_cd) %>%
    mutate(
        rank_amount = min_rank(desc(sum_amount))
        ,rank_cust = min_rank(customer_id)
        ,num_amount = n_distinct(sum_amount) # 売上金額合計の種類の数
    ) %>%
    filter(  ( (num_amount > 1) & (rank_amount == 1) ) | ( (num_amount == 1)&( rank_cust == 1 )   )  ) %>%
    select(-rank_amount,-rank_cust,-num_amount) -> df_customer_u2

「 平均点以上の人の平均点」を求めるための条件付き確率密度関数が確率密度関数であることを確認する

問題の引用元

www.youtube.com
www.youtube.com

問題の要約

  •  \rm{X} \sim N(60,10^2)
  • Q1: \rm{X}の第1四分位・第3四分位を求めよ
  • Q2:  \rm{E}[ X | X \ge 60 ] を求めよ

解(の途中まで)

Q1
  • Xは正規分布に従う。ひとまず Z = \frac{X - 60}{10}と標準化して考える。
  • 第1四分位は累積密度関数の下側25%点なので、 Z = \Phi^{-1} (0.25) = -0.674。
    • 逆標準化?して X = 10 \times (-0.674) + 10 = 53.2551
  • 第3四分位は累積密度関数の下側75%点なので、 Z = \Phi^{-1} (0.75) = 0.674。
    • 逆標準化?して X = 10 \times (0.674) + 10 = 66.7449
Q2
  • Xは正規分布に従う。ひとまず Z = \frac{X - 60}{10}と標準化して考える。
  • 標準化すると、 \rm{E}[ X | X \ge 60 ]  \rm{E}[ Z | Z \ge 0 ] と考えればよい(※同一ではないので最後に逆標準化する)。

さて、大変雑であるが…求めたいのは標準正規分布の原点から右側の期待値である。

すると、期待値は \displaystyle E[X|X\ge0 ] =  \displaystyle  \int_{-\infty}^{\infty} z \frac{ f(z|z \ge 0) }{P(z \ge 0 ) } dz    を計算すれば求まる。ここでfは標準正規分布でPは標準正規分布から計算する確率である。

ここまで(から実際の答え導出まで)の流れは上述の動画で分かりやすく解説してあるので動画を参照すれば大丈夫です。

条件付き確率密度関数確率密度関数であることの確認

ところで、 \displaystyle \frac{ f(z|z \ge 0) }{P(z \ge 0 ) }   確率密度関数なのだろうかという点に私は疑問を持ったので略証してみる。確率密度関数だったら実数の範囲で積分して1になることを示せばよい。


 \displaystyle  \int_{-\infty}^{\infty}  \frac{ f(z|z \ge 0) }{P(z \ge 0 ) } dz
 \displaystyle  = \int_{-\infty}^{\infty}  \frac{ f(z) \unicode{x1D7D9}(z \ge 0)   }{P(z \ge 0 ) } dz
 \displaystyle  = \int_{0}^{\infty}  \frac{ f(z) }{P(z \ge 0 ) } dz
 \displaystyle  = \frac{1}{ P(z \ge 0 ) }  \int_{0}^{\infty} f(z)  dz
 \displaystyle  = \frac{1}{ P(z \ge 0 ) }  P(z \ge 0 )
 \displaystyle  = 1

となるので確かに確率密度関数となっている。

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

やること

『統計学』東京図書 p.147 練習問題 問5.1の[2]、多変量正規分布の条件付き期待値・条件付き分散を解きます。

問題

(X,Y,Z)^{T}が平均ベクトル(0,0,0)^{T}、分散共分散行列が正値対称行列

 \begin{pmatrix}
1 & \rho_{xy} & \rho_{xz} \\
\rho_{xy} & 1 & \rho_{yz} \\
\rho_{xz} & \rho_{yz} & 1 
\end{pmatrix}

の時、X = xおよびZ = zを与えた下でのYの

  • 条件付き期待値
  • 条件付き分散

を求めよ。

解法

※ mathjaxだと面倒なので画像で行きます。

確率ベクトルX: p×1行列 が 平均ベクトルμ: p×1行列, 分散共分散行列Σ: p×p行列を用いた N(μ, Σ)に従っている。そして『統計学』東京図書 p.45 にあるように

f:id:axjack:20200911224418p:plain

と区分けすることを考える。

この時、『統計学』東京図書 p.46の定理2.2.6により、

f:id:axjack:20200911224838p:plain

が成り立つ。

これを踏まえると、今回の問題はX = (Y,X,Z)^{T}と順番を入れ替えて

f:id:axjack:20200911225754p:plain

と区分けすることができる。(赤線の区分けと添字の対応は以下の通り。)

f:id:axjack:20200911230603p:plain

ということで、E[Y | X = x, Z = z]を求めると、

f:id:axjack:20200911232935p:plain

となる。

一方、V[Y | X = x, Z = z]を求めると、

f:id:axjack:20200911234042p:plain

となる。

なお今回の3変量正規分布の条件付期待値と分散は、変数2つが固定されるので一次元正規分布となる。なので、ベクトル・行列の演算結果はスカラーになる。

参考

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