せいびし

専門性を持たない自分の正体を求めて

1標本の検定(ノンパラメトリック検定)のプログラム

以下の記事で作成したプログラムです。
グラフの描画のソースは省略します。
seibibibi.hatenadiary.jp

# 混合正規分布からデータのサンプリング、検定統計量の取得
sampling <- function(pi, mus, sig, N){
	# 一様乱数をサンプル
	u <- runif(N, 0, 1)
	data <- numeric(N)
	dens <- numeric(N)
	for (i in 1:N){
		for (j in 1:length(pi)){
			if (u[i] < pi[j]){
				break
			}
		}
		data[i] <- rnorm(1, mus[j], sig[j])
	}
	return(data)
}


# 0. サンプリングに関するパラメータ
pi <- c(0.2, 0.8)			# 混合比(中央値は4.3)
#pi <- c(0.5, 0.5)			# 混合比(中央値は0)
mus <- c(-5, 5)
sig <- c(sqrt(5), sqrt(5))
n = 100	# 混合正規分布からの乱数の数
N = 10000	# シミュレーション回数

# 1. データをサンプルして平均、分散、検定統計量を取得
U <- data.frame()
for (i in 1:N){
	data <- sampling(pi, mus, sig, n)
	# 検定統計量の計算
	# wilcox.testは帰無仮説はΔ=0,符号検定の帰無仮説はΔ(中央値)=0
	w_n <- wilcox.test(data)$statistic
	s_n <- SIGN.test(data, md=0)$statistic	# BSDAパッケージのダウンロードが必要
	U <- rbind(U, c(w_n, s_n))
}
colnames(U) <- c("w_n","s_n")
# 標準化
U2 <- U
exp_w_n <- n * (n+1) / 4
var_w_n <- n * (n+1) * (2*n+1) / 24
exp_s_n <- n / 2
var_s_n <- n / 4
U2$w_n <- (U$w_n-exp_w_n) / sqrt(var_w_n)
U2$s_n <- (U$s_n-exp_s_n) / sqrt(var_s_n)

# 描画
graph(U)
graph2(U2)