axjack's blog

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

一元配置分散分析を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
axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.