.Rhistory

relist(skeleton=Model$skeleton) %>% unlist() %>% Proposed2Physical(Model)
par(mfrow=c(4,4))
# Plot S-curves for the actual and MAP parameterisation
D_extent <- BinR(Dfun(Intensity, theta=Omega$theta) , Omega$zeta)
D_extent_sample <- BinR(Dfun(Intensity, theta=Omega_curr$theta) , Omega_curr$zeta)
D_MortDisp <- plnorm( Dfun(Intensity, theta=Omega$theta), Omega$Lambda1$nu, Omega$Lambda1$omega)
D_MortDisp_sample <- plnorm( Dfun(Intensity, theta=Omega_curr$theta), Omega_curr$Lambda1$nu, Omega_curr$Lambda1$omega)
D_Mort <- plnorm( Dfun(Intensity, theta=Omega$theta), Omega$Lambda2$nu, Omega$Lambda2$omega)
D_Mort_sample <- plnorm( Dfun(Intensity, theta=Omega_curr$theta), Omega_curr$Lambda2$nu, Omega_curr$Lambda2$omega)
D_Disp <- D_MortDisp - D_Mort
D_Disp_sample <- D_MortDisp_sample - D_Mort_sample
D_BD <- plnorm( Dfun(Intensity, theta=Omega$theta), Omega$Lambda3$nu, Omega$Lambda3$omega)
D_BD_sample <- plnorm( Dfun(Intensity, theta=Omega_curr$theta), Omega_curr$Lambda3$nu, Omega_curr$Lambda3$omega)
plot(Intensity, D_Mort, col='red', type='l', ylim=c(0,1)); lines(Intensity, D_Mort_sample, col='red', lty=2)
lines(Intensity, D_Disp, col='blue'); lines(Intensity, D_Disp_sample, col='blue', lty=2)
lines(Intensity, D_BD, col='pink', type='l'); lines(Intensity, D_BD_sample, col='pink', lty=2, lwd=2)
lines(Intensity, D_extent, col='green', type='l'); lines(Intensity, D_extent_sample, col='green', lty=2, lwd=2)
for(i in 1:14){
ylim=c(min(unlist(Omega)[i], output[1:it,i+1]), max(unlist(Omega)[i], output[1:it,i+1]))
plot(output[1:it,i+1], type='l', ylab='', ylim=ylim)
abline(h=unlist(Omega)[i], col='red')
}
# Save log-target and parameters
saveRDS(output,paste0(dir,"IIDIPUS_Results/output_",tag))
# Save covariance matrix
saveRDS(propCOV,paste0(dir,"IIDIPUS_Results/covariance_",tag))
saveRDS(xbar_tminus1,paste0(dir,"IIDIPUS_Results/xbar_tminus1_",tag))
print(cov2cor(propCOV))
it <- it + 1
}
return(list(PhysicalValues=output[which.max(output[,1]),2:ncol(output)] %>%
relist(skeleton=Model$skeleton) %>% unlist() %>% Proposed2Physical(Model), # MAP value
OptimisationOut=output))
}
Algorithm <- match.fun('AMCMC2')
iVals<-GetInitVals(dir,Model,AlgoParams)
# Parameterise... Here we go!
AlgoParams$AllParallel <- TRUE
AlgoParams$cores <- 6
AlgoParams$Np <- 5
AlgoParams$ABC <- 1.5
AlgoParams$itermax <- 10000
output <- Algorithm(dir=dir,
Model=Model,
iVals=iVals,
AlgoParams=AlgoParams)
pi
library(posterior)
ess(c(5,6,7,8))
ess_basic(c(5,6,7,8))
ess_basic(c(5,6,7,8,5,4))
load("/home/manderso/Downloads/DispData_EQ_V2.Rdata")
dispdata <- readRDS("/home/manderso/Downloads/DispData_EQ_V2.Rdata")
dispdata
dispdata
setGeneric("DispX", function(ODD,Omega,center, BD_params, LL,Method, epsilon=c(0.15,0.03,0.1))
standardGeneric("DispX") )
# Code that calculates/predicts the total human displacement
setMethod("DispX", "ODD", function(ODD,Omega,center, BD_params, LL=F,
Method=list(Np=20,cores=8,cap=-300), epsilon=c(0.15,0.03,0.1)
){
# Extract 0D parameters & speed up loop
Params<-FormParams(ODD,list(Np=Method$Np,center=center))
# Income distribution percentiles & extract income percentile
SincN<-seq(10,90,by = 10); Sinc<-ExtractCIndy(ODD,var = paste0("p",SincN,"p100"))
# Speed-up calculation (through accurate cpu-work distribution) to only values that are not NA
if(!LL) {notnans<-which(!(is.na(ODD$Population) | is.na(ODD$ISO3C) | is.na(ODD$GDP)))
} else notnans<-which(!(is.na(ODD$Population) | is.na(ODD$ISO3C) | is.na(ODD$GDP) |
!ODD$ISO3C%in%ODD@gmax$iso3))
BD_data_present <- ifelse(all(is.na(ODD@data$nBuildings)) , F, T)
# Calculate non-local linear predictor values
LP<-GetLP(ODD,Omega,Params,Sinc,notnans)
# Speed things up a little
hrange<-grep("hazMean",names(ODD),value = T)
# Function to predict displacement per gridpoint
CalcDam<-function(ij){
iso3c<-ODD@data$ISO3C[ij]
# Calculate local linear predictor (NOTE: is a vector due to income distribution)
locallinp<-LP$dGDP$linp[LP$dGDP$ind==LP$iGDP[ij]]*LP$Plinp[ij]*LP$linp[[iso3c]] #LOOSEEND
#locallinp<-rep(1,10) #reducing parameter space while I'm figuring out the MCMC
# Sample population per income distribution (Assumes 9 percentiles):
lPopS <- SplitSamplePop(Pop=ODD@data$Population[ij],Method$Np)
tPop <-array(0,c(3, Method$Np)) #row 1 = tDisp, #row 2 = tMort, #row 3 = tRem
tPop[3,]=colSums(lPopS)
for(h in hrange){
# for(h in c(1)){
if(is.na(ODD@data[ij,h])) next
# Resample population based on who is remaining
ind<-tPop[3,]>0
if(h!=hrange[1]) {
if(sum(ind)==0) break #if no remaining population, skip modelling
if(length(lPopS[,ind])==0) break #if no remaining population, skip modelling
#if(sum(ind)>1) sumz<-colSums(lPopS[,ind])
#else sumz<-sum(lPopS[,ind])
#lPopS[,!ind]<-0
lPopS[,ind]<-SplitSamplePop(Pop=tPop[3,ind])
}
# Sample hazard Intensity
# the uncertainty is too high... so I scale it to get some interpretable results (I know, I'm not really a statistician, I don't even have a degree, I was actually just having a look around the department when they confused me for the interviewee. I didn't have the heart to say anything. You don't hate me as much as I do)
# I_ij<-rnorm(n = Method$Np,
#             mean = ODD@data[ij,paste0("hazMean",h)],
#             sd = ODD@data[ij,paste0("hazSD",h)]/10)
I_ij<-ODD@data[ij,h]
# Separate into income distributions (as each have 10% of population, order doesn't matter)
for (s in 1:length(SincN)){
if(all(lPopS[s,]==0)) next
# Predict damage at coordinate {i,j} (vector with MC particles)
Damage <-tryCatch(fDamUnscaled(I_ij,list(I0=Params$I0, Np=Params$Np),Omega)*locallinp[s], error=function(e) NA)
if(any(is.na(Damage))) print(ij)
# Scaled damage
D_MortDisp <- plnorm(Damage, meanlog = Omega$Lambda1$nu, sdlog = Omega$Lambda1$omega)
D_Mort <- plnorm(Damage, meanlog = Omega$Lambda2$nu, sdlog = Omega$Lambda2$omega)
D_Disp<-D_MortDisp - D_Mort
D_Disp <- ifelse(D_Disp<0, 0, D_Disp)
# Accumulate the number of people displaced/deceased, but don't accumulate the remaining population
tPop[3,ind]<-0
tPop[,ind]<-tPop[,ind] + Fbdam(lPopS[s,ind],D_Disp[ind], D_Mort[ind], (1-D_Mort-D_Disp)[ind])
}
}
#ensure the total displaced, deceased or remaining does not exceed total population
tPop[tPop>ODD@data$Population[ij]]<-floor(ODD@data$Population[ij])
#if no building destruction data:
if(!BD_data_present) return(rbind(tPop[1:2,, drop=FALSE], rep(NA, Method$Np))) #return total displacement and mortality, set number of buildings destroyed to NA
#otherwise, sample from the model for the number of buildings destroyed:
#we take locallinp[5] which corresponds to locallinp for the median GDP
nBuildings = rep(ODD@data$nBuildings[ij], Method$Np)
nBD = rep(0, Method$Np)
for (h in hrange){
if(is.na(ODD@data[ij,h])) next
if(h!=hrange[1]) {
if(all(nBuildings==0)) break #if no remaining buildings, skip modelling LOOSEEND
}
I_ij<-ODD@data[ij,h]
Damage <-tryCatch(fDamUnscaled(I_ij,list(I0=Params$I0, Np=Params$Np),Omega)*locallinp[5], error=function(e) NA) #calculate unscaled damage (excluding GDP)
D_BD = plnorm(Damage, meanlog = Omega$Lambda3$nu, sdlog = Omega$Lambda3$omega)
moreBD = fBD(nBuildings, D_BD)
nBD = nBD + moreBD
nBuildings = nBuildings - moreBD
}
return(rbind(tPop[1:2,,drop=FALSE], nBD))
}
Dam<-array(0,c(nrow(ODD),Method$Np,3)) # Dam[,,1] = Displacement, Dam[,,2] = Mortality, Dam[,,3] = Buildings Destroyed
#if(Method$cores>1) { Dam[notnans,,]<-aperm(simplify2array(mclapply(X = notnans,FUN = CalcDam,mc.cores = Method$cores)), perm=c(3,2,1))
#} else  Dam[notnans,,]<- aperm(simplify2array(lapply(X = notnans,FUN = CalcDam)), perm=c(3,2,1))
Dam[notnans,,]<- aperm(simplify2array(lapply(X = notnans,FUN = CalcDam)), perm=c(3,2,1))
# return(Disp)
# If the IDMC estimate is foreseen to be a lower or upper bound, or a generally poor estimate
# for(c in ODD@gmax$iso3){
#   ind<-ODD@data$ISO3C==c  & !is.na(ODD@data$ISO3C)
#   Disp[ind,]%<>%qualifierDisp(qualifier = ODD@gmax$qualifier[ODD@gmax$iso3==c],mu = Omega$mu)
# }
funcy<-function(i,LLout=T, epsilon=AlgoParams$epsilon_min) {
tmp<-data.frame(iso3=ODD$ISO3C,IDPs=Dam[,i,1], mort=Dam[,i,2], nBD=Dam[,i,3]) %>%
group_by(iso3) %>% summarise(disp_predictor=floor(sum(IDPs,na.rm = T)),
mort_predictor=floor(sum(mort,na.rm = T)),
nBD_predictor=floor(sum(nBD,na.rm = T)),
.groups = 'drop_last')
tmp<-tmp[!is.na(tmp$iso3) & tmp$iso3%in%ODD@gmax$iso3,]
tmp%<>%merge(ODD@gmax,by="iso3")%>%arrange(desc(gmax))
#print(paste(tmp$nBD_predictor, tmp$buildDestroyed))
#print(tmp)
if(LLout) {
return(LL_IDP(tmp, epsilon))
}
return(tmp)
}
outer<-vapply(1:Method$Np,funcy,numeric(length(unique(ODD@gmax$iso3))), epsilon=epsilon)
#outer[outer<Method$cap]<-Method$cap
# Find the best fit solution
if(length(unique(ODD@gmax$iso3))>1) {
if(LL)  return(log(rowMeans(exp(outer),na.rm=T)))
MLE<-which.max(log(colSums(exp(outer),na.rm=T)))
}  else {
if(LL)  return(log(mean(exp(outer),na.rm=T)))
MLE<-which.max(log(exp(outer)))
}
if(Method$Np == 1){
MLE=1
}
# Save into ODD object
# ODD@data$Disp<-Disp[,MLE]*sum(ODD@gmax$gmax)/mean(sum(Disp[,MLE])) %>% round()
ODD@data$Disp<-Dam[,MLE,1]  #should this be named data$DispPred or something instead?
ODD@data$Mort<-Dam[,MLE,2]
ODD@data$nBD<-Dam[,MLE,3]
# I know this is kind of repeating code, but I want the other function as fast as possible
ODD@predictDisp<-funcy(MLE,LLout=F)
return(ODD)
})
unlist(c(2,3,3,2,333))
install.packages('coxed')
library(coxed)
simdata <- sim.survdata(N=20, T=20, num.data.frames=1)
simdata
?sim.survdata
res_paramW = get_param_weib(med = 1062, mu = 1134)
listCoxSim_n500_p1000 <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 500,
p = 1000,
pnonull = 20,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 0)
install.packages('surviveMS')
install.packages('survMS')
library(survMS)
res_paramW = get_param_weib(med = 1062, mu = 1134)
listCoxSim_n500_p1000 <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 500,
p = 1000,
pnonull = 20,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 0)
listCoxSim_n500_p1000
?modelSim
listCoxSim <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 50,
p = 3,
pnonull = 3,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 3)
listCoxSim
listCoxSim <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 50,
p = 3,
pnonull = 3,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 20)
listCoxSim
names(listCoxSim)
listCoxSim$Y
listCoxSim$Z
listCoxSim$model
listCoxSim$betaNorm
listCoxSim$crate_delta
listCoxSim$vecY
listCoxSim$Z
listCoxSim$Y
listCoxSim$delta
listCoxSim$crate
listCoxSim$crate_dela
listCoxSim$crate_delta
listCoxSim$hazParams
listCoxSim$hazDistr
listCoxSim$TC
listCoxSim$Z
listCoxSim$betaNorm
listCoxSim$betaNorm %*% listCoxSim$Z
listCoxSim$Z  %*% listCoxSim$betaNorm
exp(listCoxSim$Z  %*% listCoxSim$betaNorm)
hist(listCoxSim$TC)
hist(listCoxSim$TC, breaks=50)
hist(listCoxSim$TC, breaks=20)
hist(listCoxSim$TC, breaks=15)
plot(1:50, listCoxSim$TC)
plot(1:50, listCoxSim$TC, col=listCoxSim$delta)
plot(1:50, listCoxSim$TC, col=)
plot(1:50, listCoxSim$TC, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
listCoxSim <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 50,
p = 3,
pnonull = 3,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 200)
plot(1:50, listCoxSim$TC, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
listCoxSim <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 50,
p = 3,
pnonull = 3,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 0)
plot(1:50, listCoxSim$TC, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
listCoxSim$delta
listCoxSim <- modelSim(model = "cox",
matDistr = "unif",
matParam = c(-1,1),
n = 50,
p = 3,
pnonull = 3,
betaDistr = 1,
hazDistr = "weibull",
hazParams = c(res_paramW$a, res_paramW$lambda),
seed = 1,
d = 1500)
plot(1:50, listCoxSim$TC, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
listCoxSim$delta # censorship indicator
n_cov <- 3
n <- 50
beta <- c(1,2,-1)
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
covariates
beta %*% covariates
exp(covariates %*% beta)
(- log(1- runif(n, 0,1)) /exp(covariates %*% beta))
n_cov <- 3
n <- 50
beta <- c(1,2,-1)
a <- 1
lambda <- 10
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
h_t <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
T <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
Times <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
Times
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
n_cov <- 3
n <- 50
beta <- c(1,2,-1)
a <- 0.5
lambda <- 10
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
n_cov <- 3
n <- 50
beta <- c(1,2,-1)
a <- 2
lambda <- 0.1
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
Times <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
Times
plot(1:50, ST, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
plot(1:50, ST, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
hist(ST)
n_cov <- 3
n <- 50
beta <- c(1,2,-1)
a <- 2
lambda <- 0.1
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
plot(1:50, ST, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
hist(ST)
plot(1:50, ST, col=ifelse(listCoxSim$delta==T, 'red', 'black'))
hist(ST, breaks=15)
censor_times <- runif(n, 0, 20)
ST_censored <- ifelse(ST>censor_times, ST, censor_times)
ST_censored
plot(ST_censored)
points(censor_times, col='red')
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 10)
ST_censored <- ifelse(ST>censor_times, ST, censor_times)
plot(ST_censored)
points(censor_times, col='red')
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 30)
ST_censored <- ifelse(ST>censor_times, ST, censor_times)
plot(ST_censored)
points(censor_times, col='red')
plot(ST)
ST>censor_times
which(ST_censored == censor_times)
ST
ifelse(ST>censor_times, ST, censor_times)
ST>censor_times
ST
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 30)
ST_censored <- array(n)
for (ev in 1:n){
ST_censored[i] <- ifelse(ST[i]>censor_times[i], ST[i], censor_times[i])
}
plot(ST_censored)
points(censor_times, col='red')
_cov <- 3
n <- 50
beta <- c(1,2,-1)
a <- 2
lambda <- 0.1
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 30)
ST_censored <- array(n)
for (ev in 1:n){
ST_censored[i] <- ifelse(ST[i]>censor_times[i], ST[i], censor_times[i])
}
plot(ST_censored)
n_cov <- 3
n <- 50
beta <- c(1,2,-1)
a <- 2
lambda <- 0.1
covariates <- array(runif(n*n_cov,-1,1), dim=c(n, n_cov))
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 30)
ST_censored <- array(n)
for (ev in 1:n){
ST_censored[ev] <- ifelse(ST[ev]>censor_times[ev], ST[ev], censor_times[ev])
}
plot(ST_censored)
points(censor_times, col='red')
n = 1
n = 50
ev = 1
ST[ev]
censor_times[ev]
ST[ev]>censor_times[ev]
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 30)
ST_censored <- ifelse(ST<censor_times, ST, censor_times)
plot(ST_censored)
points(censor_times, col='red')
plot(seq(0,10,0.1), lambda * seq(0,10,0.1)^a)
# H_0(t) = lambda * t^a
ST <- ((- log(1- runif(n, 0,1)) /exp(covariates %*% beta)) / lambda)^(1/a)
censor_times <- runif(n, 0, 20)
ST_censored <- ifelse(ST<censor_times, ST, censor_times)
plot(ST_censored)
points(censor_times, col='red')
censor_flag <- ifelse(censor_times==ST_censored, T, F)
censor_flag
censor_times==ST_censored
censor_flag <- censor_times==ST_censored
censor_flag
install.packages(c("survival", "survminer"))
library("survival"); library("survminer")
install.packages(c("survival", "survminer"))
library("survival"); library("survminer")
?coxph
coxph(ST ~ covariates[,1])
ST
Surv(time=ST, status=censor_flag)
?Surv
Surv(time=ST)
coxph(Surv(time=ST) ~ covariates[,1])
beta
coxph(Surv(time=ST) ~ covariates[,1] + covariates[,2] + covariates[,3])
coxph(Surv(time=ST, status=ifelse(ST_censored, 0, 1)) ~ covariates[,1] + covariates[,2] + covariates[,3])
coxph(Surv(ST, ifelse(ST_censored, 0, 1)) ~ covariates[,1] + covariates[,2] + covariates[,3])
coxph(Surv(ST, event=ifelse(ST_censored, 0, 1)) ~ covariates[,1] + covariates[,2] + covariates[,3])
coxph(Surv(time=ST, event=ifelse(ST_censored, 0, 1)) ~ covariates[,1] + covariates[,2] + covariates[,3])
coxph(Surv(time=ST) ~ covariates[,1] + covariates[,2] + covariates[,3])
Inf < Inf
if (5>4) print(3)
if (5<4) print(3)
AlgoParams$smc_steps
247+579
178+166
# Where is the main folder with all the code and data
dir<-directory<-"/home/manderso/Documents/GitHub/ODDRIN/"
# Set the working directory from your environment variables
setwd(directory)
# Directory of the Data for Good data, e.g. Disaster Mapping, 4G connectivity, etc
FBdirectory<-'/home/patten/Documents/IDMC/Facebook_Data/'
# Do you want only the reduced packages or all? Choose via packred
packred<-F
source('RCode/GetEnv.R')
# Download and install the necessary packages:
source('RCode/GetODDPackages.R')
# Sourcing the data:
source('RCode/GetData.R')
# Extract model functions and priors
source('RCode/Model.R')
# Extract the model parameterisation algorithm, default = Adaptive MCMC
source('RCode/Method.R')
#Extract the functions for generating the simulations
source('RCode/Simulate.R')
#Parameterise the model and simulate the data:
Omega <- list(Lambda1 = list(nu=8,omega=4.9),
Lambda2 = list(nu= 9.2, omega=5),
Lambda3 = list(nu=7.5,omega=4.1),
Lambda4 = list(nu=8.7, omega=5),
Pdens = list(M=0.01988616, k = 6.473428),
dollar = list(M = -0.41271, k = 6.473428),
theta = list(e=0.2359788),
eps = list(eps=0.01304351))
Model$center <- simulateDataSet(100, Omega %>% addTransfParams(), Model=Model, dir = dir, outliers = FALSE)
Model$center <- simulateDataSet(100, Omega %>% addTransfParams(), Model=Model, dir = dir, outliers = FALSE)
warnings()