Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment