ICPSR国内利用協議会統計セミナー2009

講義資料

演習マニュアル

# 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