標準誤差は統計量の標準偏差であることを確かめる
やること
離散型確率分布
X | p |
---|---|
17 | 0.3 |
-1 | 0.7 |
に於いて、大きさn=100の標本を抽出する時、
- 和S の標本分布について、Sの平均、標準誤差
- 平均M の標本分布について、Mの平均、標準誤差
をそれぞれ
- 計算
- モンテカルロ シミュレーション
によりそれぞれ求める。
計算:離散型確率分布から分かること
Xの期待値・分散・標準偏差はそれぞれ、
であるから、サンプルサイズn=100における、
和S = の標本分布について
- 各は として、
- Sの期待値:=
- Sの標準誤差:=
平均M = の標本分布について
- 各は として、
- Mの期待値:=
- Mの標準誤差:=
である。
モンテカルロ シミュレーション
共通の設定
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関数の練習
(これを最初に書くべきだったのでは?)サイコロを2回投げて出た目の和が7となる確率、をRで求める。
共通
saikoro <- 1:6 # サイコロ B <- 10000 #サンプル"数"
組み合わせで割り算する素朴な求め方
# サイコロを2回投げた時の組み合わせ saikorokumi <- expand.grid(saikoro,saikoro) # 和が7になる組み合わせ ÷ サイコロを2回投げた時の組み合わせの総数 sum( apply(saikorokumi,1,sum) == 7 ) / nrow(saikorokumi)
> sum( apply(saikorokumi,1,sum) == 7 ) / nrow(saikorokumi) [1] 0.1666667
2回投げて出た目の和の分布
# サイコロを2回投げて出た目の和の分布 kiroku <- replicate(B,{ sum(sample(saikoro,2,replace=TRUE)) }) # 和が7となる確率 mean(kiroku == 7)
> mean(kiroku == 7) [1] 0.1659
2回投げて出た目の和が7となるか否かの分布
# サイコロを2回投げて出た目の和が7となるか否かの分布 kiroku2 <- replicate(B,{ sum(sample(saikoro,2,replace=TRUE)) == 7 }) # 和が7となる確率 mean(kiroku2)
> mean(kiroku2) [1] 0.1684
確率変数の取りうる値が2つの時の離散型確率分布の標準偏差を求める公式
結論
確率変数においてそれぞれ確率をとる時、標準偏差SD(X)は
となる。ここで||は絶対値の記号である。
どこで見つけたか
edXのprobabilityの教科書?であるChapter 14 Random variables | Introduction to Data Science や ここである。
計算例
離散型確率分布が以下の表、
X | p |
---|---|
17 | 0.3 |
-1 | 0.7 |
で与えられているとする。
まずは期待値・分散・標準偏差を計算し、最後に公式を用いて標準偏差を計算して2つの標準偏差が一致することを確かめる。
期待値
分散
標準偏差
分散のルートを取って求めると、
一方、冒頭の公式で求めると、
証明
確率分布は、においてそれぞれ確率であると設定する。この時、期待値は
で計算できる。
分散は、
と計算できるがここで、を消去すると、とはそれぞれ、
となる。これを分散の式に代入すると、
ゆえに、より
を得る。
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が上がっておりました。
頑張って読んでみると、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(同居家族数) )
実行結果
edXのHarvardX's Data Scienceを受講しています。
edXのHarvardX's Data Scieceとは? → HarvardX Data Science Professional Certificate | edX
とりあえずR Basicsが終わってVisualizationのIntroの途中まで来ました。2年前も同じコースを受講したのですが、その時は途中でドロップアウトしてました・・・。モチベーションとやる気があると意外とすんなりできちゃうもんですね。
感想を書いておきます。
読んでおくと良いもの
これです → Introduction to Data Science
R Basicsを受講すると良さそうな人
- 英語で説明されても頑張れる人
- 字幕もあるしスピードも変えられるので、実はそんなに身構えなくても大丈夫?でした
- Rやdata vizで使う単語を知っていると良いカモです
- Rを触ったことがある人
- RStudioをインストールしている人
- dplyrを触ったことがある人
- 周りでRを使っている人がいない人
- 独習者
- Rを体系的に学び通したい人
難しいと感じた所
- AssignmentのQuestionの英語を読み解くのが難しかったです。
よく使った関数
- class()
- str()
- names()
- match()
- which()
というかんじです。7月中にVisualizationとProbabilityを終わらせたいです。
座標軸の回転
これなら分かる最適化数学のp.48にある座標軸の回転について、回転の式を導出してみます。
x成分の計算
と、 よりとなります。
y成分の計算
とよりとなります。
まとめると
以上より、、を行列形式にまとめて、
となります。
一元配置分散分析をRで実装する
手計算でもaov
でも一元配置分散分析は出来るので、Rで実装してみようと思った次第です。F分布の累積分布を除いてほぼほぼベクトル演算?を使っています。データは水準の繰り返し数のトータルから適当に生成しています。統計学的観点はほぼ0です。
実装
# データ設定 #### ## 水準 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
参考
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)
その1*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
その2*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
その3*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
母比率の信頼区間に含まれる2次不等式を解く
母比率の信頼区間
母比率 の母集団からサイズの標本を抽出する。このとき標本割合 について、 は近似的に平均、分散 の正規分布 に従う。ただし、 である。したがって、 と標準化した は 標準正規分布 に従う。
さて、確率 から不等式 を抜き出したものに の右辺を代入すると、の%信頼区間は
と表される。例えば とすれば、である。通常ここでは未知であり不等式の左右のにあるはnが大なるときの一致性からと置き換えて、
という公式が用いられる。
2次不等式を解く
それでは、を用いず
をについて解いてみる。以下、式変形。
式変形
となる。
具体例で検証
統計学基礎のp.118の例9の数字を使って検証してみる。例9の主要な数字は、
標本サイズ n = 1200
標本比率 にて
母比率 の95%信頼区間を求める。
である。ここで、 とする。
公式を用いる
なので、 ]
となる。
2次不等式を解いた結果を用いる
なので、Rを用いて計算すると
n <- 1200
ph <- 0.054
z0 <- 1.96
A <- (-2*(n*ph) - (z0)^2)/(n+(z0)^2)
B <- (n*ph^2)/(n+z0^2)c(-(1/2)*A - sqrt( (1/4)*A ^2 - B),-(1/2)*A + sqrt( (1/4)*A ^2 - B))
> c(-(1/2)*A - sqrt( (1/4)*A ^2 - B),-(1/2)*A + sqrt( (1/4)*A ^2 - B))
[1] 0.04257642 0.06827005
より、 ]
となる。