###5章 予測分布を使って不確実性を計算



#####5.3 ラーメン屋の利益の計算
###データの入力
x <- c(94,84,108,101,103)#売上量
z <- c(244,257,242,246,251)#費用/杯
n <- length(x)#サンプル・サイズ

###事前分布のハイパーパラメータ
#ポアソン分布の事前分布
g <- 0.001;h <- 0.001#ガンマ分布
#正規分布の事前分布
m0 <- 250;b <- 0.001#正規分布
v <- 0.001;s0 <- 3#逆カイ二乗分布
#乱数の発生回数
R <- 100000

###事後分布からパラメータの乱数発生
#aの事後分布から乱数発生
a <- rgamma(R,sum(x)+g,n+h)# ガンマ分布から乱数発生
#s2の事後分布から乱数発生
v1 <- n+v
s1 <- s0+crossprod(z-mean(z))+(m0-mean(z))^2/(1/n+1/b)
s2 <- 1/rgamma(R,v1/2,s1/2)#逆カイ二乗分布（ガンマ分布を利用）から乱数発生
#mの事後分布から乱数発生
m_hat <- (b*m0+sum(z))/(n+b)
m <- rnorm(R,m_hat,sqrt(s2/(n+b)))#正規分布から乱数発生

###予測分布のシミュレーション
x_new <- rpois(R,a)#xの予測分布から乱数発生
z_new <- rnorm(R,m,sqrt(s2)) #zの予測分布から乱数発生
profit_new <- x_new*(500-z_new)#利益の予測分布

#利益の予測分布のプロット
plot(density(profit_new),xlab="利益",ylab="事後確率密度",main="利益の予測分布")

#95%信用区間はどこからどこまでか
quantile(profit_new,c(0.025,0.975))
#利益が３万円を上回る確率はどれくらいか
mean(profit_new>30000)
#利益が２万円を２日連続で下回る確率はどれくらいか
#再度計算乱数で利益の予測分布をシミュレーション
a1 <- rgamma(R,sum(x)+g,n+h)
s21 <- 1/rgamma(R,v1/2,s1/2);m1 <- rnorm(R,m_hat,sqrt(s21/(n+b)))
x_new1 <- rpois(R,a1);z_new1 <- rnorm(R,m,sqrt(s21))
profit_new1 <- x_new1*(500-z_new1)
mean(profit_new<20000&profit_new1<20000)



###5.4 クーポン付きダイレクトメールの割引率
###作業場所の指定
setwd("xxxxxx")#下のCSVデータの置き場所の指定、各自変えてください
#データの読み込み
dat <- read.csv("クーポンデータ.csv",row.names=1)

#各割引率の来店人数
n <- table(dat[,"割引率"])
#クーポンの割引率
dr <- as.numeric(names(n))
#実験クーポンの数
m <- length(dr)
#DMを出した数
N <- rep(300,m)

#事前分布のハイパーパラメータの設定
#各割引率グループで同一とする
a <- 1;b <- 1#p_jの事前分布（二項分布）
m0 <- c(0,0)#mの事前分布1（２次元正規分布）
h <- 1/10000#mの事前分布2（２次元正規分布）
v <- 5#Sの事前分布1（逆ウィシャート分布）
V0 <- diag(1/1000,2)#Sの事前分布2（逆ウィシャート分布）

#クーポンを送る人数
N_new <- 30000
#MCMCの回数
R <- 10000
#各割引率の利益の予測分布
profit <- NULL#予測分布の乱数を入れるリスト
for(j in 1:m) {#各割引率グループで繰り返す
  #割引率が同一のデータの抜き出し
  dat_j <- dat[dat[,"割引率"]==dr[j],]
  #対数化した購買金額と原価率
  Z_j <- cbind(log(dat_j[,"購買金額"]),log(dat_j[,"原価率"]))
  
  #事前の計算
  zj_bar <-colMeans(Z_j)
  V1j <- V0+cov(Z_j)*(n[j]-1)+(h*n[j])*tcrossprod(zj_bar-m0)/(h+n[j])
  m1j <- (n[j]*zj_bar+h*m0)/(n[j]+h)
  
  #予測分布のシミュレーション
  profit_j <- double(R)#乱数を入れるベクトル
  for(r in 1:R) {
    #事後分布からのサンプリング
    p_j <- rbeta(1,n[j]+1,(N[j]-n[j])+1)#反応率
    S_j <- solve(rWishart(1,v+n[j],solve(V1j))[,,1])#購買金額と原価率の分散
    m_j <- as.vector(t(chol(S_j/(h+n[j])))%*%rnorm(2)+m1j)#購買金額と原価率の平均
    #予測分布からのサンプリング
    n_j_new <- rbinom(1,N_new,p_j)#反応人数
    Z_j_new <- t(t(chol(S_j))%*%matrix(rnorm(n_j_new*2),2,n_j_new)+m_j) 
    #利益の予測分布
    mv_new <- exp(Z_j_new[,1])#購買金額
    cr_new <- exp(Z_j_new[,2])#原価率
    profit_j[r] <- sum(mv_new*(1-dr[j]-cr_new))-62*N_new
  }
  profit[[j]] <- profit_j
  cat("割引率",dr[j]*100,"%\n",sep="")
}
#予測分布の統計量
out <- NULL
for(j in 1:m) {
  out_j <- c(mean(profit[[j]]),sd(profit[[j]]),
             quantile(profit[[j]],probs=c(0.025,0.5,0.975)))
  out <- rbind(out,out_j)
}
colnames(out) <- c("事後平均","標準偏差","2.5%点","メディアン","97.5%点")
rownames(out) <- paste(dr*100,"%",sep="")
print(out)

#事後確率密度のプロット  
par(mfrow=c(2,2))#2×2に設定
for(j in 1:m) {
  plot(density((profit[[j]])/1000000),main=paste("割引率",dr[j]*100,"%",sep=""),
       xlab="利益（100万円）",ylab="事後確率密度",ylim=c(0,0.6),xlim=c(-1,10))
}
#一番利益が大きい確率を計算
Pm <- NULL#乱数列を行列化してから行で一番大きい回数を数える
for(j in 1:m) Pm <- cbind(Pm,profit[[j]])
table(apply(Pm,1,which.max))/R
