axjack's blog

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

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

やること

『統計学』東京図書 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つが固定されるので一次元正規分布となる。なので、ベクトル・行列の演算結果はスカラーになる。

参考

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

まとめ

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を機械学習する

まえがき

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
axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.