.Rhistory

counter=counter+1
natur=-C[,counter]
counter=counter+1
natur_km=C[,counter]
counter=counter+1
night=C[,counter]
counter=counter+1
sea=C[,counter]
counter=counter+1
sp=C[,counter]
counter=counter+1
it=C[,counter]
counter=counter+1
bg=C[,counter]
v1 <-  cost * (         accom * a_accom  + aisec * a_aisec	 + beach * a_beach	 + beach_km * a_beach_km	 + clear * a_clear	 + a_cost	 + crowd * a_crowd	 + cultr * a_cultr	 + cultr_km * a_cultr_km	 + days * a_days	 + natur * a_natur	 + natur_km * a_natur_km	 + night * a_night	 + sea * a_sea	 + bg * a_bg	 + it * a_it	 + sp * a_sp)
v2 <-  cost * ( asc_b + accom * b_accom	 + aisec * b_aisec	 + beach * b_beach	 + beach_km * b_beach_km	 + clear * b_clear	 + b_cost	 + crowd * b_crowd	 + cultr * b_cultr	 + cultr_km * b_cultr_km	 + days * b_days	 + natur * b_natur	 + natur_km * b_natur_km	 + night * b_night	 + sea * b_sea	 + bg * b_bg	 + it * b_it	 + sp * b_sp)
v3 <-  cost * ( asc_c + accom * c_accom	 + aisec * c_aisec	 + beach * c_beach	 + beach_km * c_beach_km	 + clear * c_clear	 + c_cost	 + crowd * c_crowd	 + cultr * c_cultr	 + cultr_km * c_cultr_km	 + days * c_days	 + natur * c_natur	 + natur_km * c_natur_km	 + night * c_night	 + sea * c_sea	 + bg * c_bg	 + it * c_it	 + sp * c_sp)
v4 <-  cost * ( asc_d + accom * d_accom	 + aisec * d_aisec	 + beach * d_beach	 + beach_km * d_beach_km	 + clear * d_clear	 + d_cost	 + crowd * d_crowd	 + cultr * d_cultr	 + cultr_km * d_cultr_km	 + days * d_days	 + natur * d_natur	 + natur_km * d_natur_km	 + night * d_night	 + sea * d_sea	 + bg * d_bg	 + it * d_it	 + sp * d_sp)
v1     <- (v1 <= -700)*-700 + ( v1 > -700 & v1 < 700 ) * v1 + (v1 >= 700)*700
v2     <- (v2 <= -700)*-700 + ( v2 > -700 & v2 < 700 ) * v2 + (v2 >= 700)*700
v3     <- (v3 <= -700)*-700 + ( v3 > -700 & v3 < 700 ) * v3 + (v3 >= 700)*700
v4     <- (v4 <= -700)*-700 + ( v4 > -700 & v4 < 700 ) * v4 + (v4 >= 700)*700
v1     <- exp(v1)
v2     <- exp(v2)
v3     <- exp(v3)
v4     <- exp(v4)
p=( choice1 * v1 + choice2 * v2 + choice3 * v3 + choice4 * v4  ) / ( v1 + v2 + v3 + v4 )
return(p)
}
doHB(likelihood, choicedata, control)
#install.packages("RSGHB")
library(RSGHB)
setwd("C:\\Work\\Papers\\Individual Models")			    # working directory
modelname <- "HB_4_WTP_trunc_normal"			# used for output
# ------------------
# DATA
# ------------------
choicedata <- read.table("data.csv",sep=",",header=T)
attach(choicedata)
choice1 = ( choice == 1 ) * 1
choice2 = ( choice == 2 ) * 1
choice3 = ( choice == 3 ) * 1
choice4 = ( choice == 4 ) * 1
Choice = choice1 * 1 + choice2 * 2 + choice3 * 3 + choice4 * 4
# ------------------
# MODEL SPECIFICATION
# ------------------
gVarNamesNormal <- c(
"ASC_B",
"ASC_C",
"ASC_D",
"ACCOM",
"AISEC",
"BEACH",
"BEACH_KM",
"CLEAR",
"COST",
"CROWD",
"CULTR",
"CULTR_KM",
"DAYS",
"NATUR",
"NATUR_KM",
"NIGHT",
"SEA",
"BG",
"IT",
"SP"
)
# For each variable, specify the distribution for its coefficient
# The options are:
# 1. normal
# 2. log-nomal
# 3. negative log-normal
# 4. normal with all values below zero massed at zero
# 5. Johnson SB with a specified min and max
# gDIST must have NIV elements
gDIST <- c(rep(1,6),4,4,4,1,4,4,4,4,4,rep(1,5))
rho <- c(
0.0043048,
0.0045383,
0.0043446,
0.0036161,
0.0033738,
0.0068626,
0.0048000,
0.00080475,
0.010786,
0.0039215,
0.0035161,
0.00017928,
0.0011455,
0.0039535,
0.00015251,
0.0035765,
0.0060286,
0.0086245,
0.0093410,
0.0088848)
svN <- rep(0,length(gVarNamesNormal))
svN <- c( 0.011572,
0.016189,
0.00014022,
-0.23862,
-0.056640,
-0.30269,
0.21136,
-0.0099214,
-2.4170,
0.11208,
-0.17610,
0.0045728,
-0.071538,
-0.095326,
0.0024947,
-0.18429,
-0.32792,
-0.19249,
-0.40145,
-0.48391)
#cbind(gVarNamesNormal, svN,rho,gDIST)
# ------------------
# ITERATION SETTINGS
# ------------------
gNCREP    <- 25000  	  # Number of iterations to use prior to convergence
gNEREP    <- 10000 	       # Number of iterations to keep for averaging after convergence has been reached
gNSKIP    <- 1			  # Number of iterations to do in between retaining draws for averaging
gINFOSKIP <- 500           # How frequently to print info about the iteration process
control <- list(
modelname=modelname,
gVarNamesNormal=gVarNamesNormal,
gDIST=gDIST,
svN=svN,
gNCREP=gNCREP,
gNEREP=gNEREP,
gNSKIP=gNSKIP,
gINFOSKIP=gINFOSKIP,
nodiagnostics=F,
rho=rho
)
# ------------------
# likelihood
# USE:     Calculates the likelihood of choice | B
#          Returns individual-specific likelihood values
# NOTES:   This is where the bulk of the computation resides so coding this efficiently
#	      is essential to reducing run time.
# ------------------
likelihood <- function(fc,C)
{
# creates vectors of individual-specific coefficients
counter =1;
asc_b=C[,counter]
counter=counter+1
asc_c=C[,counter]
counter=counter+1
asc_d=C[,counter]
counter=counter+1
accom=C[,counter]
counter=counter+1
aisec=C[,counter]
counter=counter+1
beach=C[,counter]
counter=counter+1
beach_km=C[,counter]
counter=counter+1
clear=-C[,counter]
counter=counter+1
cost=-C[,counter]
counter=counter+1
crowd=C[,counter]
counter=counter+1
cultr=-C[,counter]
counter=counter+1
cultr_km=C[,counter]
counter=counter+1
days=-C[,counter]
counter=counter+1
natur=-C[,counter]
counter=counter+1
natur_km=C[,counter]
counter=counter+1
night=C[,counter]
counter=counter+1
sea=C[,counter]
counter=counter+1
sp=C[,counter]
counter=counter+1
it=C[,counter]
counter=counter+1
bg=C[,counter]
v1 <-  cost * (         accom * a_accom  + aisec * a_aisec	 + beach * a_beach	 + beach_km * a_beach_km	 + clear * a_clear	 + a_cost	 + crowd * a_crowd	 + cultr * a_cultr	 + cultr_km * a_cultr_km	 + days * a_days	 + natur * a_natur	 + natur_km * a_natur_km	 + night * a_night	 + sea * a_sea	 + bg * a_bg	 + it * a_it	 + sp * a_sp)
v2 <-  cost * ( asc_b + accom * b_accom	 + aisec * b_aisec	 + beach * b_beach	 + beach_km * b_beach_km	 + clear * b_clear	 + b_cost	 + crowd * b_crowd	 + cultr * b_cultr	 + cultr_km * b_cultr_km	 + days * b_days	 + natur * b_natur	 + natur_km * b_natur_km	 + night * b_night	 + sea * b_sea	 + bg * b_bg	 + it * b_it	 + sp * b_sp)
v3 <-  cost * ( asc_c + accom * c_accom	 + aisec * c_aisec	 + beach * c_beach	 + beach_km * c_beach_km	 + clear * c_clear	 + c_cost	 + crowd * c_crowd	 + cultr * c_cultr	 + cultr_km * c_cultr_km	 + days * c_days	 + natur * c_natur	 + natur_km * c_natur_km	 + night * c_night	 + sea * c_sea	 + bg * c_bg	 + it * c_it	 + sp * c_sp)
v4 <-  cost * ( asc_d + accom * d_accom	 + aisec * d_aisec	 + beach * d_beach	 + beach_km * d_beach_km	 + clear * d_clear	 + d_cost	 + crowd * d_crowd	 + cultr * d_cultr	 + cultr_km * d_cultr_km	 + days * d_days	 + natur * d_natur	 + natur_km * d_natur_km	 + night * d_night	 + sea * d_sea	 + bg * d_bg	 + it * d_it	 + sp * d_sp)
v1     <- (v1 <= -700)*-700 + ( v1 > -700 & v1 < 700 ) * v1 + (v1 >= 700)*700
v2     <- (v2 <= -700)*-700 + ( v2 > -700 & v2 < 700 ) * v2 + (v2 >= 700)*700
v3     <- (v3 <= -700)*-700 + ( v3 > -700 & v3 < 700 ) * v3 + (v3 >= 700)*700
v4     <- (v4 <= -700)*-700 + ( v4 > -700 & v4 < 700 ) * v4 + (v4 >= 700)*700
v1     <- exp(v1)
v2     <- exp(v2)
v3     <- exp(v3)
v4     <- exp(v4)
p=( choice1 * v1 + choice2 * v2 + choice3 * v3 + choice4 * v4  ) / ( v1 + v2 + v3 + v4 )
return(p)
}
doHB(likelihood, choicedata, control)
library(RSGHB)
library(RSGHB)
library(RSGHB)
dnorm(10)
dnorm(-10)
dnorm(0)
dnorm(2)
dnorm(3)
dnorm(1:5)
dnorm(1:5,0,1:5)
dnorm(rep(1,5),0,1:5)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
?read.csv
nextF <- function (p, f, b, env)
{
a <- c(-0.108123785948993,-0.0821217958019051,0.0834535602658302 ,0.283944999954507,0.0424935312479905,
-2.6899917151717,	  0.449021594851894 ,-0.0135230692499273,0.136333056332813,0.263959739292714,
-0.0089962292201314,0.233720509569121)
#d <- c(0.0188030584,0.0416955760,0.0097148904,0.7044705003,0.1776460913,0.8087013895,
#       0.3843694764,0.0001350046,0.0079670233,0.3590518176,0.0002171553,0.3173359108)
d  <- rep(1e-1,12)
fnew <- f + sqrt(env$rhoF) * rnorm(env$gFIV, 0, 1)
fn <- fnew - a
fo <- f - a
pnew <- env$likelihood(fnew, b, env)
#r <- sum(log(pnew) - log(p))
r <- sum(log(pnew) + log(dnorm(fnew,a,sd=sqrt(d))) - (log(p) + log(dnorm(f,a,sd=sqrt(d)))))
ind <- (log(runif(1)) <= r)
pnew <- pnew * ind + p * (1 - ind)
fnew <- as.vector(t(fnew %*% t(ind)) + t(f %*% t(1 - ind)))
return(list(fnew, pnew))
}
library(RSGHB)
library(RSGHB)
library(RSGHB)
nextF <- function (p, f, b, env)
{
a <- c(-0.16,
-0.01,
0.17,
0.63,
0.18,
-3.76,
0.61,
-0.01,
0.19,
0.29,
-0.01,
0.5     )
#d <- c(0.0188030584,0.0416955760,0.0097148904,0.7044705003,0.1776460913,0.8087013895,
#       0.3843694764,0.0001350046,0.0079670233,0.3590518176,0.0002171553,0.3173359108)
d  <- c(0.13314425,
0.05127517,
0.06377697,
1.0201,
0.67924,
5.75926464,
0.41839361,
0.000110723,
0.040088807,
0.32848473,
0.000371945,
0.69016271)
fnew <- f + sqrt(env$rhoF) * rnorm(env$gFIV, 0, 1)
fn <- fnew - a
fo <- f - a
pnew <- env$likelihood(fnew, b, env)
#r <- sum(log(pnew) - log(p))
r <- sum(log(pnew) + log(dnorm(fnew,a,sd=sqrt(d))) - (log(p) + log(dnorm(f,a,sd=sqrt(d)))))
ind <- (log(runif(1)) <= r)
pnew <- pnew * ind + p * (1 - ind)
fnew <- as.vector(t(fnew %*% t(ind)) + t(f %*% t(1 - ind)))
return(list(fnew, pnew))
}
nextF <- function (p, f, b, env)
{
a <- c(-0.16,
-0.01,
0.17,
0.63,
0.18,
-3.76,
0.61,
-0.01,
0.19,
0.29,
-0.01,
0.5     )
#d <- c(0.0188030584,0.0416955760,0.0097148904,0.7044705003,0.1776460913,0.8087013895,
#       0.3843694764,0.0001350046,0.0079670233,0.3590518176,0.0002171553,0.3173359108)
d  <- 0.1*c(0.13314425,
0.05127517,
0.06377697,
1.0201,
0.67924,
5.75926464,
0.41839361,
0.000110723,
0.040088807,
0.32848473,
0.000371945,
0.69016271)
fnew <- f + sqrt(env$rhoF) * rnorm(env$gFIV, 0, 1)
fn <- fnew - a
fo <- f - a
pnew <- env$likelihood(fnew, b, env)
#r <- sum(log(pnew) - log(p))
r <- sum(log(pnew) + log(dnorm(fnew,a,sd=sqrt(d))) - (log(p) + log(dnorm(f,a,sd=sqrt(d)))))
ind <- (log(runif(1)) <= r)
pnew <- pnew * ind + p * (1 - ind)
fnew <- as.vector(t(fnew %*% t(ind)) + t(f %*% t(1 - ind)))
return(list(fnew, pnew))
}
library(RSGHB)
library(RSGHB)
setwd("C:\\Work\\Papers\\Individual Models\\Holdout Analysis\\E12_HO36\\conditionals")
data <- read.csv("..\\sampledobs.csv")
#data <- read.csv("..\\holdoutsample.csv")
#head(holdouts)
#parameters <- read.csv("Random.csv")
#parameters <- read.csv("MNL.csv")
#parameters <- read.csv("LCM.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv1e-3_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv1e-1_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv2e-1_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv2e-2_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv3e-2_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv4e-2_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv5e-2_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv8e-2_df1_C.csv")
#parameters <- read.csv("HB_4_WTP_normal_sample12_n741_pv1e-2_C.csv")
#parameters <- read.csv("HB_4_WTP_normal_sample12_n741_C.csv")
#parameters <- read.csv("HB_4_WTP_normal_sample12_n741_pv10_C.csv")
#parameters <- read.csv("ILM_HB_E12_HO36_pv0.1.csv")
#parameters <- read.csv("HB_ILM_PV(0.003df1).csv")
#parameters <- read.csv("HB_ILM_PV(0.003df1-mean only,pv0.1).csv")
#parameters <- read.csv("HB_ILM_PV(MSL).csv")
parameters <- read.csv("HB_ILM_PV(MSL,pv0.1CovMat).csv")
#parameters <- read.csv("Conditionals_E12_HO36.csv")
#parameters <- read.csv("HB_4_pref_normal_sample12_n741_pv1e-2_C.csv")
#parameters <- read.csv("HB_4_pref_normal_sample12_n741_C.csv")
#parameters <- read.csv("HB_4_pref_normal_sample12_n741_pv10_C.csv")
#parameters <- read.csv("HBLC_4_Pref_normal_sample12_n741_pv1e-2_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv1e-2_df1_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv1e-3_C.csv")
#parameters <- read.csv("HB_4_Pref_normal_sample12_n741_pv5e-3_df1_C.csv")
#parameters <- read.csv("ILM_ML.csv")
#parameters <- read.csv("ILM_OLS.csv")
#head(parameters)
merged.data <- merge(data,parameters,by="ID")
attach(merged.data)
# pref space
#V1 <-  (         accom * accom1     + aisec * aisec1    + cost * cost1	  + cultr * cult1	 + cultr_km * cult_km1	 + days * days1	 + natur * nat1	 + natur_km * nat_km1	 + night * night1	 )
#V2 <-  ( asc_b + accom * accom2	 + aisec * aisec2	 + cost * cost2	  + cultr * cult2	 + cultr_km * cult_km2	 + days * days2	 + natur * nat2	 + natur_km * nat_km2	 + night * night2	)
#V3 <-  ( asc_c + accom * accom3	 + aisec * aisec3	 + cost * cost3	  + cultr * cult3	 + cultr_km * cult_km3	 + days * days3	 + natur * nat3	 + natur_km * nat_km3	 + night * night3	)
#V4 <-  ( asc_d + accom * accom4	 + aisec * aisec4	 + cost * cost4	  + cultr * cult4	 + cultr_km * cult_km4	 + days * days4	 + natur * nat4	 + natur_km * nat_km4	 + night * night4	)
V1 <-  (         accom * accom1     + aisec * aisec1    + cost * cost1/1000       + cultr * cult1	 + cultr_km * cult_km1	 + days * days1	 + natur * nat1	 + natur_km * nat_km1	 + night * night1	 )
V2 <-  ( asc_b + accom * accom2	 + aisec * aisec2	 + cost * cost2/1000	  + cultr * cult2	 + cultr_km * cult_km2	 + days * days2	 + natur * nat2	 + natur_km * nat_km2	 + night * night2	)
V3 <-  ( asc_c + accom * accom3	 + aisec * aisec3	 + cost * cost3/1000	  + cultr * cult3	 + cultr_km * cult_km3	 + days * days3	 + natur * nat3	 + natur_km * nat_km3	 + night * night3	)
V4 <-  ( asc_d + accom * accom4	 + aisec * aisec4	 + cost * cost4/1000	  + cultr * cult4	 + cultr_km * cult_km4	 + days * days4	 + natur * nat4	 + natur_km * nat_km4	 + night * night4	)
# WTP space
#V1 <-  cost * (         accom * accom1   + aisec * aisec1    + cost1/1000    + cultr * cult1	 + cultr_km * cult_km1	 + days * days1	 + natur * nat1	 + natur_km * nat_km1	 + night * night1	 )
#V2 <-  cost * ( asc_b + accom * accom2	 + aisec * aisec2	 + cost2/1000	  + cultr * cult2	 + cultr_km * cult_km2	 + days * days2	 + natur * nat2	 + natur_km * nat_km2	 + night * night2	)
#V3 <-  cost * ( asc_c + accom * accom3	 + aisec * aisec3	 + cost3/1000	  + cultr * cult3	 + cultr_km * cult_km3	 + days * days3	 + natur * nat3	 + natur_km * nat_km3	 + night * night3	)
#V4 <-  cost * ( asc_d + accom * accom4	 + aisec * aisec4	 + cost4/1000	  + cultr * cult4	 + cultr_km * cult_km4	 + days * days4	 + natur * nat4	 + natur_km * nat_km4	 + night * night4	)
choice1 <- choice == 1
choice2 <- choice == 2
choice3 <- choice == 3
choice4 <- choice == 4
P <- (exp(V1)*choice1 + exp(V2)*choice2 + exp(V3)*choice3 + exp(V4)*choice4)/(exp(V1) + exp(V2) + exp(V3) + exp(V4))
RLH <- aggregate(P,by=list(ID),function(x)(prod(x))^(1/36))
summary(RLH$x)
sum(log(P))
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
checkmodel
checkModel
install.packages("compiler")
library(RSGHB)
hb
checkmodel
checkModel
library(RSGHB)
nextF
nextB
signif(0,3)
signif(0,9)
library(RSGHB)
hb
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
#}
library(RSGHB)
?chol.default
library(RSGHB)
library(RSGHB)
lcvar <- function(s, c = rep(1, ncol(s))) {
ci <- t(combn(1:length(c), 2))
cs <- matrix(c[ci], ncol = 2)
sum(c ^ 2 * diag(s)) + 2 * sum(apply(cs, MARGIN = 1, FUN = prod) * s[ci])
}
s <- matrix(c(1, -0.5, 2, -0.5, 1.3, 0.75, 2, 0.75, 0.67), nrow = 3)
c <- c(1, 0, -1/2)
lcvar(s, c)
s
ci <- t(combn(1:length(c), 2))
ci
c[ci]
diag(s)
job  <- 1:6
crop <- 1:3
prod <- c(0,5,10,20)
qual <- yield <- c(0,5,10)
window <- uptime <- c(0,10,20)
fuel   <- c(5,0,-10,-20)
acost  <- c(0,-5,-10)
price <- c(20,10,5,0,-10)
tests  <- expand.grid(job,crop,prod,qual,yield,window,uptime,fuel,acost,price)
dim(tests)
tests  <- expand.grid(prod,qual,yield,window,uptime,fuel,acost,price)
dim(tests)
setwd("\\\\i-rsg.com\\projects\\Projects1\\Deere & Company\\Enterprise Advanced Marketing\\13259 - Electrification Ecospace\\9. Analysis\\Simulator")
write.table(tests,"tests.csv",sep=",")
is.null
is.null(NULL)
is.null(1)
library(RSGHB)
library(RSGHB)
checkModel
library(RSGHB)
?write.table
?rownames)
?rownames
library(RSGHB)
library(RSGHB)
library(RSGHB)
nrow
library(RSGHB)
stop
?sotp
?stop
library(RSGHB)
library(RSGHB)
?doHB
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
detach("package:RSGHB", unload=TRUE)
tools::vignetteEngine("knitr", weave = vweave, tangle = vtangle,
pattern = "[.]Rmd$", package = "knitr")
tools::vignetteEngine("knitr", weave = vweave, tangle = vtangle,
pattern = "[.]Rmd$", package = "knitr")
Sys.getenv("R_GSCMD")
Sys.getenv("R_GSCMD") <- "C:\\Program Files\\GhostScript"
Sys.getenv("R_GSCMD")
Sys.getenv("R_GSCMD")
library(RSGHB)
trans
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)
library(RSGHB)