axjack's blog

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

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

勉強 統計

スマホで文字打ち込んでいるので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 統計

あけましておめでとうございます今年も統計学と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でやってみた

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}の第１四分位・第３四分位を求めよ
  • Q2: \rm{E}[ X | X \ge 60 ] を求めよ

解(の途中まで)

Q1
  • Xは正規分布に従う。ひとまず Z = \frac{X - 60}{10}と標準化して考える。
  • 第１四分位は累積密度関数の下側25%点なので、 Z = \Phi^{-1} (0.25) = -0.674。
    • 逆標準化？して X = 10 \times (-0.674) + 10 = 53.2551
  • 第３四分位は累積密度関数の下側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

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

参考

ryosuke-okubo.hatenablog.com

bellcurve.jp

stats.stackexchange.com

stats.stackexchange.com

math.stackexchange.com

stats.stackexchange.com

en.m.wikipedia.org

en.wikipedia.org

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

確率 線形代数 勉強 数学 統計

やること

『統計学』東京図書 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

となる。

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

参考

sweep関数の使い方を忘れがちなので自分へのメモ

R プログラミング

まとめ

sweep関数の

  • MARGIN = 1の時は
    • i 行目のそれぞれの値に対してSTATSの i番目の値をFUNと二項演算する
  • MARGIN = 2の時は
    • j 列目のそれぞれの値に対してSTATSの j番目の値をFUNと二項演算する

お気持ち

行平均の引き算のお気持ち

  • sweepで1行目の行平均を引き算したいとき、1行目の各要素から引き算する値はスカラーな行平均。
  • sweepで2行目の行平均を引き算したいとき、2行目の各要素から引き算する値はスカラーな行平均。
  • ...

となる。sweep(mat, MARGIN = 1, ...)の時、1行目ベクトルからSTATSベクトルを減算し、2行目ベクトルからSTATSベクトルを減算し・・・とはならないので要注意。 忘れたら行平均の減算、列平均の減算を思い出そう。

確認コード

MARGIN = 1の場合

> # 行列を生成し
> (m <- matrix(c(3,1,4,1,5,9,2,6,5,3,5,9),ncol=4))
     [,1] [,2] [,3] [,4]
[1,]    3    1    2    3
[2,]    1    5    6    5
[3,]    4    9    5    9
> # 行平均を確認
> rowMeans(m)
[1] 2.25 4.25 6.75
> # 行平均を各行から引く
> sweep(m, MARGIN = 1, STATS = rowMeans(m), FUN = "-")
      [,1]  [,2]  [,3] [,4]
[1,]  0.75 -1.25 -0.25 0.75
[2,] -3.25  0.75  1.75 0.75
[3,] -2.75  2.25 -1.75 2.25
>

MARGIN = 2 の場合

> # 列平均を確認
> colMeans(m)
[1] 2.666667 5.000000 4.333333 5.666667
> # 列平均を各列から引く
> sweep(m, MARGIN = 2, STATS = colMeans(m), FUN = "-")
           [,1] [,2]       [,3]       [,4]
[1,]  0.3333333   -4 -2.3333333 -2.6666667
[2,] -1.6666667    0  1.6666667 -0.6666667
[3,]  1.3333333    4  0.6666667  3.3333333
>

となる。

雑だけどirisを機械学習する

R

まえがき

Data Science: Machine Learningを受講しています。今までのR BasicsやVisualizationに比べて課題が多くてなかなか進まないです。8月中には完了させたいですね。

本題ですが、今回は「雑だけどirisを機械学習する」ということで、とりあえずcaret使って機械学習を動かしてみようと思います。理論が分からないと機械学習なんて触っちゃダメだろーと思っていましたが、まぁまずはtryですよ。

ソースコード

# テストデータで10行
# 訓練データはそれ以外
ind <- sample(iris %>% nrow(), size = 10, replace = FALSE)
iris_test  <- iris[ind, ]
iris_train <- iris[-ind,]

# trainのメソッドを欲張りセットでこしらえる
methods   <- list("lda","qda","knn","rf","nnet")

# Accuracyを出す
accuracies <- map_dbl(methods, function(me){
  trained   <- train(Species ~ ., data = iris_train, method = me )
  predicted <- predict(trained, iris_test)
  confusionMatrix(data = predicted, reference = iris_test$Species)$overall["Accuracy"]  
})

# ベクトルに名前をつけて
names(accuracies) <- methods

# 表示
accuracies

結果

2,3回ぶん回しましたがだいたい8~9割ぐらいのAccuracyが出ました。

> accuracies
 lda  qda  knn   rf nnet 
 1.0  1.0  0.8  0.8  0.9

２次元正規分布のデータから２通りの方法で回帰直線を引く

R 統計

タイトルの通りです。

  • lmで単回帰直線
  • 平均と分散共分散からE[Y|X=x] *1を計算

の２通りで回帰直線を引きます。

コード

library(MASS)
library(scales)
options(digits = 3)

# 平均
mx <- 10 
my <- 22
# 分散
Sx <- sqrt(9)
Sy <- sqrt(16)
Sxy <- sqrt(9.6)
# サンプルサイズ
n <- 10000

# 平均ベクトル
M <- c(mx, my)

# 分散共分散行列
Sig <- c(Sx^2,Sxy,Sxy,Sy^2) %>% matrix(2,2)
z <- MASS::mvrnorm(n, mu = M, Sigma = Sig )  

# dev.off()
# 散布図プロット
plot(z[,2] ~ z[,1], xlab = "x", ylab = "y",col = c(alpha('orange', 0.2)) )

# 単回帰直線(青)
lm.z <- lm(z[,2] ~ z[,1])
coef(lm.z)
abline(lm.z, col ="blue", lwd = 2, lty = 2)

# E[Y|X = x]で直線を引く(赤)
list("Intercept" = my - mx * (Sxy/(Sx*Sy))*(Sy/Sx), "Slope" = (Sxy/(Sx*Sy))*(Sy/Sx) ) %>% unlist()
abline(my - mx * (Sxy/(Sx*Sy))*(Sy/Sx), (Sxy/(Sx*Sy))*(Sy/Sx), col="red", lwd = 3, lty = 2)

結果

だいたい一致しました。青が単回帰直線、赤がE[Y|X = x]です。 f:id:axjack:20200817224258p:plain

  • lmでの係数
> coef(lm.z)
(Intercept)      z[, 1] 
     18.552       0.344
  • E[Y|X = x]での係数
> list("Intercept" = my - mx * (Sxy/(Sx*Sy))*(Sy/Sx), "Slope" = (Sxy/(Sx*Sy))*(Sy/Sx) ) %>% unlist()
Intercept     Slope 
   18.557     0.344

参考

改訂版 日本統計学会公式認定 統計検定2級対応 統計学基礎 のp.86 (2.9.18)

*1:母数から計算しているのでおそらく母回帰直線、と思われる。

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