#### 第1章 ####
#### 1.1節 ####

# ゴルフの番手，練習回数，平均，分散
PW<-rnorm(150,100,15); 
AW<-rnorm(80,80,13); 

# 度数分布表
# ランダムな正規分布に従ってデータ作成しているため、計算する毎に、表示される度数はテキストの数字と異なります。
break1<-c(0,85,95,300)
PW.tab1<-table(cut(PW,break1)); AW.tab1<-table(cut(AW,break1))
PW.tab1; AW.tab1
APW.tab1<-rbind(PW.tab1,AW.tab1)

# 表1.1
APW.tab1
APW.prob1<-APW.tab1/sum(APW.tab1)
APW.prob1

# 条件付き確率
PW.prob1<-PW.tab1/150
AW.prob1<-AW.tab1/80
PW.prob1; AW.prob1

# 図1.1
curve(dnorm(x,100,15),30,170,ylim=c(0,.06),xlab="飛距離",ylab="確率",lwd=2,cex=1.2, cex.axis=1.2, cex.lab=1.2, 
cex.main=1.2)
rect(80,-0.1,100,0.1,col="grey",border=NA,xlim=c(30,70),ylim=c(0,.06))
curve(dnorm(x,100,15),30,170,ylim=c(0,.06),xlab="飛距離",ylab="確率",lwd=2,cex=1.2, cex.axis=1.2, cex.lab=1.2, 
cex.main=1.2,add=T)
curve(dnorm(x,80,15),30,170,add=T,lwd=2,lty=2)
legend(110, 0.05,legend=c("PW","AW"),lwd=2,lty=c(1,2),cex=1.2)

# 図1.2
breaks<- seq(0,300,10)
PW.tab2<-table(cut(PW,breaks))
AW.tab2<-table(cut(AW,breaks))
barplot(PW.tab2)

# 式1.4
27/150*83/230+25/150*42/230+98/150*106/230

# 13P3
prod(1:13)/prod(1:10)
# 13C3
prod(1:13)/(prod(1:3)*prod(1:10))

#### 1.2節 ####

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

# 被説明変数もランダムなデータを作成
# 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))

# 図1.3
# 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("モデル１：Y=0.8*X+1.2","モデル２：Y=1.2*X+0.9","モデル３：Y=0.7*X+1.0","モデル４：Y=1.0*X+0.8"),
lwd=2,lty=c(1,2,3,4),cex=1.2)

# 図1.4
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)
x <- c(x1,x2,x3,x4)
y <- c(y1,y2,y3,y4)
xy.lm <- lm(y~x)
curve(xy.lm$coefficients[2]*x+xy.lm$coefficients[1],add=T,lwd=2,lty=1,xlim=c(0,21))

# 番外編
# library(MCMCpack)
# sample.gb.post1 <- MCMCregress(y0~x0,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
# 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.9, 1.1, 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.1),lwd=3,cex=1.3,cex.axis=1.5,cex.lab=1.3)
# 図1.8
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)
# 図1.9
plot(dunif,-1,2,type="l",xlab=" ",ylab=" ",lwd=3,cex.axis=1.5,cex.lab=1.4,font.lab=4)

# 図1.10
# 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)	

# 図1.11
# 試行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)
# 図1.12 (a)
# 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)
# 図1.12 (b)
# 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)))

