rfm_mcmc_exec.R

library(dplyr)
library(tidyr)
library(ggplot2)
library(rstan)
library(loo)
library(ggmcmc)

# read datasets
df <- read.csv("rfm.csv", stringsAsFactors = F)
colnames(df)[1] <- "ID"
df.spending <- read.csv("spending_mat.csv", stringsAsFactors = F)
spending.num <- read.csv("spending_freq.csv", stringsAsFactors = F, header=F)
colnames(spending.num) <- c("ID","num")

df.spending <- inner_join(df.spending, dplyr::select(df, ID), by=c("CustomerID" = "ID"))
spending.num <- inner_join(spending.num, dplyr::select(df, ID), by="ID")

# stan option settings
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# compiling
model <- stan_model("rfm_hierarchical1702.stan",
                    boost_lib = path.expand("~/boost_1_62_0/include/"),
                    eigen_lib = "/usr/local/include/eigen3")

# annualy discount = 15 %
RFM.res <- sampling(model, data=list(N=nrow(df),
                                     Time=df$Time,
                                     time=df$time,
                                     x=df$Freq,
                                     K=ncol(df.spending),
                                     NSpend=spending.num$num,
                                     Spend=as.matrix(df.spending),
                                     delta=.0027), 
                    warmup=4000, iter=6000, chain=4, 
                    pars="param3", include=F)
save(RFM.res, file=paste("result",format(Sys.Date(),"%y%m%d"), ".RData"), sep="/")
print(RFM.res, 'theta0')
print(RFM.res, 'Gamma0')
print(RFM.res, "lkj")
# サンプリングが健全か確認
data.frame(summary(RFM.res, "xi")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)
data.frame(summary(RFM.res, "tau")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)
data.frame(summary(RFM.res, "omega")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)
data.frame(summary(RFM.res, "lambda")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)
data.frame(summary(RFM.res, "mu")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)
data.frame(summary(RFM.res, "eta")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)
data.frame(summary(RFM.res, "CLV")$summary) %>% filter(Rhat >= 1.1 || n_eff/4/2000 <= .1)


ggDenseHisto <- function(result, family){
  ggmcmc(ggs(result), family=family, file=paste0("histo", gsub("\\W", "", family), ".pdf"), plot="histogram")
  ggmcmc(ggs(result), family=family, file=paste0("trace", gsub("\\W", "", family), ".pdf"), plot="traceplot")
}
ggDenseHisto(RFM.res, "theta0")
ggDenseHisto(RFM.res, "Gamma0")