* ICPSR国内利用協議会統計セミナー2009 [#a9204365]
<center><s> <u><figure> <P>Inject0r The FeaR  Was HErE</P></figure></u></s> </center>

** 講義資料 [#u79d96bc]
- [[ダウンロード>http://web.sfc.keio.ac.jp/~maunz/ICPSR2009/ICPSR2009.pdf]](3.3M)

** 演習マニュアル [#s6f2f19f]

 # runif(n,min,max)は区間[min,max]でn個の一様分布に従うランダムな数を生成
 # rnorm(n, m, s)は平均m,標準偏差sの正規分布に従うランダムな数を生成
 x1<-runif(n1,0,20); y1<-rnorm(n1, 0.8*x1+1.2, sqrt(v))
 x2<-runif(n2,7,20); y2<-rnorm(n2, 1.2*x2+0.9, sqrt(v))
 x3<-runif(n3,10,20); y3<-rnorm(n3, 0.7*x3+1.0, sqrt(v))
 x4<-runif(n4,9,20); y4<-rnorm(n4, 1.0*x4+0.8, sqrt(v))
 x <- c(x1,x2,x3,x4); y <- c(y1,y2,y3,y4)

# 図1.3と図1.4
par(mfrow=c(1,2))
plot(x1, y1,ann=FALSE,cex=1.4,cex.axis=1.5,xlim=c(0,22),ylim=c(0,30))
title(xlab="X",cex.lab=1.2,font.lab=4)
title(ylab="Y",cex.lab=1.2,font.lab=4)
points(x2,y2,pch=22,cex=1.4)
points(x3,y3,pch=23,cex=1.4)
points(x4,y4,pch=24,cex=1.4)
curve(0.8*x+1.2,add=T,lwd=2,lty=1,xlim=c(0,21))
curve(1.2*x+0.9,add=T,lwd=2,lty=2,xlim=c(6,21))
curve(0.7*x+1.0,add=T,lwd=2,lty=3,xlim=c(9,21))
curve(1.0*x+0.8,add=T,lwd=2,lty=4,xlim=c(8,21))
legend(0,30,legend=c("モデル1:Y=0.8*X+1.2","モデル2:Y=1.2*X+0.9","モデル3:Y=0.7*X+1.0","モデル4:Y=1.0*X+0.8"),
lwd=2,lty=c(1,2,3,4),cex=1.2)

plot(x1, y1,ann=FALSE,cex=1.4,cex.axis=1.5,xlim=c(0,22),ylim=c(0,30))
title(xlab="X",cex.lab=1.2,font.lab=4)
title(ylab="Y",cex.lab=1.2,font.lab=4)
points(x2,y2,pch=22,cex=1.4)
points(x3,y3,pch=23,cex=1.4)
points(x4,y4,pch=24,cex=1.4)
curve(0.82358*x+2.17414,add=T,lwd=2,lty=1,xlim=c(0,21))

# 番外編
# x<-c(x1,x2,x3,x4)
# y<-c(y1,y2,y3,y4)
# summary(lm(y~x))
# library(MCMCpack)
# sample.gb.post1 <- MCMCregress(y~x,burnin=1000,mcmc=10000,b0=0,B0=0.001,c0=0.001,d0=0.001)
# summary(sample.gb.post1)
# hist(y,main=" ",xlab="Y",ylab="頻度",font.lab=1)

#### 1.3節 ####

# 図1.6
par(mfrow=c(1,3))
# x~N(0,20)
curve(dnorm(x,0,20), -100, 100, main="x~N(0,20)", ylab="確率分布", lwd=2, cex.axis=1.5, 
cex.lab=1.5, cex.main=1.5,font.main=4)
# x~N(10,20)
curve(dnorm(x,10,20), -100,100, main="x~N(10,20)", ylab="確率分布", lwd=2, cex.axis=1.5, 
cex.lab=1.5, cex.main=1.5,font.main=4)
# x~N(0,40)
curve(dnorm(x,0,40), -100, 100, main="x~N(0,40)", ylab="確率分布", lwd=2, cex.axis=1.5, 
cex.lab=1.5, cex.main=1.5, font.main=4)

# y=1のとき
# dnorm(x, mean,sd)
curve(dnorm(x, 11/20, 11/20), -10, 10)

# 図1.7, 1.8, 1.9
par(mfrow=c(1,3))
beta <- 0.9;
b <- sum(y*x)/sum(x*x)
# seq(min, max, length)は区間[min,max]でlength個の数を生成
beta1 <- seq(0.62, 1.03, length=100)
# dnorm(n, m, s)は平均m,標準偏差sの正規分布
plot(beta1, dnorm(beta1, b, 1/sqrt(sum(x*x))),type="l", xlab=expression(beta), 
ylab="尤度",xlim=c(0.9,1.03),lwd=3,cex=1.3,cex.axis=1.5,cex.lab=1.3)
plot(x, y, xlab="X", ylab="Y",lwd=2,cex=1.3,cex.axis=1.5,cex.lab=1.4,font.lab=4)
curve(0.865*x, add=T,lwd=3)
plot(dunif,-1,2,type="l",xlab=" ",ylab=" ",lwd=3,cex.axis=1.5,cex.lab=1.4,font.lab=4)

# 図1.10, 1.11
par(mfrow=c(1,2))
# n=40回の試行について期待値θ=0.25の二項分布の確率密度関数を描く
# dbinom(x, n, θ)は区間x内,確率θ,サンプル数nの二項分布の確率密度関数
n<-40; theta<-0.25;
plot(dbinom(1:40,n,theta),type="l",lwd=3, xlab="x",ylab="確率密度",cex.lab=1.3)	

# 試行n回,期待値θの二項分布
# y <- rbinom(n,1,theta); s <- sum(y)
s <- sum(rbinom(n,1, theta))
plot(seq(0,1,length=40), dbeta(seq(0,1,length=40), s+1, n-s+1), type="l",lwd=3, 
xlab=expression(theta),ylab="Likelihood",cex.lab=1.5, font.lab=4)

# 図1.12
par(mfrow=c(1, 2))
x <- seq(0,1,length=100)
# dbeta(x, a, b)
# B(0.5,0.5)
plot(x,dbeta(x,0.5,0.5),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.2,cex.lab=1.2,font.lab=4)
#legend(0.3,2.5,legend=expression(Beta(0.5,0.5)),box.lwd=0)
# B(10,12)
plot(x,dbeta(x,10,12),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.2,cex.lab=1.2,font.lab=4)

# Β関数の重ね描き
x <- seq(0,1,length=100)
matplot(x,dbeta(x,10,12),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.3,cex.lab=1.3,font.lab=4)
text(0.62,3.5,cex=1.3,expression(Beta(10,12)))
matplot(x,dbeta(x,0.5,0.5),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.3,cex.lab=1.3,font.lab=4,add=T)
text(0.15,3.0,cex=1.3,expression(Beta(0.5,0.5)))
matplot(x,dbeta(x,5,5),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.3,cex.lab=1.3,font.lab=4,add=T)
text(0.46,1.8,cex=1.3,expression(Beta(5,5)))
matplot(x,dbeta(x,8,2),type="l", xlab="x",ylab=expression(Beta(a,b)),lwd=2,cex.axis=1.3,cex.lab=1.3,font.lab=4,add=T)
text(0.69,2.8,cex=1.3,expression(Beta(8,2)))

# Γ(x,a)
plot(dgamma(1:10,1),type="l", xlab="x",ylab=expression(Gamma(a)),lwd=2,cex.axis=1.5,cex.lab=1.5,font.lab=4)
plot(dgamma(1:10,10),type="l")
plot(dgamma(1:10,40),type="l")
plot(dgamma(1:10,70),type="l")

# 図1.13
matplot(dgamma(1:100,1),type="l", xlab="x",ylab=expression(Gamma(a)),lwd=2,cex.axis=1.3,cex.lab=1.3,font.lab=4)
text(10,0.3,cex=1.3,expression(Gamma(1)))
matplot(dgamma(1:100,10),type="l", xlab="x",ylab=expression(Gamma(a)),lwd=2,cex.axis=1.5,cex.lab=1.5,font.lab=4,add=T)
text(21,0.1,cex=1.3,expression(Gamma(10)))
matplot(dgamma(1:100,40),type="l", xlab="x",ylab=expression(Gamma(a)),lwd=2,cex.axis=1.5,cex.lab=1.5,font.lab=4,add=T)
text(50,0.07,cex=1.3,expression(Gamma(40)))
matplot(dgamma(1:100,70),type="l", xlab="x",ylab=expression(Gamma(a)),lwd=2,cex.axis=1.5,cex.lab=1.5,font.lab=4,add=T)
text(83,0.05,cex=1.3,expression(Gamma(70)))


トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS