#Scale Mixture RI

rm(list=ls())
setwd("~/GitHub/OWN_LOCAL_PATH/Section3/3.4")

library("RcppArmadillo")
library(bayesm)
source('RIDC2_SIM.R')
Rcpp::sourceCpp('BA.cpp')
Rcpp::sourceCpp('RIDCM_MH_hier_singlecmplx_parallel_scale_mix_one_inside.cpp')


P = 1000#number of individuals
N = 20  #number of choices per individual
lambda = 0.25

#attribute levels
rho_simple=c(2.5,3,3.5,4,4.5) #simple attribute levels
rho_complex=c(0,0.5,1,1.5,2)  #complex attribute levels

#generating heterogenous preferences

# settings RI model
states=rbind(rho_complex,rep(0,length(rho_complex)))
set.seed(1000)
  pricecoefs=rnorm(P,-1,0.2) #drawing price coefficients in population
  beta_price=matrix(c(rnorm(P, 2.5, 0.5), #brand coefficients
                      pricecoefs,
                      -pricecoefs), ncol= 3)
  
  max(pricecoefs) #check if all pricecoefs are negative
  
  loglambda_all=rnorm(P,-1.4,0.3)
  
  #Simulating Data - Version Feb 2022
 # set.seed(1)
  SIM_LIST<-NULL
  for (i in 1:P){
    SIM_LIST[[i]]=RIDC2_Sim(beta_price[i,], exp(loglambda_all[i]), rho_complex, rho_simple,N, complex="price")
  }
  
  #collecting Data in List form
  Data = vector(mode="list",P)
  for(i in 1:P){
    Data[[i]]=vector(mode="list",N)
  }
  
  for(j in 1:P){
    for(i in 1:N){
      Data[[j]][[i]]$X = SIM_LIST[[j]][[i]]$X
      Data[[j]][[i]]$y = SIM_LIST[[j]][[i]]$choice
      Data[[j]][[i]]$v = SIM_LIST[[j]][[i]]$curr_state
    }
  }
  
  

 
  out_RI=RIDC_MH_hier_cpp_parallel_joint(Data=Data, R=10000, lambda=lambda,
                           simp_idx=c(0,1), comp_idx=2, states=states,
                           Bbar=matrix(rep(0,2),ncol=2),
                           Areg=matrix(0.01),NUreg=4,Vreg=diag(1.5,2), 
                           propvar=diag(0.01,2),
                           beta_start=c(2.4,-0.95),
                           hier_beta_start=c(2.4,-0.95),hier_var_start=diag(2))
  
  out_scale_mix=RIDC_MH_hier_scale_mix_joint_cpp(Data=Data, R=10000,
                                          simp_idx=c(0,1), comp_idx=2, states=states,
                                          Bbar=matrix(rep(0,2),ncol=2),
                                          Areg=matrix(0.01),NUreg=4,Vreg=diag(1.5,2), 
                                          propvar=diag(0.1,2),
                                          beta_start=c(2.4,-0.95),#colMeans(beta_price[,1:2]),#c(4,-1.3,1.3),
                                          hier_beta_start=c(2.4,-0.95),hier_var_start=diag(2),
                                          loglambda_start=-1.4, propsdlambda = 0.2, loglambdafix = -1.4)
  
summary.bayesm.mat(t(out_RI$hiermean[,1001:10000]))
summary.bayesm.mat(t(out_scale_mix$hiermean[,1001:10000])) 

plot(out_scale_mix$hiermean[1,])
plot(out_RI$hiermean[1,])
plot(out_RI$hiermean[2,])

plot(colSums(out_RI$lhdraw[,1001:10000]))
logMargDenNR(colSums(out_scale_mix$lhdraw[,2001:10000]))
logMargDenNR(colSums(out_RI$lhdraw[,2001:10000]))  

hiervars_RI<-hiervars_RI_scale_mix<-matrix(nrow=10000, ncol=2)
for (i in 1:10000){
  hiervars_RI[i,]=diag(out_RI$hiervar[,,i])
  hiervars_RI_scale_mix[i,]=diag(out_scale_mix$hiervar[,,i])
  
}

summary.bayesm.mat(hiervars_RI[2001:10000,])
summary.bayesm.mat(hiervars_RI_scale_mix[2001:10000,])