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$ を行列形式にまとめて、

となります。

2020-06-16

一元配置分散分析をRで実装する

R プログラミング統計

手計算でもaovでも一元配置分散分析は出来るので、Rで実装してみようと思った次第です。F分布の累積分布を除いてほぼほぼベクトル演算？を使っています。データは水準の繰り返し数のトータルから適当に生成しています。統計学的観点はほぼ０です。

実装

# データ設定 ####
## 水準
catv  <- c( rep("A",5),rep("B",5),rep("C",10),rep("D",3))

## データの範囲
data.range <- 10:30

## データ
datav <- sample(data.range,length(catv),T)


# 計算 #### 
## 総平均
mean.all <- mean(datav)

## 水準ごとの平均
mean.cat <- tapply(datav, catv, mean)

## 水準平方和
catv.sumsq <- sum( table(catv)  *  (mean.cat - mean.all)^2 )

## 自由度: 水準の種類数 - 1
catv.df    <- length(table(catv))-1
  
## 全平方和
total.sumsq <- sum((datav - mean.all)^2)

## 残差平方和
### 残差平方和 := 全平方和 - 水準平方和
zansa.sumsq <- total.sumsq - catv.sumsq

## 自由度: 水準ごとの繰り返し数-1 をそれぞれ足す
zansa.df    <- sum(table(catv)-1)

### 別解: 残差平方和 := 要素 - 水準平均 の平方和
#abs( zansa.sumsq - sum( (datav - rep(mean.cat,table(catv)))^2 ) ) < 0.01


# 答え合わせ ####
# aov
summary( aov(datav ~ catv) )

# 自力で計算
cat("\t   ","Df","Sum Sq"," Mean Sq","F value","Pr(>F)","\n"
    ," catv","     ",catv.df,"",round(catv.sumsq,2), round(catv.sumsq/catv.df,2)
    ,"   ",round((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),3)," ",round(1-pf((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),catv.df,zansa.df ),3),"\n"
    ," Residuals",zansa.df,"",round(zansa.sumsq,2), "",round(zansa.sumsq/zansa.df,2)
    )

実行例

> # 答え合わせ ####
> # aov
> summary( aov(datav ~ catv) )
            Df Sum Sq Mean Sq F value Pr(>F)
catv         3  170.9   56.96   2.349  0.105
Residuals   19  460.8   24.25               
> 
> # 自力で計算
> cat("\t   ","Df","Sum Sq"," Mean Sq","F value","Pr(>F)","\n"
+     ," catv","     ",catv.df,"",round(catv.sumsq,2), round(catv.sumsq/catv.df,2)
+     ,"   ",round((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),3)," ",round(1-pf((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),catv.df,zansa.df ),3),"\n"
+     ," Residuals",zansa.df,"",round(zansa.sumsq,2), "",round(zansa.sumsq/zansa.df,2)
+     )
        Df Sum Sq  Mean Sq F value Pr(>F) 
  catv       3  170.89 56.96     2.349   0.105 
  Residuals 19  460.77  24.25

学び

table関数の出力結果はベクトルのようである。

確信はないのだけれどtable関数の出力はベクトルっぽく扱えます。

例

> table( c(rep("A",2),rep("B",4)) )
A B 
2 4 

> table( c(rep("A",2),rep("B",4)) ) +3 
A B 
5 7

rep関数の繰り返し回数はベクトルを渡すことができる

水準ごとの平均を水準内の繰り返し数分だけ繰り返す、といった処理をしたいときにrepとtableを組み合わせて対処することができます。

例

10を1回、20を2回、30を3回、40を4回繰り返す。

> rep(c(10,20,30,40),times = c(1,2,3,4))
 [1] 10 20 20 30 30 30 40 40 40 40

参考

2020-06-10

Rで線形計画法

lpSolve*1を使って線形計画法を解いてみた。理論は難しいから敬遠していたがソルバー使うと一瞬で解けて気持ち良い。関数lpに渡すパラメータは以下の通り。第一パラメータは"max"か"min"。詳しいことはrdocumentationなどを参照。

f_obj : ベクトル。目的関数の係数を格納する。
f_con : 行列。制約条件の係数を格納する。
f_dir : ベクトル。不等号または等号を格納する。
f_rhs : ベクトル。right hand side(右辺)、制約条件の右辺の係数(値)を格納する。

以下実例。

install and load library

install.packages('lpSolve')
library(lpSolve)

その１*2

線型計画法

f_obj <- c(1,1) #目的関数の係数
f_con <- matrix(c(1,2,2,1),ncol=2) #制約条件を行列で
f_dir <- c("<=","<=") #不等号?
f_rhs <- c(4,4) # 不等号の右側

最適解を表示

lp("max",f_obj,f_con,f_dir,f_rhs)

最適解を取る時の引数

lp("max",f_obj,f_con,f_dir,f_rhs)$solution

整数計画法

lp("max",f_obj,f_con,f_dir,f_rhs,int.vec=1:2) #int.vecはintegerなvectorのindicesを与えていると思われる。
lp("max",f_obj,f_con,f_dir,f_rhs, int.vec=1:2)$solution

0-1 整数計画法

f_obj01 <- c(16,19,23,11)
f_con01 <- matrix( c(2,3,4,5), ncol=4, nrow=1)
f_dir01 <- c("<=")
f_rhs01 <- c("6")
lp("max",f_obj01,f_con01,f_dir01,f_rhs01,all.bin=T)
lp("max",f_obj01,f_con01,f_dir01,f_rhs01,all.bin=T)$sol

その２*3

f_obj3 <- c(4,6)
f_con3 <- matrix(c(1,3,0,2,0,4),nrow=3,ncol=2)
f_dir3 <- rep("<=",3)
f_rhs3 <- c(8,12,12)
f_sol3 <- lp("max", f_obj3,f_con3,f_dir3,f_rhs3)
f_sol3
f_sol3$sol

f_obj4 <- c(20,30)
f_con4 <- matrix(c(1,3,3,2,4,1),nrow=3,ncol=2)
f_dir4 <- rep("<=",3)
f_rhs4 <- c(800,1800,1500)
f_sol4 <- lp("max",f_obj4,f_con4,"<=",f_rhs4)
f_sol4
f_sol4$sol

その３*4

f_obj5 <- c(1,1)
f_con5 <- matrix(c(3,1,1,3),nrow=2,ncol=2)
f_dir5 <- rep("<=",2)
f_rhs5 <- c(9,6)
f_sol5 <- lp("max",f_obj5,f_con5,f_dir5,f_rhs5)
f_sol5
f_sol5$solution

*1:see https://www.rdocumentation.org/packages/lpSolve/versions/5.6.15

*2:see http://s-ir.sap.hokkyodai.ac.jp/dspace/bitstream/123456789/8049/1/67-1-sizen-01.pdf

*3:see https://www.fit.ac.jp/~song/lecture/06/LP0606.pdf

*4:see https://juken-mikata.net/how-to/mathematics/linear-programming.html

コード

結果

参考