require(e1071)
require(truncdist)
require(msm)

#BOOTSTRAPPING
bootci.ext=function(temp) {
x.resample<-numeric(); xmean.resample<-numeric(); xmedian.resample<-numeric();x25.resample<-numeric();
x75.resample<-numeric();x95.resample<-numeric();x05.resample<-numeric();iqr.resample<-numeric();range.resample<-numeric();
x975.resample<-numeric();x025.resample<-numeric();
for (i in 1:100) {
for (j in 1:length(temp)) {x.resample[j]<-sample(temp,1,replace=T)}
xmean.resample[i]<-mean(x.resample, na.rm=T)
xmedian.resample[i]<-median(x.resample, na.rm=T)
x25.resample[i]<-quantile(x.resample, c(0.25), na.rm=T)
x75.resample[i]<-quantile(x.resample, c(0.75), na.rm=T)
iqr.resample[i]=x75.resample[i]-x25.resample[i]
x95.resample[i]<-quantile(x.resample, c(0.95), na.rm=T)
x975.resample[i]<-quantile(x.resample, c(0.975), na.rm=T)
x05.resample[i]<-quantile(x.resample, c(0.05), na.rm=T)
x025.resample[i]<-quantile(x.resample, c(0.025), na.rm=T)
range.resample[i]=x975.resample[i]-x025.resample[i]
}
out=list(resample=xmean.resample, resamplemedian=xmedian.resample, resample25=x25.resample, resample75=x75.resample, resampleiqr=iqr.resample, 
resamplerange=range.resample, resample25=x25.resample, resample95=x95.resample,resample05=x05.resample)
out
}

###########################################################################################################
#CORBULA AGES
###########################################################################################################
#the excel file needs to be saved as tab-delimited.txt file
setwd("C:/Data/Adriatic_sea_project/Piran/Revision")
data.raw <- read.delim("Supplement Table 3-Corbula age data.txt", header=TRUE)
finalAges=data.raw[,"Age"]
remove=which(is.na(finalAges))
data.raw=data.raw[-remove,]
finalAges=data.raw[,"Age"]

Corbula.Piran.ages=finalAges[data.raw$Station=="PiranII M53"]
Corbula.Piran.5depths=as.numeric(as.vector(data.raw[data.raw$Station=="PiranII M53",]$"X5.cm.group..cm."))
Corbula.Piran.rev5depths=Corbula.Piran.5depths
Corbula.Piran.station=data.raw[data.raw$Station=="PiranII M53","Station"]
Corbula.Piran.5depths=ifelse(Corbula.Piran.5depths<25,Corbula.Piran.5depths-2,Corbula.Piran.5depths-2.5)
Corbula.Piran2.carbonages=c(8594,8554,8750,8594,5739)
Corbula.Piran2.carbonage.depths=c(87.5,87.5,102.5,102.5,67.5)

#variance parameter for lognormal calibration model (TDK1)
Corbula.ln.d=-2.896
#variance parameter for gamma calibration model (SPK0)
Corbula.gamma.ln.d=4.954087
###########################################################################################################
#GOULDIA AGES
###########################################################################################################
#the excel file needs to be saved as tab-delimited.txt file
data.raw <- read.delim("Supplement Table 2-Gouldia age data.txt", header=TRUE)
finalAges=data.raw[,"Age"]
remove=which(is.na(finalAges))
data.raw=data.raw[-remove,]
finalAges=data.raw[,"Age"]

Gouldia.Piran2.ages=finalAges[data.raw$Station=="Piran II"]
Gouldia.Piran2.maxdepths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran II",]$"max..depth..cm."))
Gouldia.Piran2.5depths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran II",]$"X5.cm.group..cm."))
Gouldia.Piran2.5depths=ifelse(Gouldia.Piran2.5depths<25,Gouldia.Piran2.5depths-2,Gouldia.Piran2.5depths-2.5)
Gouldia.Piran2.station=data.raw[data.raw$Station=="Piran II","Station"]
Gouldia.Piran1.ages=finalAges[data.raw$Station=="Piran I"]
Gouldia.Piran1.maxdepths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"max..depth..cm."))
Gouldia.Piran1.5depths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"X5.cm.group..cm."))
Gouldia.Piran1.5depths=ifelse(Gouldia.Piran1.5depths<25,Gouldia.Piran1.5depths-2,Gouldia.Piran1.5depths-2.5)
Gouldia.Piran1.10depths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"X10.cm.group..cm."))
Gouldia.Piran1.station=data.raw[data.raw$Station=="Piran I","Station"]
Gouldia.Piran2.carbonages=as.numeric(as.vector(data.raw[data.raw$Station=="Piran II",]$"ageMedian"))
Gouldia.Piran2.carbonage.depths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran II",]$"max..depth..cm."))
Gouldia.Piran1.carbonages=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"ageMedian"))
Gouldia.Piran1.carbonage.depths=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"max..depth..cm."))
Gouldia.Piran2.carbonage.lci=as.numeric(as.vector(data.raw[data.raw$Station=="Piran II",]$"ageYng"))
Gouldia.Piran2.carbonage.uci=as.numeric(as.vector(data.raw[data.raw$Station=="Piran II",]$"ageOld"))
Gouldia.Piran1.carbonage.lci=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"ageYng"))
Gouldia.Piran1.carbonage.uci=as.numeric(as.vector(data.raw[data.raw$Station=="Piran I",]$"ageOld"))
id1=which(!is.na(Gouldia.Piran1.carbonages))
id2=which(!is.na(Gouldia.Piran2.carbonages))
Gouldia.Piran2.carbonages=Gouldia.Piran2.carbonages[id2]
Gouldia.Piran2.carbonage.lci=Gouldia.Piran2.carbonage.lci[id2]
Gouldia.Piran2.carbonage.uci=Gouldia.Piran2.carbonage.uci[id2]
Gouldia.Piran2.carbonage.depths=Gouldia.Piran2.carbonage.depths[id2]
Gouldia.Piran1.carbonages=Gouldia.Piran1.carbonages[id1]
Gouldia.Piran1.carbonage.lci=Gouldia.Piran1.carbonage.lci[id1]
Gouldia.Piran1.carbonage.uci=Gouldia.Piran1.carbonage.uci[id1]
Gouldia.Piran1.carbonage.depths=Gouldia.Piran1.carbonage.depths[id1]

#variance parameter for gamma calibration model for Gouldia
gamma.ln.d=log(51.16437)

##################################################################################################################################
#ARCA AND OSTREA AGES (CALIBRATED 2013 AD)
##################################################################################################################################
Ostrea.Piran1.age=c(6777, 6813,7582, 6317, 5918)
Ostrea.Piran1.age.lci=c(6613, 6661,7472, 6183, 5758)
Ostrea.Piran1.age.uci=c(6944, 6984,7710, 6454, 6058)
Ostrea.Piran1.depth=c(13,22.5,57.5,19,17)
Arca.Piran1.age=c(2296,6353)
Arca.Piran1.age.lci=c(2125,6229)
Arca.Piran1.age.uci=c(2434,6494)
Arca.Piran1.depth=c(11,17)
Arca.Piran2.age=c(4869,5126,3545)
Arca.Piran2.age.lci=c(4673,4936,3405)
Arca.Piran2.age.uci=c(5028,5318,3695)
Arca.Piran2.depth=c(27.5,13,3)

par(mfrow=c(2,2))
boxplot(Corbula.Piran.ages, Ostrea.Piran1.age,c(Arca.Piran1.age,Arca.Piran2.age),c(Gouldia.Piran1.ages,Gouldia.Piran2.ages), col="gray",
names=c("Corbula","Ostrea","Arca","Gouldia"), range=0, las=2)

##################################################################################################################################
#AGE BOXPLOTS
##################################################################################################################################
par(mfrow=c(1,2))
out1=boxplot(split(Gouldia.Piran1.ages,Gouldia.Piran1.5depths), at=-as.numeric(names(split(Gouldia.Piran1.ages,Gouldia.Piran1.5depths))), 
col="gray81", range=1.5, boxwex=3, horizontal=TRUE, ylab="Sediment depth (cm)", xlab="Time since 2013", main="Piran 1", xlim=c(-150,0), 
ylim=c(0,12000), cex.main=0.9, outline=F)
points(Gouldia.Piran1.carbonages, -Gouldia.Piran1.carbonage.depths, pch=21, bg="gray",cex=1.4)
points(Ostrea.Piran1.age, -Ostrea.Piran1.depth, pch=22, bg="black",cex=1.4)
points(Arca.Piran1.age, -Arca.Piran1.depth, pch=22, bg="white",cex=1.4)
segments(x0=-3000, y0=-c(12,20,65), x1 = 0, y1 = -c(12,20,65),lty=2)
rect(xleft=10800, ybottom=-150, xright=12500, ytop=-70, col="gray", border="white")
rect(xleft=10800, ybottom=-70, xright=12500, ytop=-35, col="gray", border="white")
rect(xleft=10800, ybottom=-35, xright=12500, ytop=-8, col="gray", border="white")
rect(xleft=10800, ybottom=-8, xright=12500, ytop=0, col="gray", border="white")
text(11750,-100,labels="sandy mud",srt=90, cex=0.8)
text(11750,-50,labels="sand",srt=90, cex=0.8)
text(11750,-20,labels="shell b.",srt=90, cex=0.8)
text(11750,-4,labels="sand",srt=0, cex=0.8)


out2=boxplot(split(Gouldia.Piran2.ages,Gouldia.Piran2.5depths), at=-as.numeric(names(split(Gouldia.Piran2.ages,Gouldia.Piran2.5depths))), 
col="gray81", range=1.5, boxwex=3, horizontal=TRUE, ylab="Sediment depth (cm)", xlab="Time since 2013", main="Piran 2", xlim=c(-150,0), 
ylim=c(0,12000), cex.main=0.9, outline=F)
out3=boxplot(split(Corbula.Piran.ages,Corbula.Piran.5depths), at=-as.numeric(names(split(Corbula.Piran.ages,Corbula.Piran.5depths)))+2, add=T, 
col="gray51", range=1.5, boxwex=3, horizontal=TRUE, yaxt="n", outline=F)
#legend(x="topright", pch=22, pt.bg=c("gray71","gray31"), legend=c("Gouldia","Corbula"), bty="n")
points(Gouldia.Piran2.carbonages, -Gouldia.Piran2.carbonage.depths, pch=21, bg="gray",cex=1.4)
points(Arca.Piran2.age, -Arca.Piran2.depth, pch=21, bg="white",cex=1.4)
rect(xleft=10800, ybottom=-150, xright=12500, ytop=-75, col="gray", border="white")
rect(xleft=10800, ybottom=-75, xright=12500, ytop=-35, col="gray", border="white")
rect(xleft=10800, ybottom=-35, xright=12500, ytop=-8, col="gray", border="white")
rect(xleft=10800, ybottom=-8, xright=12500, ytop=0, col="gray", border="white")
text(11750,-100,labels="sandy mud",srt=90, cex=0.8)
text(11750,-50,labels="sand",srt=90, cex=0.8)
text(11750,-20,labels="shell b.",srt=90, cex=0.8)
text(11750,-4,labels="sand",srt=0, cex=0.8)


##################################################################################################################################
#SUMMARY PER UNIT
##################################################################################################################################
par(mfrow=c(2,3))
tempage=c(Gouldia.Piran1.ages,Arca.Piran1.age,Ostrea.Piran1.age)
tempunit=c(Gouldia.Piran1.5depths, 10, 18, 14, 22.5, 57.5, 18,18)
tempunit=ifelse(tempunit<8, 1, tempunit)
tempunit=ifelse(tempunit>8 & tempunit<16, 2, tempunit)
tempunit=ifelse(tempunit>16 & tempunit<35, 3, tempunit)
tempunit=ifelse(tempunit>35 & tempunit<65, 4, tempunit)
tempunit=ifelse(tempunit>65, 5, tempunit)
out5=boxplot(rev(split(-tempage, tempunit)), range=0, outline=F, col="gray", ylab="Stratigraphic units", ylim=c(0,-12000), 
cex.main=0.9, main="Piran 1", horizontal=TRUE, yaxt="n")
axis(2, at=c(1,2,3,4,5), labels=c(5,4,3,2,1), las=2)
uci=unlist(lapply(split(tempage, tempunit), quantile, 0.975))
lci=unlist(lapply(split(tempage, tempunit), quantile, 0.025))

tempage=c(Gouldia.Piran2.ages,Corbula.Piran.ages,Arca.Piran2.age)
tempunit=c(Gouldia.Piran2.5depths,Corbula.Piran.5depths,  27.5, 14, 2)
tempunit=ifelse(tempunit<8, 1, tempunit)
tempunit=ifelse(tempunit>8 & tempunit<16, 2, tempunit)
tempunit=ifelse(tempunit>16 & tempunit<35, 3, tempunit)
tempunit=ifelse(tempunit>35 & tempunit<65, 4, tempunit)
tempunit=ifelse(tempunit>65, 5, tempunit)
out7=boxplot(rev(split(-tempage, tempunit)), range=0, outline=F, col="gray", ylab="Stratigraphic units", ylim=c(0,-12000), 
cex.main=0.9, main="Piran 2", horizontal=TRUE, yaxt="n")
axis(2, at=c(1,2,3,4,5), labels=c(5,4,3,2,1), las=2)
uci=unlist(lapply(split(tempage, tempunit), quantile, 0.975))
lci=unlist(lapply(split(tempage, tempunit), quantile, 0.025))

boxplot(-c(Corbula.Piran.ages),-c(Ostrea.Piran1.age), -c(Arca.Piran2.age,Arca.Piran1.age),-c(Gouldia.Piran1.ages,Gouldia.Piran2.ages),col="gray",range=0, 
cex.main=0.9,main="Both stations", names=c("Corbula","Ostrea","Arca","Gouldia"), horizontal=TRUE, yaxt="n", ylim=c(0,-12000))
axis(2, at=c(1,2,3,4), labels=c("Corbula","Ostrea","Arca","Gouldia"), las=2)
##################################################################################################################################
par(mfrow=c(2,2))
tempage1=c(Gouldia.Piran1.ages,Arca.Piran1.age,Ostrea.Piran1.age)
tempunit1=c(Gouldia.Piran1.5depths, 10, 18, 14, 22.5, 57.5, 18,18)
tempage2=c(Gouldia.Piran2.ages,Corbula.Piran.ages,Arca.Piran2.age)
tempunit2=c(Gouldia.Piran2.5depths,Corbula.Piran.5depths,  27.5, 14, 2)
tempage=c(tempage1,tempage2)
tempunit=c(tempunit1,tempunit2)

tempunit=ifelse(tempunit<8, 1, tempunit)
tempunit=ifelse(tempunit>8 & tempunit<16, 2, tempunit)
tempunit=ifelse(tempunit>16 & tempunit<35, 3, tempunit)
tempunit=ifelse(tempunit>35 & tempunit<65, 4, tempunit)
tempunit=ifelse(tempunit>65, 5, tempunit)
out8=boxplot(split(-tempage, tempunit), range=1.5, outline=F, col="gray", ylab="Age (years)", ylim=c(-11000,0), cex.main=0.9, main="Piran 1-2")
uci=unlist(lapply(split(tempage, tempunit), quantile, 0.975))
lci=unlist(lapply(split(tempage, tempunit), quantile, 0.025))

boxplot(c(Gouldia.Piran2.ages,Gouldia.Piran2.ages),Corbula.Piran.ages,c(Arca.Piran2.age,Arca.Piran1.age),Ostrea.Piran1.age, col="gray",
las=2, names=c("Gouldia","Corbula","Arca","Ostrea"))

##################################################################################################################################
par(mfrow=c(1,3))
par(mar=c(5,4,2,1))
plot(out4$stats[3,],-as.numeric(out4$names), xlim=c(0,12000), pch=16, ylab="Sediment depth", ylim=c(-160,0), cex.axis=1.2, cex.lab=1.2,
frame=F, main="Piran 1", xlab="Age (y)", cex.main=1.2, cex=1.4)
rect(xleft=0, ybottom=-150, xright=12000, ytop=-65, col="gray31", border="white")
rect(xleft=0, ybottom=-65, xright=12000, ytop=-35, col="gray51", border="white")
rect(xleft=0, ybottom=-35, xright=12000, ytop=-8, col="gray71", border="white")
rect(xleft=0, ybottom=-8, xright=12000, ytop=0, col="gray91", border="white")
text(11000,-100,labels="sandy mud",srt=90, cex=1)
text(11000,-50,labels="sand",srt=90, cex=1)
text(11000,-20,labels="shell b.",srt=90, cex=1)
text(11000,-4,labels="sand",srt=0, cex=1)
abline(v=c(1500,3500,6500,7500), lty=2)
segments(x0=out4$stats[2,], x1=out4$stats[4,], y0=-as.numeric(out4$names), y1=-as.numeric(out4$names))
points(out4$stats[3,],-as.numeric(out4$names), pch=16, cex=1.4)


##################################################################################################################################
par(mfrow=c(1,3))
par(mar=c(5,4,2,1))
plot(out1$stats[3,],-as.numeric(out1$names), xlim=c(0,12000), pch=16, ylab="Sediment depth", ylim=c(-160,0), cex.axis=1.2, cex.lab=1.2,
frame=F, main="Piran 1", xlab="Age (y)", cex.main=1.2, cex=1.4)
rect(xleft=0, ybottom=-150, xright=12000, ytop=-65, col="gray31", border="white")
rect(xleft=0, ybottom=-65, xright=12000, ytop=-35, col="gray51", border="white")
rect(xleft=0, ybottom=-35, xright=12000, ytop=-8, col="gray71", border="white")
rect(xleft=0, ybottom=-8, xright=12000, ytop=0, col="gray91", border="white")
text(11000,-100,labels="sandy mud",srt=90, cex=1)
text(11000,-50,labels="sand",srt=90, cex=1)
text(11000,-20,labels="shell b.",srt=90, cex=1)
text(11000,-4,labels="sand",srt=0, cex=1)
segments(x0=out1$stats[2,], x1=out1$stats[4,], y0=-as.numeric(out1$names), y1=-as.numeric(out1$names))
points(out1$stats[3,],-as.numeric(out1$names), pch=16, cex=1.4)
points(Gouldia.Piran1.carbonages, -Gouldia.Piran1.carbonage.depths, pch=21, bg="gray",cex=1.4)
points(Ostrea.Piran1.age, -Ostrea.Piran1.depth, pch=22, bg="black",cex=1.4)
points(Arca.Piran1.age, -Arca.Piran1.depth, pch=22, bg="white",cex=1.4)
abline(v=c(1500,3500,6500,7500), lty=2)

plot(out2$stats[3,],-as.numeric(out2$names), xlim=c(0,12000), pch=16, ylab="Sediment depth", ylim=c(-160,0), cex.axis=1.2, cex.lab=1.2,
frame=F, main="Piran 2", xlab="Age (y)", cex.main=1.2, cex=1.4)
rect(xleft=0, ybottom=-150, xright=12000, ytop=-65, col="gray31", border="white")
rect(xleft=0, ybottom=-65, xright=12000, ytop=-35, col="gray51", border="white")
rect(xleft=0, ybottom=-35, xright=12000, ytop=-8, col="gray71", border="white")
rect(xleft=0, ybottom=-8, xright=12000, ytop=0, col="gray91", border="white")
text(11000,-100,labels="sandy mud",srt=90, cex=1)
text(11000,-50,labels="sand",srt=90, cex=1)
text(11000,-20,labels="shell b.",srt=90, cex=1)
text(11000,-4,labels="sand",srt=0, cex=1)
segments(x0=out2$stats[2,], x1=out2$stats[4,], y0=-as.numeric(out2$names), y1=-as.numeric(out2$names))
segments(x0=out3$stats[2,], x1=out3$stats[4,], y0=-as.numeric(out3$names), y1=-as.numeric(out3$names))
points(out2$stats[3,],-as.numeric(out2$names), pch=21, bg="black", cex=1.4)
points(out3$stats[3,],-as.numeric(out3$names), pch=21, bg="gray", cex=1.4)
points(Arca.Piran2.age, -Arca.Piran2.depth, pch=22, bg="gray",cex=1.8)
abline(v=c(1500,3500,6500,7500), lty=2)



###################################
#FIGURE 05 MEDIAN AGE AND BOXPLOTS#
###################################
par(mfrow=c(1,2))
par(mar=c(10,2.5,2,1))
out1=boxplot(split(Gouldia.Piran1.ages,Gouldia.Piran1.5depths), at=-as.numeric(names(split(Gouldia.Piran1.ages,Gouldia.Piran1.5depths))), 
col="gray81", range=1.5, boxwex=3, horizontal=TRUE, ylab="Sediment depth (cm)", xlab="Calendar date since 2013 (AD)", main="Piran 1", xlim=c(-150,0), 
ylim=c(0,12000), cex.main=0.9, outline=F, las=2, cex.axis=0.7, xaxt="n", cex.lab=0.8, frame=F)
axis(1, at=seq(0,10000,by=2500), labels=seq(0,10000,by=2500), cex.axis=0.7)
segments(x0=-3000, y0=-c(8,16,35,65), x1 = 0, y1 = -c(8,16,35,65),lty=2)
rect(xleft=10800, ybottom=-150, xright=12500, ytop=-65, col="gray", border="white")
rect(xleft=10800, ybottom=-65, xright=12500, ytop=-35, col="gray", border="white")
rect(xleft=10800, ybottom=-35, xright=12500, ytop=-8, col="gray", border="white")
rect(xleft=10800, ybottom=-8, xright=12500, ytop=0, col="gray", border="white")
text(11750,-100,labels="sandy mud",srt=90, cex=0.8)
text(11750,-50,labels="sand",srt=90, cex=0.8)
text(11750,-20,labels="shell b.",srt=90, cex=0.8)
text(11750,-4,labels="sand",srt=0, cex=0.8)
points(Gouldia.Piran1.carbonages, -Gouldia.Piran1.carbonage.depths, pch=21, bg="gray81",cex=1.4)
points(Ostrea.Piran1.age, -Ostrea.Piran1.depth, pch=22, bg="black",cex=1.4)
points(Arca.Piran1.age, -Arca.Piran1.depth, pch=22, bg="white",cex=1.4)
legend(x=7000,y=-115, legend=c("Arca","Ostrea","Gouldia","Corbula"), pch=c(22,22,21,21), pt.bg=c("white","black","gray81","gray51"))

out2=boxplot(split(Gouldia.Piran2.ages,Gouldia.Piran2.5depths), at=-as.numeric(names(split(Gouldia.Piran2.ages,Gouldia.Piran2.5depths))), 
col="gray81", range=1.5, boxwex=3, horizontal=TRUE, ylab="", xlab="Calendar date since 2013 (AD)", main="Piran 2", xlim=c(-150,0), 
ylim=c(0,12000), cex.main=0.9, outline=F, las=2, cex.axis=0.7, xaxt="n", cex.lab=0.8, frame=F)
axis(1, at=seq(0,10000,by=2500), labels=seq(0,10000,by=2500), cex.axis=0.7)
segments(x0=-3000, y0=-c(8,16,35,65), x1 = 0, y1 = -c(8,16,35,65),lty=2)
out3=boxplot(split(Corbula.Piran.ages,Corbula.Piran.5depths), at=-as.numeric(names(split(Corbula.Piran.ages,Corbula.Piran.5depths)))+2, add=T, 
col="gray51", range=1.5, boxwex=3,  horizontal=TRUE, yaxt="n", outline=F, las=2, cex.axis=0.7, xaxt="n",cex.lab=0.8, frame=F)
points(Gouldia.Piran2.carbonages, -Gouldia.Piran2.carbonage.depths, pch=21, bg="gray81",cex=1.4)
points(Corbula.Piran2.carbonages, -Corbula.Piran2.carbonage.depths, pch=21, bg="gray51",cex=1.4)
points(Arca.Piran2.age, -Arca.Piran2.depth, pch=22, bg="white",cex=1.4)

rect(xleft=10800, ybottom=-150, xright=12500, ytop=-65, col="gray", border="white")
rect(xleft=10800, ybottom=-65, xright=12500, ytop=-35, col="gray", border="white")
rect(xleft=10800, ybottom=-35, xright=12500, ytop=-8, col="gray", border="white")
rect(xleft=10800, ybottom=-8, xright=12500, ytop=0, col="gray", border="white")
text(11750,-100,labels="sandy mud",srt=90, cex=0.8)
text(11750,-50,labels="sand",srt=90, cex=0.8)
text(11750,-20,labels="shell b.",srt=90, cex=0.8)
text(11750,-4,labels="sand",srt=0, cex=0.8)

################################################################################
#DEPTH CHANGES - PIRAN 1
################################################################################
age=Gouldia.Piran1.ages
Coredepth=Gouldia.Piran1.5depths
Corestation=as.vector(Gouldia.Piran1.station)

LIST=length(split(age, -Coredepth))
par(mfcol=c(LIST,2))
par(mfcol=c(11,3))
par(mar=c(1,4,1,2))
samplesizes1=numeric()
depths1=numeric()
IQRrange1=numeric()
fullrange1=numeric()
medianage1=numeric()
meanage1=numeric()
range951=numeric()
minage1=numeric()
maxage1=numeric()
skewness1=numeric()
station1=numeric()
interval1=numeric()
totalrange1=numeric()
skew1=numeric()
correction=numeric()
obs.UCI95range1=numeric(); obs.LCI95range1=numeric()
obs.UCIrange1=numeric(); obs.LCIrange1=numeric()
obs.UCImedian1=numeric(); obs.LCImedian1=numeric()
obs.UCImean1=numeric(); obs.LCImean1=numeric()
obs.UCIIQR1=numeric(); obs.LCIIQR1=numeric()
obs.UCIskewness1=numeric(); obs.LCIskewness1=numeric()

for (i in 1:(LIST-1)) {
samplesizes1[i]=length(unlist(split(age, Coredepth)[i]))
depths1[i]=names(split(age, Coredepth)[i])
interval1[i]=as.character(as.numeric(depths1[i]))
hist(unlist(split(age, Coredepth)[i]), xlab="Age (years)", ylab="",  breaks=seq(0,13000,by=500), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="gray", ylim=c(0,16),main=paste(interval1[i],"cm","(n =",samplesizes1[i],")"))
IQRrange1[i]=IQR(unlist(split(age, Coredepth)[i]))
fullrange1[i]=diff(range(unlist(split(age, Coredepth)[i])))
medianage1[i]=median(unlist(split(age, Coredepth)[i]))
meanage1[i]=mean(unlist(split(age, Coredepth)[i]))
minage1[i]=min(unlist(split(age, Coredepth)[i]))
maxage1[i]=max(unlist(split(age, Coredepth)[i]))
range951[i]=diff(quantile(unlist(split(age, Coredepth)[i]),c(0.025,0.975)))
totalrange1[i]=diff(range(unlist(split(age, Coredepth)[i])))
skew1[i]=skewness(unlist(split(age, Coredepth)[i]))
station1[i]=unlist(split(Corestation, Coredepth)[i])[1]
abline(v=median(unlist(split(age, Coredepth)[i])),lwd=2)
tempage=unlist(split(age, Coredepth)[i])

bootmedian<-numeric(); bootmean<-numeric(); bootrange<-numeric();  bootskew<-numeric(); bootIQR<-numeric();  boot95range<-numeric()
for (j in 1:1000) {
largesample<-sample(tempage,1000, replace=T)
temp=sample(largesample,length(tempage), replace=F)
bootmedian[j]=median(temp)
bootmean[j]=mean(temp)
bootskew[j]=skewness(temp)
bootrange[j]=diff(range(temp))
bootIQR[j]=IQR(temp)
temp1=quantile(temp, c(0.025))
temp2=quantile(temp, c(0.975))
boot95range[j]=temp2-temp1
}
correction[i]=mean(bootrange)-totalrange1[i]
if (correction[i] < 0) {bootrange=bootrange-correction[i]}
obs.UCI95range1[i]=quantile(boot95range,c(0.975));  obs.LCI95range1[i]=quantile(boot95range,c(0.025))
obs.UCIrange1[i]=quantile(bootrange,c(0.975)); obs.LCIrange1[i]=quantile(bootrange,c(0.025))
obs.UCImedian1[i]=quantile(bootmedian,c(0.975)); obs.LCImedian1[i]=quantile(bootmedian,c(0.025))
obs.UCImean1[i]=quantile(bootmean,c(0.975)); obs.LCImean1[i]=quantile(bootmean,c(0.025))
obs.UCIIQR1[i]=quantile(bootIQR,c(0.975));  obs.LCIIQR1[i]=quantile(bootIQR,c(0.025))
obs.UCIskewness1[i]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness1[i]=quantile(bootskew,c(0.025), na.rm=T)
print(i)
}

samplesizes1[LIST]=length(unlist(split(age, Coredepth)[LIST]))
IQRrange1[LIST]=IQR(unlist(split(age, Coredepth)[LIST]))
fullrange1[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
medianage1[LIST]=median(unlist(split(age, Coredepth)[LIST]))
meanage1[LIST]=mean(unlist(split(age, Coredepth)[LIST]))
minage1[LIST]=min(unlist(split(age, Coredepth)[LIST]))
maxage1[LIST]=max(unlist(split(age, Coredepth)[LIST]))
range951[LIST]=diff(quantile(unlist(split(age, Coredepth)[LIST]),c(0.025,0.975)))
totalrange1[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
skew1[LIST]=skewness(unlist(split(age, Coredepth)[LIST]))
depths1[LIST]=names(split(age, Coredepth)[LIST])
station1[LIST]=unlist(split( Corestation, Coredepth)[LIST])[1]
interval1[LIST]=as.character(as.numeric(depths1[LIST]))
hist(unlist(split(age, Coredepth)[LIST]), xlab="Age (years)", ylab="", breaks=seq(0,13000,by=500), cex.main=0.8, cex.axis=0.8,
col="gray", ylim=c(0,16),main=paste(interval1[LIST],"cm","(n =",samplesizes1[LIST],")"))
abline(v=median(unlist(split(age, Coredepth)[LIST])),lwd=2)

bootmedian<-numeric(); bootmean<-numeric(); bootrange<-numeric();  bootskew<-numeric(); bootIQR<-numeric();  boot95range<-numeric()
for (j in 1:1000) {
largesample<-sample(tempage,1000, replace=T)
temp=sample(largesample,length(tempage), replace=F)
bootmedian[j]=median(temp)
bootmean[j]=mean(temp)
bootskew[j]=skewness(temp)
bootrange[j]=diff(range(temp))
bootIQR[j]=IQR(temp)
temp1=quantile(temp, c(0.025))
temp2=quantile(temp, c(0.975))
boot95range[j]=temp2-temp1
}
correction[LIST]=mean(bootrange)-totalrange1[LIST]
if (correction[LIST] < 0) {bootrange=bootrange-correction[LIST]}
obs.UCI95range1[LIST]=quantile(boot95range,c(0.975));  obs.LCI95range1[LIST]=quantile(boot95range,c(0.025))
obs.UCIrange1[LIST]=quantile(bootrange,c(0.975)); obs.LCIrange1[LIST]=quantile(bootrange,c(0.025))
obs.UCImedian1[LIST]=quantile(bootmedian,c(0.975)); obs.LCImedian1[LIST]=quantile(bootmedian,c(0.025))
obs.UCImean1[LIST]=quantile(bootmean,c(0.975)); obs.LCImean1[LIST]=quantile(bootmean,c(0.025))
obs.UCIIQR1[LIST]=quantile(bootIQR,c(0.975));  obs.LCIIQR1[LIST]=quantile(bootIQR,c(0.025))
obs.UCIskewness1[LIST]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness1[LIST]=quantile(bootskew,c(0.025), na.rm=T)

plot.new()
plot.new()

#out1=cbind(rep("Gouldia minima",length(station)),station,interval,depths,samplesizes, minage, medianage,meanage,maxage,IQRrange,range95,fullrange)
#colnames(out1)=c("Species","Station","Depth interval (cm)","Max  depth (cm)", "N", "Min age","Median age","Mean age","Max age","IQR range","95% range","Range")

################################################################################
#DEPTH CHANGES - PIRAN 2
################################################################################
age=Gouldia.Piran2.ages
Coredepth=Gouldia.Piran2.5depths
Corestation=as.vector(Gouldia.Piran2.station)

LIST=length(split(age, -Coredepth))
samplesizes2=numeric()
depths2=numeric()
IQRrange2=numeric()
fullrange2=numeric()
medianage2=numeric()
meanage2=numeric()
range952=numeric()
minage2=numeric()
maxage2=numeric()
skewness2=numeric()
station2=numeric()
interval2=numeric()
totalrange2=numeric()
skew2=numeric()
correction=numeric()
obs.UCI95range2=numeric(); obs.LCI95range2=numeric()
obs.UCIrange2=numeric(); obs.LCIrange2=numeric()
obs.UCImedian2=numeric(); obs.LCImedian2=numeric()
obs.UCImean2=numeric(); obs.LCImean2=numeric()
obs.UCIIQR2=numeric(); obs.LCIIQR2=numeric()
obs.UCIskewness2=numeric(); obs.LCIskewness2=numeric()

for (i in 1:(LIST-1)) {
samplesizes2[i]=length(unlist(split(age, Coredepth)[i]))
depths2[i]=names(split(age, Coredepth)[i])
interval2[i]=as.character(as.numeric(depths2[i]))
hist(unlist(split(age, Coredepth)[i]), xlab="Age (years)", ylab="", breaks=seq(0,13000,by=500), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="gray", ylim=c(0,16),main=paste(interval2[i],"cm","(n =",samplesizes2[i],")"))
IQRrange2[i]=IQR(unlist(split(age, Coredepth)[i]))
fullrange2[i]=diff(range(unlist(split(age, Coredepth)[i])))
medianage2[i]=median(unlist(split(age, Coredepth)[i]))
meanage2[i]=mean(unlist(split(age, Coredepth)[i]))
minage2[i]=min(unlist(split(age, Coredepth)[i]))
maxage2[i]=max(unlist(split(age, Coredepth)[i]))
range952[i]=diff(quantile(unlist(split(age, Coredepth)[i]),c(0.025,0.975)))
totalrange2[i]=diff(range(unlist(split(age, Coredepth)[i])))
skew2[i]=skewness(unlist(split(age, Coredepth)[i]))
station2[i]=unlist(split(Corestation, Coredepth)[i])[1]
abline(v=median(unlist(split(age, Coredepth)[i])),lwd=2)
tempage=unlist(split(age, Coredepth)[i])

bootmedian<-numeric(); bootmean<-numeric(); bootrange<-numeric();  bootskew<-numeric(); bootIQR<-numeric();  boot95range<-numeric()
for (j in 1:1000) {
largesample<-sample(tempage,1000, replace=T)
temp=sample(largesample,length(tempage), replace=F)
bootmedian[j]=median(temp)
bootmean[j]=mean(temp)
bootskew[j]=skewness(temp)
bootrange[j]=diff(range(temp))
bootIQR[j]=IQR(temp)
temp1=quantile(temp, c(0.025))
temp2=quantile(temp, c(0.975))
boot95range[j]=temp2-temp1
}
correction[i]=mean(bootrange)-totalrange2[i]
if (correction[i] < 0) {bootrange=bootrange-correction[i]}
obs.UCI95range2[i]=quantile(boot95range,c(0.975));  obs.LCI95range2[i]=quantile(boot95range,c(0.025))
obs.UCIrange2[i]=quantile(bootrange,c(0.975)); obs.LCIrange2[i]=quantile(bootrange,c(0.025))
obs.UCImedian2[i]=quantile(bootmedian,c(0.975)); obs.LCImedian2[i]=quantile(bootmedian,c(0.025))
obs.UCImean2[i]=quantile(bootmean,c(0.975)); obs.LCImean2[i]=quantile(bootmean,c(0.025))
obs.UCIIQR2[i]=quantile(bootIQR,c(0.975));  obs.LCIIQR2[i]=quantile(bootIQR,c(0.025))
obs.UCIskewness2[i]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness2[i]=quantile(bootskew,c(0.025), na.rm=T)
print(i)
}

samplesizes2[LIST]=length(unlist(split(age, Coredepth)[LIST]))
IQRrange2[LIST]=IQR(unlist(split(age, Coredepth)[LIST]))
fullrange2[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
medianage2[LIST]=median(unlist(split(age, Coredepth)[LIST]))
meanage2[LIST]=median(unlist(split(age, Coredepth)[LIST]))
minage2[LIST]=min(unlist(split(age, Coredepth)[LIST]))
maxage2[LIST]=max(unlist(split(age, Coredepth)[LIST]))
range952[LIST]=diff(quantile(unlist(split(age, Coredepth)[LIST]),c(0.025,0.975)))
totalrange2[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
skew2[LIST]=skewness(unlist(split(age, Coredepth)[LIST]))
depths2[LIST]=names(split(age, Coredepth)[LIST])
station2[LIST]=unlist(split( Corestation, Coredepth)[LIST])[1]
interval2[LIST]=as.character(paste(as.numeric(depths2[LIST])-10,"-",depths2[LIST], sep=""))
interval2[LIST]=as.character(as.numeric(depths2[LIST]))
hist(unlist(split(age, Coredepth)[LIST]), xlab="Age (years)", ylab="", breaks=seq(0,13000,by=500), cex.main=0.8, cex.axis=0.8,
col="gray", ylim=c(0,16),main=paste(interval2[LIST],"cm","(n =",samplesizes2[LIST],")"))
abline(v=median(unlist(split(age, Coredepth)[LIST])),lwd=2)

bootmedian<-numeric(); bootmean<-numeric(); bootrange<-numeric();  bootskew<-numeric(); bootIQR<-numeric();  boot95range<-numeric()
for (j in 1:1000) {
largesample<-sample(tempage,1000, replace=T)
temp=sample(largesample,length(tempage), replace=F)
bootmedian[j]=median(temp)
bootmean[j]=mean(temp)
bootskew[j]=skewness(temp)
bootrange[j]=diff(range(temp))
bootIQR[j]=IQR(temp)
temp1=quantile(temp, c(0.025))
temp2=quantile(temp, c(0.975))
boot95range[j]=temp2-temp1
}
correction[LIST]=mean(bootrange)-totalrange2[LIST]
if (correction[LIST] < 0) {bootrange=bootrange-correction[LIST]}
obs.UCI95range2[LIST]=quantile(boot95range,c(0.975));  obs.LCI95range2[LIST]=quantile(boot95range,c(0.025))
obs.UCIrange2[LIST]=quantile(bootrange,c(0.975)); obs.LCIrange2[LIST]=quantile(bootrange,c(0.025))
obs.UCImedian2[LIST]=quantile(bootmedian,c(0.975)); obs.LCImedian2[LIST]=quantile(bootmedian,c(0.025))
obs.UCImean2[LIST]=quantile(bootmean,c(0.975)); obs.LCImean2[LIST]=quantile(bootmean,c(0.025))
obs.UCIIQR2[LIST]=quantile(bootIQR,c(0.975));  obs.LCIIQR2[LIST]=quantile(bootIQR,c(0.025))
obs.UCIskewness2[LIST]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness2[LIST]=quantile(bootskew,c(0.025), na.rm=T)

###############
#CORBULA PIRAN#
###############
age=Corbula.Piran.ages
Coredepth=Corbula.Piran.5depths
Corestation=as.vector(Gouldia.Piran2.station)

LIST=length(split(age, -Coredepth))
#par(mfcol=c(LIST,2))
#par(mar=c(1,5,1,5))
samplesizes3=numeric()
depths3=numeric()
IQRrange3=numeric()
fullrange3=numeric()
medianage3=numeric()
meanage3=numeric()
range953=numeric()
minage3=numeric()
maxage3=numeric()
skewness3=numeric()
station3=numeric()
interval3=numeric()
totalrange3=numeric()
skew3=numeric()
correction=numeric()
obs.UCI95range3=numeric(); obs.LCI95range3=numeric()
obs.UCIrange3=numeric(); obs.LCIrange3=numeric()
obs.UCImedian3=numeric(); obs.LCImedian3=numeric()
obs.UCImean3=numeric(); obs.LCImean3=numeric()
obs.UCIIQR3=numeric(); obs.LCIIQR3=numeric()
obs.UCIskewness3=numeric(); obs.LCIskewness3=numeric()

for (i in 1:(LIST-1)) {
samplesizes3[i]=length(unlist(split(age, Coredepth)[i]))
depths3[i]=names(split(age, Coredepth)[i])
interval3[i]=as.character(as.numeric(depths3[i]))
hist(unlist(split(age, Coredepth)[i]), xlab="Age (years)", ylab="", breaks=seq(0,14000,by=500), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="white", ylim=c(0,16),main=paste(interval3[i],"cm","(n =",samplesizes3[i],")"))
IQRrange3[i]=IQR(unlist(split(age, Coredepth)[i]))
fullrange3[i]=diff(range(unlist(split(age, Coredepth)[i])))
medianage3[i]=median(unlist(split(age, Coredepth)[i]))
meanage3[i]=mean(unlist(split(age, Coredepth)[i]))
minage3[i]=min(unlist(split(age, Coredepth)[i]))
maxage3[i]=max(unlist(split(age, Coredepth)[i]))
range953[i]=diff(quantile(unlist(split(age, Coredepth)[i]),c(0.025,0.975)))
totalrange3[i]=diff(range(unlist(split(age, Coredepth)[i])))
skew3[i]=skewness(unlist(split(age, Coredepth)[i]))
#station3[i]=unlist(split(Corestation, Coredepth)[i])[1]
abline(v=median(unlist(split(age, Coredepth)[i])),lwd=2)
tempage=unlist(split(age, Coredepth)[i])

bootmedian<-numeric(); bootmean<-numeric(); bootrange<-numeric();  bootskew<-numeric(); bootIQR<-numeric();  boot95range<-numeric()
for (j in 1:1000) {
largesample<-sample(tempage,1000, replace=T)
temp=sample(largesample,length(tempage), replace=F)
bootmedian[j]=median(temp)
bootmean[j]=mean(temp)
bootskew[j]=skewness(temp)
bootrange[j]=diff(range(temp))
bootIQR[j]=IQR(temp)
temp1=quantile(temp, c(0.025))
temp2=quantile(temp, c(0.975))
boot95range[j]=temp2-temp1
}
correction[i]=mean(bootrange)-totalrange3[i]
if (correction[i] < 0) {bootrange=bootrange-correction[i]}
obs.UCI95range3[i]=quantile(boot95range,c(0.975));  obs.LCI95range3[i]=quantile(boot95range,c(0.025))
obs.UCIrange3[i]=quantile(bootrange,c(0.975)); obs.LCIrange3[i]=quantile(bootrange,c(0.025))
obs.UCImedian3[i]=quantile(bootmedian,c(0.975)); obs.LCImedian3[i]=quantile(bootmedian,c(0.025))
obs.UCImean3[i]=quantile(bootmean,c(0.975)); obs.LCImean3[i]=quantile(bootmean,c(0.025))
obs.UCIIQR3[i]=quantile(bootIQR,c(0.975));  obs.LCIIQR3[i]=quantile(bootIQR,c(0.025))
obs.UCIskewness3[i]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness3[i]=quantile(bootskew,c(0.025), na.rm=T)
print(i)
}

samplesizes3[LIST]=length(unlist(split(age, Coredepth)[LIST]))
IQRrange3[LIST]=IQR(unlist(split(age, Coredepth)[LIST]))
fullrange3[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
medianage3[LIST]=median(unlist(split(age, Coredepth)[LIST]))
meanage3[LIST]=median(unlist(split(age, Coredepth)[LIST]))
minage3[LIST]=min(unlist(split(age, Coredepth)[LIST]))
maxage3[LIST]=max(unlist(split(age, Coredepth)[LIST]))
range953[LIST]=diff(quantile(unlist(split(age, Coredepth)[LIST]),c(0.025,0.975)))
totalrange3[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
skew3[LIST]=skewness(unlist(split(age, Coredepth)[LIST]))
depths3[LIST]=names(split(age, Coredepth)[LIST])
#station3[LIST]=unlist(split( Corestation, Coredepth)[LIST])[1]
interval3[LIST]=as.character(as.numeric(depths3[LIST]))
hist(unlist(split(age, Coredepth)[LIST]), xlab="Age (years)", ylab="", breaks=seq(0,12000,by=500), cex.main=0.8, cex.axis=0.8,
col="white", ylim=c(0,16),main=paste(interval3[LIST],"cm","(n =",samplesizes3[LIST],")"))
abline(v=median(unlist(split(age, Coredepth)[LIST])),lwd=2)

bootmedian<-numeric(); bootmean<-numeric(); bootrange<-numeric();  bootskew<-numeric(); bootIQR<-numeric();  boot95range<-numeric()
for (j in 1:1000) {
largesample<-sample(tempage,1000, replace=T)
temp=sample(largesample,length(tempage), replace=F)
bootmedian[j]=median(temp)
bootmean[j]=mean(temp)
bootskew[j]=skewness(temp)
bootrange[j]=diff(range(temp))
bootIQR[j]=IQR(temp)
temp1=quantile(temp, c(0.025))
temp2=quantile(temp, c(0.975))
boot95range[j]=temp2-temp1
}
correction[LIST]=mean(bootrange)-totalrange3[LIST]
if (correction[LIST] < 0) {bootrange=bootrange-correction[LIST]}
obs.UCI95range3[LIST]=quantile(boot95range,c(0.975));  obs.LCI95range3[LIST]=quantile(boot95range,c(0.025))
obs.UCIrange3[LIST]=quantile(bootrange,c(0.975)); obs.LCIrange3[LIST]=quantile(bootrange,c(0.025))
obs.UCImedian3[LIST]=quantile(bootmedian,c(0.975)); obs.LCImedian3[LIST]=quantile(bootmedian,c(0.025))
obs.UCImean3[LIST]=quantile(bootmean,c(0.975)); obs.LCImean3[LIST]=quantile(bootmean,c(0.025))
obs.UCIIQR3[LIST]=quantile(bootIQR,c(0.975));  obs.LCIIQR3[LIST]=quantile(bootIQR,c(0.025))
obs.UCIskewness3[LIST]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness3[LIST]=quantile(bootskew,c(0.025), na.rm=T)

################################################################################
#DOWNCORE TRENDS IN RAW TIME AVERAGING
################################################################################
depth1=as.numeric(depths1)
depth2=as.numeric(depths2)
depth3=as.numeric(depths3)

par(mfrow=c(2,2))
plot(range951, -depth1,xlim=c(1,10000),log="x", xlab="Time averaging (years)", ylab="Sediment depth (cm)", main="Piran 1", cex.main=0.9)
segments(x0=obs.UCI95range1,x1=obs.LCI95range1,y0=-depth1, y1=-depth1)
points(range951, -depth1, pch=21, bg="gray")
segments(x0=obs.UCIIQR1,x1=obs.LCIIQR1,y0=-depth1, y1=-depth1)
points(IQRrange1, -depth1, pch=16)
plot(skew1, -depth1,xlim=c(-3,3), xlab="Skewness", ylab="Sediment depth (cm)", pch=16, main="Piran 1", cex.main=0.9)
segments(x0=obs.UCIskewness1,x1=obs.LCIskewness1,y0=-depth1, y1=-depth1)

plot(range952, -depth2,xlim=c(1,10000),log="x", xlab="Time averaging (years)", ylab="Sediment depth (cm)", main="Piran 2", cex.main=0.9)
segments(x0=obs.UCI95range2,x1=obs.LCI95range2,y0=-depth2, y1=-depth2)
points(range952, -depth2, pch=21, bg="gray")
segments(x0=obs.UCIIQR2,x1=obs.LCIIQR2,y0=-depth2, y1=-depth2)
points(IQRrange2, -depth2, pch=16)
plot(skew2, -depth2,xlim=c(-3,3), xlab="Skewness", ylab="Sediment depth (cm)", pch=16, main="Piran 2", cex.main=0.9)
segments(x0=obs.UCIskewness2,x1=obs.LCIskewness2,y0=-depth2, y1=-depth2)


##########################################################################################################
#FIGURE SURFACE AND CORE AGE-FREQUENCY DISTRIBUTIONS AT PIRAN
##########################################################################################################
par(mfrow=c(3,3))
par(mar=c(4,4,2,1))
####################
#TOP CORE
####################
BINS=100
out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<4 & Gouldia.Piran1.ages<3000], ylab="N specimens", col="gray", breaks=seq(0,3000,by=BINS), xlab="", main="Top-core 12 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<4]), col="gray31")
out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<4 & Gouldia.Piran2.ages<3000], ylab="", col="gray", breaks=seq(0,3000,by=BINS), xlab="", main="Top-core 12 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<4]), col="gray31")
out=hist(Corbula.Piran.ages[Corbula.Piran.5depths<4 & Corbula.Piran.ages<3000], ylab="", col="gray", breaks=seq(0,3000,by=BINS), xlab="", main="Top-core 12 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths<4]), col="gray31")

out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<12 & Gouldia.Piran1.ages<3000], ylab="Number of specimens", col="gray", breaks=seq(0,3000,by=BINS), xlab="", main="Top-core 12 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<12]), col="gray31")
out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<12 & Gouldia.Piran2.ages<3000], ylab="", col="gray", breaks=seq(0,3000,by=BINS), xlab="", main="Top-core 12 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<12]), col="gray31")
out=hist(Corbula.Piran.ages[Corbula.Piran.5depths<12 & Corbula.Piran.ages<3000], ylab="", col="gray", breaks=seq(0,3000,by=BINS), xlab="", main="Top-core 12 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths<12]), col="gray31")

####################
#WHOLE CORE
####################
hist(Gouldia.Piran1.ages, breaks=seq(0,13000,by=200), 
xlab="Postmortem age (years)", col="gray", main="Piran 1")
hist(Gouldia.Piran2.ages, breaks=seq(0,13000,by=200),
xlab="Postmortem age (years)", col="gray", main="Piran 2")
hist(Corbula.Piran.ages, breaks=seq(0,13000,by=200),
xlab="Postmortem age (years)", col="gray", main="Piran 2")


##########################################################################################################
#FIGURE AGE-FREQUENCY DISTRIBUTIONS IN INTERVALS AT PIRAN
##########################################################################################################
par(mfrow=c(5,3))
par(mar=c(2,4,2,1))
####################
#TOP CORE
####################
BINS=200
out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8], ylab="N specimens", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="Top-core 8 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8]), col="gray31")
skewness(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8])
IQR(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8])
diff(quantile(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8], c(0.025,0.975)))
mlv(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8], method = "hrm", bw=0.1)
which(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<8]<112)

out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="Top-core 8 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8]), col="gray31")
skewness(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8])
IQR(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8])
diff(quantile(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8], c(0.025,0.975)))
mlv(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8], method = "hrm", bw=0.1)
which(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<8]<112)

out=hist(Corbula.Piran.ages[Corbula.Piran.5depths<8], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="Top-core 8 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths<8]), col="gray31")
skewness(Corbula.Piran.ages[Corbula.Piran.5depths<8])
IQR(Corbula.Piran.ages[Corbula.Piran.5depths<8])
diff(quantile(Corbula.Piran.ages[Corbula.Piran.5depths<8], c(0.025,0.975)))
mlv(Corbula.Piran.ages[Corbula.Piran.5depths<8], method = "hrm", bw=0.1)
length(which(Corbula.Piran.ages[Corbula.Piran.5depths<8]<112))/length(Corbula.Piran.ages[Corbula.Piran.5depths<8])

####################
#SHELL BED
####################
out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>8 & Gouldia.Piran1.5depths<16], ylab="N specimens", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="shell bed-8-16 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>8 & Gouldia.Piran1.5depths<16]), col="gray31")
skewness(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>8 & Gouldia.Piran1.5depths<16])
IQR(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>8 & Gouldia.Piran1.5depths<16])
diff(quantile(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>8 & Gouldia.Piran1.5depths<16], c(0.025,0.975)))

out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>8 & Gouldia.Piran2.5depths<16], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="shell bed-8-16 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>8 & Gouldia.Piran2.5depths<16]), col="gray31")
skewness(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>8 & Gouldia.Piran2.5depths<16])
IQR(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>8 & Gouldia.Piran2.5depths<16])
diff(quantile(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>8 & Gouldia.Piran2.5depths<16], c(0.025,0.975)))

out=hist(Corbula.Piran.ages[Corbula.Piran.5depths>8 & Corbula.Piran.5depths<16], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="shell bed-8-16 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths>8 & Corbula.Piran.5depths<16]), col="gray31")
skewness(Corbula.Piran.ages[Corbula.Piran.5depths>8 & Corbula.Piran.5depths<16])
IQR(Corbula.Piran.ages[Corbula.Piran.5depths>8 & Corbula.Piran.5depths<16])
diff(quantile(Corbula.Piran.ages[Corbula.Piran.5depths>8 & Corbula.Piran.5depths<16], c(0.025,0.975)))

####################
#SHELL BED
####################
out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>16 & Gouldia.Piran1.5depths<30], ylab="N specimens", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="shell bed-16-35 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>16 & Gouldia.Piran1.5depths<30]), col="gray31")
skewness(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>16 & Gouldia.Piran1.5depths<30])
IQR(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>16 & Gouldia.Piran1.5depths<30])
diff(quantile(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>16 & Gouldia.Piran1.5depths<30], c(0.025,0.975)))

out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>16 & Gouldia.Piran2.5depths<30], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="shell bed-16-35 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>16 & Gouldia.Piran2.5depths<30]), col="gray31")
skewness(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>16 & Gouldia.Piran2.5depths<30])
IQR(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>16 & Gouldia.Piran2.5depths<30])
diff(quantile(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>16 & Gouldia.Piran2.5depths<30], c(0.025,0.975)))

out=hist(Corbula.Piran.ages[Corbula.Piran.5depths>16 & Corbula.Piran.5depths<30], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", main="shell bed-16-35 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths>6 & Corbula.Piran.5depths<30]), col="gray31")
skewness(Corbula.Piran.ages[Corbula.Piran.5depths>16 & Corbula.Piran.5depths<30])
IQR(Corbula.Piran.ages[Corbula.Piran.5depths>16 & Corbula.Piran.5depths<30])
diff(quantile(Corbula.Piran.ages[Corbula.Piran.5depths>16 & Corbula.Piran.5depths<30], c(0.025,0.975)))

####################
#EARLY HST
####################
out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>30 & Gouldia.Piran1.5depths<65], ylab="N specimens", col="gray", breaks=seq(0,13000,by=BINS), xlab="", 
main="Muddy sand-35-65 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>30 & Gouldia.Piran1.5depths<65]), col="gray31")
skewness(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>30 & Gouldia.Piran1.5depths<65])
IQR(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>30 & Gouldia.Piran1.5depths<65])
diff(quantile(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>30 & Gouldia.Piran1.5depths<65], c(0.025,0.975)))

out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>30 & Gouldia.Piran2.5depths<65], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", 
main="Muddy sand-35-65 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>30 & Gouldia.Piran2.5depths<65]), col="gray31")
skewness(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>30 & Gouldia.Piran2.5depths<65])
IQR(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>30 & Gouldia.Piran2.5depths<65])
diff(quantile(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>30 & Gouldia.Piran2.5depths<65], c(0.025,0.975)))

out=hist(Corbula.Piran.ages[Corbula.Piran.5depths>30 & Corbula.Piran.5depths<65], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", 
main="Muddy sand-35-65 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths>30 & Corbula.Piran.5depths<65]), col="gray31")
skewness(Corbula.Piran.ages[Corbula.Piran.5depths>30 & Corbula.Piran.5depths<65])
IQR(Corbula.Piran.ages[Corbula.Piran.5depths>30 & Corbula.Piran.5depths<65])
diff(quantile(Corbula.Piran.ages[Corbula.Piran.5depths>30 & Corbula.Piran.5depths<65], c(0.025,0.975)))

####################
#TST
####################
out=hist(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>65], ylab="N specimens", col="gray", breaks=seq(0,13000,by=BINS), xlab="", 
main="Sandy mud-65-150 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>65]), col="gray31")
skewness(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>65])
IQR(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>65])
diff(quantile(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>65], c(0.025,0.975)))

out=hist(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>65], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", 
main="Sandy mud-65-150 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>65]), col="gray31")
skewness(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>65])
IQR(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>65])
diff(quantile(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>65], c(0.025,0.975)))

out=hist(Corbula.Piran.ages[Corbula.Piran.5depths>65], ylab="", col="gray", breaks=seq(0,13000,by=BINS), xlab="", 
main="Sandy mud-65-150 cm", ylim=c(0,20), cex.main=1)
abline(v=median(Corbula.Piran.ages[Corbula.Piran.5depths>65]), col="gray31")
skewness(Corbula.Piran.ages[Corbula.Piran.5depths>65])
IQR(Corbula.Piran.ages[Corbula.Piran.5depths>65])
diff(quantile(Corbula.Piran.ages[Corbula.Piran.5depths>65], c(0.025,0.975)))


###############################################################################################################
#TIME AVERAGING PER UNIT
###############################################################################################################
####################
#GOULDIA AT PIRAN 1#
####################
Gouldia.Piran1.layers=Gouldia.Piran1.5depths
Gouldia.Piran1.layers[Gouldia.Piran1.5depths<8]=1
Gouldia.Piran1.layers[Gouldia.Piran1.5depths>8 & Gouldia.Piran1.5depths<16]=2
Gouldia.Piran1.layers[Gouldia.Piran1.5depths>16 & Gouldia.Piran1.5depths<35]=3
Gouldia.Piran1.layers[Gouldia.Piran1.5depths>35 & Gouldia.Piran1.5depths<65]=4
Gouldia.Piran1.layers[Gouldia.Piran1.5depths>65]=5
Gouldia.Piran1.IQR.layers=tapply(Gouldia.Piran1.ages,Gouldia.Piran1.layers, IQR)
Gouldia.Piran1.95.layers=tapply(Gouldia.Piran1.ages,Gouldia.Piran1.layers, function(x) {diff(quantile(x,c(0.025, 0.975)))})
Gouldia.Piran1.IQR.CI.layers=array(NA, dim=c(5,3))
Gouldia.Piran1.IQR.CI.layers[1,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==1])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran1.IQR.CI.layers[2,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==2])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran1.IQR.CI.layers[3,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==3])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran1.IQR.CI.layers[4,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==4])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran1.IQR.CI.layers[5,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==5])$resampleiqr, c(0.025,0.5,0.975))

Gouldia.Piran1.95.CI.layers=array(NA, dim=c(5,3))
Gouldia.Piran1.95.CI.layers[1,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==1])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran1.95.CI.layers[2,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==2])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran1.95.CI.layers[3,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==3])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran1.95.CI.layers[4,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==4])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran1.95.CI.layers[5,]=quantile(bootci.ext(Gouldia.Piran1.ages[Gouldia.Piran1.layers==5])$resamplerange, c(0.025,0.5,0.975))

SIM.meaniqr.x.age=array(NA, dim=c(5,100));  SIM.mean95.x.age=array(NA, dim=c(5,100))
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:5) {
if (i == 1) {ages=Gouldia.Piran1.ages[Gouldia.Piran1.layers==1]}
if (i == 2) {ages=Gouldia.Piran1.ages[Gouldia.Piran1.layers==2]}
if (i == 3) {ages=Gouldia.Piran1.ages[Gouldia.Piran1.layers==3]}
if (i == 4) {ages=Gouldia.Piran1.ages[Gouldia.Piran1.layers==4]}
if (i == 5) {ages=Gouldia.Piran1.ages[Gouldia.Piran1.layers==5]}
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
temp=rgamma(100,shape=ages[x]/exp(gamma.ln.d),scale=exp(gamma.ln.d))
temp=rtrunc(100, spec="gamma", shape=ages[x]/exp(gamma.ln.d),scale=exp(gamma.ln.d),a=0, b=10000)
x.age[x]=IQR(temp)
x.age.95[x]=diff(quantile(temp, c(0.025,0.975)))
}
meaniqr.x.age[i]=mean(x.age); mean95.x.age[i]=mean(x.age.95);
}
SIM.meaniqr.x.age[,j]=meaniqr.x.age;  SIM.mean95.x.age[,j]=mean95.x.age
}

SIM.mean.95=apply(SIM.mean95.x.age, 1, mean)
SIM.mean.95LCI=apply(SIM.mean95.x.age, 1, quantile, 0.025)
SIM.mean.95UCI=apply(SIM.mean95.x.age, 1,  quantile, 0.975)
SIM.mean.iqr=apply(SIM.meaniqr.x.age, 1, mean)
SIM.mean.iqrLCI=apply(SIM.meaniqr.x.age, 1, quantile, 0.025)
SIM.mean.iqrUCI=apply(SIM.meaniqr.x.age, 1, quantile, 0.975)

Gouldia.Piran1.corr.iqr.unit=Gouldia.Piran1.IQR.layers-SIM.mean.iqr
Gouldia.Piran1.corr.iqr.lci.unit=Gouldia.Piran1.IQR.CI.layers[,1]-SIM.mean.iqr
Gouldia.Piran1.corr.iqr.uci.unit=Gouldia.Piran1.IQR.CI.layers[,3]-SIM.mean.iqr
Gouldia.Piran1.corr.95.unit=Gouldia.Piran1.95.layers-SIM.mean.95


####################
#GOULDIA AT PIRAN 2#
####################
Gouldia.Piran2.layers=Gouldia.Piran2.5depths
Gouldia.Piran2.layers[Gouldia.Piran2.5depths<8]=1
Gouldia.Piran2.layers[Gouldia.Piran2.5depths>8 & Gouldia.Piran2.5depths<16]=2
Gouldia.Piran2.layers[Gouldia.Piran2.5depths>16 & Gouldia.Piran2.5depths<35]=3
Gouldia.Piran2.layers[Gouldia.Piran2.5depths>35 & Gouldia.Piran2.5depths<65]=4
Gouldia.Piran2.layers[Gouldia.Piran2.5depths>65]=5
Gouldia.Piran2.IQR.layers=tapply(Gouldia.Piran2.ages,Gouldia.Piran2.layers, IQR)
Gouldia.Piran2.95.layers=tapply(Gouldia.Piran2.ages,Gouldia.Piran2.layers, function(x) {diff(quantile(x,c(0.025, 0.975)))})

Gouldia.Piran2.IQR.CI.layers=array(NA, dim=c(5,3))
Gouldia.Piran2.IQR.CI.layers[1,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==1])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran2.IQR.CI.layers[2,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==2])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran2.IQR.CI.layers[3,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==3])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran2.IQR.CI.layers[4,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==4])$resampleiqr, c(0.025,0.5,0.975))
Gouldia.Piran2.IQR.CI.layers[5,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==5])$resampleiqr, c(0.025,0.5,0.975))

Gouldia.Piran2.95.CI.layers=array(NA, dim=c(5,3))
Gouldia.Piran2.95.CI.layers[1,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==1])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran2.95.CI.layers[2,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==2])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran2.95.CI.layers[3,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==3])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran2.95.CI.layers[4,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==4])$resamplerange, c(0.025,0.5,0.975))
Gouldia.Piran2.95.CI.layers[5,]=quantile(bootci.ext(Gouldia.Piran2.ages[Gouldia.Piran2.layers==5])$resamplerange, c(0.025,0.5,0.975))

SIM.meaniqr.x.age=array(NA, dim=c(5,100));  SIM.mean95.x.age=array(NA, dim=c(5,100))
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:5) {
if (i == 1) {ages=Gouldia.Piran2.ages[Gouldia.Piran2.layers==1]}
if (i == 2) {ages=Gouldia.Piran2.ages[Gouldia.Piran2.layers==2]}
if (i == 3) {ages=Gouldia.Piran2.ages[Gouldia.Piran2.layers==3]}
if (i == 4) {ages=Gouldia.Piran2.ages[Gouldia.Piran2.layers==4]}
if (i == 5) {ages=Gouldia.Piran2.ages[Gouldia.Piran2.layers==5]}
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
temp=rgamma(100,shape=ages[x]/exp(gamma.ln.d),scale=exp(gamma.ln.d))
temp=rtrunc(100, spec="gamma", shape=ages[x]/exp(gamma.ln.d),scale=exp(gamma.ln.d),a=0, b=10000)
x.age[x]=IQR(temp)
x.age.95[x]=diff(quantile(temp, c(0.025,0.975)))
}
meaniqr.x.age[i]=mean(x.age); mean95.x.age[i]=mean(x.age.95);
}
SIM.meaniqr.x.age[,j]=meaniqr.x.age;  SIM.mean95.x.age[,j]=mean95.x.age
}

SIM.mean.95=apply(SIM.mean95.x.age, 1, mean)
SIM.mean.iqr=apply(SIM.meaniqr.x.age, 1, mean)
Gouldia.Piran2.corr.95.unit=Gouldia.Piran2.95.layers-SIM.mean.95

Gouldia.Piran2.corr.iqr.unit=Gouldia.Piran2.IQR.layers-SIM.mean.iqr
Gouldia.Piran2.corr.iqr.lci.unit=Gouldia.Piran2.IQR.CI.layers[,1]-SIM.mean.iqr
Gouldia.Piran2.corr.iqr.uci.unit=Gouldia.Piran2.IQR.CI.layers[,3]-SIM.mean.iqr

##################
#CORBULA AT PIRAN#
##################
Corbula.Piran.layers=Corbula.Piran.5depths
Corbula.Piran.layers[Corbula.Piran.5depths<8]=1
Corbula.Piran.layers[Corbula.Piran.5depths>8 & Corbula.Piran.5depths<16]=2
Corbula.Piran.layers[Corbula.Piran.5depths>16 & Corbula.Piran.5depths<35]=3
Corbula.Piran.layers[Corbula.Piran.5depths>35 & Corbula.Piran.5depths<65]=4
Corbula.Piran.layers[Corbula.Piran.5depths>65]=5
Corbula.Piran.IQR.layers=tapply(Corbula.Piran.ages,Corbula.Piran.layers, IQR)
Corbula.Piran.95.layers=tapply(Corbula.Piran.ages,Corbula.Piran.layers, function(x) {diff(quantile(x,c(0.025, 0.975)))})

Corbula.Piran.IQR.CI.layers=array(NA, dim=c(5,3))
Corbula.Piran.IQR.CI.layers[1,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==1])$resampleiqr, c(0.025,0.5,0.975))
Corbula.Piran.IQR.CI.layers[2,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==2])$resampleiqr, c(0.025,0.5,0.975))
Corbula.Piran.IQR.CI.layers[3,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==3])$resampleiqr, c(0.025,0.5,0.975))
Corbula.Piran.IQR.CI.layers[4,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==4])$resampleiqr, c(0.025,0.5,0.975))
Corbula.Piran.IQR.CI.layers[5,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==4])$resampleiqr, c(0.025,0.5,0.975))

Corbula.Piran.95.CI.layers=array(NA, dim=c(5,3))
Corbula.Piran.95.CI.layers[1,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==1])$resamplerange, c(0.025,0.5,0.975))
Corbula.Piran.95.CI.layers[2,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==2])$resamplerange, c(0.025,0.5,0.975))
Corbula.Piran.95.CI.layers[3,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==3])$resamplerange, c(0.025,0.5,0.975))
Corbula.Piran.95.CI.layers[4,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==4])$resamplerange, c(0.025,0.5,0.975))
Corbula.Piran.95.CI.layers[5,]=quantile(bootci.ext(Corbula.Piran.ages[Corbula.Piran.layers==5])$resamplerange, c(0.025,0.5,0.975))

SIM.meaniqr.x.age=array(NA, dim=c(5,100));  SIM.mean95.x.age=array(NA, dim=c(5,100))
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:5) {
if (i == 1) {ages=Corbula.Piran.ages[Corbula.Piran.layers==1]}
if (i == 2) {ages=Corbula.Piran.ages[Corbula.Piran.layers==2]}
if (i == 3) {ages=Corbula.Piran.ages[Corbula.Piran.layers==3]}
if (i == 4) {ages=Corbula.Piran.ages[Corbula.Piran.layers==4]}
if (i == 5) {ages=Corbula.Piran.ages[Corbula.Piran.layers==5]}
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
temp=rlnorm(100,log(ages[x]),sdlog=sqrt(exp(Corbula.ln.d)))
temp=rtrunc(100, spec="lnorm", mean=log(ages[x]),sdlog=sqrt(exp(Corbula.ln.d)),a=0, b=10000)

x.age[x]=IQR(temp)
x.age.95[x]=diff(quantile(temp, c(0.025,0.975)))
}
meaniqr.x.age[i]=mean(x.age); mean95.x.age[i]=mean(x.age.95);
}
SIM.meaniqr.x.age[,j]=meaniqr.x.age;  SIM.mean95.x.age[,j]=mean95.x.age
}

SIM.mean.95=apply(SIM.mean95.x.age, 1, mean)
SIM.mean.iqr=apply(SIM.meaniqr.x.age, 1, mean)
Corbula.corr.95.unit=Corbula.Piran.95.layers-SIM.mean.95
Corbula.corr.95.lci.unit=Corbula.Piran.95.CI.layers[,1]-SIM.mean.95
Corbula.corr.95.uci.unit=Corbula.Piran.95.CI.layers[,3]-SIM.mean.95

Corbula.corr.iqr.unit=Corbula.Piran.IQR.layers-SIM.mean.iqr
Corbula.corr.iqr.lci.unit=Corbula.Piran.IQR.CI.layers[,1]-SIM.mean.iqr
Corbula.corr.iqr.uci.unit=Corbula.Piran.IQR.CI.layers[,3]-SIM.mean.iqr

##########################################################################
#FIGURE DOWNCORE TIME AVERAGING CORRECTED PER UNIT
##########################################################################
par(mfrow=c(2,2))
plot(Gouldia.Piran1.corr.iqr.unit, c(-4,-12, -25, -50,-110), pch=16, xlim=c(0,3500), 
ylab="Sediment depth (cm)", frame=F, xlab="Time averaging (y)-interquartile range", main="Piran 1", ylim=c(-140,0))
rect(xleft=0, ybottom=-35, xright=3500, ytop=-8, col="gray", border=NA)
abline(h=c(-65,-35,-16,-8), lty=2)
lines(Gouldia.Piran1.corr.iqr.unit, c(-4,-12, -25, -50,-110), lty=3)
segments(x0=Gouldia.Piran1.corr.iqr.lci.unit,x1=Gouldia.Piran1.corr.iqr.uci.unit,y0=c(-4,-12, -25, -50,-110),y1=c(-4,-12, -25, -50,-110), lty=1)
points(Gouldia.Piran1.corr.iqr.unit, c(-4,-12, -25, -50,-110), pch=21, bg="gray", cex=1.4)

plot(Corbula.corr.iqr.unit, c(-4,-12, -25, -50,-110), pch=16, xlim=c(0,3500), 
ylab="Sediment depth (cm)", frame=F, xlab="Time averaging (y)-interquartile range", main="Piran 2", type="n", ylim=c(-140,0))
rect(xleft=0, ybottom=-35, xright=3500, ytop=-8, col="gray", border=NA)
abline(h=c(-65,-35,-16,-8), lty=2)
lines(Corbula.corr.iqr.unit, c(-4,-12, -25, -50,-110)-2, lty=3)
segments(x0=Corbula.corr.iqr.lci.unit,x1=Corbula.corr.iqr.uci.unit,y0=c(-4,-12, -25, -50,-110)-2,y1=c(-4,-12, -25, -50,-110)-2, lty=1)
points(Corbula.corr.iqr.unit, c(-4,-12, -25, -50,-110)-2, pch=21, bg="black", cex=1.4)
lines(Gouldia.Piran2.corr.iqr.unit, c(-4,-12, -25, -50,-110), lty=3)
segments(x0=Gouldia.Piran2.corr.iqr.lci.unit,x1=Gouldia.Piran2.corr.iqr.uci.unit,y0=c(-4,-12, -25, -50,-110),y1=c(-4,-12, -25, -50,-110), lty=1)
points(Gouldia.Piran2.corr.iqr.unit, c(-4,-12, -25, -50,-110), pch=21, bg="gray", cex=1.4)
legend(x="bottomright", legend=c("Gouldia","Corbula"), pch=21, pt.bg=c("gray","black"), bty="n", cex=1.2)

##########################################################################
#STRATIGRAPHIC DISORDER
##########################################################################
par(mfrow=c(2,3))
plot(Gouldia.Piran1.ages,-Gouldia.Piran1.5depths,pch=16)
cor.test(Gouldia.Piran1.ages,-Gouldia.Piran1.5depths, method="s")
cor.test(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<60],-Gouldia.Piran1.5depths[Gouldia.Piran1.5depths<60], method="s")
cor.test(Gouldia.Piran1.ages[Gouldia.Piran1.5depths>60],-Gouldia.Piran1.5depths[Gouldia.Piran1.5depths>60], method="s")
plot(Gouldia.Piran2.ages,-Gouldia.Piran2.5depths,pch=16)
cor.test(Gouldia.Piran2.ages,-Gouldia.Piran2.5depths, method="s")
cor.test(Gouldia.Piran2.ages[Gouldia.Piran2.5depths<60],-Gouldia.Piran2.5depths[Gouldia.Piran2.5depths<60], method="s")
cor.test(Gouldia.Piran2.ages[Gouldia.Piran2.5depths>60],-Gouldia.Piran2.5depths[Gouldia.Piran2.5depths>60], method="s")
plot(Corbula.Piran.ages,-Corbula.Piran.5depths,pch=16)
cor.test(Corbula.Piran.ages,-Corbula.Piran.5depths, method="s")
cor.test(Corbula.Piran.ages[Corbula.Piran.5depths<60],-Corbula.Piran.5depths[Corbula.Piran.5depths<60], method="s")
cor.test(Corbula.Piran.ages[Corbula.Piran.5depths>60],-Corbula.Piran.5depths[Corbula.Piran.5depths>60], method="s")

################################################################################################
#SURFACE DISTRIBUTIONS
################################################################################################
#DENSITY FOR RANDOM 2-PHASE EXPONENTIAL MODEL
################################################################################################
FitExpMix=function(X,niters=10000,a=2,b=0.001,rmean=mean(X)){
N=length(X)
d=function(x){
        r1*exp(-r1*x)*beta+r2*exp(-r2*x)*(1-beta)
        }
beta=0.5;
r1=1;
r2=1;
best=c(beta,r1,r2);
dbest=sum(log(d(X)))
cluster=rbinom(N,1,beta);
c=which(cluster==1)

for(iter in 1:niters){
        beta=rbeta(1,a+length(c),a+N-length(c))
        r1=rgamma(1,b+length(c),rate=rmean*b+sum(X[c]));
        r2=rgamma(1,b+N-length(c),rate=rmean*b+sum(X[-c]))      
        L1=beta*r1*exp(-r1*X)
        L2=(1-beta)*r2*exp(-r2*X)
        cluster=rbinom(N,1,L1/(L1+L2));
        c=which(cluster==1)
        dnew=sum(log(d(X)));
        if(dnew>dbest){
                dbest=dnew;best=c(beta,r1,r2)
        }
}
beta=best[1];r1=best[2];r2=best[3]
if(r1>r2){beta=1-beta;r1temp=r1;r1=r2;r2=r1temp}
out=list();out$beta=beta;out$r1=r1;out$r2=r2
out$alpha=(beta*r1)/(beta*r1+(1-beta)*r2)
out$tau=out$alpha*(r2-r1);out$lambda1=r2-out$tau;out$lambda2=r1
out$lik=dbest
out
}

#Surface
surfaceAFD.Piran1.12=Gouldia.Piran1.ages[Gouldia.Piran1.5depths<12]
surfaceAFD.Piran2.12=Gouldia.Piran2.ages[Gouldia.Piran2.5depths<12]
surfaceAFD.Piran.12=c(Gouldia.Piran1.ages[Gouldia.Piran1.5depths<12],Gouldia.Piran2.ages[Gouldia.Piran2.5depths<12])
surfaceAFD.Corbula.Piran.12=c(Corbula.Piran.ages[Corbula.Piran.5depths<12],Corbula.Piran.ages[Corbula.Piran.5depths<12])
par(mfrow=c(2,2))
hist(surfaceAFD.Piran.12, col="gray",breaks=seq(0,8000,by=100), xlab="Postmortem age (years)", ylim=c(0,25), main="Surface 12 cm-Piran", cex.main=0.9)
hist(surfaceAFD.Piran1.12, col="gray",breaks=seq(0,8000,by=100), xlab="Postmortem age (years)", ylim=c(0,25), main="Surface-12 cm-Piran 1", cex.main=0.9)
hist(surfaceAFD.Piran2.12, col="gray",breaks=seq(0,8000,by=100), xlab="Postmortem age (years)", ylim=c(0,25), main="Surface-12 cm-Piran 2", cex.main=0.9)
hist(surfaceAFD.Corbula.Piran.12, col="gray",breaks=seq(0,8000,by=100), xlab="Postmortem age (years)", ylim=c(0,25), main="Surface-C. gibba-12 cm-Piran 2", cex.main=0.9)
age=surfaceAFD.Corbula.Piran.12

####################
#SIMPLE EXPONENTIAL#
####################
ages=age
onephaserate<-1/mean(age)
lkexp=function(rho) {-length(age)*log(rho) + rho*sum(age)}
like<--(sum(log(onephaserate)-onephaserate*age))	#negative loglikelihood for exponential
AIC=2*1-(-2*like)
AICc=AIC+((2*1*(1+1))/(length(age)-1-1))
onephaseAIC<-AICc
#################################################
#RANDOM 2-PHASE EXPONENTIAL
#################################################
X=age
source("C:/Data/Codes/function_random-time_2-phase_exponential.txt")
out<-FitExpMix(X)
randomtwophaserate1<-out$r1
randomtwophaserate2<-out$r2
randomtwophasealpha<-out$alpha
randomtwophasebeta<-out$beta
lik<--(out$lik)
likrandom=lik
rAIC=2*3-(-2*lik)
randomAIC=rAIC+((2*3*(3+1))/(length(age)-3-1))
tau<-randomtwophasealpha*(randomtwophaserate2-randomtwophaserate1)
lambda1<-randomtwophaserate1*randomtwophasealpha+randomtwophaserate2*(1-randomtwophasealpha)
lambda2<-randomtwophaserate1
lambda1
tau
lambda2

age=surfaceAFD.Corbula.Piran.12
breaks=seq(0,12000,by=100)
hdata<-hist(age, breaks=breaks, cex.main=0.9, col="gray",main="", xlab="Time-since-death (years)", ylab="Number of specimens", ylim=c(0,25), xlim=c(0,10000))
times=hdata$mids-hdata$mids[1]
predictRANDOM<-randomtwophaserate1*exp(-randomtwophaserate1*times)*randomtwophasebeta+randomtwophaserate2*exp(-randomtwophaserate2*times)*(1-randomtwophasebeta)
lines(times, predictRANDOM/sum(predictRANDOM)*length(age), col="gray",lty=1, lwd=2)

#########################################################################################################################
#UNMIXING
#########################################################################################################################
#######################################################################################
#PIRAN 1 - CORE AGE DISTRIBUTION ACCOUNTING FOR ABUNDANCES AND INTERPOLATING UNDATED INCREMENTS 
#######################################################################################
#vector with absolute abundance of Gouldia minima in 4 cm increments, the first object is the uppermost core increment,
#the last is the lowermost increment
N.full.totalsizePiran1=c(160,266,481,576,507,557,363,301,184,146,121,104,52,62,34,29,22,23,20,13,25,17,7,19,11,8,8,10,41,6)
Piran1.list1=split(Gouldia.Piran1.ages,Gouldia.Piran1.5depths)
Full.Piran1bootout.50=array(NA, dim=c(100,300))
Full.Piran1bootout.200=array(NA, dim=c(100,75))
for (j in 1:100) {
Piran1bootages=numeric()
Piran1.list1=split(Gouldia.Piran1.ages,Gouldia.Piran1.5depths)
Piran1.list1$"2"=Piran1.list1$"2"
Piran1.list1$"6"=c(Piran1.list1$"2",Piran1.list1$"10")
Piran1.list1$"10"=Piran1.list1$"10"
Piran1.list1$"14"=Piran1.list1$"14"
Piran1.list1$"18"=rtnorm(100,mean=mean(c(Piran1.list1$"22.5",Piran1.list1$"14")), sd=sd(c(Piran1.list1$"22.5",Piran1.list1$"14")), lower=0, upper=11000)
Piran1.list1$"22.5"=Piran1.list1$"22.5"
Piran1.list1$"27.5"=rtnorm(100,mean=mean(c(Piran1.list1$"22.5",Piran1.list1$"42.5")), sd=sd(c(Piran1.list1$"22.5",Piran1.list1$"42.5")), lower=0, upper=11000)
Piran1.list1$"32.5"=rtnorm(100,mean=mean(c(Piran1.list1$"22.5",Piran1.list1$"42.5")), sd=sd(c(Piran1.list1$"22.5",Piran1.list1$"42.5")), lower=0, upper=11000)
Piran1.list1$"37.5"=rtnorm(100,mean=mean(c(Piran1.list1$"22.5",Piran1.list1$"42.5")), sd=sd(c(Piran1.list1$"22.5",Piran1.list1$"42.5")), lower=0, upper=11000)
Piran1.list1$"42.5"=Piran1.list1$"42.5"
Piran1.list1$"47.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"42.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"42.5")), lower=0, upper=11000)
Piran1.list1$"52.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"42.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"42.5")), lower=0, upper=11000)
Piran1.list1$"57.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"42.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"42.5")), lower=0, upper=11000)
Piran1.list1$"62.5"=Piran1.list1$"62.5"
Piran1.list1$"67.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"72.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"77.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"82.5"=rtnorm(100,mean=mean(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"62.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"87.5"=Piran1.list1$"87.5"
Piran1.list1$"92.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"97.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"102.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"107.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"87.5")), lower=0, upper=11000)
Piran1.list1$"112.5"=Piran1.list1$"112.5"
Piran1.list1$"117.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), lower=0, upper=11000)
Piran1.list1$"122.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), lower=0, upper=11000)
Piran1.list1$"127.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), lower=0, upper=11000)
Piran1.list1$"132.5"=rtnorm(100,mean=mean(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), sd=sd(c(Piran1.list1$"112.5",Piran1.list1$"137.5")), lower=0, upper=11000)
Piran1.list1$"137.5"=Piran1.list1$"137.5"
temp=names(Piran1.list1)
temp=as.numeric(temp)
ord=order(temp)
Piran1.list1=Piran1.list1[ord]
for (i in 1:length(Piran1.list1)) {
temp=unlist(Piran1.list1[i])
bootages=sample(temp,N.full.totalsizePiran1[i],replace=T)
Piran1bootages=append(Piran1bootages,bootages,length(Piran1bootages))
}
for (i in 1:length(Piran1bootages)) {
Piran1bootages[i]=rtrunc(1,spec="gamma", shape=Piran1bootages[i]/exp(gamma.ln.d),scale=exp(gamma.ln.d),a=0, b=10000)
}
out2=hist(Piran1bootages, breaks=c(0,seq(13,15000,by=200))) 
out4=hist(Piran1bootages, breaks=c(0,seq(13,15000,by=50))) 
Full.Piran1bootout.200[j,]=out2$counts
Full.Piran1bootout.50[j,]=out4$counts
}

#######################################################################################
#PIRAN 2 - CORE AGE DISTRIBUTION ACCOUNTING FOR ABUNDANCES AND INTERPOLATING UNDATED INCREMENTS 
#######################################################################################
#vector with absolute abundance of Gouldia minima in 4 cm increments, the first object is the uppermost core increment,
#the last is the lowermost increment
N.full.totalsizePiran2=c(385,489,1158,692,876,648,768,304,316,164,197,132,98,122,50,33,31,31,22,23,9,17,11,12,15,16,23,29,19,10,25,16)
Piran2.list1=split(Gouldia.Piran2.ages,Gouldia.Piran2.5depths)
Full.Piran2bootout.200=array(NA, dim=c(100,75))
Full.Piran2bootout.50=array(NA, dim=c(100,300))
for (j in 1:100) {
Piran2bootages=numeric()
Piran2.list1=split(Gouldia.Piran2.ages,Gouldia.Piran2.5depths)
Piran2.list1$"2"=Piran2.list1$"2"
Piran2.list1$"6"=c(Piran2.list1$"2",Piran2.list1$"10")
Piran2.list1$"10"=Piran2.list1$"10"
Piran2.list1$"14"=Piran2.list1$"14"
Piran2.list1$"18"=Piran2.list1$"18"
Piran2.list1$"22.5"=rtnorm(100,mean=mean(c(Piran2.list1$"18",Piran2.list1$"27.5")), sd=sd(c(Piran2.list1$"18",Piran2.list1$"27.5")), lower=0, upper=11000)
Piran2.list1$"27.5"=Piran2.list1$"27.5"
Piran2.list1$"32.5"=rtnorm(100,mean=mean(c(Piran2.list1$"27.5",Piran2.list1$"47.5")), sd=sd(c(Piran2.list1$"27.5",Piran2.list1$"47.5")), lower=0, upper=11100)
Piran2.list1$"37.5"=rtnorm(100,mean=mean(c(Piran2.list1$"27.5",Piran2.list1$"47.5")), sd=sd(c(Piran2.list1$"27.5",Piran2.list1$"47.5")), lower=0, upper=11100)
Piran2.list1$"42.5"=rtnorm(100,mean=mean(c(Piran2.list1$"27.5",Piran2.list1$"47.5")), sd=sd(c(Piran2.list1$"27.5",Piran2.list1$"47.5")), lower=0, upper=11000)
Piran2.list1$"47.5"=Piran2.list1$"47.5"
Piran2.list1$"52.5"=rtnorm(100,mean=mean(c(Piran2.list1$"67.5",Piran2.list1$"47.5")), sd=sd(c(Piran2.list1$"67.5",Piran2.list1$"47.5")), lower=0, upper=11100)
Piran2.list1$"57.5"=rtnorm(100,mean=mean(c(Piran2.list1$"67.5",Piran2.list1$"47.5")), sd=sd(c(Piran2.list1$"67.5",Piran2.list1$"47.5")), lower=0, upper=11100)
Piran2.list1$"62.5"=rtnorm(100,mean=mean(c(Piran2.list1$"67.5",Piran2.list1$"47.5")), sd=sd(c(Piran2.list1$"67.5",Piran2.list1$"47.5")), lower=0, upper=11100)
Piran2.list1$"67.5"=Piran2.list1$"67.5"
Piran2.list1$"72.5"=rtnorm(100,mean=mean(c(Piran2.list1$"67.5",Piran2.list1$"87.5")), sd=sd(c(Piran2.list1$"67.5",Piran2.list1$"87.5")), lower=0, upper=11100)
Piran2.list1$"77.5"=rtnorm(100,mean=mean(c(Piran2.list1$"67.5",Piran2.list1$"87.5")), sd=sd(c(Piran2.list1$"67.5",Piran2.list1$"87.5")), lower=0, upper=11100)
Piran2.list1$"82.5"=rtnorm(100,mean=mean(c(Piran2.list1$"67.5",Piran2.list1$"87.5")), sd=sd(c(Piran2.list1$"67.5",Piran2.list1$"87.5")), lower=0, upper=11100)
Piran2.list1$"87.5"=Piran2.list1$"87.5"
Piran2.list1$"92.5"=rtnorm(100,mean=mean(c(Piran2.list1$"102.5",Piran2.list1$"87.5")), sd=sd(c(Piran2.list1$"102.5",Piran2.list1$"87.5")), lower=0, upper=11100)
Piran2.list1$"97.5"=rtnorm(100,mean=mean(c(Piran2.list1$"102.5",Piran2.list1$"87.5")), sd=sd(c(Piran2.list1$"102.5",Piran2.list1$"87.5")), lower=0, upper=11100)
Piran2.list1$"102.5"=Piran2.list1$"102.5"
Piran2.list1$"107.5"=rtnorm(100,mean=mean(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), sd=sd(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), lower=0, upper=11100)
Piran2.list1$"112.5"=rtnorm(100,mean=mean(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), sd=sd(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), lower=0, upper=11100)
Piran2.list1$"117.5"=rtnorm(100,mean=mean(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), sd=sd(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), lower=0, upper=11100)
Piran2.list1$"122.5"=rtnorm(100,mean=mean(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), sd=sd(c(Piran2.list1$"102.5",Piran2.list1$"127.5")), lower=0, upper=11100)
Piran2.list1$"127.5"=Piran2.list1$"127.5"
Piran2.list1$"132.5"=rtnorm(100,mean=mean(c(Piran2.list1$"127.5",Piran2.list1$"147.5")), sd=sd(c(Piran2.list1$"127.5",Piran2.list1$"147.5")), lower=0, upper=11100)
Piran2.list1$"137.5"=rtnorm(100,mean=mean(c(Piran2.list1$"127.5",Piran2.list1$"147.5")), sd=sd(c(Piran2.list1$"127.5",Piran2.list1$"147.5")), lower=0, upper=11100)
Piran2.list1$"142.5"=rtnorm(100,mean=mean(c(Piran2.list1$"127.5",Piran2.list1$"147.5")), sd=sd(c(Piran2.list1$"127.5",Piran2.list1$"147.5")), lower=0, upper=11100)
Piran2.list1$"147.5"=Piran2.list1$"147.5"
temp=names(Piran2.list1)
temp=as.numeric(temp)
ord=order(temp)
Piran2.list1=Piran2.list1[ord]
for (i in 1:length(Piran2.list1)) {
temp=unlist(Piran2.list1[i])
bootages=sample(temp,N.full.totalsizePiran2[i],replace=T)
Piran2bootages=append(Piran2bootages,bootages,length(Piran2bootages))
}
for (i in 1:length(Piran2bootages)) {
Piran2bootages[i]=rtrunc(1,spec="gamma", shape=Piran2bootages[i]/exp(gamma.ln.d),scale=exp(gamma.ln.d),a=0, b=10000)
}
out2=hist(Piran2bootages, breaks=c(0,seq(13,15000,by=200))) 
out4=hist(Piran2bootages, breaks=c(0,seq(13,15000,by=50))) 
Full.Piran2bootout.200[j,]=out2$counts
Full.Piran2bootout.50[j,]=out4$counts
}


#######################################################################################
#PIRAN 1 AND 2 - CORE AGE DISTRIBUTION ACCOUNTING FOR ABUNDANCES AND INTERPOLATING UNDATED INCREMENTS 
#######################################################################################
#vector with absolute abundance of Gouldia minima in 4 cm increments - both Piran cores pooled
N.full.totalsizePiran12=c(N.full.totalsizePiran1,0,0)+N.full.totalsizePiran2
#temporary object, will be rewritten
Piran12.list1=Piran2.list1
Full.Piran12bootout.200=array(NA, dim=c(100,75))
Full.Piran12bootout.50=array(NA, dim=c(100,300))
for (j in 1:100) {
Piran12bootages=numeric()
Piran12.list1$"2"=c(Piran1.list1$"2",Piran2.list1$"2")
Piran12.list1$"10"=c(Piran1.list1$"10",Piran2.list1$"10")
Piran12.list1$"6"=c(Piran1.list1$"2",Piran2.list1$"2",Piran1.list1$"10",Piran2.list1$"10")
Piran12.list1$"14"=c(Piran1.list1$"14",Piran2.list1$"14")
Piran12.list1$"18"=Piran2.list1$"18"
Piran12.list1$"22.5"=Piran1.list1$"22.5"
Piran12.list1$"27.5"=Piran2.list1$"27.5"
Piran12.list1$"42.5"=Piran1.list1$"42.5"
Piran12.list1$"47.5"=Piran2.list1$"47.5"
Piran12.list1$"62.5"=Piran1.list1$"62.5"
Piran12.list1$"67.5"=Piran2.list1$"67.5"
Piran12.list1$"87.5"=c(Piran1.list1$"87.5",Piran2.list1$"87.5")
Piran12.list1$"102.5"=Piran2.list1$"102.5"
Piran12.list1$"112.5"=Piran1.list1$"112.5"
Piran12.list1$"127.5"=Piran2.list1$"127.5"
Piran12.list1$"137.5"=Piran1.list1$"137.5"
Piran12.list1$"147.5"=Piran2.list1$"147.5"
Piran12.list1$"32.5"=rtnorm(100,mean=mean(c(Piran12.list1$"27.5",Piran12.list1$"42.5")), sd=sd(c(Piran12.list1$"27.5",Piran12.list1$"42.5")), lower=0, upper=11000)
Piran12.list1$"37.5"=rtnorm(100,mean=mean(c(Piran12.list1$"27.5",Piran12.list1$"42.5")), sd=sd(c(Piran12.list1$"27.5",Piran12.list1$"42.5")), lower=0, upper=11000)
Piran12.list1$"52.5"=rtnorm(100,mean=mean(c(Piran12.list1$"62.5",Piran12.list1$"47.5")), sd=sd(c(Piran12.list1$"62.5",Piran12.list1$"47.5")), lower=0, upper=11000)
Piran12.list1$"57.5"=rtnorm(100,mean=mean(c(Piran12.list1$"62.5",Piran12.list1$"47.5")), sd=sd(c(Piran12.list1$"62.5",Piran12.list1$"47.5")), lower=0, upper=11000)
Piran12.list1$"72.5"=rtnorm(100,mean=mean(c(Piran12.list1$"67.5",Piran12.list1$"87.5")), sd=sd(c(Piran12.list1$"67.5",Piran12.list1$"87.5")), lower=0, upper=11000)
Piran12.list1$"77.5"=rtnorm(100,mean=mean(c(Piran12.list1$"67.5",Piran12.list1$"87.5")), sd=sd(c(Piran12.list1$"67.5",Piran12.list1$"87.5")), lower=0, upper=11000)
Piran12.list1$"82.5"=rtnorm(100,mean=mean(c(Piran12.list1$"67.5",Piran12.list1$"87.5")), sd=sd(c(Piran12.list1$"67.5",Piran12.list1$"87.5")), lower=0, upper=11000)
Piran12.list1$"92.5"=rtnorm(100,mean=mean(c(Piran12.list1$"102.5",Piran12.list1$"87.5")), sd=sd(c(Piran12.list1$"102.5",Piran12.list1$"87.5")), lower=0, upper=11000)
Piran12.list1$"97.5"=rtnorm(100,mean=mean(c(Piran12.list1$"102.5",Piran12.list1$"87.5")), sd=sd(c(Piran12.list1$"102.5",Piran12.list1$"87.5")), lower=0, upper=11000)
Piran12.list1$"107.5"=rtnorm(100,mean=mean(c(Piran12.list1$"102.5",Piran12.list1$"112.5")), sd=sd(c(Piran12.list1$"102.5",Piran12.list1$"112.5")), lower=0, upper=11000)
Piran12.list1$"117.5"=rtnorm(100,mean=mean(c(Piran12.list1$"112.5",Piran12.list1$"127.5")), sd=sd(c(Piran12.list1$"112.5",Piran12.list1$"127.5")), lower=0, upper=11000)
Piran12.list1$"122.5"=rtnorm(100,mean=mean(c(Piran12.list1$"112.5",Piran12.list1$"127.5")), sd=sd(c(Piran12.list1$"112.5",Piran12.list1$"127.5")), lower=0, upper=11000)
Piran12.list1$"132.5"=rtnorm(100,mean=mean(c(Piran12.list1$"127.5",Piran12.list1$"137.5")), sd=sd(c(Piran12.list1$"127.5",Piran12.list1$"137.5")), lower=0, upper=11000)
Piran12.list1$"142.5"=rtnorm(100,mean=mean(c(Piran12.list1$"137.5",Piran12.list1$"147.5")), sd=sd(c(Piran12.list1$"137.5",Piran2.list1$"147.5")), lower=0, upper=11000)
temp=names(Piran12.list1)
temp=as.numeric(temp)
ord=order(temp)
Piran12.list1=Piran12.list1[ord]
for (i in 1:length(Piran12.list1)) {
temp=unlist(Piran12.list1[i])
bootages=sample(temp,N.full.totalsizePiran12[i],replace=T)
Piran12bootages=append(Piran12bootages,bootages,length(Piran12bootages))
}
for (i in 1:length(Piran12bootages)) {
Piran12bootages[i]=rtrunc(1,spec="gamma", shape=Piran12bootages[i]/exp(gamma.ln.d),scale=exp(gamma.ln.d),a=0, b=10000)
}
out2=hist(Piran12bootages, breaks=c(0,seq(13,15000,by=200))) 
out4=hist(Piran12bootages, breaks=c(0,seq(13,15000,by=50))) 
Full.Piran12bootout.200[j,]=out2$counts
Full.Piran12bootout.50[j,]=out4$counts
}


#######################################################################################################################
#CORBULA PIRAN 2 - CORE AGE DISTRIBUTION ACCOUNTING FOR ABUNDANCES AND INTERPOLATING UNDATED INCREMENTS 
#######################################################################################################################
#vector with absolute abundance of Corbula gibba in 4 cm increments at Piran 2, the first object is the uppermost core increment,
#the last is the lowermost increment
full.totalsizePiran.Corbula=c(138,168,304,180,236,192,376,168,170,110,107,71,89,96,65,55,61,52,41,38,33,52,60,49,51,41,79,63,70,43,48)
Full.Piranbootout.200=array(NA, dim=c(100,70))
Full.Piranbootout.50=array(NA, dim=c(100,280))

for (j in 1:100) {
Piranbootages=numeric()
Piran.list1=split(Corbula.Piran.ages,Corbula.Piran.5depths)
Piran.list1$"2"=Piran.list1$"2"
Piran.list1$"6"=c(Piran.list1$"2",Piran.list1$"10")
Piran.list1$"10"=Piran.list1$"10"
Piran.list1$"14"=rtnorm(100,mean=mean(c(Piran.list1$"10",Piran.list1$"27.5")), sd=sd(c(Piran.list1$"10",Piran.list1$"27.5")), lower=0, upper=11000)
Piran.list1$"18"=rtnorm(100,mean=mean(c(Piran.list1$"10",Piran.list1$"27.5")), sd=sd(c(Piran.list1$"10",Piran.list1$"27.5")), lower=0, upper=11000)
Piran.list1$"22.5"=rtnorm(100,mean=mean(c(Piran.list1$"10",Piran.list1$"27.5")), sd=sd(c(Piran.list1$"10",Piran.list1$"27.5")), lower=0, upper=11000)
Piran.list1$"27.5"=Piran.list1$"27.5"
Piran.list1$"32.5"=rtnorm(100,mean=mean(c(Piran.list1$"27.5",Piran.list1$"47.5")), sd=sd(c(Piran.list1$"27.5",Piran.list1$"47.5")), lower=0, upper=11000)
Piran.list1$"37.5"=rtnorm(100,mean=mean(c(Piran.list1$"27.5",Piran.list1$"47.5")), sd=sd(c(Piran.list1$"27.5",Piran.list1$"47.5")), lower=0, upper=11000)
Piran.list1$"42.5"=rtnorm(100,mean=mean(c(Piran.list1$"27.5",Piran.list1$"47.5")), sd=sd(c(Piran.list1$"27.5",Piran.list1$"47.5")), lower=0, upper=11000)
Piran.list1$"47.5"=Piran.list1$"47.5"
Piran.list1$"52.5"=rtnorm(100,mean=mean(c(Piran.list1$"47.5",Piran.list1$"67.5")), sd=sd(c(Piran.list1$"47.5",Piran.list1$"67.5")), lower=0, upper=11000)
Piran.list1$"57.5"=rtnorm(100,mean=mean(c(Piran.list1$"47.5",Piran.list1$"67.5")), sd=sd(c(Piran.list1$"47.5",Piran.list1$"67.5")), lower=0, upper=11000)
Piran.list1$"62.5"=rtnorm(100,mean=mean(c(Piran.list1$"47.5",Piran.list1$"67.5")), sd=sd(c(Piran.list1$"47.5",Piran.list1$"67.5")), lower=0, upper=11000)
Piran.list1$"67.5"=Piran.list1$"67.5"
Piran.list1$"72.5"=rtnorm(100,mean=mean(c(Piran.list1$"67.5",Piran.list1$"87.5")), sd=sd(c(Piran.list1$"67.5",Piran.list1$"87.5")), lower=0, upper=11000)
Piran.list1$"77.5"=rtnorm(100,mean=mean(c(Piran.list1$"67.5",Piran.list1$"87.5")), sd=sd(c(Piran.list1$"67.5",Piran.list1$"87.5")), lower=0, upper=11000)
Piran.list1$"82.5"=rtnorm(100,mean=mean(c(Piran.list1$"67.5",Piran.list1$"87.5")), sd=sd(c(Piran.list1$"67.5",Piran.list1$"87.5")), lower=0, upper=11000)
Piran.list1$"87.5"=Piran.list1$"87.5"
Piran.list1$"92.5"=rtnorm(100,mean=mean(c(Piran.list1$"87.5",Piran.list1$"102.5")), sd=sd(c(Piran.list1$"87.5",Piran.list1$"102.5")), lower=0, upper=11000)
Piran.list1$"97.5"=rtnorm(100,mean=mean(c(Piran.list1$"87.5",Piran.list1$"102.5")), sd=sd(c(Piran.list1$"87.5",Piran.list1$"102.5")), lower=0, upper=11000)
Piran.list1$"102.5"=Piran.list1$"102.5"
Piran.list1$"107.5"=rtnorm(100,mean=mean(c(Piran.list1$"102.5",Piran.list1$"127.5")), sd=sd(c(Piran.list1$"102.5",Piran.list1$"127.5")), lower=0, upper=11000)
Piran.list1$"112.5"=rtnorm(100,mean=mean(c(Piran.list1$"102.5",Piran.list1$"127.5")), sd=sd(c(Piran.list1$"102.5",Piran.list1$"127.5")), lower=0, upper=11000)
Piran.list1$"117.5"=rtnorm(100,mean=mean(c(Piran.list1$"102.5",Piran.list1$"127.5")), sd=sd(c(Piran.list1$"102.5",Piran.list1$"127.5")), lower=0, upper=11000)
Piran.list1$"122.5"=rtnorm(100,mean=mean(c(Piran.list1$"102.5",Piran.list1$"127.5")), sd=sd(c(Piran.list1$"102.5",Piran.list1$"127.5")), lower=0, upper=11000)
Piran.list1$"127.5"=Piran.list1$"127.5"
Piran.list1$"132.5"=rtnorm(100,mean=mean(c(Piran.list1$"127.5",Piran.list1$"147.5")), sd=sd(c(Piran.list1$"127.5",Piran.list1$"147.5")), lower=0, upper=11000)
Piran.list1$"137.5"=rtnorm(100,mean=mean(c(Piran.list1$"127.5",Piran.list1$"147.5")), sd=sd(c(Piran.list1$"127.5",Piran.list1$"147.5")), lower=0, upper=11000)
Piran.list1$"142.5"=rtnorm(100,mean=mean(c(Piran.list1$"127.5",Piran.list1$"147.5")), sd=sd(c(Piran.list1$"127.5",Piran.list1$"147.5")), lower=0, upper=11000)
Piran.list1$"147.5"=Piran.list1$"147.5"
temp=names(Piran.list1)
temp=as.numeric(temp)
ord=order(temp)
Piran.list1=Piran.list1[ord]
for (i in 1:length(Piran.list1)) {
temp=unlist(Piran.list1[i])
bootages=sample(temp,full.totalsizePiran.Corbula[i],replace=T)
Piranbootages=append(Piranbootages,bootages,length(Piranbootages))
}
for (i in 1:length(Piranbootages)) {
Piranbootages[i]=rtrunc(1, spec="lnorm", mean=log(Piranbootages[i]),sdlog=sqrt(exp(Corbula.ln.d)),a=0, b=11000)
}
out2=hist(Piranbootages[Piranbootages<13800], breaks=c(0,seq(13,14000,by=200))) 
out4=hist(Piranbootages[Piranbootages<13750], breaks=c(0,seq(13,14000,by=50))) 
Full.Piranbootout.200[j,]=out2$counts
Full.Piranbootout.50[j,]=out4$counts
}


##############################################################################################################################
#UNMIXED RECORDS - MEANS AND QUANTILES FOR 200-YEAR COHORTS
##############################################################################################################################
Piran1.means=apply(Full.Piran1bootout.200,MARGIN=2,mean)
Piran1.025=apply(Full.Piran1bootout.200,MARGIN=2,quantile, 0.025)
Piran1.975=apply(Full.Piran1bootout.200,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=200))
Piran2.means=apply(Full.Piran2bootout.200,MARGIN=2,mean)
Piran2.025=apply(Full.Piran2bootout.200,MARGIN=2,quantile, 0.025)
Piran2.975=apply(Full.Piran2bootout.200,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=200))
Piran.means=apply(Full.Piranbootout.200,MARGIN=2,mean)
Piran.025=apply(Full.Piranbootout.200,MARGIN=2,quantile, 0.025)
Piran.975=apply(Full.Piranbootout.200,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=200))
Piran12.means=apply(Full.Piran12bootout.200,MARGIN=2,mean)
Piran12.025=apply(Full.Piran12bootout.200,MARGIN=2,quantile, 0.025)
Piran12.975=apply(Full.Piran12bootout.200,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=200))
#######################################################################################################
#UNMIXED RECORDS - 200 YEAR COHORTS, INCLUDES LOSS IN THE TAPHONOMIC ACTIVE ZONE
#######################################################################################################
temp.lambda1=lambda1
temp.tau=tau
#lambda 2 is set to zero because no or negligible burial is assumed below the core
temp.lambda2=0
alpha=temp.tau/(temp.tau+temp.lambda1-temp.lambda2)

VECTOR=BREAKS[1:100]
survival.twophase=(alpha*exp(-temp.lambda2*VECTOR)+(1-alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
survival.onephase=exp(-onephaserate*VECTOR)
Piran1.exp.input=Piran1.means[2:51]/survival.onephase
Piran1.exp2.input=Piran1.means[2:51]/survival.twophase
Piran1.exp2.input.lci=Piran1.025[2:51]/survival.twophase
Piran1.exp2.input.uci=Piran1.975[2:51]/survival.twophase

Piran2.exp.input=Piran2.means[2:51]/survival.onephase
Piran2.exp2.input=Piran2.means[2:51]/survival.twophase
Piran2.exp2.input.lci=Piran2.025[2:51]/survival.twophase
Piran2.exp2.input.uci=Piran2.975[2:51]/survival.twophase

Piran2.exp2.input=Piran2.means[2:51]/survival.twophase
Corbula.Piran2.exp.input=Piran.means[2:51]/survival.onephase
Corbula.Piran2.exp2.input=Piran.means[2:51]/survival.twophase
Corbula.Piran2.exp2.input.lci=Piran.025[2:51]/survival.twophase
Corbula.Piran2.exp2.input.uci=Piran.975[2:51]/survival.twophase



#######################################################################################################
#FIGURE 11 IN THE MANUSCRIPT
#######################################################################################################

par(mfrow=c(2,3))
par(mar=c(4,4,2,1))
plot(Piran1.means[2:50], 2013-BREAKS[2:50], xlab="Abundance/200 y",  frame=F,cex.axis=1.2,cex.lab=1.2,
type="n", ylab="Year (BC/AD)",pch=16, log="", xlim=c(0,500), ylim=c(-8000,2000), main="", cex.main=0.9)
lines(Piran1.means[2:50], 2013-BREAKS[2:50], lwd=3, col="gray")
abline(h=c(1800,-2500,-4500),lty=2)
par(new=T)
plot(Piran1.exp2.input[2:50], 2013-BREAKS[2:50], xlab="",  frame=F, xaxt="n",cex.axis=1.2,cex.lab=1.2,
type="l", lwd=3, ylab="",pch=16, log="",  ylim=c(-8000,2000), xlim=c(0,1500),main="", cex.main=0.9)
lines(Piran1.exp2.input.lci[2:50], 2013-BREAKS[2:50], lty=2)
lines(Piran1.exp2.input.uci[2:50], 2013-BREAKS[2:50], lty=2)
axis(3, at=c(0,500,1000,1500), labels=c(0,500,1000,1500), cex.axis=1.2, cex.lab=1.2)

plot(Piran2.means[2:50], 2013-BREAKS[2:50], xlab="Abundance/200 y",  frame=F,cex.axis=1.2,cex.lab=1.2,
type="n", ylab="",pch=16, log="", xlim=c(0,500), ylim=c(-8000,2000), main="", cex.main=0.9)
lines(Piran2.means[2:50], 2013-BREAKS[2:50], lwd=3, col="gray")
abline(h=c(1800,-2500,-4500),lty=2)
par(new=T)
plot(Piran2.exp2.input[2:50], 2013-BREAKS[2:50], xlab="",  frame=F, xaxt="n",cex.axis=1.2,cex.lab=1.2,
type="l", lwd=3, ylab="",pch=16, log="",  ylim=c(-8000,2000), xlim=c(0,1500),main="", cex.main=0.9)
lines(Piran2.exp2.input.lci[2:50], 2013-BREAKS[2:50], lty=2)
lines(Piran2.exp2.input.uci[2:50], 2013-BREAKS[2:50], lty=2)
axis(3, at=c(0,500,1000,1500), labels=c(0,500,1000,1500), cex.axis=1.2, cex.lab=1.2)

plot(Piran.means[2:50], 2013-BREAKS[2:50], xlab="Abundance/200 y",  frame=F,cex.axis=1.2,cex.lab=1.2,
type="n", ylab="",pch=16, log="", xlim=c(0,500), ylim=c(-8000,2000), main="", cex.main=0.9)
lines(Piran.means[2:50], 2013-BREAKS[2:50], lwd=3, col="gray")
legend(x="bottomright", col=c("gray","black"), lwd=3, legend=c("with loss in TAZ","without loss in TAZ"), bty="n")
abline(h=c(1800,-2500,-4500),lty=2)
par(new=T)
plot(Corbula.Piran2.exp2.input[2:50], 2013-BREAKS[2:50], xlab="",  frame=F, xaxt="n",cex.axis=1.2,cex.lab=1.2,
type="l", lwd=3, ylab="",pch=16, log="",  ylim=c(-8000,2000), xlim=c(0,1500),main="", cex.main=0.9)
lines(Corbula.Piran2.exp2.input.lci[2:50], 2013-BREAKS[2:50], lty=2)
lines(Corbula.Piran2.exp2.input.uci[2:50], 2013-BREAKS[2:50], lty=2)
axis(3, at=c(0,500,1000,1500), labels=c(0,500,1000,1500),cex.axis=1.2,cex.lab=1.2)

#######################################################################################################
#UNMIXED 50 YEAR COHORTS, NO SHELL LOSS IN TAZ
#######################################################################################################
Piran1.means=apply(Full.Piran1bootout.50,MARGIN=2,mean)
Piran1.025=apply(Full.Piran1bootout.50,MARGIN=2,quantile, 0.025)
Piran1.975=apply(Full.Piran1bootout.50,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=50))

Piran2.means=apply(Full.Piran2bootout.50,MARGIN=2,mean)
Piran2.025=apply(Full.Piran2bootout.50,MARGIN=2,quantile, 0.025)
Piran2.975=apply(Full.Piran2bootout.50,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=50))

Piran12.means=apply(Full.Piran12bootout.50,MARGIN=2,mean)
Piran12.025=apply(Full.Piran12bootout.50,MARGIN=2,quantile, 0.025)
Piran12.975=apply(Full.Piran12bootout.50,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=50))

Piran.means=apply(Full.Piranbootout.50,MARGIN=2,mean)
Piran.025=apply(Full.Piranbootout.50,MARGIN=2,quantile, 0.025)
Piran.975=apply(Full.Piranbootout.50,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,13600,by=50))

#######################################################################################################
#UNMIXED 50 YEAR COHORTS, ADDING SHELL LOSS IN TAZ
#######################################################################################################

VECTOR=BREAKS[1:250]
survival.twophase=(alpha*exp(-temp.lambda2*VECTOR)+(1-alpha)*exp(-(temp.tau+temp.lambda1)*VECTOR))
survival.onephase=exp(-onephaserate*VECTOR)
Piran1.exp.input=Piran1.means[2:251]/survival.onephase
Piran1.exp2.input=Piran1.means[2:251]/survival.twophase
Piran1.exp2.input.lci=Piran1.025[2:251]/survival.twophase
Piran1.exp2.input.uci=Piran1.975[2:251]/survival.twophase

Piran2.exp.input=Piran2.means[2:251]/survival.onephase
Piran2.exp2.input=Piran2.means[2:251]/survival.twophase
Piran2.exp2.input.lci=Piran2.025[2:251]/survival.twophase
Piran2.exp2.input.uci=Piran2.975[2:251]/survival.twophase

Piran12.exp.input=Piran12.means[2:251]/survival.onephase
Piran12.exp2.input=Piran12.means[2:251]/survival.twophase
Piran12.exp2.input.lci=Piran12.025[2:251]/survival.twophase
Piran12.exp2.input.uci=Piran12.975[2:251]/survival.twophase

Corbula.Piran2.exp.input=Piran.means[2:251]/survival.onephase
Corbula.Piran2.exp2.input=Piran.means[2:251]/survival.twophase
Corbula.Piran2.exp2.input.lci=Piran.025[2:251]/survival.twophase
Corbula.Piran2.exp2.input.uci=Piran.975[2:251]/survival.twophase

#######################################################################################################
#BOTTOM ROW FOR TAZ AND 50 YR COHORTS POOLED SUMMARY FOR BOTH CORES FOR GOULDIA
#######################################################################################################
age=surfaceAFD.Corbula.Piran.12
breaks=seq(0,12000,by=100)
hdata<-hist(age, breaks=breaks, cex.main=1, col="gray", xlab="Time-since-death (years)", 
ylab="Number of specimens", ylim=c(0,25), main="Loss in TAZ", xlim=c(0,8000), cex.lab=1.2, cex.axis=1.2)
times=hdata$mids-hdata$mids[1]
predictRANDOM<-randomtwophaserate1*exp(-randomtwophaserate1*times)*randomtwophasebeta+randomtwophaserate2*exp(-randomtwophaserate2*times)*(1-randomtwophasebeta)
lines(times, predictRANDOM/sum(predictRANDOM)*length(age), col="gray",lty=1, lwd=3)

plot(Piran1.means[2:40], 2013-BREAKS[2:40], xlab="Predicted  abundance/50y", cex.lab=1.2, cex.axis=1.2, 
type="n", ylab="Year (AD)",pch=16, log="", xlim=c(0,300), ylim=c(500,2000), main="", cex.main=1, frame=F)
lines(Piran12.means[2:40], 2013-BREAKS[2:40], lwd=3, col="gray")
lines(Piran12.025[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="gray")
lines(Piran12.975[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="gray")
abline(h=c(1800,1900), lty=3)
text(250,1750, labels="1800 AD")
par(new=T)
plot(Corbula.Piran2.exp2.input[2:50], 2013-BREAKS[2:50], xlab="",  frame=F, xaxt="n",cex.axis=1.2,cex.lab=1.2,
type="n", lwd=3, ylab="",pch=16, log="",  xlim=c(0,600), ylim=c(500,2000), main="", cex.main=0.9)
axis(3, at=c(0,250,500,1000,1500), labels=c(0,250,500,1000,1500),cex.axis=1.2,cex.lab=1.2)
lines(Piran12.exp2.input[2:40], 2013-BREAKS[2:40], lwd=3, col="black")
lines(Piran12.exp2.input.lci[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="black")
lines(Piran12.exp2.input.uci[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="black")

plot(Piran.means[2:40], 2013-BREAKS[2:40], xlab="Predicted abundance/50y", frame=F,
type="n", ylab="",pch=16, log="",  ylim=c(500,2000), xlim=c(0,150), main="", xaxt="n",yaxt="n",
cex.main=1, cex.axis=1.2, cex.lab=1.2)
axis(1, at=c(0,50,100,150,200,250,300), labels=c(0,50,100,150,200,250,300), cex.axis=1.2)
axis(2, at=c(0,500,1000,1500,2000), labels=c(0,500,1000,1500,2000), cex.axis=1.2)
lines(Corbula.Piran2.exp2.input[2:40], 2013-BREAKS[2:40], lwd=3, col="black")
lines(Corbula.Piran2.exp2.input.lci[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="black")
lines(Corbula.Piran2.exp2.input.uci[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="black")
abline(h=c(1800,1900), lty=3)
text(100,1750, labels="1800 AD")
par(new=T)
plot(Corbula.Piran2.exp2.input[2:50], 2013-BREAKS[2:50], xlab="",  frame=F, xaxt="n",cex.axis=1.2,cex.lab=1.2,
type="n", lwd=3, ylab="",pch=16, log="",  xlim=c(0,50), ylim=c(500,2000), main="", cex.main=0.9)
lines(Piran.means[2:40], 2013-BREAKS[2:40], lwd=3, col="gray")
lines(Piran.025[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="gray")
lines(Piran.975[2:40], 2013-BREAKS[2:40], lwd=1, lty=2, col="gray")
axis(3, at=c(0,25,50,75,100), labels=c(0,25,50,75,100),cex.axis=1.2,cex.lab=1.2)