2020-06-16

一元配置分散分析を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

2020-05-17

母比率の信頼区間に含まれる２次不等式を解く

統計調べ物 R

母比率の信頼区間
２次不等式を解く
- 式変形
具体例で検証
- 公式を用いる
- ２次不等式を解いた結果を用いる
参考

母比率の信頼区間

母比率 $p$ の母集団からサイズ $n$ の標本を抽出する。このとき標本割合 $\hat{p} = x / n$ について、 $\hat{p}$ は近似的に平均 $p$ 、分散 $pq /n$ の正規分布 $N(p, pq /n )$ に従う。ただし、 $q = 1 - p$ である。したがって、 $z = \frac{ \hat{p} - p } { \sqrt{pq /n } }$ と標準化した $z$ は標準正規分布 $N(0,1)$ に従う。

さて、確率 $\rm{P}( -z_0 < z \lt z_0 ) = 1 - \alpha$ から不等式 $-z_0 < z < z_0$ を抜き出したものに $z = \frac{ \hat{p} - p } { \sqrt{pq /n } }$ の右辺を代入すると、 $p$ の $100\alpha/2$ %信頼区間は

$\displaystyle \hat{p} - z_0 \sqrt{ \frac{pq}{n} } \lt p \lt \hat{p} + z_0 \sqrt{ \frac{pq}{n} }$

と表される。例えば $\alpha = 0.10$ とすれば、 $\hat{p} - 1.64 \sqrt{ \frac{pq}{n} } \lt p \lt \hat{p} + 1.64 \sqrt{ \frac{pq}{n} }$ である。通常ここで $p$ は未知であり不等式の左右の $\sqrt{ \, \, }$ にある $p$ はnが大なるとき $\hat{p}$ の一致性から $\frac{pq}{n} = \frac{ \hat{p}(1-\hat{p} ) }{n}$ と置き換えて、

$\displaystyle \hat{p} - z_0 \sqrt{ \frac{ \hat{p}(1-\hat{p} )}{n} } \lt p \lt \hat{p} + z_0 \sqrt{ \frac{\hat{p}(1-\hat{p}) }{n} }$

という公式が用いられる。

２次不等式を解く

それでは、 $\frac{pq}{n} = \frac{ \hat{p}(1-\hat{p} ) }{n}$ を用いず
$\hat{p} - z_0 \sqrt{ \frac{pq}{n} } \lt p \lt \hat{p} + z_0 \sqrt{ \frac{pq}{n} }$
を $p$ について解いてみる。以下、式変形。

式変形

$\hat{p} - z_0 \sqrt{ \frac{pq}{n} } \lt p \lt \hat{p} + z_0 \sqrt{ \frac{pq}{n} }$
$\iff - z_0 \sqrt{ \frac{pq}{n} } \lt p - \hat{p} \lt + z_0 \sqrt{ \frac{pq}{n} }$
$\iff | p - \hat{p} | \lt z_0 \sqrt{ \frac{pq}{n} }$
$\iff | p - \hat{p} |^2 \lt \Bigr( z_0 \sqrt{ \frac{pq}{n} } \Bigl)^2$
$\iff p^2 - 2p\hat{p} + (\hat{p})^2 - \Bigr( z_0 \sqrt{ \frac{pq}{n} } \Bigl)^2 < 0$
$\iff p^2 - 2p\hat{p} + (\hat{p})^2 - (z_0)^2 \frac{pq}{n} < 0$
$\iff np^2 - 2np\hat{p} + n(\hat{p})^2 - (z_0)^2 pq < 0$
$\iff np^2 - 2np\hat{p} + n(\hat{p})^2 - (z_0)^2 p(1-p) < 0$
$\iff np^2 - 2np\hat{p} + n(\hat{p})^2 - (z_0)^2 p + (z_0)^2p^2 < 0$
$\iff \Bigl(n + (z_0)^2\Bigr) p^2 + \Bigl( -2n\hat{p} - (z_0)^2 \Bigr)p + n(\hat{p})^2 < 0$
$\iff p^2 + \Bigl( \frac{ -2n\hat{p} - (z_0)^2 } { n + (z_0)^2 } \Bigr)p + \frac{ n(\hat{p})^2 }{{ n + (z_0)^2 } } < 0$
$\iff \Bigl( p + \frac{1}{2} ( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} ) \Bigr)^2 - \Bigl( \frac{1}{2} ( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} ) \Bigr)^2 + \frac{ n(\hat{p})^2 }{{ n + (z_0)^2 } } < 0$
$\iff \Bigl( p + \frac{1}{2} ( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} ) \Bigr)^2 < \frac{1}{4} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr)^2 - \frac{ n(\hat{p})^2 }{{ n + (z_0)^2 } }$
$\iff p + \frac{1}{2} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr) < \pm \sqrt{ \frac{1}{4} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr)^2 - \frac{ n(\hat{p})^2 }{{ n + (z_0)^2 } } }$
$\iff p< - \frac{1}{2} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr) \pm \sqrt{ \frac{1}{4} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr)^2 - \frac{ n(\hat{p})^2 }{{ n + (z_0)^2 } } }$

となる。

具体例で検証

統計学基礎のp.118の例９の数字を使って検証してみる。例９の主要な数字は、

標本サイズ n = 1200
標本比率 $\hat{p} = 0.054$ にて
母比率 $p$ の95%信頼区間を求める。

である。ここで、 $z_0 = 1.96$ とする。

公式を用いる

$\hat{p} - z_0 \sqrt{ \frac{ \hat{p}(1-\hat{p}) }{n} } \lt p \lt \hat{p} + z_0 \sqrt{ \frac{ \hat{p}(1-\hat{p}) }{n} }$
なので、 $0.054 \pm 1.96 \sqrt{ \frac{0.054(1-0.054)}{1200} } = [0.04121184, 0.06678816$ ]
となる。

２次不等式を解いた結果を用いる

$p< - \frac{1}{2} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr) \pm \sqrt{ \frac{1}{4} \Bigl( \frac{-2n\hat{p}-(z_0)^2}{ n + (z_0)^2} \Bigr)^2 - \frac{ n(\hat{p})^2 }{{ n + (z_0)^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

より、 $[0.04257642,0.06827005$ ]
となる。

参考

統計学の基礎(9.26)

2020-05-05

２次元正規分布の確率密度関数

勉強統計

東京大学出版会『統計学入門』の p.145 図7.5 ２次元正規分布について、与えられたパラメータから楕円群(等高線の式)を導出する。

パラメータ

$\displaystyle \mu = \left[ \begin{array} {c} 55.24 \\ 34.97 \end{array} \right]$ , $\displaystyle \Sigma = \left[ \begin{array} {cc} 210.54 && 126.99\\126.99 && 119.68 \end{array} \right]$

二次元正規分布の確率密度関数

$\displaystyle f( \rm x) = \frac{1}{\sqrt{(2\pi)^k \rm det(\Sigma) }} \rm exp \left( \frac{1}{2} ( \rm x - \rm \mu )'\Sigma^{-1} (\rm x - \rm \mu) \right)$

代入

$\det (\Sigma)と\Sigma^{-1}$ を計算すると、

$\det \Sigma = 9070.9671$

$\sqrt{(2\pi)^k \rm det(\Sigma) } = 238.6748277$

$\Sigma^{-1} = \left[ \begin{array} {cc} 0.013193742 && -0.013999609 \\-0.013999609 && 0.023210314 \end{array} \right]$

となるので、確率密度関数に代入するとexpの中身は $\rm x = \left( \begin{array}{c} x \\ y \end{array} \right)$ として

$\displaystyle - \frac{1}{2} \left( 0.013193742(x - 55.24)^2 + 2(-0.013999609)(x-55.24)(y-34.97) + 0.023210134(y-34.97)^2 \right)$

となる。ここでexpの中身を適当な定数c と置き、-200倍すると、

$\displaystyle 1.3193742(x - 55.24)^2 -2.7999218(x-55.24)(y-34.97) + 2.3210134(y-34.97)^2$

となり、p.145にある楕円群の式とほぼ等しくなる。

2020-05-01

AB = E ならば BA = Eをランクを用いて示す

線形代数

チャート式シリーズ　大学教養　線形代数を買ってから、行列の構造、特に「ランク」の理解が深まりました。ということでタイトルの命題の証明です。

AB = E ならば BA = E

AもBもn次正方行列、Eはn次単位行列とします。
ここで、rankA = r ≦ n , rankB = s ≦ n とし、rankE = nです。

まず、基本行列Pを左からAに掛けた結果をPA = X とすると、
rankPA = rankA = r = rankX です。また、rankP = nです。

次に、PA = Xを用い、P(AB) = PE = P および (PA)B = XB を得、総じて XB = P を得ます。

さて、XB = P のランクについて考えると、右辺：rankP = nなので rankXB = nを得ます。
rankXB = n ということは、rankX = nかつrankB = n ということが分かります。

ここで、XはAをPを用いて階段化した行列であり、
なおかつrankX = nよりXは正則であることから、
X = Eを得ます。したがって、XB = P ⇔ EB = P ⇔ B = P です。

以上より、仮定はAB = Eだったので
AB = E ⇔ AP = E ⇔ P = A‪⁻¹‬ですので、
BA = PA = A‪⁻¹‬A = E となります。■

参考

線形代数Ｉ/行列の階数 - 武内＠筑波大

2020-04-29

行簡約行列をRで

R 線形代数

pracmaのrrefを使って行簡約行列を出してみます。これで適当な行列を手計算で簡約化して答え合わせができますね。

ソース

# install.packages('pracma')
# library('pracma')

## 行ベクトルを4本
a1 <- c(1,2,0)
a2 <- c(2,4,0)
a3 <- c(0,1,3)
a4 <- c(1,3,3)

## 行列にする
M <- matrix(c(a1,a2,a3,a4),nrow = 3, ncol=4)

## 行列を表示
print(M)

## 行簡約行列を表示 
rref(M)

## 行列の階数を表示
Rank(M)

実行結果

> a1 <- c(1,2,0)
> a2 <- c(2,4,0)
> a3 <- c(0,1,3)
> a4 <- c(1,3,3)
> M <- matrix(c(a1,a2,a3,a4),nrow = 3, ncol=4)
> print(M)
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    1
[2,]    2    4    1    3
[3,]    0    0    3    3
> Rank(M)
[1] 2
> rref(M)
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    1
[2,]    0    0    1    1
[3,]    0    0    0    0
>

2020-04-17

update文

SQL

とある１レコードの値で本番のデータを全件更新してしまうというヤラカシをやってしまったので、次回は気をつける意を込めた記録です。 SQLはcompile sql server onlineを使用します。

仕様

tbl1.val1をtbl2.val2へと更新する
更新の条件
- tbl1.idとtbl2.idを結合し、tbl2.val2 = 'foo'の時だけ更新する
- tbl2.val2 <> 'foo'の時は更新しない

tbl1

f:id:axjack:20200417235435p:plain

tbl2

f:id:axjack:20200417235458p:plain

ということでこの場合、id = {10,30,80}が更新対象となります。

問題のコード(ヤラカシあり)

-- データ作成 ここから
create table tbl1  (
 id int primary key
  , val1 nvarchar(10) null
)
;

insert into tbl1 values
  (10,'a')
  ,(20,'b')
  ,(30,'c')
  ,(40,'d')
  ,(50,'e')
  ,(60,'f')
  ,(70,'g')
  ,(80,'h')
  ,(90,'i')
;

create table tbl2 (
 id int primary key
  , val2 nvarchar(10) null
  , constraint fk_id foreign key(id) references  tbl1
)
;

insert into tbl2 values
  (10,'foo')
  ,(30,'foo')
  ,(50,'baz')
  ,(60,'baz')
  ,(80,'foo')
;
-- データ作成 ここまで



--  一時テーブルに tbl1 のレコードを格納
select * into #t1 from tbl1
;


-- tbl1 更新前
select * from tbl1
;

-- 更新処理
update tbl1
set val1 = tbl2.val2
from #t1 
  inner join tbl2 
    on #t1.id = tbl2.id 
    and tbl2.val2 = 'foo'
;


-- tbl1 更新後
select * from tbl1
;

実行結果

上：更新前のtbl1, 下：更新後のtbl1なのですが、見ての通り残念な感じに更新されてしまいました。。

f:id:axjack:20200417234647p:plain

better なコード

★の箇所のように、updateで更新対象となるテーブルをinner join してあげるとinner join で抽出される結果通りに更新されます。

update tbl1
set val1 = tbl2.val2
from #t1 
  inner join tbl2 
    on #t1.id = tbl2.id 
    and tbl2.val2 = 'foo'
  inner join tbl1 on tbl1.id = #t1.id -- ★
;

better なコードの実行結果

仕様を満たした更新結果となりました。 f:id:axjack:20200417235104p:plain

まなび

更新対象のテーブルをちゃんとjoin しましょう。

2020-04-05

カーネルは線型空間

線形代数

示したいこと

線型空間V, W に対して線型写像 $\displaystyle T:V \rightarrow W$ が与えられているとする。
この時、カーネル $\displaystyle \rm Ker \ T = \{ x \in V | T(x) = 0_w \in W \}$ は線型空間であることを示す。

証明

$\displaystyle \forall x, y \in \rm Ker \ T$ および $\displaystyle \forall a, b \in \mathbb{R}$ をそれぞれ取る。
すると、 $\displaystyle ax+by \in V$ より、 $\displaystyle T(ax+by) = T(ax) + T(by) = aT(x) + bT(y) = a0_w + b0_w = 0_w + 0_w = 0_w \in W$
となるので、Ker T の任意の元はベクトル和とスカラー倍に関して閉じていることが分かった。よってKer T は線型空間であることが示せた。