二標本コルモゴロフ・スミルノフ検定 Last modified: Aug 20, 2009
目的
二標本コルモゴロフ・スミルノフ検定を行う
注: R には ks.test という関数が用意されている(原データを用いる)
使用法
ks2(obs1, obs2)
plot.ks2(obj, method=c("barplot", "polygon"), name=c(obj$name1, obj$name2), col=c("grey30", "grey90"),
pch=c(19, 2), col2=c(1,4), lty=c(1,2), xlab="", ylab="パーセント", x="topright", y=NULL, ...)
引数
obs1 第一群の度数分布
obs2 第二群の度数分布
obj ks2 関数が返すオブジェクト
method=c("barplot", "polygon") barplot か度数多角形か
name=c(obj$name1, obj$name2) 凡例に表示する各群の名前
col=c("grey30", "grey90") barplot の色
pch=c(19, 2) 度数多角形のマーク
col2=c(1,4) 度数多角形の線の色
lty=c(1,2) 度数多角形の線の種類
xlab="", ylab="パーセント" 軸の名称
x="topright", y=NULL legend の位置
... barplot, matplot への引数
ソース
インストールは,以下の 1 行をコピーし,R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/ks2.R", encoding="euc-jp")
# 二標本コルモゴロフ・スミルノフ検定を行う
ks2 <- function( obs1, # 第一群の度数分布
obs2) # 第二群の度数分布
{
stopifnot(length(obs1) == length(obs2)) # 要素数は同じでなくてはならない
name1 <- deparse(substitute(obs1))
name2 <- deparse(substitute(obs2))
data.name <- paste(name1, "and", name2)
method <- "二標本コルモゴロフ・スミルノフ検定"
n1 <- sum(obs1) # 第一群のサンプルサイズ
n2 <- sum(obs2) # 第二群のサンプルサイズ
cum1 <- cumsum(obs1)/n1 # 第一群の累積相対度数
cum2 <- cumsum(obs2)/n2 # 第二群の累積相対度数
D <- max(abs(cum1-cum2)) # 差の絶対値の最大値
chi <- 4*D^2*n1*n2/(n1+n2) # 検定統計量
df <- 2
p <- min(1, pchisq(chi, df, lower.tail=FALSE)*2) # P 値
stat <- c(D=D, "X-squared"=chi)
names(df) <- "df"
return(structure(list(statistic=stat, parameter=df, p.value=p,
method=method, data.name=data.name,
obs1=obs1, obs2=obs2, name1=name1, name2=name2),
class=c("htest", "ks2"))) # 結果をまとめて返す
}
# plot メソッド(群別の度数分布図を描く)
plot.ks2 <- function( obj, # ks2 関数が返すオブジェクト
method=c("barplot", "polygon"), # barplot か度数多角形か
name=c(obj$name1, obj$name2), # 各群の名前
col=c("grey30", "grey90"), # barplot の色
pch=c(19, 2), # 度数多角形のマーク
col2=c(1,4), # 度数多角形の線の色
lty=c(1,2), # 度数多角形の線の種類
xlab="", ylab="パーセント", # 軸の名称
x="topright", y=NULL, # legend の位置
...) # barplot, matplot への引数
{
d <- rbind(obj$obs1, obj$obs2)
d <- d/rowSums(d)*100
nc <- ncol(d)
if (match.arg(method) == "barplot") {
barplot(d, beside=TRUE, col=col, ylab=ylab, ...)
legend(x, y, legend=name, fill=col)
axis(1, at=1:nc*3-1, labels=1:nc, pos=0)
}
else {
matplot(t(d), type="l", col=col2, lty=lty, xlab=xlab, ylab=ylab, xaxt="n", bty="n", ...)
points(rep(1:nc, each=2), d, col=rep(col2, nc), pch=rep(pch, nc))
legend(x, y, legend=name, col=col2, lty=lty, pch=pch)
axis(1, at=1:nc, labels=1:nc, pos=0)
}
}
使用例
> obs1 <- c(1, 2, 1, 3, 2, 3, 3, 2, 7, 11, 10, 9, 13, 13, 22, 17, 23, 20, 17, 14, 13, 5, 2, 1, 1, 1)
> obs2 <- c(0, 1, 2, 2, 4, 5, 6, 8, 10, 7, 16, 17, 17, 13, 19, 13, 18, 10, 4, 6, 6, 5, 1, 3, 0, 0)
> (a <- ks2(obs1, obs2))
二標本コルモゴロフ・スミルノフ検定
data: obs1 and obs2
D = 0.1892, X-squared = 14.5969, df = 2, p-value = 0.001353
> plot(a) # barplot による度数分布図
> plot(a, method="polygon") # 度数多角形による度数分布図
解説ページ
直前のページへ戻る
E-mail to Shigenobu AOKI