Gedevan-Aleksizde

# --------------------------------
# calculate DIC* from stanfit
#
# * Spiegelhalter et al. (2002)
# --------------------------------

DIC <- function(stanfit, df.input, dev ){
  # stanfit: stanfit object
  # df.input: input data.frame
  # dev: function that calculate the dev. of post.mean; 

Gedevan-Aleksizde

require(stats4) # 3.3.1

nloglik_lnormal <- function(lmu, lsigma){
  # df columns: b.u(pper border), b.l(ower border), n(umber of households)
  return(-sum(with(temp,
                   n*log(plnorm(b.u, meanlog=lmu, sdlog = exp(lsigma))
                         - plnorm(b.l, meanlog=lmu, sdlog = exp(lsigma))
                   )
                   -log(n))
  ) - sum(log(1:sum(df$n)))

Gedevan-Aleksizde

glmnet.mock <- function(formula=y~., family="gaussian", data, ...){
  fam.link <- strsplit(family, "\\(|\\)")[[1]]
  family <- fam.link[1]
  if(length(fam.link) >=2) link <- fam.link[2]
  else link <- "identity"
  return(glmnet(x=model.matrix(formula, data),
                y=get(link)(model.response(model.frame(formula, data))),
                family = family,
                ...)
        )

Gedevan-Aleksizde

library(ggmcmc) # ver. 1.1
library(data.table)
ggs2 <- function(result, pattern=NULL, invert=F, ...){
  if(class(result)[1] =="stanfit" ){
    nThin       <- result@sim$thin
    nChains     <- result@sim$chains
    nBurnin     <- result@sim$warmup2[1]
    nIterations  <- result@sim$n_save[1] - nBurnin
    result <- as.data.frame(result) %>% 
      mutate(Chain=rep(1:nChains, nIterations), Iteration=rep(1:nIterations, nChains)) %>% 

Gedevan-Aleksizde

/* -----------------------------------------------------------------------------
Estimate latent variable hierarchical Bayes RFM model 

stan ver. 2.14

CREATED: by ill-identified, 18/APR/2016 
REVISED: 25/FEB/2017, REVISED
------------------------------------------------------------------------------- */

data {

Gedevan-Aleksizde

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"

Gedevan-Aleksizde

# multivariate normal Kalman filter
require(dplyr)
require(tidyr)
require(ggplot2)
require(animation)

# ARIMA(1,1) + 線形トレンド の乱数生成
N      <- 50
phi1   <- .5
theta1 <- .2

Gedevan-Aleksizde

require(mlogit) # 0.2-4

as.mldata <- function(data){
  # convert HC dataset
  # The alternatives are
  # 1. Gas central heat with cooling (gcc)
  # 2. Electric central resistence heat with cooling (ecc)
  # 3. Electric room resistence heat with cooling (erc)
  # 4. Electric heat pump, which provides cooling also (hpc)
  # 5. Gas central heat without cooling (gc)

Gedevan-Aleksizde

# ------ common part ----

require(ggplot2)
require(dplyr)
require(tidyr)
require(dlm) # 1.1-4
require(KFAS) # 1.2.9
require(bsts) # 0.7.1

# calculate the reference price

Gedevan-Aleksizde

source("common.R", encoding = "utf-8")

# ---- dlm ------

# model 3
res <- data.frame()
for( a in seq(from=.1, to=.95, by=.05) ){
  RP <- calc_RP(RP = df$RP095, p = df$AP, a = a)
  z <- calc_Z(RP = RP, p = df$AP)
  z1 <- z$z1