《R语言与统计分析》 11章代码

#例题 11.3.1

library(dplyr)
library(reshape2)
library(ggplot2)

x <- seq(0, 1, 0.01)
n <- 5
y <- 3
alpha <- c(0.5, 0.5, 0.5, 1, 1, 1, 1.5, 1.5, 1.5)
beta <- c(0.5, 1, 1.5, 0.5, 1, 1.5, 0.5, 1, 1.5)
z.post <- as.data.frame(matrix(0, nrow = length(x), ncol = length(alpha)))
z.prior <- as.data.frame(matrix(0, nrow = length(x), ncol = length(alpha)))
for (i in c(1:9)) {
  z.post[, i] <- dbeta(x, alpha[i] + y, n - y + beta[i])
  z.prior[, i] <- dbeta(x, alpha[i], beta[i])
}

z.post <- select(z.post, post.V1 = V1, post.V2 = V2, post.V3 = V3, post.V4 = V4,
                 post.V5 = V5, post.V6 =V6, post.V7 = V7, post.V8 = V8, post.V9 = V9)
z.prior <- select(z.prior, prior.V1 = V1, prior.V2 = V2, prior.V3 = V3,
                  prior.V4 = V4, prior.V5 = V5, prior.V6 = V6, prior.V7 = V7, prior.V8 = V8,
                  prior.V9 = V9)

z.wide <- cbind(x,z.post,z.prior)
z.long1 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V1","prior.V1"),variable.name = "bayes",value.name = "value")%>%mutate(group="V1")
z.long2 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V2","prior.V2"),variable.name = "bayes",value.name = "value")%>%mutate(group="V2")
z.long3 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V3","prior.V3"),variable.name = "bayes",value.name = "value")%>%mutate(group="V3")
z.long4 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V4","prior.V4"),variable.name = "bayes",value.name = "value")%>%mutate(group="V4")
z.long5 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V5","prior.V5"),variable.name = "bayes",value.name = "value")%>%mutate(group="V5")
z.long6 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V6","prior.V6"),variable.name = "bayes",value.name = "value")%>%mutate(group="V6")
z.long7 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V7","prior.V7"),variable.name = "bayes",value.name = "value")%>%mutate(group="V7")
z.long8 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V8","prior.V8"),variable.name = "bayes",value.name = "value")%>%mutate(group="V8")
z.long9 <- melt(z.wide,id.vars=("x"),measure.vars = c("post.V9","prior.V9"),variable.name = "bayes",value.name = "value")%>%mutate(group="V9")

z <- rbind(z.long1,z.long2,z.long3,z.long4,z.long5,z.long6,z.long7,z.long8,z.long9)


ggplot(z)+
  geom_line(aes(x,value,colour = factor(bayes)))+
  facet_wrap(~ group)
  #facet_grid(. ~ group)

ggplot(z.wide, aes(x)) +
  geom_line(aes(y = post.V1, colour = "post"), size = 0.5) +
  geom_line(aes(y = prior.V1, colour = "prior"), size = 0.5)

postmean<-(alpha[1]+y)/(alpha[1]+y+n-y+beta[1])


#例题 11.3.2
x=seq(0,1,0.01)
n=5
y <- c(0,1,2,3,4,5)
alpha <- 1
beta <- 1
z <- as.data.frame(matrix(0,nrow = length(x),ncol=length(y)))
for(i in c(1:6)){
  z[,i] <- dbeta(x,alpha+y[i],n-y[i]+beta)
}

matplot(x, z, ylim=c(0,4), xlab="x", ylab="",
          col=1:6, type="l", lwd=2)
text(0.20,3,"y=0")
text(0.25,2.5,"y=1")
text(0.4,2.2,"y=2")
text(0.55,2.2,"y=3")
text(0.75,2.5,"y=4")
text(0.95,3,"y=5")
text(0.08,1.1,"Prior")