二標本コルモゴロフ・スミルノフ検定     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 による度数分布図

fig

> plot(a, method="polygon")                  # 度数多角形による度数分布図

fig

・ 解説ページ


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh