################################################################################################# #GENERATE TRANSITION MATRIX FOR TWO LAYERS (SURFACE AND SUBSURFACE) WITH DIAGENETIC STABILIZATION ################################################################################################# if (model=="low-resolution" & stabilization==TRUE) { Amat<-function(par1,n){ par=pmax(pmin(par1, 100),-100) Qmat<-matrix(data=0,nrow=n,ncol=n) ################################################### #disintegration rate in surface layer 1 - TAZ ################################################### Qmat[1,4]<-exp(par[1]) ################################################### #disintegration rate in subsurface layer 2 - SZ ################################################### Qmat[2,4]<-exp(par[2]) ################################################### #latent disintegration rate in TAZ for shells exhumed from SZ #################################################### Qmat[3,4]<-exp(par[5]) ################################################### #burial from surface layer (TAZ) to to subsurface SZ ################################################### Qmat[1,2]<-exp(par[3]) ################################################### #exhumation from SZ to TAZ ################################################### Qmat[2,3]<-exp(par[4]) ################################################### #burial from TAZ for exhumed shells back to SZ ################################################### Qmat[3,2]<-exp(par[3]) Qmat[4,1]<-Qmat[4,2]<-Qmat[4,3]<-Qmat[4,4]<-0 Qmat[1,1]<- -sum(Qmat[1,-1]) Qmat[2,2]<- -sum(Qmat[2,-2]) Qmat[3,3]<- -sum(Qmat[3,-3]) Qmat[4,4]<- -sum(Qmat[4,-4]) Qmat } n=4 par<-rep(0,5) Amat(par1=par,n) init=par } ################################################################################################# #GENERATE TRANSITION MATRIX FOR TWO LAYERS WITHOUT DIAGENETIC STABILIZATION IN THE SUBSURFACE LAYER ################################################################################################# if (model=="low-resolution" & stabilization==FALSE) { Amat<-function(par1,n){ par=pmax(pmin(par1, 100),-100) Qmat<-matrix(data=0,nrow=n,ncol=n) ################################################### #disintegration rate in layer 1 - TAZ ################################################### Qmat[1,3]<-exp(par[1]) ################################################### #disintegration rate in layer 2 - SZ ################################################### Qmat[2,3]<-exp(par[2]) ################################################### #burial from 1 to 2 ################################################### Qmat[1,2]<-exp(par[3]) ################################################### #exhumation from layer 2 to latent 1 ################################################### Qmat[2,3]<-exp(par[4]) Qmat[3,1]<-Qmat[3,2]<-Qmat[3,3]<-0 Qmat[1,1]<- -sum(Qmat[1,-1]) Qmat[2,2]<- -sum(Qmat[2,-2]) Qmat[3,3]<- -sum(Qmat[3,-3]) Qmat } n=3 par<-rep(0,4) Amat(par1=par,n) init=par } ################################################################################################################################# #Estimate maximum-likelihood parameters, given the definition of the transition-rate matrix with Amat function ################################################################################################################################# myest<-function(par,Amat,n,data,model,stabilization){ #compute negative log likelihood for a specific parameter and model #function is used in the optim or optimParallel function mylik<-function(par,Amat,n,data,model, stabilization){ require(expm) #Amat functions are defined below Qmat<-Amat(par,n) ##likelihood - defined separately for each model if (model=="TAZ-1-1 layer-sequestration") { myf<-function(v){ s<-c(1,rep(0,n-1))%*%expm(as.numeric(v[2])*Qmat) (s[1,1]+s[1,2])*as.numeric(v[1]==1) } } if (model=="low-resolution" & stabilization==TRUE) { myf<-function(v){ s<-c(1,rep(0,n-1))%*%expm(as.numeric(v[2])*Qmat) #in the scenario with stabilization, densitities for the upper increment (TAZ) with shells that were not sequestered into SZ #and densities for shells that were exhumed from the SZ to the TAZ are added s[1,v[1]]+s[1,3]*as.numeric(v[1]==1) } } if (model=="low-resolution" & stabilization==FALSE) { myf<-function(v){ s<-c(1,rep(0,n-1))%*%expm(as.numeric(v[2])*Qmat) s[1,1]*as.numeric(v[1]==1)+s[1,2]*as.numeric(v[1]==2) } } 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,data=data,model=model,stabilization=stabilization, 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,data=data,stabilization=stabilization, 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) #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 } #################################################################################################### #final function that involves parameter estimation /myest function/ #################################################################################################### onesim<-function(data){ est=myest(par,Amat,n,data,model, stabilization) model=model par<-est$Parameters #transition matrix generated by Amat functions below (depending on whether model has diagenetic stabilization or not) Qmat<-est$Transition_matrix #transition matrix generated by Amat functions below (depending on whether #age range over which age distributions are modeled, in years/1000, here from 1 year to 55,000 years 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 or increment) #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){ require(expm) 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 visualizations 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 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 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 }