SedimentDNA_MIII<-function(d,PrPfp=c(1,1),PrPfn0=c(1,1),Prr_decay=1000,PrPA=0.5,lambda_c=1000,lambda_e=1000,niter=1000,samplefreq=1,burnin=0){ #d= matrix containing sample depth (in first column, called "depth") #followed by parameter on Bernoulli prior on DNA amplification success. One column for each primer. #PrPfp= 2-element vector containing the parameters of the beta prior on the false positive rate #PrPfn0= 2-element vector containing the parameters of the beta prior on the false negative rate at zero depth #Prr_decay=parameter of the exponential prior on the rate of DNA decay with depth. #PrPA= prior on the probability that dna is present. Either a scalar (which then applies at each sample depth), or #a vector with separate values for all sample depths, or a two-element vector providing values for the #upper and lower samples, from which values for the samples in between are calculated by linear interpolation. #lambda_c= the (Dirac delta prior on) colonization rate #lambda_e= the (Dirac delta prior on) extinction rate #niter= (integer)number of MCMC iterations #samplefreq= (integer) frequency to retain MCMC samples #burnin= number of iterations to discard as burn-in. the MCMC chain will consist of (niter-burnin)/samplefreq samples. GetLik<-function(obs,PA,Pfn,Pfp){ #likelihood function L<-matrix(nrow=length(PA), ncol=length(Pfp))#empty matrix to store likelihoods for (j in 1:length(Pfp)) {#one primer pair at a time IS LOOP NECESSARY HERE? ind<-PA==1 & (obs[,j]==0) #false negatives L[ind,j]<-Pfn[ind,j]#[ind] ind<-PA==1 & (obs[,j]==1) #not false negatives L[ind,j]=1-Pfn[ind,j] ind<-PA==0 & (obs[,j]==0) #not false positives L[ind,j]=1-Pfp[j]#[ind] ind<-PA==0 & (obs[,j]==1)#false positives L[ind,j]=Pfp[j] }#end of for-loop GetLik<-log(L) }#end of function GetLik Getdi<-function(d,dmax,PrPA,lambda_c,lambda_e){#Sample vector of colonization/extinction depths di<-matrix(,nrow=0,ncol=2)#initialize matrix with depths and events PA<-!logical(length=length(d$depth))#vector to contain presence/absence status at each sampled depth iPA<-(runif(1)0) { k<-k+1 if(cstat==0) {lambda<-1/lambda_c} else {lambda<-1/lambda_e} if(k==1 && iPA!=cPA){ cdepth<-dmax+log(1-runif(1)*(1-exp(-dmax*lambda)))/lambda } else { cdepth<-cdepth-rexp(1,lambda) #depth of next event without conditioning } if (cdepth<0) {break} #top of sediment core reached else { cstat<-!cstat di<-rbind(di,cbind(cdepth,cstat))#0=extinction, 1=colonization PA[d$depthburnin & ((iter-burnin) %% samplefreq)==0){ #this is just to store the values in the output variable ind=(iter-burnin)/samplefreq mcmc$L[ind]<-sum(sum(L)) mcmc$LPA[ind]<-sum(LPA) mcmc$PA[ind,]<-PA mcmc$d_ce[[ind]]<-di#the length of di may vary between iterations mcmc$Pfp[ind,]<-Pfp mcmc$Pfn0[ind,]<-Pfn0 mcmc$r_decay[ind,]<-r_decay }#end if }#end of mcmc iteration loop SedimentDNA_MIII<-mcmc #assign output variable }#end of main function