##  01  ##
xy <- sapply(seq(10000),function(x){(rnorm(1,68,13)-rnorm(1,41,21))})
	#「2つの異なる正規分布に従う群から1つずつ数値を取り出して差を求める」試行を独立に10000回繰り返す
	#rnorm;正規分布から乱数を抽出
hist(xy)
	#横軸に差の値、縦軸に生起回数のヒストグラム
plot(density(xy))
	#確率密度を推定してグラフに表示

##  02  ##
m <- sapply(seq(100),function(x){mean(runif(10000))})	#runif;一様分布から乱数を抽出
m
	#n=10000の標本を100回抽出してそれぞれの標本平均を求める

mean(m)
var(m)
	#その標本平均の平均値と分散

##  03  ##
sd(m)
(sqrt(1/12))/(sqrt(na))
	#標本平均の分布の標準偏差と標本分布から算出される標準誤差
	#sqrt;平方根


##  04  ##

na <- 10000
a <- sapply(seq(1),function(x){runif(na)})
a
	#n=10000の標本を1回抽出する

mean(a)
sd(a)
	#{a}の標本平均と標準偏差

aa <- (a-mean(a))
(sum(aa^2))/(na)+(1/12*1/na)
(sum(aa^2))/(na-1)
	#母分散の推定値
1/12
	#母集団分散の理論値

sda <- sqrt((sum(aa^2))/(na-1))
sda/sqrt(na)
	#標準誤差の推定値

##  05  ##

za <- aa/sd(a)
za
mean(za)
sd(za)
	#{a}を標準化し、平均値、SDを計算

zm <- scale(m)	#scale;標準化
par(mfrow=c(1,2))	#グラフエリアを1行2列に分割
plot(density(m))
plot(density(zm))
	#{m}を使って標準化した分布の比較

##  06  ##

integrate(dnorm, 1.65, Inf)	#integrate;積分, Inf;無限大
integrate(dnorm, -Inf, -1.65)
	#標準正規分布で片側確率を算出


##  07  ##
###あとで訂正###
b <- sapply(seq(1),function(x){(rnorm(100, 0, 1.0))})
b
mean(b)
c <- sapply(seq(1),function(x){(rnorm(100, -0.6, 1.0))})
c
mean(c)
	#bは区間(0, 1.0)より抽出、cは区間(0.6, 1.0)より抽出

zb <- (mean(b)-0.5)/sqrt((1/12)/100)
zb
zc <- (mean(c)-0.5)/sqrt((1/12)/100)
zc
	#それぞれの群を標準化

azb <- abs(zb)
azc <- abs(zc)
integrate(dnorm, azb, Inf)
integrate(dnorm, azc, Inf)
	#積分してz>=|z0|となる確率密度を計算

##  08  ##

k <- sapply(seq(100000),function(x){(((rnorm(1,5,2)-5))/2)^2})
par(mfrow=c(1,2))
hist(k)
plot(density(k))
	#自由度1のカイ2乗分布

k2 <- sapply(seq(100000),function(x){(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,15,2)-15))/2)^2})
par(mfrow=c(1,2))
hist(k2)
plot(density(k2))
	#自由度2のカイ2乗分布

k5 <- sapply(seq(100000),function(x){(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2
})
par(mfrow=c(1,2))
hist(k5)
plot(density(k5), xlim=c(0,30), ylim=c(0,1.0))
	#自由度5のカイ2乗分布

k10 <- sapply(seq(100000),function(x){(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2+
(((rnorm(1,5,2)-5))/2)^2
})
par(mfrow=c(1,2))
hist(k10)
plot(density(k10), xlim=c(0,30), ylim=c(0,1.0))
	#自由度10のカイ2乗分布


par(mfrow=c(1,1))
plot(density(k))
lines(density(k2), col=2)
lines(density(k5), col=3)
lines(density(k10), col=4)
	#まとめてグラフ内で比較

##  09  ##

tx1 <-sapply(seq(100),function(x){rnorm(2)})
t1 <- list()
for(i in 1:100){
t1[[i]] <- (mean(tx1[,i])-0)/sqrt(sd(tx1[,i])^2/2)
}
t1.t <- unlist(t1)
plot(density(t1.t))
hist(t1.t)
	#自由度1のt分布

tx5 <-sapply(seq(100),function(x){rnorm(6)})
t5 <- list()
for(i in 1:100){
t5[[i]] <- (mean(tx5[,i])-0)/sqrt(sd(tx5[,i])^2/6)
}
t5.t <- unlist(t5)
plot(density(t5.t))
hist(t5.t)
	#自由度5のt分布

tx100 <-sapply(seq(100),function(x){rnorm(101)})
t100 <- list()
for(i in 1:100){
t100[[i]] <- (mean(tx100[,i])-0)/sqrt(sd(tx100[,i])^2/101)
}
t100.t <- unlist(t100)
plot(density(t100.t))
hist(t100.t)
	#自由度100のt分布

tx10000 <-sapply(seq(100),function(x){rnorm(10001)})
t10000 <- list()
for(i in 1:100){
t10000[[i]] <- (mean(tx10000[,i])-0)/sqrt(sd(tx10000[,i])^2/10001)
}
t10000.t <- unlist(t10000)
plot(density(t10000.t))
hist(t10000.t)
	#自由度10000のt分布

plot(density(t10000.t))
lines(density(t1.t), col=2)
lines(density(t5.t), col=3)
lines(density(t100.t), col=4)
	#まとめてグラフで比較

##  10  ##
kukan1x <-sapply(seq(100),function(x){rnorm(2)})
	#n=2の標本を100回抽出して95%信頼区間を求める
kukan1l <- list()
kukan1h <- list()
for (i in 1:100){
kukan1l[[i]] <- ((mean(kukan1x[,i])-abs(qt(0.025, df=2)))*(sd(kukan1x[,i])/sqrt(2)))
}&{
kukan1h[[i]] <- ((mean(kukan1x[,i])+abs(qt(0.025, df=2)))*(sd(kukan1x[,i])/sqrt(2)))
}	#qt(a, df=f);自由度fのt分布における100*a%点の理論値
kukan1 <- cbind(low=kukan1l, high=kukan1h)	#信頼下限をlow,信頼上限をhighとする

kukan1t <- list()
for (i in 1:100){
if (kukan1[i,1] < 0 & 0 < kukan1[i,2]){
kukan1t[[i]] <- "true"
}else{
kukan1t[[i]] <- "false"
}
}	#区間内に母平均(0)が入っていればtrue,入っていなければfalseを返す
kukan1t <- unlist(kukan1t)
kukan1t


kukan100x <-sapply(seq(100),function(x){rnorm(101)})
	#n=101の標本を100回抽出して95%信頼区間を求める
kukan100l <- list()
kukan100h <- list()
for (i in 1:100){
kukan100l[[i]] <- ((mean(kukan100x[,i])-abs(qt(0.025, df=2)))*(sd(kukan100x[,i])/sqrt(101)))
}&{
kukan100h[[i]] <- ((mean(kukan100x[,i])+abs(qt(0.025, df=2)))*(sd(kukan100x[,i])/sqrt(101)))
}
kukan100 <- cbind(low=kukan100l, high=kukan100h)

kukan100t <- list()
for (i in 1:100){
if (kukan100[i,1] < 0 & 0 < kukan100[i,2]){
kukan100t[[i]] <- "true"
}else{
kukan100t[[i]] <- "false"
}
}
kukan100t <- unlist(kukan100t)
kukan100t



##  作図用  ##

plot(dunif,-1,2)

n <- sapply(seq(1),function(x){runif(20)})
mean(n)
nn <- (n-mean(n))
nnn <- sum(nn^2)
nnn/19
var(n)


par(mfrow=c(2,2))
hist(sapply(seq(5),function(x){mean(runif(5))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="5回の場合")
hist(sapply(seq(10),function(x){mean(runif(5))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="10回の場合")
hist(sapply(seq(50),function(x){mean(runif(5))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="50回の場合")
hist(sapply(seq(100),function(x){mean(runif(5))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="100回の場合")

par(mfrow=c(2,2))
hist(sapply(seq(100),function(x){mean(runif(3))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="n=3の場合")
hist(sapply(seq(100),function(x){mean(runif(5))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="n=5の場合")
hist(sapply(seq(100),function(x){mean(runif(10))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="n=10の場合")
hist(sapply(seq(100),function(x){mean(runif(100))}),xlim=c(0,1),
xlab="meanx", ylab="出現回数", main="n=100の場合")