#obs_ssKansas7.txt is the observed summary statistic  vector calculated  from Oaks et al. (Dryad doi:10.5061/dryad.5s07m)
# using the msbayes function 
#$obsSumStats.pl -T obs_ssKansas7.table batch_table1_prior.txt > obs_ssKansas7.txt
#obs_ssKansas7HEADLESS.txt is the same observed summary statistic  vector without the header and with a dummy first column


rm(list=ls(all=TRUE))
source("loc2plot_d.r")
source("make_pd2.r")
library(locfit)
library(nnet)
library(abc)
SIMS<- scan(("Posterior_M12345678_Z.ZZ"))
sims<-matrix(SIMS,ncol=467,byrow=T)
OBS<- scan(("obs_ssKansas7HEADLESS.txt"),nlines=1)
obs<-matrix(OBS,ncol=467,byrow=T)
Psi0.out <- postpr(obs[1,c(50:115,138:159)], sims[,2], sims[,c(50:115,138:159)], tol=1.0, method="mnlogistic", kernel="epanechnikov", sizenet = 5,maxit = 500) 
Psi0_out_summary <- summary(Psi0.out, unadj = FALSE, intvl = .95, print = TRUE, digits = max(3, getOption("digits")-3),)

#1000 below refers to the number of prior samples a particular posterior is sampling from (ie number of accepted prior simulations)

Z1noLL8 <- try(abc(obs[1,c(50:115,138:159)], sims[,5], sims[,c(50:115,138:159)], tol=1.0,method="loclinear"),T)
for (i in 1:1000){
if (Z1noLL8$adj.values[i,1]<0.0)
{Z1noLL8$adj.values[i,1]=0.0} else {Z1noLL8$adj.values[i,1]=Z1noLL8$adj.values[i,1]}
}
TDnoLLRmode8 <- loc1stats(Z1noLL8$adj.values, prob=0.95)[1]
TDnoLLRmode8raw <- loc1stats(sims[,5], prob=0.95)[1]

Z1noLL9 <- try(abc(obs[1,c(50:115,138:159)], sims[,4], sims[,c(50:115,138:159)], tol=1.0,method="loclinear"),T)

for (i in 1:1000){
if (Z1noLL9$adj.values[i,1]<0.0)
{Z1noLL9$adj.values[i,1]=0.0} else {Z1noLL9$adj.values[i,1]=Z1noLL9$adj.values[i,1]}
}
TDnoLLRmode9 <- loc1stats(Z1noLL9$adj.values, prob=0.95)[1]


Omega<-Z1noLL8$adj.values
Et<-Z1noLL9$adj.values


#local linear transformed 
postscript("E.t_Omega_cond_Mk.ps")
plot(Omega[,1],Et[,1],xlim=c(0,1.0),ylim=c(0,1.0),lty=2,lwd=0.5,xlab="V(t)/E(t)",ylab="E(t)")
dev.off()



#simple rejection 
postscript("E.t_Omega_cond_Mk_raw.ps")
plot(sims[,5],sims[,4],xlim=c(0,1.0),ylim=c(0,1.0),lty=2,lwd=0.5,xlab="V(t)/E(t)",ylab="E(t)")
dev.off()