################################################################################################################################### #GENERATE MATRIX WITHOUT DIAGENETIC STABILIZATION AND WITHOUT ANY SHIFT IN SEDIMENTATION ################################################################################################################################### if (stabilization==FALSE & SED.SHIFT==0) { par<-par1 Amat<-function(par,n,INCREMENTS, TAZ, SED.SHIFT){ Qmat<-matrix(data=0,nrow=n,ncol=n) SED.SHIFT=SED.SHIFT TAZ=TAZ DEPTH=INCREMENTS dis.params=INCREMENTS mix.params=INCREMENTS-1 bur.param=dis.params+mix.params+1 #################################### #disintegration in individual increments #################################### for (i in 1:(INCREMENTS-1)) { Qmat[i,(INCREMENTS+1)]<-exp(par[i]) } #################################### #loss from the last increment (by disintegration and/or burial below it) #################################### Qmat[INCREMENTS,(INCREMENTS+1)]<-exp(par[INCREMENTS])+exp(par[bur.param]) #loss in TAZ + burial #################################### #exhumation by burrowers #################################### if (DEPTH>1) { for (ii in 2:DEPTH) { Qmat[ii,1]<-exp(par[(INCREMENTS-1)+ii]) } } #################################### #burial by mixing/burrowers #################################### for (ii in 2:DEPTH) { Qmat[1,ii]<-exp(par[(INCREMENTS-1)+ii]) } #################################### #burial by sedimentation #################################### for (ii in 1:(INCREMENTS-1)) { Qmat[ii,ii+1]<-exp(par[bur.param]) } #################################### #burial by sedimentation and burial by mixing overwrites the transition from the first to the second increment #################################### if (DEPTH>1) Qmat[1,2]<-exp(par[dis.params+mix.params+1])+exp(par[INCREMENTS+1]) #mixing+burial for (ii in 1:INCREMENTS) { Qmat[ii,ii]<- -sum(Qmat[ii,-ii]) } Qmat } n<-INCREMENTS+1; par<-par; init=par; } ################################################################################################################################### #GENERATE MATRIX WITH DIAGENETIC STABILIZATION AND WITHOUT ANY SHIFT IN SEDIMENTAT ACCUMULATION ################################################################################################################################### if (stabilization==TRUE & SED.SHIFT==0) { par<-par1 Amat<-function(par,n,INCREMENTS, TAZ, SED.SHIFT){ Qmat<-matrix(data=0,nrow=n,ncol=n) for (i in 1:(INCREMENTS-1)) { Qmat[i,(INCREMENTS+TAZ+1)]<-exp(par[i]) #disintegrations } SED.SHIFT=SED.SHIFT TAZ=TAZ DEPTH=INCREMENTS dis.params=INCREMENTS mix.params=INCREMENTS-1 bur.param=dis.params+mix.params+1 seq.param=bur.param+1 #################################### #if TAZ=1, no exhumation from 2nd increment to TAZ, but from the latent increment to TAZ instead #exhumation (avoid exhumation from SZ) #################################### if (DEPTH>1) { for (ii in 2:(TAZ)) { Qmat[ii,1]<-exp(par[(INCREMENTS-1)+ii]) } for (ii in (TAZ+1):INCREMENTS) { Qmat[ii,INCREMENTS+1]<-exp(par[(INCREMENTS-1)+ii]) } } #################################### #burial by mixing in the ML #################################### for (ii in 2:DEPTH) { Qmat[1,ii]<-exp(par[(INCREMENTS-1)+ii]) } #################################### #burial #################################### for (ii in 1:(INCREMENTS-1)) { Qmat[ii,ii+1]<-exp(par[bur.param]) } #################################### #loss from the last core increment #################################### Qmat[INCREMENTS,(INCREMENTS+1+TAZ)]<-exp(par[bur.param])+exp(par[INCREMENTS]) #################################### #first burial and burial by mixing in the first cell overwrites 1-2 transition #################################### Qmat[1,2]<-exp(par[bur.param])+exp(par[INCREMENTS+1]) #mixing+burial ############################################################################################### #LATENT PART APPLIES TO SHELLS THAT WERE SEQUESTERED AND EXHUMED BACK INTO THE TAZ ############################################################################################### #exhumation part from SZ to TAZ #################################### Qmat[TAZ+1,TAZ]<-0 #################################### #diagenetic stabilization in latent increments #################################### for (ii in 1:TAZ) { Qmat[(INCREMENTS+ii),(INCREMENTS+TAZ+1)]<-exp(par[seq.param]) } #################################### #new exhumation among latent increments #################################### if (TAZ>1) { for (ii in 1:(TAZ-1)) { if (ii=DEPTH) Qmat[(INCREMENTS+ii+1),(INCREMENTS+1)]<-0 } #################################### #new burial by mixing among latent increments #################################### for (ii in 1:(TAZ-1)) { if (ii1) Qmat[(INCREMENTS+ii),(INCREMENTS+ii+1)]<-exp(par[bur.param]) } #################################### #loss from the last TAZ increment into the top SZ increment #################################### Qmat[INCREMENTS+TAZ,(1+TAZ)]<-exp(par[bur.param]) } ############################################ for (ii in 1:(INCREMENTS+TAZ)) { Qmat[ii,ii]<- -sum(Qmat[ii,-ii]) } Qmat } n<-INCREMENTS+1+TAZ; par<-par; init=par; } ################################################################################################################################### #GENERATE MATRIX WITHOUT DIAGENETIC STABILIZATION AND WITH ONE SHIFT IN SEDIMENTATION ################################################################################################################################### if (stabilization==FALSE & SED.SHIFT>0) { par<-par1 Amat<-function(par,n,INCREMENTS, TAZ, SED.SHIFT){ Qmat<-matrix(data=0,nrow=n,ncol=n) SED.SHIFT=SED.SHIFT TAZ=TAZ DEPTH=INCREMENTS dis.params=INCREMENTS mix.params=INCREMENTS-1 bur.param.1=dis.params+mix.params+1 bur.param.2=bur.param.1+1 ######################################## #disintegration ######################################## for (i in 1:(INCREMENTS-1)) { Qmat[i,(INCREMENTS+1)]<-exp(par[i]) #loss in TAZ } ######################################## #loss from the last increment - ######################################## Qmat[INCREMENTS,(INCREMENTS+1)]<-exp(par[INCREMENTS])+exp(par[bur.param.2]) #loss + burial ######################################## #exhumation ######################################## if (DEPTH>1) { for (ii in 2:DEPTH) { Qmat[ii,1]<-exp(par[(INCREMENTS-1)+ii]) } } ######################################## #burial by mixing in the ML ######################################## for (ii in 2:DEPTH) { Qmat[1,ii]<-exp(par[(INCREMENTS-1)+ii]) } ######################################## #burial ######################################## for (ii in 1:(INCREMENTS-1)) { if (ii=SED.SHIFT) Qmat[ii,ii+1]<-exp(par[bur.param.2]) } ######################################## #first burial and burial by mixing in the first cell overwrites 1-2 transition ######################################## if (DEPTH>1) Qmat[1,2]<-exp(par[bur.param.1])+exp(par[INCREMENTS+1]) #mixing+burial for (ii in 1:INCREMENTS) { Qmat[ii,ii]<- -sum(Qmat[ii,-ii]) } Qmat } n<-INCREMENTS+1; par<-par; init=par; } ################################################################################################################################### #GENERATE MATRIX WITH DIAGENETIC STABILIZATION AND WITH SHIFT IN SEDIMENTATION ################################################################################################################################### if (stabilization==TRUE & SED.SHIFT>0) { par<-par1 Amat<-function(par,n,INCREMENTS, TAZ,SED.SHIFT){ Qmat<-matrix(data=0,nrow=n,ncol=n) for (i in 1:(INCREMENTS-1)) { Qmat[i,(INCREMENTS+TAZ+1)]<-exp(par[i]) #disintegrations } SED.SHIFT=SED.SHIFT TAZ=TAZ DEPTH=INCREMENTS dis.params=INCREMENTS mix.params=INCREMENTS-1 bur.param.1=dis.params+mix.params+1 bur.param.2=bur.param.1+1 seq.param=bur.param.1+2 ######################################## #if TAZ=1, no exhumation from 2nd increment to TAZ, but from the latent increment to TAZ instead ######################################## if (DEPTH>1) { for (ii in 2:(TAZ)) { Qmat[ii,1]<-exp(par[(INCREMENTS-1)+ii]) } for (ii in (TAZ+1):INCREMENTS) { Qmat[ii,INCREMENTS+1]<-exp(par[(INCREMENTS-1)+ii]) } } ######################################## #burial by mixing in the ML ######################################## for (ii in 2:DEPTH) { Qmat[1,ii]<-exp(par[INCREMENTS+ii]) } ######################################## #burial ######################################## for (ii in 1:(INCREMENTS-1)) { if (ii=SED.SHIFT) Qmat[ii,ii+1]<-exp(par[bur.param.2]) } ######################################## #loss from the last core increment ######################################## Qmat[INCREMENTS,(INCREMENTS+1+TAZ)]<-exp(par[bur.param.2])+exp(par[INCREMENTS]) ######################################## #first burial and burial by mixing in the first cell overwrites 1-2 transition ######################################## Qmat[1,2]<-exp(par[bur.param.1])+exp(par[INCREMENTS+1]) #mixing+burial ############################################################################################### #LATENT PART APPLIES TO SHELLS THAT WERE SEQUESTERED AND EXHUMED BACK INTO THE TAZ ############################################################################################### #exhumation part from SZ to TAZ Qmat[TAZ+1,TAZ]<-0 #diagenetic stabilization in latent increments for (ii in 1:TAZ) { Qmat[(INCREMENTS+ii),(INCREMENTS+TAZ+1)]<-exp(par[seq.param]) } #exhumation among latent increments if (TAZ>1) { for (ii in 1:(TAZ-1)) { if (ii=DEPTH) Qmat[(INCREMENTS+ii+1),(INCREMENTS+1)]<-0 } #burial by mixing among latent increments for (ii in 1:(TAZ-1)) { if (ii1 & ii1 & ii>=SED.SHIFT) Qmat[(INCREMENTS+ii),(INCREMENTS+ii+1)]<-exp(par[bur.param.2]) } #loss from the last TAZ increment by burial Qmat[INCREMENTS+TAZ,(1+TAZ)]<-exp(par[bur.param.2]) } ############################################ for (ii in 1:(INCREMENTS+TAZ)) { Qmat[ii,ii]<- -sum(Qmat[ii,-ii]) } Qmat } n<-INCREMENTS+1+TAZ; par<-par; init=par; } ############################################################################################### #Find maximum-likelihood estimates of parameters ############################################################################################### myest<-function(par,Amat,n,data,INCREMENTS, TAZ, SED.SHIFT, model, stabilization){ #compute negative log likelihood for a specific parameter and model #using optim or optimParallel functions below mylik<-function(par,Amat,n,data,INCREMENTS, TAZ, SED.SHIFT, model, stabilization){ require(expm) Qmat<-Amat(par,n, INCREMENTS, TAZ, SED.SHIFT) ######################### #likelihood functions ######################### if (model=="high-resolution" & stabilization==FALSE) { myf<-function(v){ s<-c(1,rep(0,n-1))%*%expm(as.numeric(v[2])*Qmat) s[1,v[1]] } } if (model=="high-resolution" & stabilization==TRUE) { myf<-function(v){ s<-c(1,rep(0,n-1))%*%expm(as.numeric(v[2])*Qmat) TOTAL=INCREMENTS+TAZ sum.temp=numeric() #densities for TAZ increments - sum of densities for shells that were not sequestered from the underlying SZ and sequestered shells for (xx in 1:TAZ) { sum.temp[xx]=(s[1,xx]+s[1,xx+INCREMENTS])*as.numeric(v[1]==xx) } #densities for SZ increments - the same process as in the model without stabilization for (xx in (TAZ+1):INCREMENTS) { sum.temp[xx]=s[1,xx]*as.numeric(v[1]==xx) } s[1,1:INCREMENTS]=sum.temp } } myff<-function(t){ c(1,rep(0,n-1))%*%expm(t*Qmat)%*%c(rep(1,n-1),0) } tseq<-exp(seq(-6.5,4,by=0.01)) term1=apply(data,1,myf) term2=sum(sapply(tseq,myff)*(tseq-c(0,tseq[-length(tseq)]))) term1=ifelse(term1<10e-100, 10e-100, term1) term2=ifelse(term2<10e-100, 10e-100, term2) #the signs are opposite as one is the numerator and the other is the denominator. -sum(log(term1))+dim(data)[1]*log(term2) } ################################################################################################################################# #Use quasi-newton method to minimize the negative log likelihood (i.e. maximize likelihood) if (PARALLEL==TRUE) { cl <- makeCluster(4) # set the number of processor cores setDefaultCluster(cl=cl) # set 'cl' as default cluster fit1<-optimParallel(par,mylik, Amat=Amat,n=n,INCREMENTS=INCREMENTS, TAZ=TAZ, SED.SHIFT=SED.SHIFT, stabilization=stabilization, data=data,model=model,method="L-BFGS-B", lower = log(1/100000*1000), upper = log(10*1000)) stopCluster(cl) } if (PARALLEL==FALSE) { fit1<-optim(par,mylik, Amat=Amat,model=model,n=n,INCREMENTS=INCREMENTS, TAZ=TAZ, SED.SHIFT=SED.SHIFT, stabilization=stabilization, data=data,method="L-BFGS-B", lower = log(1/100000*1000), upper = log(10*1000)) } par.est<-fit1$par #estimated transition matrix Qmat.est<-Amat(par.est,n, INCREMENTS, TAZ, SED.SHIFT) #maximized likelihood, i.e., adding minus sign just reversed the minus sign used in minimization lik.est<--fit1$value result<-list("Parameters"=par.est,"Transition_matrix"=Qmat.est,"Log_likelihood"=lik.est) result } #################################################################################################### #input = data from the function upload age data #onesim is a function that includes the myest function that finds maximum-likelihood parameters #################################################################################################### onesim<-function(data){ est=myest(par,Amat,n,data,INCREMENTS, TAZ, SED.SHIFT, model, stabilization) model=model par<-est$Parameters Qmat<-est$Transition_matrix tseq<-exp(seq(-6.5,4,by=0.01)) #each age interval is multiplied by the probability of shell preservation in the given state (layer) #and the resulting square matrix is exponentiated, this is performed for all xx states, so "res" is array with xx rows and xx ages myf<-function(t){ c(1,rep(0,n-1))%*%expm(t*Qmat) } BIN=500 res<-sapply(tseq,myf,simplify=TRUE) scale<-nrow(data)/sum(res[-n,]) #age distribution output binned to 500 years, useful for fit visualization non.log.tseq<-seq(50,50000,by=BIN)/1000 non.log.res<-sapply(non.log.tseq,myf,simplify=TRUE) non.log.scale<-nrow(data)/sum(non.log.res[-n,]) #age distribution output binned to 50 years, useful for fit visualization non.log.tseq.bin50<-seq(25,50000,by=50)/1000 non.log.res.bin50<-sapply(non.log.tseq.bin50,myf,simplify=TRUE) non.log.scale.bin50<-nrow(data)/sum(non.log.res.bin50[-n,]) #age distribution output binned to 10 years, useful for fit visualization non.log.tseq.bin10<-seq(0.1,50000,by=10)/1000 non.log.res.bin10<-sapply(non.log.tseq.bin10,myf,simplify=TRUE) non.log.scale.bin10<-nrow(data)/sum(non.log.res.bin10[-n,]) out <- list(pars = c(est$Parameters), res=res, scale=scale, Qmat=Qmat, tseq=tseq, non.log.res=non.log.res, non.log.scale=non.log.scale, non.log.tseq=non.log.tseq, non.log.res.bin50=non.log.res.bin50, non.log.scale.bin50=non.log.scale.bin50, non.log.tseq.bin50=non.log.tseq.bin50, non.log.res.bin10=non.log.res.bin10, non.log.scale.bin10=non.log.scale.bin10, non.log.tseq.bin10=non.log.tseq.bin10, loglik=est$"Log_likelihood") class(out) <- "onesim" out }