#Estimating change in lambda

rm(list=ls())
setwd("~/GitHub/OWN_LOCAL_PATH/Appendix/A_1_Delta_Lambda")

library("RcppArmadillo")
library(bayesm)
source('RIDCM2_SIM.R')
Rcpp::sourceCpp('BA_single_stable_vectorised.cpp')
Rcpp::sourceCpp('RIDCM_MH_singlecmplx_joint_lambda_delta.cpp')
Rcpp::sourceCpp('RIDCM_MH_singlecomplex_cpp_joint.cpp')



set.seed(66)
N = 1000#number of individuals
lambda_low = 0.1
lambda_high = 0.3
#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
betas=c(2.5,-1,1)

SIM_LIST1<-SIM_LIST2<-NULL
set.seed(66)
for (i in 1:N){
  SIM_LIST1[[i]]=RIDCM2_Sim(betas, lambda_low, rho_complex, rho_simple,1, comp_idx=3, brand="TRUE")
  SIM_LIST2[[i]]=RIDCM2_Sim(betas, lambda_high, rho_complex, rho_simple,1, comp_idx=3, brand="TRUE")
  
}



Data_low<-Data_high<-vector(mode="list",N)
for(i in 1:N){
  Data_low[[i]]$X = SIM_LIST1[[i]][[1]]$X
  Data_low[[i]]$y = SIM_LIST1[[i]][[1]]$choice
  Data_low[[i]]$v = SIM_LIST1[[i]][[1]]$curr_state
  Data_high[[i]]$X = SIM_LIST2[[i]][[1]]$X
  Data_high[[i]]$y = SIM_LIST2[[i]][[1]]$choice
  Data_high[[i]]$v = SIM_LIST2[[i]][[1]]$curr_state
} 

Data=c(Data_low, Data_high)

states=rbind(rho_complex,rep(0,length(rho_complex)))

out=RIDCM_MH_singlecomplex_joint_cpp_lambda(Data1=Data_high,Data2=Data_low,R=10000,
                                       simp_idx=c(0,1),comp_idx=c(2), #cpp indexing
                                       states=states,
                                       betavar=diag(0.01,2),beta_start=c(2.4,-0.9),
                                       priomean=c(0,0),priovar=diag(100,2),
                                       lambdastd = 0.01, loglambda_start =log(lambda_low) , lambda_fix = lambda_high,
                                       loglambdaprio = 0, loglambdasd = 5)

out_one_lambda=RIDCM_MH_singlecomplex_cpp_joint(Data=Data, R=10000, lambda=0.2,
                                                simp_idx=c(0,1),comp_idx=c(2), #cpp indexing
                                                states=states,
                                                stepsizes=diag(0.01,2),beta_start=c(2.4,-0.9),
                                                priomean=c(0,0),priovar=diag(100,2))


windows()
par(mfrow=c(2,3))
plot(out$betadraw[1,], type="l")
plot(out$betadraw[2,], type="l")
plot((out$loglambda), type="l")
plot(out_one_lambda$betadraw[1,], type="l")
plot(out_one_lambda$betadraw[2,], type="l")

#Table 14 - Results
#LMD
logMargDenNR(out$lhdraw[5001:10000])
logMargDenNR(out_one_lambda$lhdraw[5001:10000])

#Betas
summary.bayesm.mat(t(rbind(out$betadraw[,5000:10000], exp(out$loglambda[5000:10000]))))
summary.bayesm.mat(t(out_one_lambda$betadraw[,5001:10000]))