メモ

• 上の記事でF分布についてメモした
• 少しANOVAを。こちらも参考

n.iter <- 1000
ret <- rep(0,n.iter)
n.g <- 2
n.s <- 200
for(i in 1:n.iter){
	X <- rnorm(n.s)
	P <- sample(1:n.g,n.s,replace=TRUE)
	s1 <- which(P==1)
	s2 <- which(P==2)
	ret[i] <- t.test(X[s1],X[s2])$p.value
}
plot(sort(ret))

ret <- ret.2 <- rep(0,n.iter)
n.g <- 2
n.s <- 200
for(i in 1:n.iter){
	X <- rnorm(n.s)
	P <- sample(1:n.g,n.s,replace=TRUE)
	s1 <- which(P==1)
	s2 <- which(P==2)
	ret[i] <- t.test(X[s1],X[s2])$p.value

	tmp <- lm(X~P)
	ret.2[i] <- anova(tmp)[[5]][1]

}
plot(ret,ret.2)

ret <- ret.2 <- ret.3 <- rep(0,n.iter)
n.g <- 2
n.s <- 200
for(i in 1:n.iter){
	X <- rnorm(n.s)
	P <- sample(1:n.g,n.s,replace=TRUE)
	P <- sort(P)
	s1 <- which(P==1)
	s2 <- which(P==2)
	ret[i] <- t.test(X[s1],X[s2])$p.value

	tmp <- lm(X~P)
	ret.2[i] <- anova(tmp)[[5]][1]

	D <- as.matrix(dist(X))
	D.sq <- D^2
	sum.same <- sum(D.sq[s1,s1])/2/length(s1) + sum(D.sq[s2,s2])/2/length(s2)
	sum.not.same <- sum(D.sq)/2/length(P) - sum.same
	tmp.F <- (sum.same/(n.s-n.g))/(sum.not.same/(n.g-1))
	ret.3[i] <- pf(tmp.F,n.s-n.g,n.g-1)
}
plot(ret.2,ret.3)

ret <- ret.2 <- rep(0,n.iter)
n.g <- 5
n.s <- 200
for(i in 1:n.iter){
	X <- rnorm(n.s)
	P <- sample(1:n.g,n.s,replace=TRUE)
	tmp <- lm(X~as.factor(P))
	ret[i] <- anova(tmp)[[5]][1]
	D <- as.matrix(dist(X))
	D.sq <- D^2
	sum.same <- 0
	for(j in 1:n.g){
		tmp.s <- which(P==j)
		sum.same <- sum.same + sum(D.sq[tmp.s,tmp.s])/2/length(tmp.s)
	}
	sum.not.same <- sum(D.sq)/2/length(P) - sum.same
	tmp.F <- (sum.same/(n.s-n.g))/(sum.not.same/(n.g-1))
	ret.2[i] <- pf(tmp.F,n.s-n.g,n.g-1)

}
par(mfcol=c(1,2))
plot(sort(ret))
plot(ret,ret.2)
par(mfcol=c(1,1))

ret <- ret.2 <- rep(0,n.iter)
n.g <- 2
n.s <- 200
for(i in 1:n.iter){
	X <- matrix(rnorm(n.s*n.g),ncol=n.g)
	ret[i] <- t.test(X[,1],X[,2])$p.value
	ret.2[i] <- t.test(X[,1],X[,2],paired =TRUE)$p.value
	
}
plot(ret,ret.2)

ret <- ret.2 <- rep(0,n.iter)
n.g <- 2
n.s <- 200
for(i in 1:n.iter){
	X <- matrix(0,n.s,n.g)
	X[,1] <- rnorm(n.s)
	X[,2] <- X[,1] + rnorm(n.s)
	ret[i] <- t.test(X[,1],X[,2])$p.value
	ret.2[i] <- t.test(X[,1],X[,2],paired =TRUE)$p.value
	
}
par(mfcol=c(2,3))
boxplot(X)
matplot(t(X),type="l")
plot(sort(ret))
plot(sort(ret.2))
plot(ret,ret.2)
par(mfcol=c(1,1))