�@.�@

�uR�Ŋw�ԓ��v��́v�f�ڃv���O����

�@.�@
�@�@3��

�@�@4��

�@�@5��

�@�@6��

�@�@7��

�@�@8��

�@�@9��

�@�@10��


�uR�Ŋw�ԓ��v��́v�Ɍf�ڂ����u�v���O������v�u���s��v��R�̃v���O�������f�ڂ��܂����B�R�s�[�y�[�X�g��R�̃G�f�B�^�Ɏ�荞�ނ��Ƃ��o���܂��B


3��

# �v���O������P�F�W�{�����ƏW�v
population <- c(rep("�x��",40000), rep("�s�x��",60000))
z <- sapply(rep(1000,1000), function(x) {
		ss = sample(population, x)
		sum(ss == "�x��")})
hist(z)
c(min(z),max(z))

# �v���O������Q�F�W�{�����ƏW�v�i�䗦�j
z <- sapply(rep(1000,1000), function(x) {
		ss = sample(population, x)
		mean(ss == "�x��")})
mean(abs(z-0.4) < 0.02)

# �v���O������R�F�W�{�����ƏW�v�i�W�c�I�����j
z <- sapply(rep(100,2000), function(x) {
		ss = sample(population, x)
		mean(ss == "�x��")})
c(mean(z), sd(z))
hist(z, breaks=20, freq=F)
curve(dnorm(x,mean(z),sd(z)), add=T)

# �v���O������S�F�W�{�����i���K��W�c�j
n <- 10000; m <- 510; s <- 5;
w <- rnorm(n, m, s);
hist(w, breaks=20, freq=F)
curve(dnorm(x,m,s), add=T)
c(mean(w),sd(w))


# �v���O������T�F�W�{�����A�����Ώې��̈Ⴂ
par(mfrow=c(2,1))
n <- 20; m <- 510; s <- 5; N <- 1000;
ww <- sapply(rep(n,N), function(x) mean(rnorm(x, m, s)))
hist(ww, freq=F)
curve(dnorm(x,mean(ww),sd(ww)), add=T)
n <- 100
ww <- sapply(rep(n,N), function(x) mean(rnorm(x, m, s)))
hist(ww, freq=F)
curve(dnorm(x,mean(ww),sd(ww)), add=T)

4��

# �v���O������F2�ϗʐ��K���z�̖��x�֐�
x <- y <- seq(-3, 3, length=60);
rho <- 0.8
norm3D <- function(x, y)
exp(-(x^2-2*rho*x*y+y^2)/2/(1-rho^2))/2/pi/sqrt(1-rho^2)
z <- outer(x, y, norm3D)
persp(x, y, z, theta=30, phi=30, expand=0.5, ticktype="detailed")

# �v���O������F�h���A�u�����v���X�̒藝
n <- 50; p <- 0.3;
plot(0:n, dbinom(0:n, n, p), type="h")
curve(dnorm(x,n*p,sqrt(n*p*(1-p))), add=T)

5��

# �v���O������F�J�C2�敪�z�Ɋ֘A����W�{����
m <- 10; sd <- 2; n <- 5; N <- 1000
z <- sapply(1:N, function(x){var(rnorm(n,m,sd))*(n-1)/sd/sd})
hist(z, freq=F, main=paste("���R�x",n-1))
curve(dchisq(x,n-1), add=T)
curve(dchisq(x,n), add=T, col=2)

# �v���O������F�X�`���[�f���g�̂����z�Ɋ֘A����W�{����
tSampling <- function(n, N) {
	z <- sapply(1:N, 
		function(x){x=rnorm(n); mean(x)/(sd(x)/sqrt(n))})
	za <- floor(min(z)); zb = ceiling(max(z))
	hist(z, breaks=za:zb, freq=F, main=paste("���R�x�F",n-1))
	curve(dt(x,n-1), za, zb, add=T)
	curve(dnorm(x), za, zb, add=T, col="red")
}
tSampling(5, 10000)

# �v���O������FF���z�Ɋ֘A����W�{����
m <- 50; n <- 40;
w <- sapply(1:1000, function(x) {
za <- rnorm(m); zb <- rnorm(n,100,1); var(za)/var(zb)})
hist(w, freq=F,main=paste("m,n=",m,n))
curve(df(x,m-1,n-1),add=T)

6��

# �v���O������F��Ԑ���̐M����
coverage <- function(n = 10, a = 0.05, m = 1000) {
	tn <- qt(1-a/2,n-1) / sqrt(n)
	sampling <- function(n) {y = rnorm(n);
		ifelse(abs(mean(y)) < tn*sd(y), 1, 0)
	}
	mean(sapply(rep(n,m), sampling)
}
sapply(rep(n,10), coverage)
sapply(rep(n,10), coverage, 0.01, 10000)

7��

# ���s��F������̗�
z <- c(507.6, 519.2, 518.9,08.7, 518, 512.2, 514.8, 511.2, 509.6, 512)
t.test(z, nu=S10, alternative="greater")

8��

# �v���O������
normapr <- function(sd1=1, n1=10, sd2=2, n2=20) {
	x1 <- rnorm(n1, 0, sd1)
	x2 <- rnorm(n2, 0, sd2)
	return((mean(x1)-mean(x2)) / sqrt(var(x1)/n1 + var(x2)/n2))
}
z <- sapply(1:10000, normapr)
hist(z, ylim=c(0,0.4), freq=F)
curve(dnorm(x), add=T)

# ���s��
zz <- rnorm(20, 170, 10)
ww <- rnorm(20, 160, 5)
var.test(zz, ww)

# �v���O������
data <- c(3.5, 4.9, 4.7, 4.2, 3.9, 6.3, 5.1, 3.5, 5.3, 3.8, 
	3.9, 5.7, 7.2, 6.3, 4.4, 4.4, 4.1, 3.7, 5.0, 3.9)
qqnorm(data, main="Q-Q �v���b�g")
qqline(data)
pd <- pnorm(data, mean(data), sd(data))
qqplot((c(1:20)-0.5)/20, pd, main="P-P �v���b�g")
abline(0,1)

# �v���O������
normalPlot <- function() {
	plot(c(-3,3.5), c(-3,3), type="n", axes=F, 
		xlab="", ylab="", main="���K�m����")
	xa <- c(0.01,0.05,0.1*c(1:9),0.95,0.99)
	axis(2,qnorm(xa), pnorm(qnorm(xa)), label=F)
	xb <- c(0.01, 0.1, 0.5, 0.9, 0.99)
	axis(2,qnorm(xb),pnorm(qnorm(xb)))
	sapply(xa, function(x) 
		lines(c(-3,3),c(qnorm(x),qnorm(x)), lty=3))
	sapply(-3:3, function(x) lines(c(x,x), c(-2.5,2.5), lty=2))
	sapply(c(-2,-1,1,2), function(x) lines(c(-3,3),c(x,x),lty=2))
	lines(c(-3,3), c(0,0))
	text(3,2, expression(paste(mu,+2,sigma)), pos=4)
	text(3,1, expression(mu+sigma), pos=4)
	text(3,0, expression(mu), pos=4)
	text(3,-1, expression(mu-sigma), pos=4)
	text(3,-2, expression(paste(mu,-2,sigma)), pos=4)
}
normalPlot()

# �v���O������
#### �f�[�^
ht <- c(166.2,169.2,170.3,169.5,179.6,162.7,175.3,173.1,166.8,173.9,
	187.3,176.4,171.4,176.3,163.8,168.3,164.8,175.2,168.4,157.7)
n <- length(ht)
#### Q-Q �v���b�g
qd <- qnorm((c(1:n)-0.5)/n, mean(ht), sd(ht))
plot(sort(ht), qd, main="Q-Q �v���b�g")
abline(0,1)

# �v���O������
ks. test (ht, "pnornrt , mean(ht) , sd(ht) )
ks.test(ht, "pnorm", nean(ht), sd(ht), alternati.ve=,,g,')
ks.test(ht, "p[orn", nean(ht), sd(ht), alternative=rI,r)

# �v���O������
z <- c(45,15,25,15)
chisq.test(c(45,15,25,15), p=c(0.38,0.22,0.31,0.09))

# �v���O������
m <- 20; n <- length(rr) ## rr �ɂ͎��v���f�[�^���L������Ă���
class <- qnorm(c(0:m)/m) ## �S�̂�m ����
dosu <- table(cut((rr-mean(rr))/sd(rr), class))
chisq.test(dosu, p=rep(1/m,m))

# ���s��
count <- c(22,63,49,35,21,10)
m <- sum(count*(0:5))/sum(count)
chisq.test(count,p=c(dPois(0:4,m),ppois(4,m,lower.tail=F)))
pchisq(3.4311, 4, lower.tail=F)
qchisq(0.05, 4, lower.tail=F)

# �v���O������
n <- 20; r <- 0.5; N <- 1000;
rr <- sapply(rep(n,N), function(n) {x <- rnorm(n);
cor(x,r * x + sqrt(1-r^2) * rnorm(n))})
hist(rr, freq=F)

9��

# �v���O������
accident <- read.table("clipboard",header=T,row.names=1)
attach(accident)
names(accident)
plot(accident)

# ���s��
lm(���̌���~�ۗL�ԗ���)
abline(lm(���̌���~�ۗL�ԗ���))

# �v���O������
ahat <- model$coefficients[1] ## �ؕЂ̐���l
bhat <- model$coefficients[2] ## ��A�W���̐���l
�덷<- ���̌���- (ahat + bhat * �ۗL�ԗ���)
sgm <- sqrt(var(�덷)*(n-1)/(n-2)) ## �s�Ε��U�̕�����
for(i in 1:100) {
	y <- sapply(1:n, 
		function(m) rnorm(1,ahat+bhat*�ۗL�ԗ���[m],sgm))
	abline(lm(y~�ۗL�ԗ���), lty=3) ## ��A����������
}

# �v���O������
aa <- c(); bb <- c();
for(i in 1:1000) {
	y <- sapply(1:n, 
		function(m) rnorm(1,ahat+bhat*�ۗL�ԗ���[m],sgm)) 
				## �ϑ��f�[�^
	md <- lm(y~�ۗL�ԗ���) ## ��A���f���̌v�Z
	aa <- c(aa, md$coefficients[1]) ## �ؕЂ̐���l���L�^
	bb <- c(bb, md$coefficients[2]) ## ��A�W���̐���l���L�^
}
par(mfrow=c(2,1))
hist(aa, breaks=20, freq=F, main="�ؕ�") 
				## �ؕЂ̐���l�̃q�X�g�O����
curve(dnorm(x,ahat,sgm * sqrt((xv+xm^2) / (n*xv))), add=T, col=2) 
				## ���K���z�Ȑ�
hist(bb, breaks=20, freq=F, main="��A�W��") 
				## ��A�W���̐���l�̃q�X�g�O����
curve(dnorm(x,bhat,sgm * sqrt(1 / (n*xv))), add=T, col=2) 
				## ���K���z�Ȑ�

# �v���O������
civ <- function(x, ap) qt(1-ap/2,n-2)*sgm*sqrt((xv+(x-xm)^2)/(n*xv))
curve(ahat+bhat*x + civ(x,0.2), add=T, col=2)
curve(ahat+bhat*x - civ(x,0.2), add=T, col=2)
curve(ahat+bhat*x + civ(x,0.05), add=T, col=2)
curve(ahat+bhat*x - civ(x,0.05), add=T, col=2)

10��

# �v���O������
za <- c(99.8,104.6,101.7,100.0,104.0,101.8,101.0,105.0)
zb <- c(99.2,104.7,104.0, 101.8, 100.5, 101.8, 101.0, 103.1)
zc <- c(104.7, 101.6, 100.9, 101.5, 102.7, 101.4, 102.6, 99.7)
zd <- c(106.2, 101.2, 105.8, 104.1, 106.6, 107.9, 103.6, 105.9)
data <- c(za, zb, zc, zd)
level <- c(rep(1,8), rep(2,8), rep(3,8), rep(4,8))
plot(level,data,xaxp=c(1,4,3),xlab="��ƈ�")
points(1:4,c(mean(za),mean(zb),mean(zc),mean(zd)),pch=16)
lines(1:4,c(mean(za),mean(zb),mean(zc),mean(zd)))
abline(h=mean(data), lty=2)
legend(2,108,c("�•ʕ���","������"),lty=c(1,2),cex=0.8)

# �v���O������
level <- factor(level)
summary(aov(data ~ level))

# �v���O������
data <- c(za, zb, zc)
level <- factor(c(rep("A",8), rep("B",8), rep("C",8)))
summary(aov(data ~ level))

# ���s��
Z <- c(99.8,104.6,101.7,99.2,104.7,104.0,100.7,107.6,100.9,100.2,
				110.2,105.8)
X <- factor(c(rep("1",3),rep("2",3),rep("3",3),rep("4",3)))
Y <- factor(rep(c("1m","2a","3e"),4))
summary(aov(Z ~ X + Y))

# �v���O������
bartlett = function(s1=1, s2=1, s3=1.1, n1=10, n2=10, n3=10) {
	m1 <- 0; z1 <- rnorm(n1,m1,s1); v1 <- var(z1)
	m2 <- 1; z2 <- rnorm(n2,m2,s2); v2 <- var(z2)
	m3 <- 2; z3 <- rnorm(n3,m3,s3); v3 <- var(z3)
	fe <- n1 + n2 + n3 - 3
	ve <- (v1*(n1-1) + v2*(n2-1) + v3*(n3-1)) / fe
	c <- 1 + (1/(n1-1)+1/(n2-1)+1/(n3-1) - 1/fe)/3/2
	(fe*log(ve) - (n1-1)*log(v1)-(n2-1)*log(v2)-(n3-1)*log(v3))/c
}
nn <- 5000; zz <- numeric(nn)
for(i in 1:nn) zz[i] <- bartlett(1,1,1.1,10,10,10)
mean(zz > qchisq(0.05,2,lower.tail=F)) # ����5% �𒴂���䗦
hist(pchisq(zz,2,lower.tail=F))
dosu <- table(cut(pchisq(zz,2,lower.tail=F),c(0:10)/10))
chisq.test(dosu, p=rep(0.1,10))