2020-08-28

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
>

となる。

2020-08-24

雑だけど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

2020-08-17

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

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:母数から計算しているのでおそらく母回帰直線、と思われる。

2020-08-04

標準誤差は統計量の標準偏差であることを確かめる

R 確率統計

やること

離散型確率分布

X	p
17	0.3
-1	0.7

に於いて、大きさn=100の標本を抽出する時、

和S の標本分布について、Sの平均、標準誤差
平均M の標本分布について、Mの平均、標準誤差

をそれぞれ

計算
モンテカルロシミュレーション

によりそれぞれ求める。

計算：離散型確率分布から分かること

Xの期待値・分散・標準偏差はそれぞれ、

$E(X) = \sum X_i p_i = 17\times0.3 + (-1) \times 0.7 = 4.4$
$V(X) = {SD(X)}^{2} \approx {8.249}^{2} \approx 68.04$
$SD(X) = |X_1 - X_2|\sqrt{p(1-p)} = |17 - (-1) |\sqrt{0.3\times0.7} \approx 8.249$

であるから、サンプルサイズn=100における、

和S = $\sum{X_i}$ の標本分布について
- 各 $X_i$ は $X_i \sim X$ として、
- Sの期待値:= $E(S) = E(\sum X_i ) = E(nX_i) = nE(X) = 100 \times 4.4 = 440$
- Sの標準誤差:= $SE(S) = \sqrt{ V(S) } = \sqrt{ V(\sum{X_i} ) } = \sqrt{ \sum{V(X_i)} } = \sqrt{ \sum{V(X)} } = \sqrt{ nV(X)} = \sqrt{n}\sqrt{V(X)} = \sqrt{n} SD(X) = 10\times8.249 = 82.49$
平均M = $\frac{1}{n}\sum{X_i}$ の標本分布について
- 各 $X_i$ は $X_i \sim X$ として、
- Mの期待値:= $E(M) = E(\frac{1}{n}\sum X_i ) = E(X) = 4.4$
- Mの標準誤差:= $SE(M) = \sqrt{ V(M) } = \sqrt{ V( \frac{1}{n} \sum{X_i} ) } = \sqrt{ \frac{1}{n^{2} } \sum{ V(X_i) } } = \frac{1}{n} \sqrt{ \sum{ V(X)} } = \frac{1}{n} \sqrt{n}\sqrt{V(X)} = \frac{1}{\sqrt{n}} SD(X) = \frac{8.249}{\sqrt{100}} = 0.8249$

である。

モンテカルロシミュレーション

共通の設定

X <- c(17,-1) #確率変数
p <- c(0.3,0.7) #確率
n <- 100 #サンプルサイズ
B <- 10000 #サンプル"数"

和の分布

# 和の分布 ####
S <- replicate(B,{
  sum(sample(X,size=n,replace=TRUE,prob=p))
})

# 結果
c(平均値 = mean(S), 標準誤差 = sd(S))

> c(平均値 = mean(S), 標準誤差 = sd(S))
   平均値  標準誤差 
438.75620  82.56416

平均値の分布

# 平均値の分布
M <- replicate(B,{
  mean(sample(X,size=n,replace=TRUE,prob=p)) 
})

# 結果
c(平均値 = mean(M), 標準誤差 = sd(M) )

> c(平均値 = mean(M), 標準誤差 = sd(M) )
   平均値  標準誤差 
4.4000180 0.8259942

おわりに

やってることは例によってChapter 13 Probability | Introduction to Data Scienceで学んだことです。標準誤差は統計量の標準偏差であるってのは数式だけだといまいちピンとこないのですが、このようにRでシミュレーションするとだいたい同じになるので腑に落ちました。replicate関数の第一引数にサンプル"数"を、第二引数には統計量(の式)を記述することによって、統計量の分布をシミュレーションすることができます。便利です。

おまけ: replicate関数の練習

(これを最初に書くべきだったのでは？)サイコロを２回投げて出た目の和が7となる確率、をRで求める。

共通

saikoro <- 1:6 # サイコロ
B <- 10000 #サンプル"数"

組み合わせで割り算する素朴な求め方

# サイコロを２回投げた時の組み合わせ
saikorokumi <- expand.grid(saikoro,saikoro)

# 和が7になる組み合わせ ÷ サイコロを２回投げた時の組み合わせの総数
sum( apply(saikorokumi,1,sum) == 7 ) / nrow(saikorokumi)

> sum( apply(saikorokumi,1,sum) == 7 ) / nrow(saikorokumi)
[1] 0.1666667

２回投げて出た目の和の分布

# サイコロを２回投げて出た目の和の分布
kiroku <- replicate(B,{
  sum(sample(saikoro,2,replace=TRUE))
})
# 和が7となる確率
mean(kiroku == 7)

> mean(kiroku == 7)
[1] 0.1659

２回投げて出た目の和が７となるか否かの分布

# サイコロを２回投げて出た目の和が７となるか否かの分布
kiroku2 <- replicate(B,{
  sum(sample(saikoro,2,replace=TRUE)) == 7 
})
# 和が7となる確率
mean(kiroku2)

> mean(kiroku2)
[1] 0.1684

2020-08-03

確率変数の取りうる値が２つの時の離散型確率分布の標準偏差を求める公式

統計確率

結論

確率変数 $X = {X_1, X_2}$ においてそれぞれ確率 $p,1-p$ をとる時、標準偏差SD(X)は

$\displaystyle SD(X) = |X_1 - X_2| \sqrt{ p(1-p) }$

となる。ここで||は絶対値の記号である。

どこで見つけたか

edXのprobabilityの教科書？であるChapter 14 Random variables | Introduction to Data Science やここである。

f:id:axjack:20200803234358p:plain *1 f:id:axjack:20200803234436p:plain *2

計算例

離散型確率分布が以下の表、

X	p
17	0.3
-1	0.7

で与えられているとする。

まずは期待値・分散・標準偏差を計算し、最後に公式を用いて標準偏差を計算して２つの標準偏差が一致することを確かめる。

期待値

$\mu = E(X) = \sum(X = X_i)p_i = 17\times0.3 + (-1)\times0.7 = 4.4$

分散

$V(X) = E({( X-\mu) }^2 ) \\ = \sum {( X_i - \mu)}^2p_i = {(17 - 4.4)}^2\times 0.3 + {(-1 - 4.4)}^2\times0.7 \\ = 68.04$

標準偏差

分散のルートを取って求めると、

$SD(X) = \sqrt{ V(X) } = \sqrt{ 68.04} \approx 8.248636$

一方、冒頭の公式で求めると、

$SD(X) = |X_1 - X_2 |\sqrt{ p(1-p) } = | 17 - (-1) | \sqrt{ 0.7 \times 0.3 } \approx 8.248636$

証明

確率分布は、 $X = {X_1, X_2}$ においてそれぞれ確率 $p,1-p$ であると設定する。この時、期待値は

$\mu = E(X) = \sum(X = X_i)p_i = X_1 p + X_2 (1-p)$

で計算できる。

分散は、 $V(X) = E({( X-\mu) }^2 ) = \sum {( X_i - \mu)}^2p_i = {( X_1 - \mu )}^2p+ {(X_2 - \mu)}^2(1-p)$

と計算できるがここで、 $\mu$ を消去すると、 $X_1 - \mu$ と $X_2 - \mu$ はそれぞれ、

$X_1 - \mu = X_1 - (X_1p + X_2(1-p) ) = (X_1 - X_2)(1-p)$
$X_2 - \mu = X_2 - (X_1p + X_2(1-p) ) = (X_2 - X_1)p$

となる。これを分散の式に代入すると、

$V(X) = {( X_1 - \mu )}^2p+ {(X_2 - \mu)}^2(1-p) \\ = {\bigl( (X_1 - X_2)(1-p) \bigr)}^2 p + { \bigl( (X_2 - X_1)p \bigr) }^2(1-p) \\ = { (X_1 - X_2) }^{2} {(1-p)}^{2}p+{(X_2 - X_1)}^{2}p^{2}(1-p) \\ = p(1-p)\bigl( {(X_1 - X_2)}^{2}(1-p) + {( X_2 - X_1)}^{2} p \bigr) \\ = p(1-p)\bigl( {(X_1 - X_2)}^{2}(1-p) + { \bigl( (-1)(X_1 - X_2)\bigr) }^{2} p \bigr) \\ = p(1-p)\bigl( {(X_1 - X_2)}^{2}(1-p+p) \bigr) \\ = p(1-p){(X_1 - X_2)}^2$

ゆえに、 $V(X) = {(X_1 - X_2)}^{2}p(1-p)$ より

$SD(X) = \sqrt{ V(X) } = \sqrt{ {(X_1 - X_2)}^{2}p(1-p) } = |X_1-X_2|\sqrt{ p(1-p) }$ を得る。

*1:https://rafalab.github.io/dsbook/random-variables.htmlより引用

*2:http://www.stat.ucla.edu/~vlew/stat10/archival/SP01/lectures/sp011013.pdfより引用

2020-07-26

Rで文字列→数値に変換する際、NAs introduced by coercionが出て困った時のtips

次のような変換を考えます。

"1名" → 1
"2名" → 2
"なし" → 0
"調査中" → NA

これをdplyrのパイプの中でmutate( case_when(...) )を駆使して実行していたのですが、エラーとなってしまいました。

データフレーム（見栄えのためtibble）

mydf <- tibble( 同居家族 = c("1名","2名","なし","調査中") )

ダメな例

mydf %>% mutate(
  同居家族数 = case_when( 
    同居家族 == "なし" ~ 0
    , 同居家族 == "調査中" ~ NA_integer_
    , str_detect(同居家族,"名") ~ as.numeric( str_remove(同居家族,"名") )
    )
  )

上を実行すると次のようなErrorが出ます。

Error: Problem with `mutate()` input `同居家族数`.
x must be an integer vector, not a double vector.
ℹ Input `同居家族数` is `case_when(...)`.
Run `rlang::last_error()` to see where the error occurred.
In addition: Warning message:
In eval_tidy(pair$rhs, env = default_env) : NAs introduced by coercion

ググってみると、なんとなく似たようなissueが上がっておりました。

github.com

頑張って読んでみると、case_when で分岐してもas.numeric( str_remove(同居家族,"名") )が実行されるのはstr_detect(同居家族,"名")の時のみではないっぽいことがわかりました。ということで結局、as.numeric()するときに数値に変換できない文字列が入っていたためエラーになっていた、というわけです。

改善策

改善策としては、面倒臭がらず列を数字に変換できる文字列に一旦変換したのちに、改めてas.numeric()してあげれば良いわけです。

mydf %>% 
  mutate(
    同居家族数 = case_when( 
      同居家族 == "なし" ~ "0"
      , 同居家族 == "調査中" ~ NA_character_
      , str_detect(同居家族,"名") ~  str_remove(同居家族, "名")
    ),
    同居家族数 = as.numeric(同居家族数)
  )

実行結果

f:id:axjack:20200726145758p:plain

2020-07-19

edXのHarvardX's Data Scienceを受講しています。

edXのHarvardX's Data Scieceとは？ → HarvardX Data Science Professional Certificate | edX

とりあえずR Basicsが終わってVisualizationのIntroの途中まで来ました。２年前も同じコースを受講したのですが、その時は途中でドロップアウトしてました・・・。モチベーションとやる気があると意外とすんなりできちゃうもんですね。

感想を書いておきます。

読んでおくと良いもの

これです → Introduction to Data Science

R Basicsを受講すると良さそうな人

英語で説明されても頑張れる人
- 字幕もあるしスピードも変えられるので、実はそんなに身構えなくても大丈夫？でした
- Rやdata vizで使う単語を知っていると良いカモです
Rを触ったことがある人
RStudioをインストールしている人
dplyrを触ったことがある人
周りでRを使っている人がいない人
独習者
Rを体系的に学び通したい人

難しいと感じた所

AssignmentのQuestionの英語を読み解くのが難しかったです。

よく使った関数

class()
str()
names()
match()
which()

というかんじです。7月中にVisualizationとProbabilityを終わらせたいです。

2020-07-06

座標軸の回転

勉強数学

これなら分かる最適化数学のp.48にある座標軸の回転について、回転の式を導出してみます。 f:id:axjack:20200706211716p:plain

x成分の計算

f:id:axjack:20200706211844p:plain

$\rm \overline{Oa} = x'cos\theta$ と $\rm \overline{ax} = y'sin\theta$ 、 $\rm \overline{Ox} =\overline{Oa} - \overline{ ax }$ より $\rm x = x'cos\theta - y'sin\theta$ となります。

y成分の計算

f:id:axjack:20200706212253p:plain

$\rm \overline{Ob} = y'cos\theta$ と $\rm \overline{yb} = \overline{Pb'};\overline{Ox'}=\overline{y'P}; \overline{yb} = x'sin\theta$ より $\rm y = x'sin\theta + y'cos\theta$ となります。

まとめると

以上より、 $\rm x = x'cos\theta - y'sin\theta$ 、 $\rm y = x'sin\theta + y'cos\theta$ を行列形式にまとめて、

となります。

まとめ