せいびし

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

2標本の検定・尺度母数の推定(ノンパラメトリック検定)

以下の記事で作成したプログラムです。グラフの部分は省略しました。
seibibibi.hatenadiary.jp

# seibisi
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.3, 0.7)
mus1 <- c(10, 20)
mus2 <- c(10, 20)
sig1 <- c(sqrt(5), sqrt(5))
sig2 <- c(sqrt(10), sqrt(10))
m <- 100
n <- 50
N <- 10000
# 1. データをサンプルして検定統計量の計算
MO1 <- data.frame()
MO2 <- data.frame()
F <- data.frame()
for (i in 1:N){
	data1 <- sampling(pi, mus1, sig1, m)
	data2 <- sampling(pi, mus2, sig2, n)
	MO1 <- rbind(MO1, mood.test(data1, data2)$statistic)
	MO2 <- rbind(MO2, mood.test(data2, data1)$statistic)
	F <- rbind(F, var.test(data1, data2)$statistic)
}
MO <- data.frame(MO1, MO2)
colnames(MO) <- c("MO1","MO2")
colnames(F) <- "F"
data <- data.frame(data1, data2)
# 2. 棄却域に落ちた個数
p_u <- qf(df1=100-1, df2=50-1, p=0.05, lower.tail=FALSE)
F_num <- length(F[F$F > p_u,])
p_l <- qnorm(p=0.025, mean=0, sd=1, lower.tail=TRUE)
p_u <- qnorm(p=0.025, mean=0, sd=1, lower.tail=FALSE)
MO_num <- length(MO$MO1[MO$MO1 < p_l]) + length(MO$MO1[MO$MO1 > p_u])
F_num
MO_num