require(e1071)
require(psych)
require(nlme)
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
}

#DENSITY FOR WEIBULL MODEL
weibull=function(x) {
        age=age
        x[1]=exp(x[1])
        x[2]=exp(x[2])          
        C<-(1/x[1])*gamma(1+1/x[2])
        lik<-(1/C)*exp(-(x[1]*age)^x[2])                
        lik=log(lik)
return(-sum(lik))
}

#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
}
###########################################################################################################
data.raw <- read.delim("Supplementary table 2.txt", header=TRUE)

###########################################################################################################
finalAges=data.raw[,"Age"]
remove=which(is.na(finalAges))
data.raw=data.raw[-remove,]
finalAges=data.raw[,"Age"]

Corbula.Panzano.ages=finalAges[data.raw$Station=="Pan M28"]
Corbula.Panzano.5depths=as.numeric(as.vector(data.raw[data.raw$Station=="Pan M28",]$"X5.cm.group..cm."))
#this object is used for some figures to pool shells from adjacent increments with smaller sample sizes
Corbula.Panzano.rev5depths=ifelse(Corbula.Panzano.5depths==70,65,Corbula.Panzano.5depths)
Corbula.Panzano.rev5depths=ifelse(Corbula.Panzano.rev5depths==90,95,Corbula.Panzano.rev5depths)
Corbula.Panzano.rev5depths=ifelse(Corbula.Panzano.5depths<25,Corbula.Panzano.rev5depths-2,Corbula.Panzano.rev5depths)
Corbula.Panzano.rev5depths=ifelse(Corbula.Panzano.5depths>20,Corbula.Panzano.rev5depths-2.5,Corbula.Panzano.rev5depths)
Corbula.Panzano.station=data.raw[data.raw$Station=="Pan M28","Station"]
Corbula.Panzano.maxdepths=as.numeric(as.vector(data.raw[data.raw$Station=="Pan M28",]$"max..depth..cm."))

Corbula.Po3.ages=finalAges[data.raw$Station=="Po3 M13"]
Corbula.Po3.5depths=as.numeric(as.vector(data.raw[data.raw$Station=="Po3 M13",]$"X5.cm.group..cm."))
Corbula.Po3.rev5depths=ifelse(Corbula.Po3.5depths==135,140,Corbula.Po3.5depths)
Corbula.Po3.rev5depths=ifelse(Corbula.Po3.rev5depths==110,115,Corbula.Po3.rev5depths)
Corbula.Po3.5depths=ifelse(Corbula.Po3.5depths<25,Corbula.Po3.5depths-2,Corbula.Po3.5depths)
Corbula.Po3.5depths=ifelse(Corbula.Po3.5depths>20,Corbula.Po3.5depths-2.5,Corbula.Po3.5depths)
#this object is used for some figures to pool shells from adjacent increments with smaller sample sizes
Corbula.Po3.station=data.raw[data.raw$Station=="Po3 M13","Station"]

Corbula.Po4.ages=finalAges[data.raw$Station=="Po4 M21"]
Corbula.Po4.5depths=as.numeric(as.vector(data.raw[data.raw$Station=="Po4 M21",]$"X5.cm.group..cm."))
Corbula.Po4.rev5depths=ifelse(Corbula.Po4.5depths==130,135,Corbula.Po4.5depths)
Corbula.Po4.5depths=ifelse(Corbula.Po4.5depths<25,Corbula.Po4.5depths-2,Corbula.Po4.5depths)
Corbula.Po4.5depths=ifelse(Corbula.Po4.5depths>20,Corbula.Po4.5depths-2.5,Corbula.Po4.5depths)
#this object is used for some figures to pool shells from adjacent increments with smaller sample sizes
Corbula.Po4.station=data.raw[data.raw$Station=="Po4 M21","Station"]

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"]

###########################################################################################
#variance parameter for lognormal calibration model (TDK1)
ln.d=-2.896
#variance parameter for gamma calibration model (SPK0)
gamma.ln.d=4.954087
###########################################################################################


################################################################################
#DOWNCORE CHANGES IN DISTRIBUTIONS AT PANZANO
################################################################################
#par(new=T)
age=Corbula.Panzano.ages
Coredepth=Corbula.Panzano.rev5depths
Corestation=as.vector(Corbula.Panzano.station)


LIST=length(split(age, -Coredepth))
par(mfcol=c(LIST,2))
par(mar=c(1,5,1,5))
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,850,by=25), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="gray", ylim=c(0,25),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]))
tempage=unlist(split(age, Coredepth)[LIST])
hist(tempage[tempage<1000], xlab="Age (years)", ylab="", breaks=seq(0,850,by=25), cex.main=0.8, cex.axis=0.8,
col="gray", ylim=c(0,25),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)



################################################################################
#DEPTH CHANGES PO3
################################################################################
age=Corbula.Po3.ages
Coredepth=Corbula.Po3.5depths
Corestation=as.vector(Corbula.Po3.station)

LIST=length(split(age, -Coredepth))
par(mfcol=c(LIST,3))
par(mar=c(1,3,1,3))
#par(mar=c(2,15,1,15))
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(paste(as.numeric(depths[i])-10,"-",depths[i], sep=""))
interval2[i]=as.character(as.numeric(depths2[i]))
hist(unlist(split(age, Coredepth)[i]), xlab="Age (years)", ylab="", breaks=seq(0,250,by=10), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="gray", ylim=c(0,20),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)
}

tempage=unlist(split(age, Coredepth)[LIST])
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(as.numeric(depths2[LIST]))
hist(unlist(split(age, Coredepth)[LIST]), xlab="Age (years)", ylab="", breaks=seq(0,250,by=10), cex.main=0.8, cex.axis=0.8,
col="gray", ylim=c(0,20),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)

#######################################################################################################
#DEPTH CHANGES PO4
#######################################################################################################
age=Corbula.Po4.ages
Coredepth=Corbula.Po4.5depths
Corestation=as.vector(Corbula.Po4.station)

LIST=length(split(age, -Coredepth))
#par(mfrow=c(10,1))
#par(mar=c(2,15,1,15))
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(paste(as.numeric(depths[i])-10,"-",depths[i], sep=""))
interval3[i]=as.character(as.numeric(depths3[i]))
hist(unlist(split(age, Coredepth)[i]), xlab="Age (years)", ylab="", breaks=seq(0,250,by=10), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="gray", ylim=c(0,20),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)
}

tempage=unlist(split(age, Coredepth)[LIST])
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,260,by=10), cex.main=0.8, cex.axis=0.8,
col="gray", ylim=c(0,20),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)


#######################################################################################################
#DEPTH CHANGES PIRAN
#######################################################################################################
age=Corbula.Piran.ages
Coredepth=Corbula.Piran.5depths
Corestation=as.vector(Corbula.Piran.station)

LIST=length(split(age, -Coredepth))
par(mfrow=c(9,1))
par(mar=c(2,15,1,15))
samplesizes4=numeric()
depths4=numeric()
IQRrange4=numeric()
fullrange4=numeric()
medianage4=numeric()
meanage4=numeric()
range954=numeric()
minage4=numeric()
maxage4=numeric()
skewness4=numeric()
station4=numeric()
interval4=numeric()
totalrange4=numeric()
skew4=numeric()
correction=numeric()
obs.UCI95range4=numeric(); obs.LCI95range4=numeric()
obs.UCIrange4=numeric(); obs.LCIrange4=numeric()
obs.UCImedian4=numeric(); obs.LCImedian4=numeric()
obs.UCImean4=numeric(); obs.LCImean4=numeric()
obs.UCIIQR4=numeric(); obs.LCIIQR4=numeric()
obs.UCIskewness4=numeric(); obs.LCIskewness4=numeric()

for (i in 1:(LIST-1)) {
samplesizes4[i]=length(unlist(split(age, Coredepth)[i]))
depths4[i]=names(split(age, Coredepth)[i])
#interval4[i]=as.character(paste(as.numeric(depths[i])-10,"-",depths[i], sep=""))
interval4[i]=as.character(as.numeric(depths4[i]))
hist(unlist(split(age, Coredepth)[i]), xlab="Age (years)", ylab="", breaks=seq(0,16000,by=500), cex.main=0.8, cex.axis=0.8, xaxt="n",
col="gray", ylim=c(0,20),main=paste(interval4[i],"cm","(n =",samplesizes4[i],")"))
IQRrange4[i]=IQR(unlist(split(age, Coredepth)[i]))
fullrange4[i]=diff(range(unlist(split(age, Coredepth)[i])))
medianage4[i]=median(unlist(split(age, Coredepth)[i]))
meanage4[i]=mean(unlist(split(age, Coredepth)[i]))
minage4[i]=min(unlist(split(age, Coredepth)[i]))
maxage4[i]=max(unlist(split(age, Coredepth)[i]))
range954[i]=diff(quantile(unlist(split(age, Coredepth)[i]),c(0.025,0.975)))
totalrange4[i]=diff(range(unlist(split(age, Coredepth)[i])))
skew4[i]=skewness(unlist(split(age, Coredepth)[i]))
station4[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)-totalrange4[i]
if (correction[i] < 0) {bootrange=bootrange-correction[i]}
obs.UCI95range4[i]=quantile(boot95range,c(0.975));  obs.LCI95range4[i]=quantile(boot95range,c(0.025))
obs.UCIrange4[i]=quantile(bootrange,c(0.975)); obs.LCIrange4[i]=quantile(bootrange,c(0.025))
obs.UCImedian4[i]=quantile(bootmedian,c(0.975)); obs.LCImedian4[i]=quantile(bootmedian,c(0.025))
obs.UCImean4[i]=quantile(bootmean,c(0.975)); obs.LCImean4[i]=quantile(bootmean,c(0.025))
obs.UCIIQR4[i]=quantile(bootIQR,c(0.975));  obs.LCIIQR4[i]=quantile(bootIQR,c(0.025))
obs.UCIskewness4[i]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness4[i]=quantile(bootskew,c(0.025), na.rm=T)
print(i)
}

tempage=unlist(split(age, Coredepth)[LIST])
samplesizes4[LIST]=length(unlist(split(age, Coredepth)[LIST]))
IQRrange4[LIST]=IQR(unlist(split(age, Coredepth)[LIST]))
fullrange4[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
medianage4[LIST]=median(unlist(split(age, Coredepth)[LIST]))
meanage4[LIST]=median(unlist(split(age, Coredepth)[LIST]))
minage4[LIST]=min(unlist(split(age, Coredepth)[LIST]))
maxage4[LIST]=max(unlist(split(age, Coredepth)[LIST]))
range954[LIST]=diff(quantile(unlist(split(age, Coredepth)[LIST]),c(0.025,0.975)))
totalrange4[LIST]=diff(range(unlist(split(age, Coredepth)[LIST])))
skew4[LIST]=skewness(unlist(split(age, Coredepth)[LIST]))
depths4[LIST]=names(split(age, Coredepth)[LIST])
station4[LIST]=unlist(split( Corestation, Coredepth)[LIST])[1]
interval4[LIST]=as.character(as.numeric(depths4[LIST]))
hist(unlist(split(age, Coredepth)[LIST]), xlab="Age (years)", ylab="", breaks=seq(0,15000,by=500), cex.main=0.8, cex.axis=0.8,
col="gray", ylim=c(0,20),main=paste(interval4[LIST],"cm","(n =",samplesizes4[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)-totalrange4[LIST]
if (correction[LIST] < 0) {bootrange=bootrange-correction[LIST]}
obs.UCI95range4[LIST]=quantile(boot95range,c(0.975));  obs.LCI95range4[LIST]=quantile(boot95range,c(0.025))
obs.UCIrange4[LIST]=quantile(bootrange,c(0.975)); obs.LCIrange4[LIST]=quantile(bootrange,c(0.025))
obs.UCImedian4[LIST]=quantile(bootmedian,c(0.975)); obs.LCImedian4[LIST]=quantile(bootmedian,c(0.025))
obs.UCImean4[LIST]=quantile(bootmean,c(0.975)); obs.LCImean4[LIST]=quantile(bootmean,c(0.025))
obs.UCIIQR4[LIST]=quantile(bootIQR,c(0.975));  obs.LCIIQR4[LIST]=quantile(bootIQR,c(0.025))
obs.UCIskewness4[LIST]=quantile(bootskew,c(0.975), na.rm=T); obs.LCIskewness4[LIST]=quantile(bootskew,c(0.025), na.rm=T)


###############################################################################################
#FIGURE PO AND PANZANO AGE SUMMARY
###############################################################################################
par(mfrow=c(2,3))
par(mar=c(4,4,2,3))
out2=boxplot(split(Corbula.Po3.ages,Corbula.Po3.5depths), at=-as.numeric(names(split(Corbula.Po3.ages,Corbula.Po3.5depths))), 
col="gray", range=0, boxwex=3, horizontal=TRUE, ylim=c(0,250), ylab="Sediment depth (cm)", xlim=c(-155,0),
xlab="Time since 2013 (y)", main="Po 3", cex.main=1, las=2, yaxt="n")
axis(4, at=-as.numeric(depths2), labels=round(2013-medianage2), las=1)
axis(2, at=-seq(0,150,by=25), labels=seq(0,150,by=25), las=2)
rect(xleft=0, ybottom=-25, xright=250, ytop=0,  col = "gray41", border = NA)
rect(xleft=0, ybottom=-90, xright=250, ytop=-25,  col = "gray51", border = NA)
rect(xleft=0, ybottom=-155, xright=250, ytop=-85,  col = "gray71", border = NA)
out2=boxplot(split(Corbula.Po3.ages,Corbula.Po3.5depths), at=-as.numeric(names(split(Corbula.Po3.ages,Corbula.Po3.5depths))), 
col="gray", range=0, boxwex=3, horizontal=TRUE, ylim=c(0,250), add=T,  xlim=c(-155,0), las=2, yaxt="n")
cbind(2013-out2$stats[2,], 2013-out2$stats[3,], 2013-out2$stats[4,])
text(190,-5, labels="Post-2000")
text(190,-30, labels="Late 20th c.")
text(190,-90, labels="Early 20th c.")

out3=boxplot(split(Corbula.Po4.ages,Corbula.Po4.5depths), at=-as.numeric(names(split(Corbula.Po4.ages,Corbula.Po4.5depths))), 
col="gray", range=0, boxwex=3, horizontal=TRUE, ylim=c(0,250), ylab="", xlim=c(-155,0), 
xlab="Time since 2013 (y)", main="Po 4", cex.main=1, las=2, yaxt="n")
axis(4, at=-as.numeric(depths3), labels=round(2013-medianage3), las=1)
rect(xleft=0, ybottom=-25, xright=250, ytop=0,  col = "gray41", border = NA)
rect(xleft=0, ybottom=-85, xright=250, ytop=-25,  col = "gray51", border = NA)
rect(xleft=0, ybottom=-155, xright=250, ytop=-85,  col = "gray71", border = NA)
out3=boxplot(split(Corbula.Po4.ages,Corbula.Po4.5depths), at=-as.numeric(names(split(Corbula.Po4.ages,Corbula.Po4.5depths))), 
col="gray", range=0, boxwex=3, horizontal=TRUE, ylim=c(0,250), add=T,  xlim=c(-155,0), las=2, yaxt="n")
cbind(2013-out3$stats[2,], 2013-out3$stats[3,], 2013-out3$stats[4,])
text(190,-5, labels="Post-2000")
text(190,-30, labels="Late 20th c.")
text(190,-90, labels="Early 20th c.")

out1=boxplot(split(Corbula.Panzano.ages,Corbula.Panzano.rev5depths), at=-as.numeric(names(split(Corbula.Panzano.ages,Corbula.Panzano.rev5depths))), 
col="gray", range=0, boxwex=3, horizontal=TRUE, ylim=c(0,500), ylab="", xlim=c(-155,0), 
xlab="Time since 2013 (y)", main="Panzano", cex.main=1, las=2, yaxt="n")
axis(4, at=-as.numeric(depths1), labels=round(2013-medianage1), las=1)
rect(xleft=0, ybottom=-4, xright=500, ytop=0,  col = "gray41", border = NA)
rect(xleft=0, ybottom=-16, xright=500, ytop=-4,  col = "gray51", border = NA)
rect(xleft=0, ybottom=-25, xright=500, ytop=-16,  col = "gray71", border = NA)
rect(xleft=0, ybottom=-45, xright=500, ytop=-25,  col = "gray81", border = NA)
rect(xleft=0, ybottom=-75, xright=500, ytop=-45,  col = "gray91", border = NA)
out1=boxplot(split(Corbula.Panzano.ages,Corbula.Panzano.rev5depths), at=-as.numeric(names(split(Corbula.Panzano.ages,Corbula.Panzano.rev5depths))), 
col="gray", range=0, boxwex=3, horizontal=TRUE, ylim=c(0,250), ylab="", add=T, xlim=c(-155,0), las=2, yaxt="n")
text(400,-10, labels="L. 20th c.")
text(400,-20, labels="E. 20th c.")
text(400,-35, labels="L. 19th c.")
text(400,-70, labels="E. 19th c.")


###############################################################################################################
#TIME AVERAGING PER UNIT (UNITS 1-3) 
###############################################################################################################
######
#PO 3#
######
Corbula.Po3.layers=Corbula.Po3.rev5depths
Corbula.Po3.layers[Corbula.Po3.rev5depths<21]=1
Corbula.Po3.layers[Corbula.Po3.rev5depths>21 & Corbula.Po3.rev5depths<86]=2
Corbula.Po3.layers[Corbula.Po3.rev5depths>86]=3
Po3.layer2.age=Corbula.Po3.ages[Corbula.Po3.ages>13 & Corbula.Po3.layers==2]
Po3.layer3.age=Corbula.Po3.ages[Corbula.Po3.ages>43 & Corbula.Po3.layers==3]
tapply(Corbula.Po3.ages,Corbula.Po3.layers, IQR)
tapply(Corbula.Po3.ages,Corbula.Po3.layers, function(x) {diff(quantile(x,c(0.025, 0.975)))})
tapply(Corbula.Po3.ages[Corbula.Po3.ages>13],Corbula.Po3.layers[Corbula.Po3.ages>13], IQR)[2]
tapply(Corbula.Po3.ages[Corbula.Po3.ages>43],Corbula.Po3.layers[Corbula.Po3.ages>43], IQR)[3]
Po3.layer2.95range=diff(unlist(tapply(Corbula.Po3.ages[Corbula.Po3.ages>13],Corbula.Po3.layers[Corbula.Po3.ages>13], quantile, c(0.025,0.975))[2]))
Po3.layer3.95range=diff(unlist(tapply(Corbula.Po3.ages[Corbula.Po3.ages>43],Corbula.Po3.layers[Corbula.Po3.ages>43], quantile, c(0.025,0.975))[3]))
Po3.layer2.iqr=tapply(Corbula.Po3.ages[Corbula.Po3.ages>13],Corbula.Po3.layers[Corbula.Po3.ages>13], IQR)[2]
Po3.layer3.iqr=tapply(Corbula.Po3.ages[Corbula.Po3.ages>43],Corbula.Po3.layers[Corbula.Po3.ages>43], IQR)[3]

IQR(Corbula.Po3.ages[Corbula.Po3.layers==2])
IQR(Corbula.Po3.ages[Corbula.Po3.layers==3])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.layers==2], c(0.025,0.975)))
diff(quantile(Corbula.Po3.ages[Corbula.Po3.layers==3], c(0.025,0.975)))

Po3.layer2.iqr=IQR(Corbula.Po3.ages[Corbula.Po3.ages>13 & Corbula.Po3.layers==2])
Po3.layer3.iqr=IQR(Corbula.Po3.ages[Corbula.Po3.ages>43 & Corbula.Po3.layers==3])
quantile(bootci.ext(Po3.layer2.age)$resampleiqr, c(0.025,0.5,0.975))
quantile(bootci.ext(Po3.layer3.age)$resampleiqr, c(0.025,0.5,0.975))
quantile(bootci.ext(Po3.layer2.age)$resamplerange, c(0.025,0.5,0.975))
quantile(bootci.ext(Po3.layer3.age)$resamplerange, c(0.025,0.5,0.975))

######
#PO 4#
######
Corbula.Po4.layers=Corbula.Po4.rev5depths
Corbula.Po4.layers[Corbula.Po4.rev5depths<21]=1
Corbula.Po4.layers[Corbula.Po4.rev5depths>21 & Corbula.Po4.rev5depths<86]=2
Corbula.Po4.layers[Corbula.Po4.rev5depths>86]=3
Po4.layer2.age=Corbula.Po4.ages[Corbula.Po4.ages>13 & Corbula.Po4.layers==2]
Po4.layer3.age=Corbula.Po4.ages[Corbula.Po4.ages>43 & Corbula.Po4.layers==3]
tapply(Corbula.Po4.ages,Corbula.Po4.layers, IQR)
tapply(Corbula.Po4.ages,Corbula.Po4.layers, function(x) {diff(quantile(x,c(0.025, 0.975)))})
tapply(Corbula.Po4.ages[Corbula.Po4.ages>13],Corbula.Po4.layers[Corbula.Po4.ages>13], IQR)[2]
tapply(Corbula.Po4.ages[Corbula.Po4.ages>43],Corbula.Po4.layers[Corbula.Po4.ages>43], IQR)[2]
diff(unlist(tapply(Corbula.Po4.ages[Corbula.Po4.ages>13],Corbula.Po4.layers[Corbula.Po4.ages>13], quantile, c(0.025,0.975))[2]))
diff(unlist(tapply(Corbula.Po4.ages[Corbula.Po4.ages>43],Corbula.Po4.layers[Corbula.Po4.ages>43], quantile, c(0.025,0.975))[2]))

Corbula.Panzano.layers=Corbula.Panzano.5depths
Corbula.Panzano.layers[Corbula.Panzano.5depths<14]=1
Corbula.Panzano.layers[Corbula.Panzano.5depths>14 & Corbula.Panzano.rev5depths<70]=2
Corbula.Panzano.layers[Corbula.Panzano.5depths>70]=3
tapply(Corbula.Panzano.ages,Corbula.Panzano.layers, IQR)
tapply(Corbula.Panzano.ages,Corbula.Panzano.layers, function(x) {diff(quantile(x,c(0.025, 0.975)))})
tapply(Corbula.Panzano.ages[Corbula.Panzano.ages>13],Corbula.Panzano.layers[Corbula.Panzano.ages>13], IQR)
tapply(Corbula.Panzano.ages[Corbula.Panzano.ages>43],Corbula.Panzano.layers[Corbula.Panzano.ages>43], function(x) {diff(quantile(x,c(0.025, 0.975)))})

Po4.layer2.iqr=IQR(Corbula.Po4.ages[Corbula.Po4.ages>13 & Corbula.Po4.layers==2])
Po4.layer3.iqr=IQR(Corbula.Po4.ages[Corbula.Po4.ages>43 & Corbula.Po4.layers==3])
quantile(bootci.ext(Po4.layer2.age)$resampleiqr, c(0.025,0.975))
quantile(bootci.ext(Po4.layer3.age)$resampleiqr, c(0.025,0.975))
quantile(bootci.ext(Po4.layer2.age)$resamplerange, c(0.025,0.975))
quantile(bootci.ext(Po4.layer3.age)$resamplerange, c(0.025,0.975))


###############################################################################################################
#TIME AVERAGING PER INCREMENT, TRUNCATING POSTDEPOSITIONAL COHORTS
###############################################################################################################
#PANZANO
###############################################################################################################
IQR1=round(tapply(Corbula.Panzano.ages,Corbula.Panzano.rev5depths,IQR))
mean1=round(tapply(Corbula.Panzano.ages,Corbula.Panzano.rev5depths,mean))
median1=round(tapply(Corbula.Panzano.ages,Corbula.Panzano.rev5depths,median))
LCI1=round(tapply(Corbula.Panzano.ages,Corbula.Panzano.rev5depths,quantile,c(0.025)))
UCI1=round(tapply(Corbula.Panzano.ages,Corbula.Panzano.rev5depths,quantile,c(0.975)))
pre.depth1=(tapply(Corbula.Panzano.rev5depths,Corbula.Panzano.rev5depths,function (x) x[1]))

preemption=c(rep(13,3),rep(43,4),rep(43,4))
list1=split(Corbula.Panzano.ages, Corbula.Panzano.rev5depths)
pre.range951=rep(NA,length(list1))
pre.IQR1=rep(NA,length(list1))
for (i in 1:length(list1)) {
temp=unlist(list1[i])
temp=temp[temp>preemption[i]]
pre.range951[i]=diff(quantile(temp,c(0.025,0.975)))
pre.IQR1[i]=IQR(temp)
}

SIM.meaniqr.x.age=array(NA, dim=c(length(list1),100)); SIM.mean95.x.age=array(NA, dim=c(length(list1),100)); 
for (j in 1:100) {
meaniqr.x.age=numeric();  mean95.x.age=numeric(); 
for (i in 1:length(list1)) {
ages=unlist(list1[i])
ages=ages[ages>preemption[i]]
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), 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)
pre.cor.range951=pre.range951-SIM.mean.95
pre.cor.range951.LCI=pre.range951-apply(SIM.mean95.x.age,1,quantile, 0.025)
pre.cor.range951.UCI=pre.range951-apply(SIM.mean95.x.age,1,quantile, 0.975)

pre.cor.IQR1=pre.IQR1-SIM.mean.iqr
pre.cor.IQR1.LCI=pre.IQR1-apply(SIM.meaniqr.x.age,1,quantile, 0.025)
pre.cor.IQR1.UCI=pre.IQR1-apply(SIM.meaniqr.x.age,1,quantile, 0.975)


###############################################################################################################
#PO 3
###############################################################################################################
IQR2=round(tapply(Corbula.Po3.ages,Corbula.Po3.rev5depths,IQR))
mean2=round(tapply(Corbula.Po3.ages,Corbula.Po3.rev5depths,mean))
median2=round(tapply(Corbula.Po3.ages,Corbula.Po3.rev5depths,median))
LCI2=round(tapply(Corbula.Po3.ages,Corbula.Po3.rev5depths,quantile,c(0.025)))
UCI2=round(tapply(Corbula.Po3.ages,Corbula.Po3.rev5depths,quantile,c(0.975)))
pre.depth2=round(tapply(Corbula.Po3.rev5depths,Corbula.Po3.rev5depths,function (x) x[1]))

preemption=c(rep(0,5),rep(13,4),rep(43,3))
list1=split(Corbula.Po3.ages, Corbula.Po3.rev5depths)
pre.range952=rep(NA,length(list1))
pre.IQR2=rep(NA,length(list1))
for (i in 1:length(list1)) {
temp=unlist(list1[i])
temp=temp[temp>preemption[i]]
pre.range952[i]=diff(quantile(temp,c(0.025,0.975)))
pre.IQR2[i]=IQR(temp)
}

SIM.meaniqr.x.age=array(NA, dim=c(length(list1),100)); SIM.medianiqr.x.age=array(NA, dim=c(length(list1),100));  SIM.mean95.x.age=array(NA, dim=c(length(list1),100)); 
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:length(list1)) {
ages=unlist(list1[i])
ages=ages[ages>preemption[i]]
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), 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)
pre.cor.range952=pre.range952-SIM.mean.95
pre.cor.range952.LCI=pre.range952-apply(SIM.mean95.x.age,1,quantile, 0.025)
pre.cor.range952.UCI=pre.range952-apply(SIM.mean95.x.age,1,quantile, 0.975)

pre.cor.IQR2=pre.IQR2-SIM.mean.iqr
pre.cor.IQR2.LCI=pre.IQR2-apply(SIM.meaniqr.x.age,1,quantile, 0.025)
pre.cor.IQR2.UCI=pre.IQR2-apply(SIM.meaniqr.x.age,1,quantile, 0.975)

###############################################################################################################
#PO 4
###############################################################################################################
IQR3=round(tapply(Corbula.Po4.ages,Corbula.Po4.rev5depths,IQR))
mean3=round(tapply(Corbula.Po4.ages,Corbula.Po4.rev5depths,mean))
median3=round(tapply(Corbula.Po4.ages,Corbula.Po4.rev5depths,median))
LCI3=round(tapply(Corbula.Po4.ages,Corbula.Po4.rev5depths,quantile,c(0.025)))
UCI3=round(tapply(Corbula.Po4.ages,Corbula.Po4.rev5depths,quantile,c(0.975)))
pre.depth3=round(tapply(Corbula.Po4.rev5depths,Corbula.Po4.rev5depths,function (x) x[1]))

preemption=c(rep(0,5),rep(13,4),rep(43,4))
list1=split(Corbula.Po4.ages, Corbula.Po4.rev5depths)
pre.range953=rep(NA,length(list1))
pre.IQR3=rep(NA,length(list1))
for (i in 1:length(list1)) {
temp=unlist(list1[i])
temp=temp[temp>preemption[i]]
pre.range953[i]=diff(quantile(temp,c(0.025,0.975)))
pre.IQR3[i]=IQR(temp)
}

SIM.meaniqr.x.age=array(NA, dim=c(length(list1),100)); SIM.medianiqr.x.age=array(NA, dim=c(length(list1),100));  SIM.mean95.x.age=array(NA, dim=c(length(list1),100)); 
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:length(list1)) {
ages=unlist(list1[i])
ages=ages[ages>preemption[i]]
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), 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)
pre.cor.range953=pre.range953-SIM.mean.95
pre.cor.range953.LCI=pre.range953-apply(SIM.mean95.x.age,1,quantile, 0.025)
pre.cor.range953.UCI=pre.range953-apply(SIM.mean95.x.age,1,quantile, 0.975)

pre.cor.IQR3=pre.IQR3-SIM.mean.iqr
pre.cor.IQR3.LCI=pre.IQR3-apply(SIM.meaniqr.x.age,1,quantile, 0.025)
pre.cor.IQR3.UCI=pre.IQR3-apply(SIM.meaniqr.x.age,1,quantile, 0.975)


######################################################################################################################
#CORRECTON FOR TIME AVERAGING PER INCREMENT, NOT TRUNCATING POSTDEPOSITIONAL COHORTS
######################################################################################################################
#########################################################
#PANZANO
#########################################################
list1=split(Corbula.Panzano.ages, Corbula.Panzano.rev5depths)
SIM.meaniqr.x.age=array(NA, dim=c(length(list1),100)); SIM.medianiqr.x.age=array(NA, dim=c(length(list1),100));  SIM.mean95.x.age=array(NA, dim=c(length(list1),100)); 
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:length(list1)) {
ages=unlist(list1[i])
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), c(0.025,0.975)))
}
meaniqr.x.age[i]=mean(x.age); medianiqr.x.age[i]=median(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)
cor.range951=range951-SIM.mean.95
cor.range951.LCI=range951-apply(SIM.mean95.x.age,1,quantile, 0.025)
cor.range951.UCI=range951-apply(SIM.mean95.x.age,1,quantile, 0.975)
cor.IQR1=IQR1-SIM.mean.iqr
cor.IQR1.LCI=IQR1-apply(SIM.meaniqr.x.age,1,quantile, 0.025)
cor.IQR1.UCI=IQR1-apply(SIM.meaniqr.x.age,1,quantile, 0.975)

#########################################################
#P0 3
#########################################################
list1=split(Corbula.Po3.ages, Corbula.Po3.rev5depths)
out=lapply(list1, quantile, c(0.025,0.975))
t.range952=unlist(lapply(out, diff))
t.IQR2=unlist(lapply(list1, IQR))

SIM.meaniqr.x.age=array(NA, dim=c(length(list1),100)); SIM.medianiqr.x.age=array(NA, dim=c(length(list1),100));  SIM.mean95.x.age=array(NA, dim=c(length(list1),100)); 
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:length(list1)) {
ages=unlist(list1[i])
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), 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)

cor.range952=t.range952-SIM.mean.95
cor.range952.LCI=t.range952-apply(SIM.mean95.x.age,1,quantile, 0.025)
cor.range952.UCI=t.range952-apply(SIM.mean95.x.age,1,quantile, 0.975)
cor.IQR2=t.IQR2-SIM.mean.iqr
cor.IQR2.LCI=t.IQR2-apply(SIM.meaniqr.x.age,1,quantile, 0.025)
cor.IQR2.UCI=t.IQR2-apply(SIM.meaniqr.x.age,1,quantile, 0.975)

#########################################################
#P0 4
#########################################################
list1=split(Corbula.Po4.ages, Corbula.Po4.rev5depths)
out=lapply(list1, quantile, c(0.025,0.975))
t.range953=unlist(lapply(out, diff))
t.IQR3=unlist(lapply(list1, IQR))

SIM.meaniqr.x.age=array(NA, dim=c(length(list1),100));  SIM.mean95.x.age=array(NA, dim=c(length(list1),100)); 
for (j in 1:100) {
meaniqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:length(list1)) {
ages=unlist(list1[i])
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), 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)
cor.range953=t.range953-SIM.mean.95
cor.range953.LCI=t.range953-apply(SIM.mean95.x.age,1,quantile, 0.025)
cor.range953.UCI=t.range953-apply(SIM.mean95.x.age,1,quantile, 0.975)

cor.IQR3=t.IQR3-SIM.mean.iqr
cor.IQR3.LCI=t.IQR3-apply(SIM.meaniqr.x.age,1,quantile, 0.025)
cor.IQR3.UCI=t.IQR3-apply(SIM.meaniqr.x.age,1,quantile, 0.975)

##################################################################################################################################
depth1=round(tapply(Corbula.Panzano.5depths,Corbula.Panzano.rev5depths,function (x) x[1]))
depth2=(tapply(Corbula.Po3.5depths,Corbula.Po3.rev5depths,function (x) x[1]))
depth3=(tapply(Corbula.Po4.5depths,Corbula.Po4.rev5depths,function (x) x[1]))

med1=round(tapply(Corbula.Panzano.ages,Corbula.Panzano.rev5depths,median))
med2=round(tapply(Corbula.Po3.ages,Corbula.Po3.rev5depths,median))
med3=round(tapply(Corbula.Po4.ages,Corbula.Po4.rev5depths,median))

par(mfrow=c(2,2))
plot(cor.IQR2, -as.numeric(names(cor.IQR2)), pch=21, bg="gray", log="", ylab="Sediment depth",xlim=c(0.5,100), xlab="Inter-quartile age range (y)", main="Po 3")
segments(x0=cor.IQR2.LCI,x1=cor.IQR2.UCI,y0=-as.numeric(names(cor.IQR2)), y1=-as.numeric(names(cor.IQR2)))
points(cor.IQR2, -as.numeric(names(cor.IQR2)), pch=21,bg="gray")
points(IQR2, -as.numeric(names(cor.IQR2)), pch=21,bg="white")

plot(cor.IQR3, -as.numeric(names(cor.IQR3)), pch=21, bg="gray", log="", xlim=c(0.5,100), xlab="Inter-quartile age range (y)", main="Po 4")
segments(x0=cor.IQR3.LCI,x1=cor.IQR3.UCI,y0=-as.numeric(names(cor.IQR3)), y1=-as.numeric(names(cor.IQR3)))
points(cor.IQR3, -as.numeric(names(cor.IQR3)), pch=21,bg="gray")
points(IQR3, -as.numeric(names(cor.IQR3)), pch=21,bg="white")

plot(cor.IQR1, -depth1, pch=21, bg="gray", log="",  xlim=c(0.5,500), xlab="Inter-quartile age range (y)", main="Panzano")
segments(x0=cor.IQR1.LCI,x1=cor.IQR1.UCI,y0=-depth1, y1=-depth1)
points(cor.IQR1, -depth1, pch=21,bg="gray")
points(IQR1, -depth1, pch=21,bg="white")

##################################################################################################################################
#FIGURE - DOWNCORE CHANGES IN TIME AVERAGING
##################################################################################################################################
top=c(cor.range952[depth2<20],cor.range953[depth3<20])
mid=c(cor.range952[depth2>20 & depth2<90],cor.range953[depth3>20 & depth3<90])
bottom=c(cor.range952[depth2>90],cor.range953[depth3>90])

par(mfrow=c(2,4))
boxplot(cor.range952[depth2<20],cor.range952[depth2>20 & depth2<90],cor.range952[depth2>90], ylim=c(0,200))
boxplot(cor.range953[depth3<20],cor.range953[depth3>20 & depth3<90],cor.range953[depth3>90], ylim=c(0,200))
boxplot(cor.IQR2[depth2<20],cor.IQR2[depth2>20 & depth2<90],cor.IQR2[depth2>90], ylim=c(0,200))
boxplot(cor.IQR3[depth3<20],cor.IQR3[depth3>20 & depth3<90],cor.IQR3[depth3>90], ylim=c(0,200))

LOG=""
par(mfrow=c(2,3))
par(mar=c(4,4,2,1))
plot(cor.IQR2, -depth2, pch=21, bg="gray", log=LOG, ylab="Sediment depth", xlim=c(0.5,200), ylim=c(-155,0), xlab="Age range (y)", main="Po 3")
segments(x0=cor.IQR2.LCI,x1=cor.IQR2.UCI,y0=-depth2, y1=-depth2)
segments(x0=cor.range952.LCI,x1=cor.range952.UCI,y0=-depth2, y1=-depth2)
points(cor.IQR2, -depth2, pch=21,bg="gray", cex=1.2)
points(cor.range952, -depth2, pch=21,bg="black", cex=1.2)
legend(x="topright", pch=21, pt.bg=c("gray","black"), legend=c("IQR","95% range"))
plot(cor.IQR3, -depth3, pch=21, bg="gray", log=LOG, xlim=c(0.5,200), ylim=c(-155,0), ylab="", xlab="Age range (y)", main="Po 4")
segments(x0=cor.IQR3.LCI,x1=cor.IQR3.UCI,y0=-depth3, y1=-depth3)
segments(x0=cor.range953.LCI,x1=cor.range953.UCI,y0=-depth3, y1=-depth3)
points(cor.IQR3, -depth3, pch=21,bg="gray", cex=1.2)
points(cor.range953, -depth3, pch=21,bg="black", cex=1.2)
plot(cor.IQR1, -depth1, pch=21, bg="gray", log=LOG, xlim=c(0.5,300), ylim=c(-155,0), ylab="", xlab="Age range (y)", main="Panzano")
segments(x0=cor.IQR1.LCI,x1=cor.IQR1.UCI,y0=-depth1, y1=-depth1)
segments(x0=cor.range951.LCI,x1=cor.range951.UCI,y0=-depth1, y1=-depth1)
points(cor.IQR1, -depth1, pch=21,bg="gray", cex=1.2)
points(cor.range951, -depth1, pch=21,bg="black", cex=1.2)


#####################################################################################
#TOTAL ABUNDANCES OF CORBULA IN INCREMENTS, ORDERED FROM THE UPPERMOST TO LOWERMOST INCREMENT
#####################################################################################
full.totalsizePanzano=c(83, 117, 100, 134,  79,  59,  86,  74,  96, 139,  71,  50,  34,  43,  40,  27,  32,  11,  21,  26,  39,  66,  48,  49,  47,  45, 48,  76,  50,  26,  37)
full.totalsizePo3=c(4, 13, 32, 42, 43, 53, 46, 76, 30, 55, 38, 40, 49, 41, 43, 74, 44, 36, 37, 22, 43, 19, 17, 12, 13, 11,  4, 11,  9,  8,  5,  2)
full.totalsizePo4=c(15,  40,   8,  15,  23,  80, 106,  79,  84,  54, 102, 113,  28, 100,  61,  75, 103,  47,  15,  49,  20,  25,  29,  19,  24,  15, 12,  13,   3,   6,  3,   7)
full.totalsizePiran=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)

#######################################################################################
#CORE AGE DISTRIBUTION ACCOUNTING FOR ABUNDANCES AND INTERPOLATING UNDATED INCREMENTS 
#FIGURE AGE-FREQUENCY DISTRIBUTIONS ACC STRATIGRAPHIC LAYERS
#######################################################################################################

surfaceAFD.Panzano.6=c(Corbula.Panzano.ages[Corbula.Panzano.maxdepths<8])
surfaceAFD.Po3.20=c(Corbula.Po3.ages[Corbula.Po3.5depths<21])
surfaceAFD.Po4.20=c(Corbula.Po4.ages[Corbula.Po4.5depths<21])
surfaceAFD.Piran.4=Corbula.Piran.ages[Corbula.Piran.5depths==4]
###############################################
#ONE-PHASE EXPONENTIAL AND WEIBULL MODEL#
###############################################

expAIC=numeric()
weibullAIC=numeric()
likratio=numeric()
onephaserate=numeric()
weibullrate=numeric()
weibullshape=numeric()
randomAIC=numeric()
tau=numeric()
lambda1=numeric()
lambda2=numeric()
skew=numeric()
iqr.layer=numeric()
range95.layer=numeric()

par(mfrow=c(3,3))
par(mar=c(4,4,2,1))
for (i in 1:3) {
if (i==1) {age=surfaceAFD.Po3.20; thickness=20}
if (i==2) {age=surfaceAFD.Po4.20; thickness=20}
if (i==3) {age=surfaceAFD.Panzano.6; thickness=6}

skew[i]=skewness(age)
iqr.layer[i]=IQR(age)
range95.layer[i]=diff(quantile(age, c(0.025,0.975)))

onephaserate[i]<-1/mean(age)
like<--(sum(log(onephaserate[i])-onephaserate[i]*age))	
AIC=2*1-(-2*like)
AICc=AIC+((2*1*(1+1))/(length(age)-1-1))
expAIC[i]<-AICc
onephaseAIC<-AICc

start1=c(log(0.0005),log(0.1))
fit1=optim(par=start1, fn=weibull)
weibullrate[i]<-exp(fit1$par[1])
weibullshape[i]<-exp(fit1$par[2])
likweibull<-fit1$value
AIC=2*2-(-2*likweibull)
AICc=AIC+((2*2*(2+1))/(length(age)-2-1))
weibullAIC[i]<-AICc
#likelihood ratio test comparing weibull model with one-phase exponential
dev<-2*(like-likweibull)
likratio[i]<-1-pchisq(dev,1)

#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[i]=rAIC+((2*3*(3+1))/(length(age)-3-1))
tau[i]<-randomtwophasealpha*(randomtwophaserate2-randomtwophaserate1)
lambda1[i]<-randomtwophaserate1*randomtwophasealpha+randomtwophaserate2*(1-randomtwophasealpha)
lambda2[i]<-randomtwophaserate1

LAB=""
if (i==1) {LAB="Number of specimens"}

hist(age, col="white",breaks=seq(0,300,by=5), xlab="", ylab=LAB,ylim=c(0,40), main="", cex.main=0.9, bty="u")
times=seq(0,2000, by=5)+2.5
lambda=onephaserate[i]
predictSIMPLE<-dexp(times,rate=lambda)
lines(times, (predictSIMPLE/sum(predictSIMPLE))*length(age), lty=1, col="black", lwd=3)
#predictWEIBULL<-(1/((1/weibullrate[i])*gamma(1+1/weibullshape[i])))*exp(-(weibullrate[i]*times)^weibullshape[i])
#lines(times, (predictWEIBULL/sum(predictWEIBULL))*length(age), lty=1, col="gray", lwd=2)
predictRANDOM<-randomtwophaserate1*exp(-randomtwophaserate1*times)*randomtwophasebeta+randomtwophaserate2*exp(-randomtwophaserate2*times)*(1-randomtwophasebeta)
lines(times, (predictRANDOM/sum(predictRANDOM))*length(age), col="gray31", lwd=3, lty=1)
#predictRANDOMTRUNCATED2=(proportion*dnorm(times,mean=mean1,sd=sqrt(variance)))/pnorm(mean1/sqrt(variance),0,1)+((1-proportion)*dnorm(times,mean=mean2,sd=sqrt(variance)))/pnorm(mean2/sqrt(variance),0,1)
#lines(times, predictRANDOMTRUNCATED2/sum(predictRANDOMTRUNCATED2)*length(age), col="red",lty=1, lwd=3)
abline(v=median(age), col="gray31")
text(150,40,labels=paste("FML thickness =","",thickness,"","cm"))
text(150,35,labels=paste("FML shell loss rate =","",round(onephaserate[i], digits=3)))
text(150,30,labels=paste("Time to sh. loss =","",round(1/onephaserate[i], digits=1),"","yr"))
text(150,25,labels=paste("burial =","",round(thickness/(1/onephaserate[i]), digits=1),"","cm/y"))
}

cbind(expAIC,weibullAIC,randomAIC)

skew
iqr.layer
range95.layer

BINS=5
#MIDDLE INTERVALS
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90], ylab="Number of specimens", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Interval 2 - 20-85 cm", ylim=c(0,35), cex.main=1)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90]), col="gray31")
skewness(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90])
IQR(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Interval 2 - 20-85 cm", ylim=c(0,35), cex.main=1)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90]), col="gray31")
skewness(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90])
IQR(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90], c(0.025,0.975)))

out=hist(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Interval 2 - 6-14 cm", ylim=c(0,15), cex.main=1)
abline(v=median(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16]), col="gray31")
skewness(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16])
IQR(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16])
diff(quantile(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16], c(0.025,0.975)))

#LOWER INTERVALS
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths>90], ylab="Number of specimens",col="gray", breaks=seq(0,300,by=BINS), xlab="Postmortem age (y)", main="Interval 1 - 85-150 cm", ylim=c(0,35), cex.main=1)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths>90]), col="gray31")
skewness(Corbula.Po3.ages[Corbula.Po3.5depths>90])
IQR(Corbula.Po3.ages[Corbula.Po3.5depths>90])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>90], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths>90], ylab="",col="gray", breaks=seq(0,300,by=BINS), xlab="Postmortem age (y)", main="Interval 1 - 85-150 cm", ylim=c(0,35), cex.main=1)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths>90]), col="gray31")
#axis(1, at=c(13,23,33,43,53,63,73,83,93,103,113), labels=rev(seq(1700,2000,by=BINS)), las=2)
#axis(1, at=c(13,23,33,43,53,63,73,83,93,103,113), labels=rev(seq(1700,2000,by=BINS)), las=2)
skewness(Corbula.Po4.ages[Corbula.Po4.5depths>90])
IQR(Corbula.Po4.ages[Corbula.Po4.5depths>90])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>90], c(0.025,0.975)))

temp=Corbula.Panzano.ages[Corbula.Panzano.5depths>16  & Corbula.Panzano.5depths<35]
out=hist(temp[temp<300], ylab="",col="gray", breaks=seq(0,300,by=BINS),  xlab="Postmortem age (y)", main="Interval 1 - 14-35 cm", ylim=c(0,15), cex.main=1)
abline(v=median(temp), col="gray31")
#axis(1, at=c(13,23,33,43,53,63,73,83,93,103,113), labels=rev(seq(1700,2000,by=BINS)), las=2)
skewness(temp)
IQR(temp)
diff(quantile(temp, c(0.025,0.975)))

######################################################################################################################
#FIGURE - time averaging in subintervals
######################################################################################################################
skew.sublayer.Po3=numeric()
iqr.sublayer.Po3=numeric()
range95.sublayer.Po3=numeric()

par(mfrow=c(3,2))
par(mar=c(4,4,2,1))
for (i in 1:5) {
if (i==1) {age=Corbula.Po3.ages[Corbula.Po3.5depths<20]}
if (i==2) {age=Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<50]}
if (i==3) {age=Corbula.Po3.ages[Corbula.Po3.5depths>50 & Corbula.Po3.5depths<85]}
if (i==4) {age=Corbula.Po3.ages[Corbula.Po3.5depths>85 & Corbula.Po3.5depths<125]}
if (i==5) {age=Corbula.Po3.ages[Corbula.Po3.5depths>125]}
skew.sublayer.Po3[i]=skewness(age)
iqr.sublayer.Po3[i]=IQR(age)
range95.sublayer.Po3[i]=diff(quantile(age, c(0.025,0.975)))
}
skew.sublayer.Po4=numeric()
iqr.sublayer.Po4=numeric()
range95.sublayer.Po4=numeric()
for (i in 1:5) {
if (i==1) {age=Corbula.Po4.ages[Corbula.Po4.5depths<20]}
if (i==2) {age=Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<50]}
if (i==3) {age=Corbula.Po4.ages[Corbula.Po4.5depths>50 & Corbula.Po4.5depths<85]}
if (i==4) {age=Corbula.Po4.ages[Corbula.Po4.5depths>85 & Corbula.Po4.5depths<125]}
if (i==5) {age=Corbula.Po4.ages[Corbula.Po4.5depths>125]}
skew.sublayer.Po4[i]=skewness(age)
iqr.sublayer.Po4[i]=IQR(age)
range95.sublayer.Po4[i]=diff(quantile(age, c(0.025,0.975)))
}

par(mfrow=c(5,3))
BINS=5
#SUBINTERVALS 5
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths<20], ylab="N specimens", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 3 - < 20 cm (21st century)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths<20]), col="gray31")
IQR(Corbula.Po3.ages[Corbula.Po3.5depths<20])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths<20], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths<20], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 4 - < 20 cm (21st century)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths<0]), col="gray31")
IQR(Corbula.Po4.ages[Corbula.Po4.5depths<20])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths<20], c(0.025,0.975)))

plot.new()

#SUBINTERVALS 4
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<50], ylab="N specimens", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 3 - 20-50 cm (~1975-2000 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<50]), col="gray31")
IQR(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<50])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<50], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<50], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 4 - 20-50 cm (~1975-2000 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<50]), col="gray31")
IQR(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<50])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<50], c(0.025,0.975)))
plot.new()

#SUBINTERVALS 3
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths>50 & Corbula.Po3.5depths<85], ylab="N specimens", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 3 - 50-85 cm (~1950-1975 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths>50 & Corbula.Po3.5depths<85]), col="gray31")
IQR(Corbula.Po3.ages[Corbula.Po3.5depths>50 & Corbula.Po3.5depths<85])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>50 & Corbula.Po3.5depths<85], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<85], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 4 - 50-85 cm (~1950-1975 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths>50 & Corbula.Po4.5depths<85]), col="gray31")
IQR(Corbula.Po4.ages[Corbula.Po4.5depths>50 & Corbula.Po4.5depths<85])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>50 & Corbula.Po4.5depths<85], c(0.025,0.975)))
plot.new()

#SUBINTERVALS 2
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths>85 & Corbula.Po3.5depths<125], ylab="N specimens", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 3 - 85-125 cm (~1925-1950 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths>85 & Corbula.Po3.5depths<125]), col="gray31")
IQR(Corbula.Po3.ages[Corbula.Po3.5depths>85 & Corbula.Po3.5depths<125])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>85 & Corbula.Po3.5depths<125], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths>85 & Corbula.Po4.5depths<125], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 4 - 85-125 cm (~1925-1950 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths>85 & Corbula.Po4.5depths<125]), col="gray31")
IQR(Corbula.Po4.ages[Corbula.Po4.5depths>85 & Corbula.Po4.5depths<125])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>85 & Corbula.Po4.5depths<125], c(0.025,0.975)))
plot.new()

#SUBINTERVALS 1
out=hist(Corbula.Po3.ages[Corbula.Po3.5depths>125], ylab="N specimens", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 3 - 125-150 cm (~1900-1925 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po3.ages[Corbula.Po3.5depths>125]), col="gray31")
IQR(Corbula.Po3.ages[Corbula.Po3.5depths>125])
diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>125], c(0.025,0.975)))

out=hist(Corbula.Po4.ages[Corbula.Po4.5depths>125], ylab="", col="gray", breaks=seq(0,300,by=BINS), xlab="", main="Po 4 - 125-150 cm (~1900-1925 AD)", ylim=c(0,35), cex.main=0.9)
abline(v=median(Corbula.Po4.ages[Corbula.Po4.5depths>125]), col="gray31")
IQR(Corbula.Po4.ages[Corbula.Po4.5depths>125])
diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>125], c(0.025,0.975)))
plot.new()

######################################################################################################################
#TIME AVERAGING FOR THREE UNITS
######################################################################################################################
#interquartile range
temp1=IQR(surfaceAFD.Po3.20)
temp2=IQR(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90])
temp3=IQR(Corbula.Po3.ages[Corbula.Po3.5depths>90])
temp4=IQR(surfaceAFD.Po4.20)
temp5=IQR(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90])
temp6=IQR(Corbula.Po4.ages[Corbula.Po4.5depths>90])
temp7=IQR(surfaceAFD.Panzano.6)
temp8=IQR(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16])
temp9=IQR(Corbula.Panzano.ages[Corbula.Panzano.5depths>16  & Corbula.Panzano.5depths<36])
int.iqr=c(temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,temp9)

#95% range
temp1=diff(quantile(surfaceAFD.Po3.20, c(0.025,0.975)))
temp2=diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90], c(0.025,0.975)))
temp3=diff(quantile(Corbula.Po3.ages[Corbula.Po3.5depths>90], c(0.025,0.975)))
temp4=diff(quantile(surfaceAFD.Po4.20, c(0.025,0.975)))
temp5=diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90], c(0.025,0.975)))
temp6=diff(quantile(Corbula.Po4.ages[Corbula.Po4.5depths>90], c(0.025,0.975)))
temp7=diff(quantile(surfaceAFD.Panzano.6, c(0.025,0.975)))
temp8=diff(quantile(Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16], c(0.025,0.975)))
temp9=diff(quantile(Corbula.Panzano.ages[Corbula.Panzano.5depths>16  & Corbula.Panzano.5depths<36], c(0.025,0.975)))
int.range95=c(temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,temp9)

#CORRECTION FOR CALIBRATION ERROR
SIM.meaniqr.x.age=array(NA, dim=c(9,100)); 
SIM.mean95.x.age=array(NA, dim=c(9,100)); 
for (j in 1:100) {
meaniqr.x.age=numeric(); medianiqr.x.age=numeric(); mean95.x.age=numeric(); 
for (i in 1:9) {
if (i == 1) {ages=surfaceAFD.Po3.20}
if (i == 2) {ages=Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90]}
if (i == 3) {ages=Corbula.Po3.ages[Corbula.Po3.5depths>90]}
if (i == 4) {ages=surfaceAFD.Po4.20}
if (i == 5) {ages=Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90]}
if (i == 6) {ages=Corbula.Po4.ages[Corbula.Po4.5depths>90]}
if (i == 7) {ages=surfaceAFD.Panzano.6}
if (i == 8) {ages=Corbula.Panzano.ages[Corbula.Panzano.5depths>4 & Corbula.Panzano.5depths<16]}
if (i == 9) {ages=Corbula.Panzano.ages[Corbula.Panzano.5depths>16  & Corbula.Panzano.5depths<36]}
x.age=numeric()
x.age.95=numeric()
for (x in 1:length(ages)) {
x.age[x]=IQR(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))))
x.age.95[x]=diff(quantile(rlnorm(100,log(ages[x]),sdlog=sqrt(exp(ln.d))), c(0.025,0.975)))
}
meaniqr.x.age[i]=mean(x.age); medianiqr.x.age[i]=median(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)
int.cor.range95=int.range95-SIM.mean.95
int.cor.IQR=int.iqr-SIM.mean.iqr
int.cor.IQR

#######################################################################################################
#Disorder - Spearman correlation
#######################################################################################################
#DO NOT TRUNCATE YOUNG AGES
cor1=cor.test(Corbula.Po4.ages[Corbula.Po4.5depths<20],Corbula.Po4.5depths[Corbula.Po4.5depths<20], method="s")
cor2=cor.test(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90],Corbula.Po4.5depths[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90], method="s")
cor3=cor.test(Corbula.Po4.ages[Corbula.Po4.5depths>90],Corbula.Po4.5depths[Corbula.Po4.5depths>90], method="s")
cor4=cor.test(Corbula.Po3.ages[Corbula.Po3.5depths<20],Corbula.Po3.5depths[Corbula.Po3.5depths<20], method="s")
cor5=cor.test(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90],Corbula.Po3.5depths[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90], method="s")
cor6=cor.test(Corbula.Po3.ages[Corbula.Po3.5depths>90],Corbula.Po3.5depths[Corbula.Po3.5depths>90], method="s")
cor7=cor.test(Corbula.Panzano.ages[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14],Corbula.Panzano.5depths[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14], method="s")
cor8=cor.test(Corbula.Panzano.ages[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40],Corbula.Panzano.5depths[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40], method="s")
cor1conf=cor.ci(cbind(Corbula.Po4.ages[Corbula.Po4.5depths<20],Corbula.Po4.5depths[Corbula.Po4.5depths<20]), method="s")$ci[c(1,4)]
cor2conf=cor.ci(cbind(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90],Corbula.Po4.5depths[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90]), method="s")$ci[c(1,4)]
cor3conf=cor.ci(cbind(Corbula.Po4.ages[Corbula.Po4.5depths>90],Corbula.Po4.5depths[Corbula.Po4.5depths>90]), method="s")$ci[c(1,4)]
cor4conf=cor.ci(cbind(Corbula.Po3.ages[Corbula.Po3.5depths<20],Corbula.Po3.5depths[Corbula.Po3.5depths<20]), method="s")$ci[c(1,4)]
cor5conf=cor.ci(cbind(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90],Corbula.Po3.5depths[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90]), method="s")$ci[c(1,4)]
cor6conf=cor.ci(cbind(Corbula.Po3.ages[Corbula.Po3.5depths>90],Corbula.Po3.5depths[Corbula.Po3.5depths>90]), method="s")$ci[c(1,4)]
cor7conf=cor.ci(cbind(Corbula.Panzano.ages[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14],Corbula.Panzano.5depths[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14]), method="s")$ci[c(1,4)]
cor8conf=cor.ci(cbind(Corbula.Panzano.ages[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40],Corbula.Panzano.5depths[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40]), method="s")$ci[c(1,4)]

#TRUNCATE YOUNG AGES
cor1=cor.test(Corbula.Po4.ages[Corbula.Po4.5depths<20],Corbula.Po4.5depths[Corbula.Po4.5depths<20], method="s")
cor2=cor.test(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90 & Corbula.Po4.ages > 13],Corbula.Po4.5depths[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90 & Corbula.Po4.ages > 13], method="s")
cor3=cor.test(Corbula.Po4.ages[Corbula.Po4.5depths>90 & Corbula.Po4.ages > 43],Corbula.Po4.5depths[Corbula.Po4.5depths>90 & Corbula.Po4.ages > 43], method="s")
cor4=cor.test(Corbula.Po3.ages[Corbula.Po3.5depths<20],Corbula.Po3.5depths[Corbula.Po3.5depths<20], method="s")
cor5=cor.test(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90 & Corbula.Po3.ages > 13],Corbula.Po3.5depths[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90 & Corbula.Po3.ages > 13], method="s")
cor6=cor.test(Corbula.Po3.ages[Corbula.Po3.5depths>90 & Corbula.Po3.ages > 43],Corbula.Po3.5depths[Corbula.Po3.5depths>90 & Corbula.Po3.ages > 43], method="s")
cor7=cor.test(Corbula.Panzano.ages[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14],Corbula.Panzano.5depths[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14], method="s")
cor8=cor.test(Corbula.Panzano.ages[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40 & Corbula.Panzano.ages > 43],Corbula.Panzano.5depths[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40 & Corbula.Panzano.ages > 43], method="s")
cor1conf=cor.ci(cbind(Corbula.Po4.ages[Corbula.Po4.5depths<20],Corbula.Po4.5depths[Corbula.Po4.5depths<20]), method="s")$ci[c(1,4)]
cor2conf=cor.ci(cbind(Corbula.Po4.ages[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90 & Corbula.Po4.ages > 13],Corbula.Po4.5depths[Corbula.Po4.5depths>20 & Corbula.Po4.5depths<90 & Corbula.Po4.ages > 13]), method="s")$ci[c(1,4)]
cor3conf=cor.ci(cbind(Corbula.Po4.ages[Corbula.Po4.5depths>90 & Corbula.Po4.ages > 43],Corbula.Po4.5depths[Corbula.Po4.5depths>90 & Corbula.Po4.ages > 43]), method="s")$ci[c(1,4)]
cor4conf=cor.ci(cbind(Corbula.Po3.ages[Corbula.Po3.5depths<20],Corbula.Po3.5depths[Corbula.Po3.5depths<20]), method="s")$ci[c(1,4)]
cor5conf=cor.ci(cbind(Corbula.Po3.ages[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90 & Corbula.Po3.ages > 13],Corbula.Po3.5depths[Corbula.Po3.5depths>20 & Corbula.Po3.5depths<90 & Corbula.Po3.ages > 13]), method="s")$ci[c(1,4)]
cor6conf=cor.ci(cbind(Corbula.Po3.ages[Corbula.Po3.5depths>90 & Corbula.Po3.ages > 43],Corbula.Po3.5depths[Corbula.Po3.5depths>90 & Corbula.Po3.ages > 43]), method="s")$ci[c(1,4)]
cor7conf=cor.ci(cbind(Corbula.Panzano.ages[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14],Corbula.Panzano.5depths[Corbula.Panzano.5depths>0 & Corbula.Panzano.5depths<14]), method="s")$ci[c(1,4)]
cor8conf=cor.ci(cbind(Corbula.Panzano.ages[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40 & Corbula.Panzano.ages > 43],Corbula.Panzano.5depths[Corbula.Panzano.5depths>14 & Corbula.Panzano.5depths<40 & Corbula.Panzano.ages > 43]), method="s")$ci[c(1,4)]

#######################################################################################################
#FIGURE DISORDER
#######################################################################################################
par(mfrow=c(2,3))
#par(mar=c(4,2,2,0))
plot(c(1:8), c(cor1$estimate, cor2$estimate, cor3$estimate,cor4$estimate, cor5$estimate, cor6$estimate, cor7$estimate, cor8$estimate), 
pch=16, ylim=c(-0.5,1), ylab="Strat. disorder (age vs. depth correlation)", xlab="", xaxt="n", xlim=c(0.5,8.5), cex=1.4, main="", cex.main=0.95)
segments(x0=c(1:8), x1=c(1:8), y0=unlist(c(cor1conf[1],cor2conf[1],cor3conf[1],cor4conf[1],cor5conf[1],cor6conf[1],cor7conf[1],cor8conf[1])),
y1=unlist(c(cor1conf[2],cor2conf[2],cor3conf[2],cor4conf[2],cor5conf[2],cor6conf[2],cor7conf[2],cor8conf[2])), lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
abline(v=c(3.5,6.5))
abline(h=0, lty=3)
text(2,0.95, labels="Po 4",srt=0)
text(5,0.95, labels="Po 3",srt=0)
text(7.5,0.9, labels="Panzano",srt=90)

top=c(cor.range952[depth2<20],cor.range953[depth3<20])
mid=c(cor.range952[depth2>20 & depth2<90],cor.range953[depth3>20 & depth3<90])
bottom=c(cor.range952[depth2>90],cor.range953[depth3>90])

######################################################################################################################
#TIME AVERAGING FOR THREE UNITS AT THE SCALE OF INCREMENTS
######################################################################################################################
#######################################################################################################
#50% AGE RANGE (IQR)
#######################################################################################################
#######################################################################################################
#FIGURE WITHOUT AGE TRUNCATION, WITH ERROR CORRECTION
#######################################################################################################
out.iqr.err=boxplot(cor.IQR3[depth3<20],cor.IQR3[depth3>20 & depth3<90],cor.IQR3[depth3>90],
cor.IQR2[depth2<20],cor.IQR2[depth2>20 & depth2<90],cor.IQR2[depth2>90],
cor.IQR1[depth1>0 & depth1<14],cor.IQR1[depth1>14 & depth1<75], ylab="Time averaging (y)-IQR", ylim=c(0,200), col="gray", range=0, main="", cex.main=0.95)
abline(v=c(3.5,6.5),lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
text(2,190, labels="Po 4",srt=0)
text(4.5,190, labels="Po 3",srt=0)
text(7.5,185, labels="Panzano",srt=90)

#######################################################################################################
#FIGURE WITH AGE TRUNCATION
#######################################################################################################
out.iqr.tru=boxplot(pre.IQR3[pre.depth3<20],pre.IQR3[pre.depth3>20 & pre.depth3<90],pre.IQR3[pre.depth3>90],
pre.IQR2[pre.depth2<20],pre.IQR2[pre.depth2>20 & pre.depth2<90],pre.IQR2[pre.depth2>90],
pre.IQR1[pre.depth1>0 & pre.depth1<14],pre.IQR1[pre.depth1>14 & pre.depth1<75], ylab="Time averaging (y)-IQR", ylim=c(0,200), col="gray", main="", cex.main=0.95)
abline(v=c(3.5,6.5),lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
text(2,190, labels="Po 4",srt=0)
text(4.5,190, labels="Po 3",srt=0)
text(7.5,185, labels="Panzano",srt=90)

#######################################################################################################
#95% AGE RANGE
#######################################################################################################
#FIGURE TIME AVERAGING PER INCREMENT, WITHOUT AGE TRUNCATION, WITH ERROR CORRECTION
#######################################################################################################
out.95.err=boxplot(cor.range953[depth3<20],cor.range953[depth3>20 & depth3<90],cor.range953[depth3>90], 
cor.range952[depth2<20],cor.range952[depth2>20 & depth2<90],cor.range952[depth2>90],
cor.range951[depth1>0 & depth1<14],cor.range951[depth1>14 & depth1<75],  
ylab="Time averaging (y)-95% range", ylim=c(0,300), col="gray", range=0, main="95% range", cex.main=0.95)
abline(v=c(3.5,6.5),lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
text(2,295, labels="Po 4",srt=0)
text(4.5,295, labels="Po 3",srt=0)
text(7.5,280, labels="Panzano",srt=90)
par(srt=0)

#######################################################################################################
#FIGURE TIME AVERAGING PER INCREMENT
#######################################################################################################
out.95.tru=boxplot(pre.range953[pre.depth3<20],pre.range953[pre.depth3>20 & pre.depth3<90],pre.range953[pre.depth3>90], 
pre.range952[pre.depth2<20],pre.range952[pre.depth2>20 & pre.depth2<90],pre.range952[pre.depth2>90],
pre.range951[pre.depth1>0 & depth1<14],pre.range951[pre.depth1>14 & pre.depth1<75],  
ylab="Time averaging (y)-95% range", ylim=c(0,300), col="gray", main="95% range", cex.main=0.95, cex.main=0.95)
abline(v=c(3.5,6.5),lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
text(2,295, labels="Po 4",srt=0)
text(4.5,295, labels="Po 3",srt=0)
text(7.5,280, labels="Panzano",srt=90)
par(srt=0)

out.iqr.err$stats[3,]
out.iqr.tru$stats[3,]
out.95.err$stats[3,]
out.95.tru$stats[3,]


#######################################################################################################
#FIGURE TIME AVERAGING PER INCREMENT, WITH AGE TRUNCATION, WITH ERROR CORRECTION
#######################################################################################################
boxplot(pre.cor.IQR3[pre.depth3<20],pre.cor.IQR3[pre.depth3>20 & pre.depth3<90],pre.cor.IQR3[pre.depth3>90],
pre.cor.IQR2[pre.depth2<20],pre.cor.IQR2[pre.depth2>20 & pre.depth2<90],pre.cor.IQR2[pre.depth2>90],
pre.cor.IQR1[pre.depth1>0 & pre.depth1<14],pre.cor.IQR1[pre.depth1>14 & pre.depth1<75], ylab="Time averaging (y)-IQR", ylim=c(0,200), col="gray")
abline(v=c(3.5,6.5),lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
text(2,190, labels="Po 4",srt=0)
text(4.5,190, labels="Po 3",srt=0)
text(7.5,185, labels="Panzano",srt=90)

boxplot(pre.cor.range953[pre.depth3<20],pre.cor.range953[pre.depth3>20 & pre.depth3<90],pre.cor.range953[pre.depth3>90], 
pre.cor.range952[pre.depth2<20],pre.cor.range952[pre.depth2>20 & pre.depth2<90],pre.cor.range952[pre.depth2>90],
pre.cor.range951[pre.depth1>0 & depth1<14],pre.cor.range951[pre.depth1>14 & pre.depth1<75],  
ylab="Time averaging (y)-95% range", ylim=c(0,300), col="gray")
abline(v=c(3.5,6.5),lty=2)
axis(1, at=c(1:8), labels=c("Post-2000","Late 20th c.","Early 20th c.","Post-2000","Late 20th c.","Early 20th c.","Late 20th c.","19th-Early 20th c."), las=2)
text(2,295, labels="Po 4",srt=0)
text(4.5,295, labels="Po 3",srt=0)
text(7.5,280, labels="Panzano",srt=90)
par(srt=0)


#######################################################################################################
#RECONSTRUCTION OF PRODUCTION
#######################################################################################################
ERROR=FALSE
######################################################################################################
#PO3 UNMIXING
######################################################################################################
SIMULATIONS=100
N.full.totalsizePo3=full.totalsizePo3
Full.Po3bootout.2=array(NA, dim=c(SIMULATIONS,150))
Full.Po3bootout.5=array(NA, dim=c(SIMULATIONS,59))
Full.Po3bootout.10=array(NA, dim=c(SIMULATIONS,30))
Full.Po3bootout.20=array(NA, dim=c(SIMULATIONS,15))

for (j in 1:SIMULATIONS) {
Po3bootages=numeric()
Po3.list1=split(Corbula.Po3.ages,Corbula.Po3.5depths)
Po3.list1$"2"=Po3.list1$"2"
Po3.list1$"6"=Po3.list1$"6"
Po3.list1$"10"=Po3.list1$"10"
Po3.list1$"14"=Po3.list1$"14"
Po3.list1$"18"=Po3.list1$"18"
#20-30 cm
Po3.list1$"22.5"=Po3.list1$"22.5"
Po3.list1$"27.5"=rtnorm(100,mean=mean(c(Po3.list1$"22.5")), sd=sd(c(Po3.list1$"22.5")), lower=0)
#30-50 cm
Po3.list1$"37.5"=rtnorm(100,mean=mean(c(Po3.list1$"42.5")), sd=sd(c(Po3.list1$"42.5")), lower=0)
Po3.list1$"42.5"=Po3.list1$"42.5"
Po3.list1$"47.5"=rtnorm(100,mean=mean(c(Po3.list1$"42.5")), sd=sd(c(Po3.list1$"42.5")), lower=0)
Po3.list1$"32.5"=rtnorm(100,mean=mean(c(Po3.list1$"42.5")), sd=sd(c(Po3.list1$"47.5")), lower=0)
#50-85 cm
Po3.list1$"67.5"=Po3.list1$"67.5"
Po3.list1$"52.5"=rtnorm(100,mean=mean(c(Po3.list1$"67.5")), sd=sd(c(Po3.list1$"67.5")), lower=0)
Po3.list1$"57.5"=rtnorm(100,mean=mean(c(Po3.list1$"67.5")), sd=sd(c(Po3.list1$"67.5")), lower=0)
Po3.list1$"62.5"=rtnorm(100,mean=mean(c(Po3.list1$"67.5")), sd=sd(c(Po3.list1$"67.5")), lower=0)
Po3.list1$"72.5"=rtnorm(100,mean=mean(c(Po3.list1$"67.5",Po3.list1$"82.5")), sd=sd(c(Po3.list1$"67.5",Po3.list1$"82.5")), lower=0)
Po3.list1$"77.5"=rtnorm(100,mean=mean(c(Po3.list1$"67.5",Po3.list1$"82.5")), sd=sd(c(Po3.list1$"67.5",Po3.list1$"82.5")), lower=0)
Po3.list1$"82.5"=Po3.list1$"82.5"
#>85 cm
Po3.list1$"87.5"=rtnorm(100,mean=mean(c(Po3.list1$"97.5")), sd=sd(c(Po3.list1$"97.5")), lower=0)
Po3.list1$"92.5"=rtnorm(100,mean=mean(c(Po3.list1$"97.5",Po3.list1$"82.5")), sd=sd(c(Po3.list1$"97.5",Po3.list1$"82.5")), lower=0)
Po3.list1$"97.5"=Po3.list1$"97.5"
Po3.list1$"102.5"=rtnorm(100,mean=mean(c(Po3.list1$"107.5",Po3.list1$"97.5")), sd=sd(c(Po3.list1$"107.5",Po3.list1$"97.5")), lower=0)
Po3.list1$"107.5"=Po3.list1$"107.5"
Po3.list1$"112.5"=Po3.list1$"112.5"
Po3.list1$"117.5"=rtnorm(100,mean=mean(c(Po3.list1$"132.5",Po3.list1$"112.5")), sd=sd(c(Po3.list1$"132.5",Po3.list1$"112.5")), lower=0)
Po3.list1$"122.5"=rtnorm(100,mean=mean(c(Po3.list1$"132.5",Po3.list1$"112.5")), sd=sd(c(Po3.list1$"132.5",Po3.list1$"112.5")), lower=0)
Po3.list1$"127.5"=rtnorm(100,mean=mean(c(Po3.list1$"132.5",Po3.list1$"112.5")), sd=sd(c(Po3.list1$"132.5",Po3.list1$"112.5")), lower=0)
Po3.list1$"132.5"=Po3.list1$"132.5"
Po3.list1$"137.5"=Po3.list1$"137.5"


temp=names(Po3.list1)
temp=as.numeric(temp)
ord=order(temp)
Po3.list1=Po3.list1[ord]
for (i in 1:length(Po3.list1)) {
temp=unlist(Po3.list1[i])
sample.age=sample(temp,full.totalsizePo3[i],replace=T)
if (ERROR==TRUE) bootages=rlnorm(full.totalsizePo3[i],log(sample.age),sdlog=sqrt(exp(ln.d)))
if (ERROR==FALSE) bootages=sample(temp,full.totalsizePo3[i],replace=T)
Po3bootages=append(Po3bootages,bootages,length(Po3bootages))
}
out1=hist(Po3bootages[Po3bootages<298], breaks=c(0,8,seq(13,300,by=5))) 
breaks.Po3.10=2013-c(0,seq(8,303,by=10))
out2=hist(Po3bootages[Po3bootages<298], breaks=2013-breaks.Po3.10) 
out3=hist(Po3bootages[Po3bootages<293], breaks=c(0,seq(13,300,by=20))) 
out4=hist(Po3bootages[Po3bootages<298], breaks=c(0,2,4,6,8,10,12,seq(13,300,by=2))) 
Full.Po3bootout.2[j,]=out4$counts
Full.Po3bootout.5[j,]=out1$counts
Full.Po3bootout.10[j,]=out2$counts
Full.Po3bootout.20[j,]=out3$counts
}

######################################################################################################
#PO4 UNMIXING
######################################################################################################
N.full.totalsizePo4=full.totalsizePo4
Full.Po4bootout.2=array(NA, dim=c(SIMULATIONS,150))
Full.Po4bootout.5=array(NA, dim=c(SIMULATIONS,59))
Full.Po4bootout.10=array(NA, dim=c(SIMULATIONS,30))
Full.Po4bootout.20=array(NA, dim=c(SIMULATIONS,15))

for (j in 1:SIMULATIONS) {
Po4bootages=numeric()
Po4.list1=split(Corbula.Po4.ages,Corbula.Po4.5depths)
Po4.list1$"2"=Po4.list1$"2"
Po4.list1$"6"=Po4.list1$"6"
Po4.list1$"10"=Po4.list1$"10"
Po4.list1$"14"=Po4.list1$"14"
Po4.list1$"18"=Po4.list1$"18"
Po4.list1$"22.5"=rtnorm(100,mean=mean(c(Po4.list1$"32.5")), sd=sd(c(Po4.list1$"32.5")), lower=0)
Po4.list1$"27.5"=rtnorm(100,mean=mean(c(Po4.list1$"32.5")), sd=sd(c(Po4.list1$"32.5")), lower=0)
Po4.list1$"32.5"=Po4.list1$"32.5"
Po4.list1$"37.5"=rtnorm(100,mean=mean(c(Po4.list1$"32.5",Po4.list1$"47.5")), sd=sd(c(Po4.list1$"32.5",Po4.list1$"47.5")), lower=0)
Po4.list1$"42.5"=rtnorm(100,mean=mean(c(Po4.list1$"32.5",Po4.list1$"47.5")), sd=sd(c(Po4.list1$"32.5",Po4.list1$"47.5")), lower=0)
Po4.list1$"47.5"=Po4.list1$"47.5"
Po4.list1$"52.5"=rtnorm(100,mean=mean(c(Po4.list1$"62.5")), sd=sd(c(Po4.list1$"62.5")), lower=0)
Po4.list1$"57.5"=rtnorm(100,mean=mean(c(Po4.list1$"62.5")), sd=sd(c(Po4.list1$"62.5")), lower=0)
Po4.list1$"62.5"=Po4.list1$"62.5"
Po4.list1$"67.5"=rtnorm(100,mean=mean(c(Po4.list1$"77.5",Po4.list1$"62.5")), sd=sd(c(Po4.list1$"77.5",Po4.list1$"62.5")), lower=0)
Po4.list1$"72.5"=rtnorm(100,mean=mean(c(Po4.list1$"77.5",Po4.list1$"62.5")), sd=sd(c(Po4.list1$"77.5",Po4.list1$"62.5")), lower=0)
Po4.list1$"77.5"=Po4.list1$"77.5"
Po4.list1$"82.5"=rtnorm(100,mean=mean(c(Po4.list1$"77.5")), sd=sd(c(Po4.list1$"77.5")), lower=0)
Po4.list1$"87.5"=rtnorm(100,mean=mean(c(Po4.list1$"92.5")), sd=sd(c(Po4.list1$"92.5")), lower=0)
Po4.list1$"92.5"=Po4.list1$"92.5"
Po4.list1$"97.5"=rtnorm(100,mean=mean(c(Po4.list1$"107.5",Po4.list1$"92.5")), sd=sd(c(Po4.list1$"107.5",Po4.list1$"92.5")), lower=0)
Po4.list1$"102.5"=rtnorm(100,mean=mean(c(Po4.list1$"107.5",Po4.list1$"92.5")), sd=sd(c(Po4.list1$"107.5",Po4.list1$"92.5")), lower=0)
Po4.list1$"107.5"=Po4.list1$"107.5"
Po4.list1$"112.5"=rtnorm(100,mean=mean(c(Po4.list1$"107.5",Po4.list1$"127.5")), sd=sd(c(Po4.list1$"107.5",Po4.list1$"127.5")), lower=0)
Po4.list1$"117.5"=rtnorm(100,mean=mean(c(Po4.list1$"107.5",Po4.list1$"127.5")), sd=sd(c(Po4.list1$"107.5",Po4.list1$"127.5")), lower=0)
Po4.list1$"122.5"=rtnorm(100,mean=mean(c(Po4.list1$"107.5",Po4.list1$"127.5")), sd=sd(c(Po4.list1$"107.5",Po4.list1$"127.5")), lower=0)
Po4.list1$"127.5"=Po4.list1$"127.5"
Po4.list1$"132.5"=Po4.list1$"132.5"
Po4.list1$"137.5"=rtnorm(100,mean=mean(c(Po4.list1$"152.5",Po4.list1$"132.5")), sd=sd(c(Po4.list1$"152.5",Po4.list1$"132.5")), lower=0)
Po4.list1$"142.5"=rtnorm(100,mean=mean(c(Po4.list1$"152.5",Po4.list1$"132.5")), sd=sd(c(Po4.list1$"152.5",Po4.list1$"132.5")), lower=0)
Po4.list1$"147.5"=rtnorm(100,mean=mean(c(Po4.list1$"152.5",Po4.list1$"132.5")), sd=sd(c(Po4.list1$"152.5",Po4.list1$"132.5")), lower=0)
Po4.list1$"152.5"=Po4.list1$"152.5"

temp=names(Po4.list1)
temp=as.numeric(temp)
ord=order(temp)
Po4.list1=Po4.list1[ord]
for (i in 1:length(Po4.list1)) {
temp=unlist(Po4.list1[i])
sample.age=sample(temp,full.totalsizePo4[i],replace=T)
if (ERROR==TRUE) bootages=rlnorm(full.totalsizePo4[i],log(sample.age),sdlog=sqrt(exp(ln.d)))
if (ERROR==FALSE) bootages=sample(temp,full.totalsizePo4[i],replace=T)
Po4bootages=append(Po4bootages,bootages,length(Po4bootages))
}

if (max(Po4bootages)>290) {Po4bootages=Po4bootages[-which(Po4bootages>290)]}
out1=hist(Po4bootages, breaks=c(0,8,seq(13,300,by=5))) 
breaks.Po4.10=2013-c(0,seq(8,303,by=10))
out2=hist(Po4bootages, breaks=2013-breaks.Po4.10) 
out3=hist(Po4bootages, breaks=c(0,seq(13,300,by=20))) 
out4=hist(Po4bootages, breaks=c(0,2,4,6,8,10,12,seq(13,300,by=2)))
Full.Po4bootout.2[j,]=out4$counts
Full.Po4bootout.5[j,]=out1$counts
Full.Po4bootout.10[j,]=out2$counts
Full.Po4bootout.20[j,]=out3$counts
}

######################################################################################################
#PANZANO UNMIXING
######################################################################################################

N.full.totalsizePanzano=full.totalsizePanzano
Full.Panzanobootout.2=array(NA, dim=c(SIMULATIONS,1000))
Full.Panzanobootout.5=array(NA, dim=c(SIMULATIONS,399))
Full.Panzanobootout.10=array(NA, dim=c(SIMULATIONS,200))
Full.Panzanobootout.20=array(NA, dim=c(SIMULATIONS,100))
bootPanzano.lambda=numeric()
for (j in 1:SIMULATIONS) {
Panzanobootages=numeric()
Panzano.list1=split(Corbula.Panzano.ages,Corbula.Panzano.5depths)
Panzano.list1$"4"=Panzano.list1$"4"
Panzano.list1$"8"=Panzano.list1$"8"
Panzano.list1$"12"=Panzano.list1$"12"
Panzano.list1$"16"=rtnorm(100,mean=mean(c(Panzano.list1$"12",Panzano.list1$"20")), sd=sd(c(Panzano.list1$"12",Panzano.list1$"20")), lower=0)
Panzano.list1$"20"=Panzano.list1$"20"
Panzano.list1$"25"=rtnorm(100,mean=mean(c(Panzano.list1$"20",Panzano.list1$"35")), sd=sd(c(Panzano.list1$"20",Panzano.list1$"35")), lower=0)
Panzano.list1$"30"=rtnorm(100,mean=mean(c(Panzano.list1$"20",Panzano.list1$"35")), sd=sd(c(Panzano.list1$"20",Panzano.list1$"35")), lower=0)
Panzano.list1$"35"=Panzano.list1$"35"
Panzano.list1$"40"=rtnorm(100,mean=mean(c(Panzano.list1$"35",Panzano.list1$"50")), sd=sd(c(Panzano.list1$"35",Panzano.list1$"50")), lower=0)
Panzano.list1$"45"=rtnorm(100,mean=mean(c(Panzano.list1$"35",Panzano.list1$"50")), sd=sd(c(Panzano.list1$"35",Panzano.list1$"50")), lower=0)
Panzano.list1$"50"=Panzano.list1$"50"
Panzano.list1$"55"=rtnorm(100,mean=mean(c(Panzano.list1$"65",Panzano.list1$"50")), sd=sd(c(Panzano.list1$"65",Panzano.list1$"50")), lower=0)
Panzano.list1$"60"=rtnorm(100,mean=mean(c(Panzano.list1$"65",Panzano.list1$"50")), sd=sd(c(Panzano.list1$"65",Panzano.list1$"50")), lower=0)
Panzano.list1$"65"=Panzano.list1$"65"
Panzano.list1$"70"=Panzano.list1$"70"
Panzano.list1$"75"=rtnorm(100,mean=mean(c(Panzano.list1$"70",Panzano.list1$"90")), sd=sd(c(Panzano.list1$"70",Panzano.list1$"90")), lower=0)
Panzano.list1$"80"=rtnorm(100,mean=mean(c(Panzano.list1$"70",Panzano.list1$"90")), sd=sd(c(Panzano.list1$"70",Panzano.list1$"90")), lower=0)
Panzano.list1$"85"=rtnorm(100,mean=mean(c(Panzano.list1$"70",Panzano.list1$"90")), sd=sd(c(Panzano.list1$"70",Panzano.list1$"90")), lower=0)
Panzano.list1$"90"=Panzano.list1$"90"
Panzano.list1$"95"=Panzano.list1$"95"
Panzano.list1$"100"=rtnorm(100,mean=mean(c(Panzano.list1$"115",Panzano.list1$"95")), sd=sd(c(Panzano.list1$"115",Panzano.list1$"95")), lower=0)
Panzano.list1$"105"=rtnorm(100,mean=mean(c(Panzano.list1$"115",Panzano.list1$"95")), sd=sd(c(Panzano.list1$"115",Panzano.list1$"95")), lower=0)
Panzano.list1$"110"=rtnorm(100,mean=mean(c(Panzano.list1$"115",Panzano.list1$"95")), sd=sd(c(Panzano.list1$"115",Panzano.list1$"95")), lower=0)
Panzano.list1$"115"=Panzano.list1$"115"
Panzano.list1$"120"=rtnorm(100,mean=mean(c(Panzano.list1$"115",Panzano.list1$"130")), sd=sd(c(Panzano.list1$"115",Panzano.list1$"130")), lower=0)
Panzano.list1$"125"=rtnorm(100,mean=mean(c(Panzano.list1$"115",Panzano.list1$"130")), sd=sd(c(Panzano.list1$"115",Panzano.list1$"130")), lower=0)
Panzano.list1$"130"=Panzano.list1$"130"
Panzano.list1$"135"=rtnorm(100,mean=mean(c(Panzano.list1$"150",Panzano.list1$"130")), sd=sd(c(Panzano.list1$"150",Panzano.list1$"130")), lower=0)
Panzano.list1$"140"=rtnorm(100,mean=mean(c(Panzano.list1$"150",Panzano.list1$"130")), sd=sd(c(Panzano.list1$"150",Panzano.list1$"130")), lower=0)
Panzano.list1$"145"=rtnorm(100,mean=mean(c(Panzano.list1$"150",Panzano.list1$"130")), sd=sd(c(Panzano.list1$"150",Panzano.list1$"130")), lower=0)
Panzano.list1$"150"=Panzano.list1$"150"
temp=names(Panzano.list1)
temp=as.numeric(temp)
ord=order(temp)
Panzano.list1=Panzano.list1[ord]
for (i in 1:length(Panzano.list1)) {
temp=unlist(Panzano.list1[i])
sample.age=sample(temp,full.totalsizePanzano[i],replace=T)
if (ERROR==TRUE) bootages=rlnorm(full.totalsizePanzano[i],log(sample.age),sdlog=sqrt(exp(ln.d)))
if (ERROR==FALSE) bootages=sample(temp,full.totalsizePanzano[i],replace=T)
Panzanobootages=append(Panzanobootages,bootages,length(Panzanobootages))
}

bootPanzano.lambda[j]=1/mean(Panzanobootages)
if (max(Panzanobootages)>1990) {Panzanobootages=Panzanobootages[-which(Panzanobootages>1990)]}
breaks.Panzano.10=2013-c(0,seq(8,2003,by=10))
out1=hist(Panzanobootages, breaks=2013-breaks.Panzano.10) 
out2=hist(Panzanobootages, breaks=c(0,seq(13,2000,by=20))) 
out3=hist(Panzanobootages, breaks=c(0,8,seq(13,2000,by=5))) 
out4=hist(Panzanobootages, breaks=c(0,2,4,6,8,10,12,seq(13,2000,by=2))) 
Full.Panzanobootout.2[j,]=out4$counts
Full.Panzanobootout.5[j,]=out3$counts
Full.Panzanobootout.10[j,]=out1$counts
Full.Panzanobootout.20[j,]=out2$counts
}



###########################################################################################################
full.totalsizePiran.Corbula=full.totalsizePiran
Full.Piranbootout.100=array(NA, dim=c(SIMULATIONS,140))
Full.Piranbootout.200=array(NA, dim=c(SIMULATIONS,70))
Full.Piranbootout.250=array(NA, dim=c(SIMULATIONS,56))

for (j in 1:SIMULATIONS) {
Piranbootages=numeric()
Piran.list1=split(Corbula.Piran.ages,Corbula.Piran.5depths)
Piran.list1$"4"=Piran.list1$"4"
Piran.list1$"8"=rtnorm(100,mean=mean(c(Piran.list1$"4",Piran.list1$"12")), sd=sd(c(Piran.list1$"4",Piran.list1$"12")), lower=0)
Piran.list1$"12"=Piran.list1$"12"
Piran.list1$"16"=rtnorm(100,mean=mean(c(Piran.list1$"12",Piran.list1$"30")), sd=sd(c(Piran.list1$"12",Piran.list1$"30")), lower=0)
Piran.list1$"20"=rtnorm(100,mean=mean(c(Piran.list1$"12",Piran.list1$"30")), sd=sd(c(Piran.list1$"12",Piran.list1$"30")), lower=0)
Piran.list1$"25"=rtnorm(100,mean=mean(c(Piran.list1$"12",Piran.list1$"30")), sd=sd(c(Piran.list1$"12",Piran.list1$"30")), lower=0)
Piran.list1$"30"=Piran.list1$"30"
Piran.list1$"35"=rtnorm(100,mean=mean(c(Piran.list1$"30",Piran.list1$"50")), sd=sd(c(Piran.list1$"30",Piran.list1$"50")), lower=0)
Piran.list1$"40"=rtnorm(100,mean=mean(c(Piran.list1$"30",Piran.list1$"50")), sd=sd(c(Piran.list1$"30",Piran.list1$"50")), lower=0)
Piran.list1$"45"=rtnorm(100,mean=mean(c(Piran.list1$"30",Piran.list1$"50")), sd=sd(c(Piran.list1$"30",Piran.list1$"50")), lower=0)
Piran.list1$"50"=Piran.list1$"50"
Piran.list1$"55"=rtnorm(100,mean=mean(c(Piran.list1$"50",Piran.list1$"70")), sd=sd(c(Piran.list1$"50",Piran.list1$"70")), lower=0)
Piran.list1$"60"=rtnorm(100,mean=mean(c(Piran.list1$"50",Piran.list1$"70")), sd=sd(c(Piran.list1$"50",Piran.list1$"70")), lower=0)
Piran.list1$"65"=rtnorm(100,mean=mean(c(Piran.list1$"50",Piran.list1$"70")), sd=sd(c(Piran.list1$"50",Piran.list1$"70")), lower=0)
Piran.list1$"70"=Piran.list1$"70"
Piran.list1$"75"=rtnorm(100,mean=mean(c(Piran.list1$"70",Piran.list1$"90")), sd=sd(c(Piran.list1$"70",Piran.list1$"90")), lower=0)
Piran.list1$"80"=rtnorm(100,mean=mean(c(Piran.list1$"70",Piran.list1$"90")), sd=sd(c(Piran.list1$"70",Piran.list1$"90")), lower=0)
Piran.list1$"85"=rtnorm(100,mean=mean(c(Piran.list1$"70",Piran.list1$"90")), sd=sd(c(Piran.list1$"70",Piran.list1$"90")), lower=0)
Piran.list1$"90"=Piran.list1$"90"
Piran.list1$"95"=rtnorm(100,mean=mean(c(Piran.list1$"90",Piran.list1$"105")), sd=sd(c(Piran.list1$"90",Piran.list1$"105")), lower=0)
Piran.list1$"100"=rtnorm(100,mean=mean(c(Piran.list1$"90",Piran.list1$"105")), sd=sd(c(Piran.list1$"90",Piran.list1$"105")), lower=0)
Piran.list1$"105"=Piran.list1$"105"
Piran.list1$"110"=rtnorm(100,mean=mean(c(Piran.list1$"105",Piran.list1$"130")), sd=sd(c(Piran.list1$"105",Piran.list1$"130")), lower=0)
Piran.list1$"115"=rtnorm(100,mean=mean(c(Piran.list1$"105",Piran.list1$"130")), sd=sd(c(Piran.list1$"105",Piran.list1$"130")), lower=0)
Piran.list1$"120"=rtnorm(100,mean=mean(c(Piran.list1$"105",Piran.list1$"130")), sd=sd(c(Piran.list1$"105",Piran.list1$"130")), lower=0)
Piran.list1$"125"=rtnorm(100,mean=mean(c(Piran.list1$"105",Piran.list1$"130")), sd=sd(c(Piran.list1$"105",Piran.list1$"130")), lower=0)
Piran.list1$"130"=Piran.list1$"130"
Piran.list1$"135"=rtnorm(100,mean=mean(c(Piran.list1$"130",Piran.list1$"150")), sd=sd(c(Piran.list1$"130",Piran.list1$"150")), lower=0)
Piran.list1$"140"=rtnorm(100,mean=mean(c(Piran.list1$"130",Piran.list1$"150")), sd=sd(c(Piran.list1$"130",Piran.list1$"150")), lower=0)
Piran.list1$"145"=rtnorm(100,mean=mean(c(Piran.list1$"130",Piran.list1$"150")), sd=sd(c(Piran.list1$"130",Piran.list1$"150")), lower=0)
Piran.list1$"150"=Piran.list1$"150"
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])
#incorporating age error
sample.age=sample(temp,full.totalsizePiran[i],replace=T)
if (ERROR==TRUE) bootages=rlnorm(full.totalsizePiran[i],log(sample.age),sdlog=sqrt(exp(ln.d)))
if (ERROR==FALSE) bootages=sample(temp,full.totalsizePiran[i],replace=T)
Piranbootages=append(Piranbootages,bootages,length(Piranbootages))
}

out1=hist(Piranbootages[Piranbootages<13900], breaks=c(0,seq(13,14000,by=100))) 
out2=hist(Piranbootages[Piranbootages<13800], breaks=c(0,seq(13,14000,by=200))) 
out3=hist(Piranbootages[Piranbootages<13750], breaks=c(0,seq(13,14000,by=250))) 
Full.Piranbootout.100[j,]=out1$counts
Full.Piranbootout.200[j,]=out2$counts
Full.Piranbootout.250[j,]=out3$counts
}


#######################################################################################################
#SURFACE AND CORE AGE FREQUENCY DISTRIBUTION
#######################################################################################################
#SURFACE PANZANO#
#################
par(mfrow=c(2,3))
age=surfaceAFD.Panzano.6
hist(surfaceAFD.Panzano.6, col="gray",breaks=seq(0,500,by=10), xlab="Postmortem age (years)", ylab="Number of specimens",
ylim=c(0,25), main="Surface mixed layer", cex.main=0.9, bty="u")
lambda<-1/mean(age); 
times=seq(0,200, by=10)+5
predictSIMPLE<-dexp(times,rate=lambda)
lines(times, (predictSIMPLE/sum(predictSIMPLE))*length(age), lty=1, col="black", lwd=3)

#RAW CORE NOT INTERPOLATED, NO SAMPLE SIZE CORRECTION
hist(Corbula.Panzano.ages[Corbula.Panzano.ages<500 & Corbula.Panzano.5depths>4], breaks=seq(0,500,by=10), cex.main=0.9, col="gray",
ylim=c(0,25), main="Historical layers", xlab="Postmortem age (years)", ylab="Number of specimens")

hist(Corbula.Panzano.ages[Corbula.Panzano.ages<500], breaks=seq(0,500,by=10), cex.main=0.9, col="gray",
ylim=c(0,25), main="Empirical pooling", xlab="Postmortem age (years)", ylab="Number of specimens")
median(Corbula.Panzano.ages)


#######################################################################################################
#COMPUTE NUMBER OF SHELLS PER 5 OR 10-YEAR COHORT
#######################################################################################################
Po3.means.5=apply(Full.Po3bootout.5,MARGIN=2,mean)
Po3.025.5=apply(Full.Po3bootout.5,MARGIN=2,quantile, 0.025)
Po3.975.5=apply(Full.Po3bootout.5,MARGIN=2,quantile, 0.975)
names(Po3.means.5)=c(0,8,seq(13,295,by=5))
time.Po3.5=2013-c(0,8,seq(13,295,by=5))

Po4.means.5=apply(Full.Po4bootout.5,MARGIN=2,mean)
Po4.025.5=apply(Full.Po4bootout.5,MARGIN=2,quantile, 0.025)
Po4.975.5=apply(Full.Po4bootout.5,MARGIN=2,quantile, 0.975)
time.Po4.5=2013-c(0,8,seq(13,295,by=5))

Panzano.means.5=apply(Full.Panzanobootout.5,MARGIN=2,mean)
Panzano.025.5=apply(Full.Panzanobootout.5,MARGIN=2,quantile, 0.025)
Panzano.975.5=apply(Full.Panzanobootout.5,MARGIN=2,quantile, 0.975)
time.Panzano.5=2013-c(0,8,seq(13,1995,by=5))


time.Po3.10=breaks.Po3.10[-length(breaks.Po3.10)]-5
Po3.means.10=apply(Full.Po3bootout.10,MARGIN=2,mean)
Po3.025.10=apply(Full.Po3bootout.10,MARGIN=2,quantile, 0.025)
Po3.975.10=apply(Full.Po3bootout.10,MARGIN=2,quantile, 0.975)

time.Po4.10=breaks.Po4.10[-length(breaks.Po4.10)]-5
Po4.means.10=apply(Full.Po4bootout.10,MARGIN=2,mean)
Po4.025.10=apply(Full.Po4bootout.10,MARGIN=2,quantile, 0.025)
Po4.975.10=apply(Full.Po4bootout.10,MARGIN=2,quantile, 0.975)

time.Panzano.10=breaks.Panzano.10[-length(breaks.Panzano.10)]-5
Panzano.means.10=apply(Full.Panzanobootout.10,MARGIN=2,mean)
Panzano.025.10=apply(Full.Panzanobootout.10,MARGIN=2,quantile, 0.025)
Panzano.975.10=apply(Full.Panzanobootout.10,MARGIN=2,quantile, 0.975)


#######################################################################################################
#20th CENTURY RECONSTRUCTION AND PO AND PANZANO
#######################################################################################################
par(mfrow=c(2,2))
plot(Po3.means.10[1:12],time.Po3.10[1:12], type="l", ylab="Calendar age", xlab="Corbula gibba abundance/10 years", 
xlim=c(0,300), main="Po", cex.main=0.9)
segments(x0=Po3.025.10[1:12],x1=Po3.975.10[1:12],y0=time.Po3.10[1:12], y1=time.Po3.10[1:12],lty=1)
lines(Po3.means.10[1:12],time.Po3.10[1:12],lwd=2)
#plot(Po4.means.10[1:12],time.Po4.10[1:12], type="l", ylab="Calendar age", xlab="Corbula gibba abundance/10 years")
segments(x0=Po4.025.10[1:12],x1=Po4.975.10[1:12],y0=time.Po4.10[1:12], y1=time.Po4.10[1:12],lty=1,col="gray")
lines(Po4.means.10[1:12],time.Po4.10[1:12],lwd=2,col="gray")
points(Po3.means.10[1:12],time.Po3.10[1:12],pch=21,bg="black", cex=1.4)
points(Po4.means.10[1:12],time.Po4.10[1:12],pch=21,bg="gray", cex=1.4)
plot(Panzano.means.10[1:12],time.Panzano.10[1:12], type="l", ylab="Calendar age", xlab="Corbula gibba abundance/10 years", 
xlim=c(0,150), main="Panzano",cex.main=0.9)
segments(x0=Panzano.025.10[1:12],x1=Panzano.975.10[1:12],y0=time.Panzano.10[1:12], y1=time.Panzano.10[1:12],lty=1)
lines(Panzano.means.10[1:12],time.Panzano.10[1:12],lwd=2)
points(Panzano.means.10[1:12],time.Panzano.10[1:12],pch=16, cex=1.4)


#######################################################################################################
#SHELL LOSS RATES
#######################################################################################################
#max mean age at Po3 - 106 years, 1/106 = 0.0094
#max mean age at Po4 - 112 years, 1/112  = 0.0089
#for Panzano, in Tomasovych et al. 2017 in Geology paper, we have the inverse of the median age - 1/365 = 0.0028

RATE.Po3=0.0094  
RATE.Po4=0.0089
RATE.Panzano=0.0028

#######################################################################################################
survival.onephase.Panzano=exp(-RATE.Panzano*(2013-time.Panzano.10))
exp.input.Panzano=Panzano.means.10/survival.onephase.Panzano
survival.onephase.Po3=exp(-RATE.Po3*(2013-time.Po3.10))
exp.input.Po3=Po3.means.10/survival.onephase.Po3
survival.onephase.Po4=exp(-RATE.Po4*(2013-time.Po4.10))
exp.input.Po4=Po4.means.10/survival.onephase.Po4

survival.onephase.Panzano.5=exp(-RATE.Panzano*(2013-time.Panzano.5))
exp.input.Panzano.5=Panzano.means.5/survival.onephase.Panzano.5
survival.onephase.Po3.5=exp(-RATE.Po3*(2013-time.Po3.5))
exp.input.Po3.5=Po3.means.5/survival.onephase.Po3.5
survival.onephase.Po4.5=exp(-RATE.Po4*(2013-time.Po4.5))
exp.input.Po4.5=Po4.means.5/survival.onephase.Po4.5

#######################################################################################################
#INCORPORATING THE EFFECTS OF LOSS AND BURIAL
#######################################################################################################

############################################################################################################################
#FIGURE SURFACE AFD
############################################################################################################################
par(mfrow=c(2,3))
par(mar=c(4,3,2,1))
hist(surfaceAFD.Po3.20, col="gray",breaks=seq(0,210,by=5), xlab="Postmortem age (y)", ylab="Number of specimens", ylim=c(0,50), main="Surface 20 cm-Po 3", cex.main=0.9)
hist(surfaceAFD.Po4.20, col="gray",breaks=seq(0,210,by=5), xlab="Postmortem age (y)", ylab="", ylim=c(0,50), main="Surface-20 cm-Po 4", cex.main=0.9)
hist(surfaceAFD.Panzano.6, col="gray",breaks=seq(0,210,by=5), ylab="", xlab="Postmortem age (y)", ylim=c(0,50), main="Surface 6 cm-Panzano", cex.main=0.9)

#######################################################################################################
#ACCOUNTING FOR MIXING AND BURIAL
#######################################################################################################
Panzano.time=2013-time.Panzano.10
survival.onephase=exp(-RATE.Panzano*Panzano.time)
Panzano.exp.input=Panzano.means.10/survival.onephase
Panzano.exp.inputUCI=Panzano.025.10/survival.onephase
Panzano.exp.inputLCI=Panzano.975.10/survival.onephase

#time.Po3.10
Po3.time=2013-time.Po3.10
survival.onephase=exp(-RATE.Po3*Po3.time)
Po3.exp.input=Po3.means.10/survival.onephase
Po3.exp.inputUCI=Po3.025.10/survival.onephase
Po3.exp.inputLCI=Po3.975.10/survival.onephase

#time.Po4.10
Po4.time=2013-time.Po4.10
survival.onephase=exp(-RATE.Po4*Po4.time)
Po4.exp.input=Po4.means.10/survival.onephase
Po4.exp.inputUCI=Po4.025.10/survival.onephase
Po4.exp.inputLCI=Po4.975.10/survival.onephase

#time.Panzano.5
Panzano.time.5=c(0,seq(8,1995,by=5))
survival.onephase.5=exp(-RATE.Panzano*Panzano.time.5)
Panzano.exp.input.5=Panzano.means.5/survival.onephase.5
Panzano.exp.inputUCI.5=Panzano.025.5/survival.onephase.5
Panzano.exp.inputLCI.5=Panzano.975.5/survival.onephase.5

#time.Po3.5
Po3.time.5=c(0,seq(8,295,by=5))
survival.onephase.5=exp(-RATE.Po3*Po3.time.5)
Po3.exp.input.5=Po3.means.5/survival.onephase.5
Po3.exp.inputUCI.5=Po3.025.5/survival.onephase.5
Po3.exp.inputLCI.5=Po3.975.5/survival.onephase.5

#time.Po4.5
Po4.time.5=c(0,seq(8,295,by=5))
survival.onephase.5=exp(-RATE.Po4*Po4.time.5)
Po4.exp.input.5=Po4.means.5/survival.onephase.5
Po4.exp.inputUCI.5=Po4.025.5/survival.onephase.5
Po4.exp.inputLCI.5=Po4.975.5/survival.onephase.5

#######################################################################################################
#FIGURE 18-20TH CENTURY
#######################################################################################################
#FACTOR IN 10-Y COHORTS IS 12.5 BECAUSE SURFACE AREA OF TWO CORES IS 0.04 M2, 
#so 1/0.04 = 25, AND 25/2 = 12.5 BECAUSE TWO LIFESPANS PER 10-YEAR COHORT
#if residence times or time to burial of  sedimentary particles is exponentially distributed,
#inter-quartile age range of particles in the surface mixed layer is equal to the mean particle age.
#therefore, the IQR can be compared against the estimate of temporal resolution that would be expected 
#in the absence of any mixing or in the absence of mixing below the FML.

#######################################################################################################
#FIGURE PALEOBIOLOGY 20TH CENTURY -  10-YEAR COHORTS
#######################################################################################################

CEX.AXIS=0.9
LOG=""
par(mfrow=c(2,3))
par(mar=c(6,4,2,1))
FACTOR.5=25/(10/5)
FACTOR.2=25/(10/2)

labelat=-c(13,33,53,73,93); labels=2013+labelat
par(mfrow=c(2,3))
par(mar=c(6,4,3,1))

plot(Po3.exp.input[1:12]*FACTOR.5, 2013-Po3.time[1:12], cex.main=0.95, xlab="Standing abundance (m2)", col="black", xlim=c(1,5500), ylim=c(1900,2015), ylab="Year", main="Po 3")
lines(Po3.exp.input[1:12]*FACTOR.5, 2013-Po3.time[1:12], col="black")
lines(Po3.exp.input[1:12]*FACTOR.2, 2013-Po3.time[1:12], col="gray")
segments(x0=Po3.exp.inputLCI[1:12]*FACTOR.2,x1=Po3.exp.inputUCI[1:12]*FACTOR.2,y0=2013-Po3.time[1:12], y1=2013-Po3.time[1:12])
segments(x0=Po3.exp.inputLCI[1:12]*FACTOR.5,x1=Po3.exp.inputUCI[1:12]*FACTOR.5,y0=2013-Po3.time[1:12], y1=2013-Po3.time[1:12])
points(Po3.exp.input[1:12]*FACTOR.5, 2013-Po3.time[1:12], pch=21, bg="black", cex=1.1)
points(Po3.exp.input[1:12]*FACTOR.2, 2013-Po3.time[1:12], pch=21, bg="gray", cex=1.1)
legend(x="bottomright", pch=c(21,21), pt.bg=c("gray","black"), legend=c("lifespan-2y","lifespan-5y"), cex=1.1, bty="n")

plot(Po4.exp.input[1:12]*FACTOR.5, 2013-Po4.time[1:12], cex.main=0.95, xlab="Standing abundance (m2)", col="black", xlim=c(1,5500), ylim=c(1900,2015), ylab="Year", main="Po 4")
lines(Po4.exp.input[1:12]*FACTOR.5, 2013-Po4.time[1:12], col="black")
lines(Po4.exp.input[1:12]*FACTOR.2, 2013-Po4.time[1:12], col="gray")
segments(x0=Po4.exp.inputLCI[1:12]*FACTOR.2,x1=Po4.exp.inputUCI[1:12]*FACTOR.2,y0=2013-Po4.time[1:12], y1=2013-Po4.time[1:12])
segments(x0=Po4.exp.inputLCI[1:12]*FACTOR.5,x1=Po4.exp.inputUCI[1:12]*FACTOR.5,y0=2013-Po4.time[1:12], y1=2013-Po4.time[1:12])
points(Po4.exp.input[1:12]*FACTOR.5, 2013-Po4.time[1:12], pch=22, bg="black", cex=1.1)
points(Po4.exp.input[1:12]*FACTOR.2, 2013-Po4.time[1:12], pch=22, bg="gray", cex=1.1)
legend(x="bottomright", pch=c(22,22), pt.bg=c("gray","black"), legend=c("lifespan-2y","lifespan-5y"), cex=1.1, bty="n")

plot(Panzano.exp.input[1:50]*FACTOR.5, 2013-Panzano.time[1:50], cex.main=0.95, xlab="Standing abundance (m2)", col="black", xlim=c(1,2000), ylim=c(1900,2015), ylab="Year", main="Panzano")
lines(Panzano.exp.input[1:50]*FACTOR.5, 2013-Panzano.time[1:50], col="black")
lines(Panzano.exp.input[1:50]*FACTOR.2, 2013-Panzano.time[1:50], col="gray")
segments(x0=Panzano.exp.inputLCI[1:50]*FACTOR.2,x1=Panzano.exp.inputUCI[1:50]*FACTOR.2,y0=2013-Panzano.time[1:50], y1=2013-Panzano.time[1:50])
segments(x0=Panzano.exp.inputLCI[1:50]*FACTOR.5,x1=Panzano.exp.inputUCI[1:50]*FACTOR.5,y0=2013-Panzano.time[1:50], y1=2013-Panzano.time[1:50])
points(Panzano.exp.input[1:50]*FACTOR.5, 2013-Panzano.time[1:50], pch=21, bg="black", cex=1.1)
points(Panzano.exp.input[1:50]*FACTOR.2, 2013-Panzano.time[1:50], pch=21, bg="gray", cex=1.1)
legend(x="bottomright", pch=21, pt.bg=c("gray","black"), legend=c("lifespan-2y","lifespan-5y"), cex=1.1, bty="n")

#######################################################################################################
#500 YEARS AT PANZANO
#######################################################################################################
Full.Panzano.means=apply(Full.Panzanobootout.10,MARGIN=2,mean)
Full.Panzano.025=apply(Full.Panzanobootout.10,MARGIN=2,quantile, 0.025)
Full.Panzano.975=apply(Full.Panzanobootout.10,MARGIN=2,quantile, 0.975)
BREAKS=c(0,seq(13,2000,by=10))
Corbula.decadal.scale=2013-BREAKS
plot(Panzano.exp.input[1:50]*FACTOR.5, 2013-BREAKS[1:50], xlab="Pred. abundance (10 years)", type="n", ylab="Time (years)",pch=16, log="", xlim=c(0,2000))
axis(4,at=c(2000,1900,1800,1700,1600),labels=rep("",5))
#rect(xleft=0, ybottom=1900,xright=2000, ytop=2000, col="gray", border=NA)
#rect(xleft=0, ybottom=1700,xright=2000, ytop=1800, col="gray", border=NA)
#rect(xleft=0, ybottom=1500,xright=2000, ytop=1600, col="gray", border=NA)
rect(xleft=0, ybottom=1770,xright=2000, ytop=1820, col="gray75", border=NA)
rect(xleft=0, ybottom=1870,xright=2000, ytop=1900, col="gray75", border=NA)
rect(xleft=0, ybottom=1950,xright=2000, ytop=1990, col="gray75", border=NA)
lines(Panzano.exp.input[1:50]*FACTOR.5, 2013-BREAKS[1:50], lwd=2)
segments(x0=Panzano.exp.inputLCI[1:50]*FACTOR.5,x1=Panzano.exp.inputUCI[1:50]*FACTOR.5,y0=2013-BREAKS[1:50],y1=2013-BREAKS[1:50])