################################################################################################################ ###### Determinisitic Version of Lotka Volterra Model ###### ################################################################################################################ ### Attempt to assess internal and external model discrepancy using determinisitc model and match to real world data ### that has been simulated from the more advanced stochastic version of the model, with unknown ICs and time varying ### rate parameters. Observed data saved in "Possible_Observed_Data_V1.Rdata", generated from "Stochastic_LotVol_Real_World.R" library(deSolve) PDF <- 1 # 1 if plotting PDFs, 0 otherwise ### Define the Lotka Volterra Predator-Prey model Lotka_Volterra <- function(t, y, parms) { with(as.list(parms),{ dPrey = th1 * y["Prey"] - th2 * y["Prey"] * y["Pred"] dPred = th2 * y["Prey"] * y["Pred"] - th3 * y["Pred"] res <- c(dPrey, dPred) list(res) }) } ### set rate constants, output times and initial conditions ### parms = c( th1 = 1, th2 = 0.00044, th3 = 1.8 ) ### Start values for output y containing prey and pred numbers ### ystart <- c(Prey = 2000, Pred = 800 ) ### Define time period of interest ### times <- seq(0, 10, length=1001) ### run the lsoda solver ### out <- as.matrix(lsoda(ystart, times, Lotka_Volterra, parms, maxsteps=20000)) ## solve the differential eqtn if(PDF) pdf("plots_health/A1_LV_single_run.pdf",height=6,width=8) plot(out[,"time"],out[,"Prey"],lwd=2,col=6,ty="l",ylim=c(min(out[,-1]),max(out[,-1])),log="",xlab="time(t)",ylab="Pred/Prey Number") lines(out[,"time"],out[,"Pred"],lwd=2,col="orange") legend("topright",legend=c("prey","pred"),col=c(6,"orange"),lty=1,lwd=2) if(PDF) dev.off() ### load observed data from Stochastic_LotVol_Real_World.R" ### #load("Observed Data/Possible_Observed_Data_V1.Rdata") # plot_pred_prey <- function(out){ # par(mfrow=c(2,1),mar=c(4,4,1,1),cex=0.75) # plot(out[,"time"],out[,"Prey"],lwd=2,col=6,ty="l",ylim=c(min(out[,-1]),max(out[,-1]))*1.1,log="",xlab="time(t)",ylab="Pred/Prey Number") # points(t_obs,y_obs[,"Prey"],xlim=c(0,10),ylim=c(0,max(y_obs)),pch=16,col=4,cex=0.9) # arrows(t_obs,y_obs[,"Prey"] + 3*obs_sd,t_obs,y_obs[,"Prey"] - 3*obs_sd,angle=90,code=3,col=4,length=0.04,lwd=2.5) # plot(out[,"time"],out[,"Pred"],lwd=2,col="orange",ty="l",ylim=c(min(out[,-1]),max(out[,-1]))*1.1,log="",xlab="time(t)",ylab="Pred/Prey Number") # points(t_obs,y_obs[,"Pred"],xlim=c(0,10),ylim=c(0,max(y_obs)),pch=16,col=2,cex=0.9) # arrows(t_obs,y_obs[,"Pred"] + 3*obs_sd,t_obs,y_obs[,"Pred"] - 3*obs_sd,angle=90,code=3,col=2,length=0.04,lwd=2.5) # } # # if(PDF) pdf("plots_health/B1_LV_single_run_with_obs_data.pdf",height=8,width=8) # plot_pred_prey(out=out) # if(PDF) dev.off() ### plot repetition runs for next few sections ### plot_pred_prey_array <- function(out,out_wx=NULL,lw1=0.8,lw2=1.6,int_md=0,plotprey=1,plotpred=1,xylog="",set_ylim_obs=NULL,good_id=NULL,time_cut=5,obs_col=1,wv1_col=c(6,"orange")){ nreps <- dim(out)[3] if(is.null(set_ylim_obs)) ylim1 <- c(min(out[,-1,]),max(out[,-1,])) else ylim1 <- c(min(y_obs)/set_ylim_obs,max(y_obs)*set_ylim_obs) if(plotprey==1 & plotpred==1) par(mfrow=c(2,1),mar=c(4,4,1,1),cex=0.75) if(plotprey==1){ plot(out[,"time",1],out[,"Prey",1],lwd=2,col=wv1_col[1],ty="n",ylim=ylim1,xlab="time(t)",ylab="Prey Number",log=xylog) for(i in 1:nreps) lines(out[,"time",i],out[,"Prey",i],lwd=lw1,col=wv1_col[1]) if(!is.null(good_id)) for(i in good_id) lines(out[,"time",i],out[,"Prey",i],lwd=lw2,col="green") if(!is.null(out_wx)) for(i in 1:dim(out_wx)[3]) lines(out_wx[,"time",i],out_wx[,"Prey",i],lwd=lw2,col="purple") points(t_obs,y_obs[,"Prey"],xlim=c(0,10),ylim=c(0,max(y_obs)),pch=16,col=obs_col,cex=0.9) if(int_md==0) arrows(t_obs,y_obs[,"Prey"] + 3*obs_sd,t_obs,y_obs[,"Prey"] - 3*obs_sd,angle=90,code=3,col=obs_col,length=0.04,lwd=2.5) if(int_md==1) arrows(t_obs,y_obs[,"Prey"] + 3*obs_int_md_prey,t_obs,y_obs[,"Prey"] - 3*obs_int_md_prey, angle=90,code=3,col=obs_col,length=0.04,lwd=2.5) if(is.null(good_id)) legend("topleft",legend=c("Wave 1 runs","Observed Data"),col=c(6,1),lty=c(1,NA),pch=c(NA,21),lwd=1.6) else legend("topleft",legend=c("Wave 1 runs","Wave 2 runs","Wave 3 runs","Observed Data"),col=c(wv1_col[1],"green","purple",1),lty=c(1,1,1,NA),pt.bg=c(NA,NA,NA,obs_col),pch=c(NA,NA,NA,21),lwd=1.6) if(!is.null(time_cut)) abline(v=time_cut,lty=2) } if(plotpred==1){ plot(out[,"time",1],out[,"Pred",1],lwd=2,col=wv1_col[2],ty="n",ylim=ylim1,xlab="time(t)",ylab="Predator Number",log=xylog) for(i in 1:nreps) lines(out[,"time",i],out[,"Pred",i],lwd=lw1,col=wv1_col[2]) if(!is.null(good_id)) for(i in good_id) lines(out[,"time",i],out[,"Pred",i],lwd=lw2,col="green") if(!is.null(out_wx)) for(i in 1:dim(out_wx)[3]) lines(out_wx[,"time",i],out_wx[,"Pred",i],lwd=lw2,col="purple") points(t_obs,y_obs[,"Pred"],xlim=c(0,10),ylim=c(0,max(y_obs)),pch=16,col=obs_col,cex=0.9) if(int_md==0) arrows(t_obs,y_obs[,"Pred"] + 3*obs_sd,t_obs,y_obs[,"Pred"] - 3*obs_sd,angle=90,code=3,col=obs_col,length=0.04,lwd=2.5) if(int_md==1) arrows(t_obs,y_obs[,"Pred"] + 3*obs_int_md_pred,t_obs,y_obs[,"Pred"] - 3*obs_int_md_pred, angle=90,code=3,col=obs_col,length=0.04,lwd=2.5) if(is.null(good_id)) legend("topleft",legend=c("Wave 1 runs","Observed Data"),col=c(wv1_col[2],1),lty=c(1,NA),pch=c(NA,21),lwd=1.6) else legend("topleft",legend=c("Wave 1 runs","Wave 2 runs","Wave 3 runs","Observed Data"),col=c(wv1_col[2],"green","purple",1),lty=c(1,1,1,NA),pt.bg=c(NA,NA,NA,obs_col),pch=c(NA,NA,NA,21),lwd=1.6) if(!is.null(time_cut)) abline(v=time_cut,lty=2) } } ###################################################################################### ### Select set of wave 2 runs to use for Model Discrepancy Assessment ### library(tgp) ### generate Latin hypercube design for wave 1 runs ### parms = c( th1 = 1, th2 = 0.00044, th3 = 1.8 ) nruns <- 750 set.seed(6) x1 <- lhs(nruns,matrix(c(parms*0.7,parms*1.3),3,2) ) colnames(x1) <- paste("th",1:3,sep="") ### run the lsoda solver i.e. generate wave 1 runs from LV model ### out_w1 <- array(0,c(length(times),3,nruns),dimnames=list( NULL,c("time","Prey","Pred"), NULL)) for(i in 1:nruns){ out_w1[,,i] <- as.matrix(lsoda(ystart, times, Lotka_Volterra, parms=x1[i,], maxsteps=20000)) # solve the diff eqtn cat(i,"\n") } save(out_w1,x1,nruns,file="runs/Wave1_runs.Rdata") # load("runs/Wave1_runs.Rdata") ### Select good runs based on max implausibility ### # load("Observed Data/Possible_Observed_Data_V1.Rdata") al_obs <- match(t_obs,times) # indexes of observation data locations in time prey_w1 <- out_w1[al_obs,2,] # prey output at obs time pts pred_w1 <- out_w1[al_obs,3,] # pred output at obs time pts MD_approx_w1 <- 0.15*y_obs # 15 percent cautious initial wave 1 MD assessment prey_imp <- abs(prey_w1 - y_obs[,"Prey"])/sqrt(obs_sd^2 + MD_approx_w1[,"Prey"]^2) # implausibility for prey pred_imp <- abs(pred_w1 - y_obs[,"Pred"])/sqrt(obs_sd^2 + MD_approx_w1[,"Pred"]^2) # implausibility for pred # find max implausibility from t=0 to t=5 secs only (only using first 10 outputs from each of prey and pred) max_imp <- apply(rbind(prey_imp[1:10,],pred_imp[1:10,]), 2, max) ngood <- sum(max_imp<3) # impose Imp < 3 cut ngood # number of "acceptable" runs (wave 2 green runs in Figure 2) al_good <- max_imp < max_imp[order(max_imp)][ngood+1] # indexes of good runs x1_good <- x1[al_good,] # x inputs of good wave 2 green runs out_w1_good <- out_w1[,,al_good] # f outputs of good wave 2 green runs x1_good_best_ind <- which.min(max_imp[al_good]) # index of best run so far (wave 2 test pt: red pt in Fig 2) ### XXX plot wave 1 and 2 runs so far (this plot is superseded by later plot with wave 3) if(PDF) png("plots_health/J1_LV_wave1and2_runs_pred_prey_outputs.png",width=7.5*0.8,height=8*0.8,pointsize=12,res=72*2+40,units="in") plot_pred_prey_array(out=out_w1,lw1=0.6,xylog="",good_id=which(al_good==1),wv1_col=c("grey60","grey60")) if(PDF) dev.off() ### XXX plot wave 1 and 2 run input locations so far (this plot is superseded by later plot with wave 3) if(PDF) pdf("plots_health/K1_LV_wave2_nimp_inputs.pdf",height=8*0.8,width=8.5*0.8) xran <- rbind(parms*0.7,parms*1.3) lab1 <- c(expression(x[1]),expression(x[2]),expression(x[3])) col1 <- rep("green",ngood) col1[x1_good_best_ind] <- "red" pairs(rbind(x1,x1_good,xran),pch=c(rep(16,nruns),rep(21,ngood),16,16),col=c(rep("grey80",nruns),rep(1,ngood),0,0),bg=c(rep("grey80",nruns),rep(col1,ngood),0,0),labels=lab1) if(PDF) dev.off() save(out_w1_good,x1_good,ngood,x1_good_best_ind,file="runs/Wave1_goodruns.Rdata") # load("runs/Wave1_goodruns.Rdata") ### end select set of OK wave 2 runs to use for Model Discrepancy Assessment ### ###################################################################################### ########################################################################## ### calculate internal model discrepancy covariance matrix due to ICs ### ### run the lsoda solver i.e. the LV model for each perturbation ### set.seed(1) nreps <- 50 # no. of perturbations d_i out_ic <- array(0,c(length(times),3,nreps),dimnames=list( NULL,c("time","Prey","Pred"), NULL)) for(i in 1:nreps){ ystart <- c(Prey = rnorm(1,2000,40),Pred = rnorm(1,800,40)) # random initial conditions out_ic[,,i] <- as.matrix(lsoda(ystart, times, Lotka_Volterra, parms=x1_good[x1_good_best_ind,], maxsteps=20000)) # solve the differential eqtn cat(i,"\n") } ### choose which times of interest (choose all times here for plotting purposes) ### al_obs <- 1:length(times) out_at_obs_ic <- rbind(out_ic[al_obs,2,],out_ic[al_obs,3,]) ### calculate cor and cov matrices of Internal MD induced by the IC uncertainty ### md_cor_ic <- cor(t(out_at_obs_ic)) md_cov_ic <- cov(t(out_at_obs_ic)) ### end calculate model discrepancy covariance matrix due to ICs ### ######################################################################## ###################################################################################### ### calculate internal model discrepancy covariance matrix due to time varying constants ### library(MASS) set.seed(2) ### draw a few parameters for the time varying function g(t) ### a <- runif(1,min=5-0.5,max=5+0.5) b <- runif(1,min=0.6,max=0.9) c <- runif(15,0.02,0.05) ### Define function (g(t) in paper) to give time varying rate inputs ### func_t <- function(t){ 1 + c*cos( 2*pi*(t-a)/max(times) * b) } ### The time varying function used for the real system was a quadratic: ### func_t_reality <- function(t){ 1 - ((t-5)^2/25 - 0.5) * 2 * 0.05 } ### plot draws of time varying function: deliberately different from real system ### f_t_mat <- matrix(0,length(times),length(c)) for(i in 1:length(times)) f_t_mat[i,] <- func_t(times[i]) plot(times,f_t_mat[,1],ty="n",ylim=range(func_t_reality(times),f_t_mat) ) for(j in 1:length(c)) lines(times,f_t_mat[,j]) lines(times,func_t_reality(times),col=2) ### Define the Lotka Volterra Predator-Prey model with time varying rate constants ### Lotka_Volterra_tvar <- function(t, y, parms) { with(as.list(parms),{ dPrey = func_t(t)[1] * th1 * y["Prey"] - func_t(t)[2] * th2 * y["Prey"] * y["Pred"] dPred = func_t(t)[2] * th2 * y["Prey"] * y["Pred"] - func_t(t)[3] * th3 * y["Pred"] res <- c(dPrey, dPred) list(res) }) } ### run the lsoda solver i.e. the LV model for each perturbation ### nreps <- 50 out_vx <- array(0,c(length(times),3,nreps),dimnames=list( NULL,c("time","Prey","Pred"), NULL)) for(i in 1:nreps){ ### draw different time varying parameters of g(t) = func_t function ### a <- runif(1,min=5-0.5,max=5+0.5) b <- runif(1,min=0.6,max=0.9) c <- rep(runif(1,0.02,0.05),3) out_vx[,,i] <- as.matrix(lsoda(ystart, times, Lotka_Volterra_tvar, parms=x1_good[x1_good_best_ind,], maxsteps=20000)) # solve the diff eqtn cat(i,"\n") } ### choose which times of interest (choose all times here for plotting purposes) ### al_obs <- 1:length(times) out_at_obs_vx <- rbind(out_vx[al_obs,2,],out_vx[al_obs,3,]) ### calculate cor and cov matrices of Internal MD induced by the time varying inputs ### md_cor_vx <- cor(t(out_at_obs_vx)) md_cov_vx <- cov(t(out_at_obs_vx)) ### end calculate model discrepancy covariance matrix due time varying inputs ### ################################################################################# ############################################################################################### ### Now plot both IC (initial condition) and vx (time varying x inputs) components together ### if(PDF) pdf("plots_health/I1_LV_single_run_MD_contributions.pdf",height=8,width=8.5) par(mfrow=c(2,1),mar=c(4,4,3,3),cex=0.75) lwd1 <- 1.4 pp_ind <- rbind(1:1001,1002:2002) o_ind <- rbind(seq(50,1000,50),seq(1050,2000,50)) for(j in 1:2){ md_cov_tot_approx <- md_cov_vx + md_cov_ic # total internal discrepancy plot(seq(0,10,len=1001),sqrt(diag(md_cov_vx))[pp_ind[j,]],ty="l",col=2,ylim=range(sqrt(c(diag(md_cov_tot_approx)[pp_ind[j,]],diag(md_cov_vx)))), lwd=lwd1,ylab="SD of MD contribution",xlab="Time",main=paste("Contributions to MD:",c("Prey","Pred")[j])) lines(seq(0,10,len=1001),sqrt(diag(md_cov_ic))[pp_ind[j,]],ty="l",col=4,lwd=lwd1) lines(seq(0,10,len=1001),sqrt(diag(md_cov_tot_approx))[pp_ind[j,]],ty="l",col=1,lwd=lwd1) lines(seq(0,10,len=1001),out[,j+1]*0.15,lty=2) # conservative 15% MD lines(seq(0,10,len=1001),sqrt(diag(md_cov_tot_approx)[pp_ind[j,]]+(out[,j+1]*0.02)^2),ty="l",col=5,lwd=lwd1) lines(seq(0,10,len=1001),out[,j+1]*0.02,col=3,lwd=lwd1) # nominal 2% external MD legend("topleft",legend=c("Initial Condition Int MD","Time-varying x Int MD","Total Int MD","Ext MD","Total MD","Conservative MD"),col=c(4,2,1,3,5,1),lty=c(rep(1,5),2),cex=1.1,lwd=lwd1) } if(PDF) dev.off() ### end calculate model discrepancy covariance matrix due to time varying constant ### ########################################################################################## ########################################################################################################### ### calculate full internal model discrepancy covariance matrix due to ICS and varying x's over time ### library(MASS) # load("runs/Wave1_goodruns.Rdata") ### run the lsoda solver i.e. the LV model for each perturbation, now perturbing both IC and time varying x ### set.seed(1) nreps <- 200 out_full <- array(0,c(length(times),3,nreps,ngood),dimnames=list( NULL,c("time","Prey","Pred"), NULL,NULL)) for(j in 1:ngood){ for(i in 1:nreps){ ystart <- c(Prey = rnorm(1,2000,40),Pred = rnorm(1,800,40)) # random initial conditions a <- runif(1,min=5-0.5,max=5+0.5) # varying inputs x over time parameters b <- runif(1,min=0.6,max=0.9) c <- rep(runif(1,0.02,0.05),3) out_full[,,i,j] <- as.matrix(lsoda(ystart, times, Lotka_Volterra_tvar, parms=x1_good[j,], maxsteps=20000)) cat(i," ",j,"\n") } } save(out_full,nreps,ngood,file=paste("runs/Model_Disc_Assess_full_array_",nreps,"_reps.Rdata",sep="")) load(paste("runs/Model_Disc_Assess_full_array_",nreps,"_reps.Rdata",sep="")) ### find full mod dis cov matrix for each of the ngood runs ### al_obs <- 1:length(times) ### set up matrices to store Inter MD SD and Biases in ### md_full <- matrix(0,2*length(al_obs),ngood) bias_full_prey <- matrix(0,length(al_obs),ngood) bias_full_pred <- matrix(0,length(al_obs),ngood) ### loop over wave 2 acceptable runs and calculate Int MD SD and biases ### for(j in 1:ngood){ out_at_obs_full <- rbind(out_full[al_obs,2,,j],out_full[al_obs,3,,j]) # SDs (prey and pred stacked) md_cov_full <- cov(t(out_at_obs_full)) md_full[,j] <- sqrt(diag(md_cov_full)) # Biases bias_full_prey[,j] <- out_w1_good[,2,j] - apply(out_full[al_obs,2,,j],1,mean) bias_full_pred[,j] <- out_w1_good[,3,j] - apply(out_full[al_obs,3,,j],1,mean) cat(j,"\n") } ### plot full correlation structure of single good wave 2 test point ### j <- x1_good_best_ind out_at_obs_full <- rbind(out_full[al_obs,2,,j],out_full[al_obs,3,,j]) md_cor_full <- cor(t(out_at_obs_full)) if(PDF) png(file="plots_health/M0_LV_single_run_combined_IC_VX_CorMat.png",type="cairo",width=1750,height=1600,res=200) filled.contour(md_cor_full[seq(1,2*length(times),10),rev(seq(1,2*length(times),10))],main="Correlation Structure of Internal MD: Prey and Pred.", plot.axes = { axis(1, at=seq(0, 1, len = 11), labels=c(0,rep(seq(2,10,2),2)) ) axis(2, at=seq(0, 1, len = 11), labels=rev(c(0,rep(seq(2,10,2),2))) ) abline(h=0.5-0.5/201,lwd=2);abline(v=0.5+0.5/201,lwd=2)}, xlab="Prey Pred.", ylab="Pred. Prey", ) if(PDF) dev.off() ### plot SD Int MD for each wave 2 acceptable run ### if(PDF) pdf("plots_health/M1_LV_30pts_IntMD_contributions.pdf",height=7.5,width=6.4) par(mfrow=c(2,1),mar=c(4,4,3,3),cex=0.75) lwd1 <- 1.4 pp_ind <- rbind(1:1001,1002:2002) for(j in 1:2){ plot(seq(0,10,len=1001),md_full[,1][pp_ind[j,]],ty="n",ylim=c(0,max(md_full)),xlab="Time",ylab="SD of Int MD",main=paste("SD of Internal MD for 60 Input Points:",c("Prey","Pred")[j])) abline(h=0) for(k in 1:ngood) lines(seq(0,10,len=1001),md_full[,k][pp_ind[j,]]) } if(PDF) dev.off() ### plot Bias Int MD for each wave 2 acceptable run ### if(PDF) pdf("plots_health/N1_LV_30pts_IntMD_Bias.pdf",height=7.5,width=6.4) par(mfrow=c(2,1),mar=c(4,4,3,3),cex=0.75) lwd1 <- 1.4 plot(seq(0,10,len=1001),bias_full_prey[,1],ty="n",ylim=range(bias_full_prey),xlab="Time",ylab="SD of Int MD",main="Bias of Internal MD for 60 Input Points: Prey") abline(h=0) for(k in 1:ngood) lines(seq(0,10,len=1001),bias_full_prey[,k],col=4) plot(seq(0,10,len=1001),bias_full_pred[,1],ty="n",ylim=range(bias_full_pred),xlab="Time",ylab="SD of Int MD",main="Bias of Internal MD for 60 Input Points: Pred") abline(h=0) for(k in 1:ngood) lines(seq(0,10,len=1001),bias_full_pred[,k],col=4) if(PDF) dev.off() ########################################################################################################### ### Now emulate Int MD Var and Bias at the 40 observed points ### ### time pts at observations ### al_obs <- match(t_obs,times) ### create int md sd at t obs for prey and pred ### md_int_sd_tobs_prey <- md_full[al_obs,] dim(md_int_sd_tobs_prey) md_int_sd_tobs_pred <- md_full[al_obs+1001,] dim(md_int_sd_tobs_pred) ### create int md bias at t obs ### md_int_bias_tobs_prey <- bias_full_prey[al_obs,] dim(md_int_bias_tobs_prey) md_int_bias_tobs_pred <- bias_full_pred[al_obs,] dim(md_int_bias_tobs_pred) ### set up simple lm emulators for int md sd and bias ### lm_list_int_sd_prey <- vector(mode="list") lm_list_int_sd_pred <- vector(mode="list") lm_list_int_bias_prey <- vector(mode="list") lm_list_int_bias_pred <- vector(mode="list") for(i in 1:nrow(md_int_sd_tobs_prey)){ lm_list_int_sd_prey[[i]] <- lm(y~th1+th2+th3,data=data.frame(y=md_int_sd_tobs_prey[i,],x1_good)) lm_list_int_sd_pred[[i]] <- lm(y~th1+th2+th3,data=data.frame(y=md_int_sd_tobs_pred[i,],x1_good)) lm_list_int_bias_prey[[i]] <- lm(y~th1+th2+th3,data=data.frame(y=md_int_bias_tobs_prey[i,],x1_good)) lm_list_int_bias_pred[[i]] <- lm(y~th1+th2+th3,data=data.frame(y=md_int_bias_tobs_pred[i,],x1_good)) } ### function to propose new input points, evaluate model (would emulate for slower model), ### emulate int MD SD and bias and hence compute implausibility ### (both wave 1 imp: using cautious 15% and wave 2 imp: using refined emulated Int MD) ### returns inputs pts that have Wave 1 Imp < 3 i.e. are non-implausible after first wave. wave2_imp <- function(nruns_w2=1000){ ### propose pts using Latin Hypercube with reduced ranges to improve acceptance rate ### xp_w2 <- lhs(nruns_w2,matrix(c(0.8,1.1, 0.00035,0.000572, 1.5,2.2),3,2,byrow=TRUE) ) colnames(xp_w2) <- paste("th",1:3,sep="") ### setupi matrices to store Int MD SD and bias for proposed runs ### em_int_sd_prey <- matrix(0,ncol=nrow(xp_w2),nrow=length(lm_list_int_sd_prey)) em_int_sd_pred <- em_int_sd_prey em_int_bias_prey <- em_int_sd_prey em_int_bias_pred <- em_int_sd_prey ### eval emulators for Int MD SD and bias at inputs xp_w2 ### for(i in 1:length(lm_list_int_sd_prey)){ em_int_sd_prey[i,] <- predict(lm_list_int_sd_prey[[i]],newdata = as.data.frame(xp_w2)) em_int_sd_pred[i,] <- predict(lm_list_int_sd_pred[[i]],newdata = as.data.frame(xp_w2)) em_int_bias_prey[i,] <- predict(lm_list_int_bias_prey[[i]],newdata = as.data.frame(xp_w2)) em_int_bias_pred[i,] <- predict(lm_list_int_bias_pred[[i]],newdata = as.data.frame(xp_w2)) } ### Evaluate model at inputs xp_w2 (would emulate this for slower model) ### out_w2 <- array(0,c(length(times),3,nrow(xp_w2)),dimnames=list( NULL,c("time","Prey","Pred"), NULL)) for(i in 1:nruns_w2){ out_w2[,,i] <- as.matrix(lsoda(ystart, times, Lotka_Volterra, parms=xp_w2[i,], maxsteps=20000)) # solve the diff eqtn } ### restrict outputs to times pts where observations are ### al_obs <- match(t_obs,times) prey_w2 <- out_w2[al_obs,2,] pred_w2 <- out_w2[al_obs,3,] ### evaluate wave 1 implausibility for pred and prey outputs, using 15% cautious MD ### MD_approx_w1 <- 0.15*y_obs # 15 percent approx wave 1 MD prey_imp_w1 <- abs(prey_w2 - y_obs[,"Prey"])/sqrt(obs_sd^2 + MD_approx_w1[,"Prey"]^2) pred_imp_w1 <- abs(pred_w2 - y_obs[,"Pred"])/sqrt(obs_sd^2 + MD_approx_w1[,"Pred"]^2) ### evaluate wave 2 implausibility for pred and prey outputs, using refined Int SD and bias and Ext MD ### MD_ext_sd <- 0.02*y_obs # 2% of obs for Ext MD prey_imp_w2 <- abs(prey_w2 - em_int_bias_prey - y_obs[,"Prey"])/sqrt(obs_sd^2 + MD_ext_sd[,"Prey"]^2 + em_int_sd_prey^2) pred_imp_w2 <- abs(pred_w2 - em_int_bias_pred - y_obs[,"Pred"])/sqrt(obs_sd^2 + MD_ext_sd[,"Pred"]^2 + em_int_sd_pred^2) ### construct max implausibility for wave 1 and wave 2 only using first 10 outputs up to t <= 5 ### max_imp_w1 <- apply( rbind(prey_imp_w1[1:10,],pred_imp_w1[1:10,]), 2, max) # only match to first 10 outputs i.e. up to t=5 max_imp_w2 <- apply( rbind(prey_imp_w2[1:10,],pred_imp_w2[1:10,]), 2, max) # only match to first 10 outputs i.e. up to t=5 ### Impose Imp<3 cut on wave 1 imps and return runs that satisfy this along with their wave 1 and 2 imps ### good_w1 <- max_imp_w1 < 3 cbind(xp_w2,max_imp_w1,max_imp_w2)[good_w1,] } ### Simple propose and reject step to gather points that are decent after wave 2 i.e. have Imp<5 at wave 2 ### set.seed(2) xp_w2_good <- matrix(0,nrow=0,ncol=5) for(k in 1:10){ xp_temp <- wave2_imp() # propose points xp_w2_good <- rbind(xp_w2_good,xp_temp[ xp_temp[,5]<5, ]) # only keep points with wave 2 Imp < 5 cat(k,nrow(xp_w2_good),"\n") } dim(xp_w2_good) hist(xp_w2_good[,5]) save(xp_w2_good,file="runs/Wave3_purple_pts.Rdata") # load("runs/Wave3_purple_pts.Rdata") ### store non-implausible input points that have wave 2 Imp < 3, i.e. the wave 3 purple points ### xp_w2_good_plot <- xp_w2_good[xp_w2_good[,5]<3.0,1:3] ngood_w2 <- nrow(xp_w2_good_plot) ngood_w2 # number of non-implausible wave 3 points ### run model for the good runs to obtain output for plots ### out_w2_good_plot <- array(0,c(length(times),3,nrow(xp_w2_good_plot)),dimnames=list( NULL,c("time","Prey","Pred"), NULL)) for(i in 1:ngood_w2){ out_w2_good_plot[,,i] <- as.matrix(lsoda(ystart, times, Lotka_Volterra, parms=xp_w2_good_plot[i,], maxsteps=20000)) } ### now plot wave 1, 2 and 3 runs on same output plot ### if(PDF) png("plots_health/J2_LV_wave1and2_runs_pred_prey_outputs_grey.png",width=7.5*0.8,height=8*0.8,pointsize=12,res=72*2+40,units="in") plot_pred_prey_array(out=out_w1,out_wx=out_w2_good_plot,lw1=0.6,xylog="",good_id=which(al_good==1),obs_col="yellow",wv1_col=c("grey70","grey70")) if(PDF) dev.off() ### plot wave 1, 2 and 3 inputs points on same pairs plot ### if(PDF) pdf("plots_health/P1_LV_wave2_nimp_inputs_Int_Ext_MD.pdf",height=8*0.8,width=8.5*0.8) xran <- rbind(parms*0.7,parms*1.3) lab1 <- c(expression(x[1]),expression(x[2]),expression(x[3])) col1 <- rep("green",ngood) col1[x1_good_best_ind] <- "red" col2 <- "purple" pairs(rbind(x1,x1_good,xp_w2_good_plot,x1_good[x1_good_best_ind,],xran), pch=c(rep(16,nruns),rep(21,ngood),rep(21,ngood_w2),21,16,16), col=c(rep("grey80",nruns),rep(1,ngood),rep(1,ngood_w2),1,0,0), bg=c(rep("grey80",nruns),col1,rep(col2,ngood_w2),"red",0,0), cex=c(rep(1,nruns),rep(1.4,ngood),rep(1.2,ngood_w2),1.4,1,1),labels=lab1) if(PDF) dev.off() ######################################################################################################################## ###### Maximum peak analysis for section 6(c) of paper ###### library(DescTools) ### simple function to draw prediction ellispses ### plot_ellipse <- function(mu,covar,col1="green",bor1=1){ m <- covar eig <- eigen(m) eig.val.sq <- sqrt(eig$values) eig.vec <- eig$vectors DrawEllipse(x=mu[1], y=mu[2], radius.x=3*eig.val.sq[1], radius.y=3*eig.val.sq[2], rot=acos(eig.vec[1,1]), border=bor1, lwd=0.5,col=col1) } ### load max peak positions for obs data if needed, created in Stochastic_Lot_vol_Real_World.R### # load("Observed Data/Max_peaks_Obs_and_ObsErrors_V1.Rdata") ### peaks for obs data ### Forecast function which gives forecast of second peak time and magnitude for both pred and prey ### while using knowledge of the first peak. The model discrepancy epsilon is implicitly updated ### by knowledge of the first peak in this calculatione forecast <- function(out_w1_good, # the array of original runs out_full # the array of perturbation to original runs ){ ngood <- dim(out_w1_good)[3] # number of original runs ### create objects to store all core forecast quantities ### f_xstar_mat1 <- matrix(0,nrow=8,ncol=ngood) # matrix to store f(x*) at max peak imp_max_xstar <- rep(0,ngood) # vector to store imp for each run cov_z_arra1 <- array(0,dim=c(4,4,ngood)) # array for Cov[z] priors Var_yp_arra1 <- array(0,dim=c(4,4,ngood)) # array for Var_yp for each run E_yp_z_mat1 <- matrix(0,nrow=4,ncol=ngood) # array for Bayes Linear updates of E_z[y_p] Var_yp_z_arra1 <- array(0,dim=c(4,4,ngood)) # array for Bayes Linear updates of Var_z[y_p] ### loop over runs in out_w1_good (loops just for clarity, not efficiency!) ### for(g_ind in 1:ngood){ ### find index, time and height of two peaks for the actual green run g_ind ### f_prey_T1_ind <- which.max(out_w1_good[1:500,2,g_ind]) f_pred_T1_ind <- which.max(out_w1_good[1:500,3,g_ind]) f_prey_T2_ind <- 500 + which.max(out_w1_good[501:1000,2,g_ind]) f_pred_T2_ind <- 500 + which.max(out_w1_good[501:1000,3,g_ind]) f_prey_T1 <- times[f_prey_T1_ind] f_pred_T1 <- times[f_pred_T1_ind] f_prey_T2 <- times[f_prey_T2_ind] f_pred_T2 <- times[f_pred_T2_ind] f_prey_Y1 <- out_w1_good[f_prey_T1_ind,2,g_ind] f_pred_Y1 <- out_w1_good[f_pred_T1_ind,3,g_ind] f_prey_Y2 <- out_w1_good[f_prey_T2_ind,2,g_ind] f_pred_Y2 <- out_w1_good[f_pred_T2_ind,3,g_ind] ### The actual output of this xstar run ### f_xstar <- c("f_prey_T1"=f_prey_T1,"f_pred_T1"=f_pred_T1,"f_prey_Y1"=f_prey_Y1,"f_pred_Y1"=f_pred_Y1, "f_prey_T2"=f_prey_T2,"f_pred_T2"=f_pred_T2,"f_prey_Y2"=f_prey_Y2,"f_pred_Y2"=f_pred_Y2) f_xstar_mat1[,g_ind] <- f_xstar ### find index, time and height of two peaks for the 200 perturbations to green run g_ind ### prey_T1_ind <- apply(out_full[1:500,2,,g_ind],c(2),which.max) pred_T1_ind <- apply(out_full[1:500,3,,g_ind],c(2),which.max) prey_T2_ind <- 500 + apply(out_full[501:1001,2,,g_ind],c(2),which.max) pred_T2_ind <- 500 + apply(out_full[501:1001,3,,g_ind],c(2),which.max) prey_T1 <- times[prey_T1_ind] pred_T1 <- times[pred_T1_ind] prey_T2 <- times[prey_T2_ind] pred_T2 <- times[pred_T2_ind] prey_Y1 <- diag(out_full[prey_T1_ind,2,,g_ind]) pred_Y1 <- diag(out_full[pred_T1_ind,3,,g_ind]) prey_Y2 <- diag(out_full[prey_T2_ind,2,,g_ind]) pred_Y2 <- diag(out_full[pred_T2_ind,3,,g_ind]) ### Assess covariance of internal and external Model Discrepancy: Cov[epsilon_I]=sigma_I^2 and Cov[epsilon_E] ### ### the 200 perturbations for internal discrepancy assessment at this xstar run ### MD_maximums <- cbind(MD_Prey.T.1=prey_T1, MD_Pred.T.1=pred_T1, MD_Prey.Y.1=prey_Y1, MD_Pred.Y.1=pred_Y1, MD_Prey.T.2=prey_T2, MD_Pred.T.2=pred_T2, MD_Prey.Y.2=prey_Y2, MD_Pred.Y.2=pred_Y2) MD_Int_Maximums <- sqrt(diag(cov(MD_maximums))) MD_Int_Maximums_Cov <- cov(MD_maximums) MD_Int_Maximums_Cor <- cor(MD_maximums) sig_tmax_ext <- 0.03 # external MD for time location of peak sig_ymax_prey_ext <- max(y_obs[,"Prey"])*0.02/3 # external MD for magnitude of peak (prey) sig_ymax_pred_ext <- max(y_obs[,"Pred"])*0.02/3 # external MD for magnitude of peak (pred) ### Construct uncorrelated External MD Covariance matrix ### MD_Ext_Maximums_Cov <- diag(c(sig_tmax_ext,sig_tmax_ext, sig_ymax_prey_ext,sig_ymax_pred_ext, sig_tmax_ext,sig_tmax_ext, sig_ymax_prey_ext,sig_ymax_pred_ext)^2) ### Construct Total MD Covariance matrix ### MD_Tot_Maximums_Cov <- MD_Int_Maximums_Cov + MD_Ext_Maximums_Cov MD_Tot_Maximums_Cor <- cov2cor(MD_Tot_Maximums_Cov) ### Assess internal MD bias term mu_I(xstar) = E[epsilon_I] (note minus sign change in definition to usual bias) ### E_epsilon_I <- apply(MD_maximums,2,mean) - f_xstar ### construct obs error 4x4 cov matrix for first peak measurements ### obs_err_Maximums_Cov <- diag( rep(c(obs_max_T_sd,obs_max_Y_sd),c(2,2))^2 ) ### construct cov matrix for history measurement first peak z and future prediction y_p and cov between ### cov_z <- MD_Tot_Maximums_Cov[1:4,1:4] + obs_err_Maximums_Cov # Cov[z] cov_yp <- MD_Tot_Maximums_Cov[5:8,5:8] # Cov[y_p] cov_yp_z <- MD_Tot_Maximums_Cov[5:8,1:4] # Cov[y_p,z] z <- obs_maximums[1:4] # the first peak obs data ### Calculate the univariate max implausibility ### imp_max_xstar[g_ind] <- max(abs(f_xstar[1:4] + E_epsilon_I[1:4] - z)/sqrt(diag(cov_z))) ### actual run values for 1st and 2nd peak ### fh_xstar <- f_xstar[1:4] fp_xstar <- f_xstar[5:8] ### update prediction y_p by measurement z ### E_yp <- fp_xstar + E_epsilon_I[5:8] # f(xstar) + E[epsilon_I_p] E_z <- fh_xstar + E_epsilon_I[1:4] # f(xstar) + E[epsilon_I_h] Var_yp <- cov_yp # rename to Var[y_p] to match BL update formula! E_yp_z <- E_yp + cov_yp_z %*% chol2inv(chol(cov_z)) %*% (z - E_z) # Full BL update of E[y_p] -> E_z[y_p] Var_yp_z <- Var_yp - cov_yp_z %*% chol2inv(chol(cov_z)) %*% t(cov_yp_z) # Full BL update of Var[y_p] -> Var_z[y_p] ### store results in arrays etc ### cov_z_arra1[,,g_ind] <- cov_z Var_yp_arra1[,,g_ind] <- Var_yp E_yp_z_mat1[,g_ind] <- c(E_yp_z) Var_yp_z_arra1[,,g_ind] <- Var_yp_z cat(g_ind,"\n") } ### return all important quantities from forecast ### return(list(f_xstar_mat1 = f_xstar_mat1, imp_max_xstar = imp_max_xstar, cov_z_arra1 = cov_z_arra1, Var_yp_arra1 = Var_yp_arra1, E_yp_z_mat1 = E_yp_z_mat1, Var_yp_z_arra1 = Var_yp_z_arra1 )) } ### perform forecast using wave 2 (green) runs ### forec_out <- forecast(out_w1_good = out_w1_good, out_full = out_full) ### define objects for further plotting ### f_xstar_mat1 = forec_out$f_xstar_mat1 imp_max_xstar = forec_out$imp_max_xstar cov_z_arra1 = forec_out$cov_z_arra1 Var_yp_arra1 = forec_out$Var_yp_arra1 E_yp_z_mat1 = forec_out$E_yp_z_mat1 Var_yp_z_arra1 = forec_out$Var_yp_z_arra1 ### plot single run forecast (first peak then second peak) ### pdf(file="plots_health/R0_prediction_plot_single_run_A1.pdf",height=5,width=9) op <- par(mfrow=c(1,2),mar=c(5.1,4.1,3.1,0.5)) for(j in 1:2){ i <- 9 preypred=1 if(preypred==1) {ft1_ind <- c(1,3); ft2_ind <- c(5,7)} if(preypred==2) {ft1_ind <- c(2,4); ft2_ind <- c(6,8)} plot(t_obs,y_obs[,preypred],xlim=rbind(c(0,3.25),c(6,8.5))[j,],ylim=c(4000,13500),pch=16,col=c(4,2)[preypred],cex=0.8, ylab="Prey",xlab="Time(t)",main=c("First Peak","Second Peak")[j]) for(k in 1:200) lines(times,out_full[,2,k,i],col=6) # the perturbations of this single wave 2 run arrows(t_obs,y_obs[,preypred] + 3*obs_sd,t_obs,y_obs[,preypred] - 3*obs_sd,angle=90,code=3,col=c(4,2)[preypred],length=0.03,lwd=3) lines(times,out_w1_good[,2,i],lwd=2) ### naive predictions of first and second peak ### plot_ellipse(mu=f_xstar_mat1[ft1_ind,i],covar=cov_z_arra1[ft1_ind,ft1_ind,i],col1=SetAlpha("green",0.6)) plot_ellipse(mu=f_xstar_mat1[ft2_ind,i],covar=Var_yp_arra1[ft1_ind,ft1_ind,i],col1=SetAlpha("green",0.6),bor1=1) ### Bayes linear update of forecast of second peak wrt first peak information ### plot_ellipse(mu=E_yp_z_mat1[ft1_ind,i],covar=Var_yp_z_arra1[ft1_ind,ft1_ind,i],col1=SetAlpha("yellow",0.6),bor1=1) ### centres of forecasts ### points(f_xstar_mat1[ft1_ind[1],i],f_xstar_mat1[ft1_ind[2],i],pch=16,cex=0.5,col="darkgreen") points(f_xstar_mat1[ft2_ind[1],i],f_xstar_mat1[ft2_ind[2],i],pch=16,cex=0.5,col="darkgreen") points(E_yp_z_mat1[ft1_ind[1],i],E_yp_z_mat1[ft1_ind[2],i],pch=16,cex=0.5,col="red") ### actual observed peak locations ### if(preypred==1) points(obs_maximums[c(1,5)],obs_maximums[c(3,7)],pch=21,bg="cyan",cex=1.1) if(preypred==2) points(obs_maximums[c(2,6)],obs_maximums[c(4,8)],pch=21,bg="cyan",cex=1.1) ### legends ### if(j==1) legend("topleft",legend=c("Direct run prediction f(x) for 1st peak", "Actual observed 1st peak","Observed time series data", "Run output f_t(x)","Int MD perturbations"), pch=c(21,21,21,NA,NA),pt.bg=c("green","cyan",4,NA,NA), lty=c(NA,NA,NA,1,1),col=c(1,1,1,1,6),lwd=c(1,1,1,2,2),cex=0.65,pt.cex=1.3) if(j==2) legend("topleft",legend=c("Direct run prediction f(x) for 2nd peak", "Updated MD prediction for 2nd peak y_p(x)", "Observed 2nd peak (prediction target)","Observed time series data (not used)", "Run output f_t(x)","Internal MD perturbations"), pch=c(21,21,21,21,NA,NA),pt.bg=c("green","yellow","cyan",4,NA,NA), lty=c(NA,NA,NA,NA,1,1),col=c(1,1,1,1,1,6),lwd=c(1,1,1,1,2,2),cex=0.65,pt.cex=1.3) } par(op) par()$mar dev.off() ### generate perturbations for the purple wave 3 runs aka xp_w2_good_plot (would normally emulate the ### ### 8x8 internal MD covariance matrix but instead just generating the perturbations for simplicity of exposition) ### set.seed(1) nreps <- 200 out_w2_good_plot_full <- array(0,c(length(times),3,nreps,ngood_w2),dimnames=list( NULL,c("time","Prey","Pred"), NULL,NULL)) for(j in 1:ngood_w2){ for(i in 1:nreps){ ystart <- c(Prey = rnorm(1,2000,40),Pred = rnorm(1,800,40)) # random initial conditions a <- runif(1,min=5-0.5,max=5+0.5) b <- runif(1,min=0.6,max=0.9) c <- rep(runif(1,0.02,0.05),3) out_w2_good_plot_full[,,i,j] <- as.matrix(lsoda(ystart, times, Lotka_Volterra_tvar, parms=xp_w2_good_plot[j,], maxsteps=20000)) cat(i," ",j,"\n") } } save(out_w2_good_plot_full,nreps,ngood_w2,xp_w2_good_plot,file=paste("runs/Model_Disc_Assess_full_array_",nreps,"_reps_w2goodplot_aka_w3.Rdata",sep="")) # load(paste("runs/Model_Disc_Assess_full_array_",nreps,"_reps_w2goodplot_aka_w3.Rdata",sep="")) ### now perform forecast using wave 3 runs i.e. xp_w2_good_plot runs ### forec_out_w3 <- forecast(out_w1_good = out_w2_good_plot, out_full = out_w2_good_plot_full) ### plot several wave 2 and wave 3 forecast ellipses on same plot ### reasonable_first_max <- (1:ngood)[imp_max_xstar