#### 第4章 ####
#### 4.1節 ####
#### 二項プロビットモデル
# データの読み込み
leisure <- read.table("leisure.csv",sep=",",header=T)
summary(leisure)

# MCMCpackの適用
library(MCMCpack)
pro.post1 <- MCMCprobit(sat2~gend+age+budget, data=leisure, 
burnin = 1000, mcmc = 10000, b0 = 0, B0 = 0.001)
summary(pro.post1)
# 図4.1
plot(pro.post1)
# bayemsの適用
library(bayesm)
Data2 <- list(X=cbind(leisure$gend, leisure$age, leisure$budget), y=leisure$sat2)
mcmc2 <- list(R=10000,keep=1)
pro.post2 <- rbprobitGibbs(Data=Data2,Mcmc=mcmc2)
summary(pro.post2$betadraw)
# 図4.2
plot(pro.post2$betadraw)

# 稼働検査期間を指定して変数を要約
summary(pro.post2$betadraw, burnin=1000)

#le.post2 <- MCMCprobit(sat2~gend+age+person+repeat+jisha+sizen+kaimono+bunka+sonota+budget, data=leisure, burnin = 1000, mcmc = 10000, b0 = 0, B0 = 0.001)
# summary(le.post2$betadraw, burnin=1000)

# WinBUGS用コード
# probit.txt
# Lancaster (2006) を加筆修正
# 定数項を用いていない
mode{
for(i in 1:n){
y[i] ~ dbin(p[i],1)
p[i] <- phi(b1*X[i,1]+b2*X[i,2]+b3*X[i,3])}
b1 ~ dnorm(0, 0.001)
b2 ~ dnorm(0, 0.001)
b3 ~ dnorm(0, 0.001)
}

# Rコード
# プロビットモデルの推定が（初期値を変えても）うまくいかない場合、
# Rのバージョンを落とすとうまく推定できることがある
library(R2WinBUGS)
X=cbind(leisure$gend, leisure$age, leisure$budget)
y=leisure$sat2
n<-nrow(X)
data <- list("n", "y", "X")

# 各初期値を定義する
#in1 <- list(b1=rnorm(1,0,.0001), b2=rnorm(1,0,.0001), b3=rnorm(1,0,.0001))
#in2 <- list(b1=rnorm(1,0,.0001), b2=rnorm(1,0,.0001), b3=rnorm(1,0,.0001))
#in3 <- list(b1=rnorm(1,0,.0001), b2=rnorm(1,0,.0001), b3=rnorm(1,0,.0001))
#inits <- list(in1,in2,in3)

# function()関数を使って初期値を定義する
#inits<-function(){
#list (
#b1=rnorm(1,0,.2),
#b2=rnorm(1,0,.2),
#b3=rnorm(1,0,.2))}  

# function()関数を使って初期値に定数を与える
#inits<-function(){
#list(b0=0,b1=1,b2=1,b3=1)}

# probit回帰モデルのパラメータを初期値に与える
# http://gking.harvard.edu/zelig/docs/
source("http://gking.harvard.edu/zelig/install.R")
library(Zelig)
pro.out <- zelig(sat2~gend+age+budget,model="probit", data=leisure)
summary(pro.out)
inits<-function(){
list (
b1=pro.out$coeff[[2]],
b2=pro.out$coeff[[3]],
b3=pro.out$coeff[[4]])}

parameters <- c("b1", "b2", "b3")
pro.post3 <- bugs(data, inits, parameters, 
model.file="C:/Program Files/R/R-2.7.2/probit.txt", debug=TRUE,
n.chains=3, n.iter=10000, n.burnin=1000,
bugs.directory = "C:/Program Files/WinBUGS14/",
working.directory = "C:/Program Files/R/R-2.7.2")
pro.post3
# 図4.3
plot(pro.post3)

# WinBUGSで「Trap」ウィンドウが表示され、「undefined real result」という
# メッセージが表示された場合、初期値をすべて違う値にするなど、
# 初期値や事前情報を改善するとよい。

#### 4.2節 ####
#### 二項ロジットモデル
# WinBUGS用コード
# logit.txt
# Lancaster (2006) を参照
# 定数項を用いていない
mode{
for(i in 1:n){
y[i] ~ dbin(p[i],1)
logit(p[i]) <- b1*X[i,1]+b2*X[i,2]+b3*X[i,3]}
b1 ~ dnorm(0, 0.001)
b2 ~ dnorm(0, 0.001)
b3 ~ dnorm(0, 0.001)
}

# Rコード
library(R2WinBUGS)
X=cbind(leisure$gend, leisure$age, leisure$budget)
y=leisure$sat2
n<-nrow(X)
data <- list("n", "y", "X")

# 初期値をランダムな正規分布に従って与える
#inits<-function(){
#list (
#b1=rnorm(1,0,.2),
#b2=rnorm(1,0,.2),
#b3=rnorm(1,0,.2))}  

library(Zelig)
log.out <- zelig(sat2~gend+age+budget-1,model="logit", data=leisure)
summary(log.out)
inits<-function(){
list (
b1=log.out$coeff[[1]],
b2=log.out$coeff[[2]],
b3=log.out$coeff[[3]])}

parameters <- c("b1", "b2", "b3")
log.post2 <- bugs(data, inits, parameters, 
model.file="C:/Program Files/R/R-2.7.2/logit.txt", debug=FALSE,
n.chains=3, n.iter=10000, n.burnin=1000,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2")
print(log.post2, digits=3)
# 図4.4
plot(log.post2)

# MCMCpackのMCMClogitを用いる方法
library(MCMCpack)
log.post1 <- MCMClogit(sat2~gend+age+budget-1, data=leisure, burnin = 1000, mcmc = 10000, b0 = 0, B0 = 0.001)
summary(log.post1)
# 図4.5
plot(log.post1)

#### 4.3節 ####
#### トビットモデル
# library(MCMCpack)
tob.post1 <- MCMCtobit(sat2~gend+age+budget-1, data=leisure, burnin = 1000, mcmc = 10000, b0 = 0, B0 = 0.001)
summary(tob.post1)
# 図4.6
plot(tob.post1)

# WinBUGS用コード
# tobit.txt 
# Lancaster(2006)を加筆修正
# ct=0.39894=1/sqrt(2*pi)
mode{
for(i in 1:n){
ones[i] <-1
ones[i] <- dbern(p[i])
term1[i] <- ct*pow(tau,1/2)*exp(-1/2*tau*pow(y[i]-mu[i],2))
term2[i] <- phi(-mu[i]*pow(tau,1/2))
p[i] <- pow(term1[i],ind[i])*pow(term2[i],1-ind[i])/10000
mu[i] <- b1*X[i,1]+b2*X[i,2]+b3*X[i,3]}
b1 ~ dnorm(0, 0.001)
b2 ~ dnorm(0, 0.001)
b3 ~ dnorm(0, 0.001)
tau ~ dgamma(0.001, 0.001)
}

#Rコード
library(R2WinBUGS)
library(survival)
library(Zelig)

X=cbind(leisure$gend, leisure$age, leisure$budget)
y=leisure$sat2
n<-nrow(X)
K <- 10000
ct <-1/sqrt(2*3.1416)
ind <- 1*(y>0)
data <- list("n", "K", "ct", "y", "ind", "X")

#result<-survreg(Surv(y, y > 0, type ='left')~X[,1]+X[,2]+X[,3], dist = 'gaussian')
# par.1 <- result$coef

#inits <- function(){
#list(
#b1=result$coef[[1]], 
#b2=result$coef[[2]],
#b3=result$coef[[3]],
#tau=dgamma(1, 0.1, 0.1))} 

tob.out <- zelig(sat2~gend+age+budget-1,model="tobit", data=leisure)
summary(tob.out)
inits<-function(){
list (
b1=tob.out$coeff[[1]],
b2=tob.out$coeff[[2]],
b3=tob.out$coeff[[3]],
tau=dgamma(1, 0.1, 0.1))}

parameters <- c("b1", "b2", "b3")

tob.post2 <- bugs(data, inits, parameters, 
model.file="C:/Program Files/R/R-2.7.2/tobit.txt", debug=TRUE,
n.chain=1, n.iter=10000, n.burnin=1000,
n.thin=1, DIC = TRUE, digits =4, 
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2")
print(tob.post2,digits=3)

# 図4.7
plot(tob.post2)

#### 4.4節 ####
#### 順序プロビットモデル
# library(MCMCpack)
oprob.post1 <- MCMCoprobit(as.factor(sat1)~gend+age+budget, data=leisure, burnin = 1000, mcmc = 10000, b0 = 0, B0 = 0.001)
summary(oprob.post1)
# 図4.8
plot(oprob.post1)

# 順序プロビット回帰モデル
oprob.out <- zelig(as.factor(sat1)~gend+age+budget,model="oprobit", data=leisure)
summary(oprob.out)

#### 4.5節 ####
#### 多項プロビットモデル
# library(MCMCpack)
# library(bayesm)

# データ読み込み
leisure2 <- read.table("leisure2.csv",sep=",",header=T) 
p<-4 			# 選択肢数
n<-nrow(leisure2) 	# サンプル数
# 選択肢に依存する属性
na <- 2			
price <- cbind(leisure2$price1,leisure2$price2,leisure2$price3,leisure2$price4)
days <- cbind(leisure2$day1,leisure2$day2,leisure2$day3,leisure2$day4)
Xa <- cbind(price, days)
# 選択肢に依存しない属性（個人属性）
nd <- 2
gend <- leisure2$gend
age <- leisure2$age
Xd <- cbind(gend,age)
# 説明変数
X <- createX(p,na=na,nd=nd,Xa=Xa,Xd=Xd,DIFF=TRUE,base=1)

data1 <- list(p=p, y=leisure2$choice, X=X)
mcmc1 <- list(R=2000, keep=1)

le2.rmnp.post1 <- rmnpGibbs(Data=data1,Mcmc=mcmc1)
# 図4.9
plot(le2.rmnp.post1$betadraw)
# 図4.10
plot(le2.rmnp.post1$sigmadraw)

# マーケティングデータをランダムに生成する
# 説明変数行列を生成
p<-4; n<-100
gend <- sample(c(0, 1), n, replace = TRUE)
age <- sample(c(0, 1), n, replace = TRUE)
price <- matrix(runif(n*p,min=0,max=2),ncol=4,nrow=n)
days <- matrix(runif(n*p,min=0,max=2),ncol=4,nrow=n)
na <- 2			
Xa <- cbind(price, days)
nd <- 2
Xd <- cbind(gend,age)
X <- createX(p=p,na=na,nd=nd,Xa=Xa,Xd=Xd,DIFF=TRUE,base=1)

# 被説明変数ベクトルを生成
beta <- c(runif(ncol(X),min=0,max=1))
sigma <- matrix(c(1,.33,.33,.33,1,.33,.33,.33,1),ncol=(p-1))
Xbeta <- X%*%beta
w <- as.vector(crossprod(chol(sigma),matrix(rnorm((p-1)*n),ncol=n)))+ Xbeta
w <- matrix(w,ncol=(p-1),byrow=TRUE)
maxw <- apply(w,1,max)
indmax <- function(x) {which(max(x)==x)}
y <- apply(w,1,indmax)
y <- ifelse(maxw < 0,p,y)

# 多項プロビット・モデルを推定
data <- list(p=p, y=y, X=X)
mcmc <- list(R=1000, keep=1)
le2.rmnp.post2=rmnpGibbs(Data=data,Mcmc=mcmc)
# 図4.11
plot(le2.rmnp.post2$betadraw)
# 図4.12
plot(le2.rmnp.post2$sigmadraw)

#### 4.6節 ####
#### 多変量プロビットモデル
leisure <- read.table("leisure.csv",sep=",",header=T)
y <- as.matrix(leisure[,7:11])
p=ncol(y); n=nrow(y)
dimnames(y)=NULL
y=as.vector(t(y))
y=as.integer(y)
I_p=diag(p)
X=rep(I_p,n)
X=matrix(X,nrow=p)
X=t(X)
R=2000
Data=list(p=p,X=X,y=y)
Mcmc=list(R=R)
set.seed(66)
le.rmvp.post1=rmvpGibbs(Data=Data,Mcmc=Mcmc)
# 図4.13
plot(le.rmvp.post1$betadraw)
# 図4.14
plot(le.rmvp.post1$sigmadraw)




