# メモ

• 上の記事で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))
```