.　

「Rで学ぶ統計解析」掲載プログラム
　.　
　　3章

　　4章

　　5章

　　6章

　　7章

　　8章

　　9章

　　10章

　　もどる

「Rで学ぶ統計解析」に掲載した「プログラム例」「実行例」のRのプログラムを掲載しました。コピーペーストでRのエディタに取り込むことが出来ます。


3章

# プログラム例１：標本調査と集計
population <- c(rep("支持",40000), rep("不支持",60000))
z <- sapply(rep(1000,1000), function(x) {
		ss = sample(population, x)
		sum(ss == "支持")})
hist(z)
c(min(z),max(z))

# プログラム例２：標本調査と集計（比率）
z <- sapply(rep(1000,1000), function(x) {
		ss = sample(population, x)
		mean(ss == "支持")})
mean(abs(z-0.4) < 0.02)

# プログラム例３：標本調査と集計（集団的性質）
z <- sapply(rep(100,2000), function(x) {
		ss = sample(population, x)
		mean(ss == "支持")})
c(mean(z), sd(z))
hist(z, breaks=20, freq=F)
curve(dnorm(x,mean(z),sd(z)), add=T)

# プログラム例４：標本調査（正規母集団）
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))


# プログラム例５：標本調査、調査対象数の違い
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章

# プログラム例：2変量正規分布の密度関数
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")

# プログラム例：ドモアブルラプラスの定理
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章

# プログラム例：カイ2乗分布に関連する標本実験
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("自由度",n-1))
curve(dchisq(x,n-1), add=T)
curve(dchisq(x,n), add=T, col=2)

# プログラム例：スチューデントのｔ分布に関連する標本実験
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("自由度：",n-1))
	curve(dt(x,n-1), za, zb, add=T)
	curve(dnorm(x), za, zb, add=T, col="red")
}
tSampling(5, 10000)

# プログラム例：F分布に関連する標本実験
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章

# プログラム例：区間推定の信頼性
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章

# 実行例：ｔ検定の例
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章

# プログラム例
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)

# 実行例
zz <- rnorm(20, 170, 10)
ww <- rnorm(20, 160, 5)
var.test(zz, ww)

# プログラム例
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 プロット")
qqline(data)
pd <- pnorm(data, mean(data), sd(data))
qqplot((c(1:20)-0.5)/20, pd, main="P-P プロット")
abline(0,1)

# プログラム例
normalPlot <- function() {
	plot(c(-3,3.5), c(-3,3), type="n", axes=F, 
		xlab="", ylab="", main="正規確率紙")
	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()

# プログラム例
#### データ
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 プロット
qd <- qnorm((c(1:n)-0.5)/n, mean(ht), sd(ht))
plot(sort(ht), qd, main="Q-Q プロット")
abline(0,1)

# プログラム例
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)

# プログラム例
z <- c(45,15,25,15)
chisq.test(c(45,15,25,15), p=c(0.38,0.22,0.31,0.09))

# プログラム例
m <- 20; n <- length(rr) ## rr には収益率データが記憶されている
class <- qnorm(c(0:m)/m) ## 全体をm 分割
dosu <- table(cut((rr-mean(rr))/sd(rr), class))
chisq.test(dosu, p=rep(1/m,m))

# 実行例
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)

# プログラム例
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章

# プログラム例
accident <- read.table("clipboard",header=T,row.names=1)
attach(accident)
names(accident)
plot(accident)

# 実行例
lm(事故件数~保有車両数)
abline(lm(事故件数~保有車両数))

# プログラム例
ahat <- model$coefficients[1] ## 切片の推定値
bhat <- model$coefficients[2] ## 回帰係数の推定値
誤差<- 事故件数- (ahat + bhat * 保有車両数)
sgm <- sqrt(var(誤差)*(n-1)/(n-2)) ## 不偏分散の平方根
for(i in 1:100) {
	y <- sapply(1:n, 
		function(m) rnorm(1,ahat+bhat*保有車両数[m],sgm))
	abline(lm(y~保有車両数), lty=3) ## 回帰直線を引く
}

# プログラム例
aa <- c(); bb <- c();
for(i in 1:1000) {
	y <- sapply(1:n, 
		function(m) rnorm(1,ahat+bhat*保有車両数[m],sgm)) 
				## 観測データ
	md <- lm(y~保有車両数) ## 回帰モデルの計算
	aa <- c(aa, md$coefficients[1]) ## 切片の推定値を記録
	bb <- c(bb, md$coefficients[2]) ## 回帰係数の推定値を記録
}
par(mfrow=c(2,1))
hist(aa, breaks=20, freq=F, main="切片") 
				## 切片の推定値のヒストグラム
curve(dnorm(x,ahat,sgm * sqrt((xv+xm^2) / (n*xv))), add=T, col=2) 
				## 正規分布曲線
hist(bb, breaks=20, freq=F, main="回帰係数") 
				## 回帰係数の推定値のヒストグラム
curve(dnorm(x,bhat,sgm * sqrt(1 / (n*xv))), add=T, col=2) 
				## 正規分布曲線

# プログラム例
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章

# プログラム例
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)

# プログラム例
level <- factor(level)
summary(aov(data ~ level))

# プログラム例
data <- c(za, zb, zc)
level <- factor(c(rep("A",8), rep("B",8), rep("C",8)))
summary(aov(data ~ level))

# 実行例
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))

# プログラム例
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))