Here is the function called ancestor.descendant.gaps which estimates joint gap between ancestor and descendant when we have pessimistic stratigraphical scenario when ancestor is immediately succeeded by the descendant (either by gradualism or punctuation). The gap between the two species depends on following parameters – probability of sampling of ancestor (p1), it’s descendant (p2), average duration of an ancestor, and it’s descendant in number of samples (n1, n2), and the average duration of the interval between samples k, and standard deviation of the sampling interval (sd). Number of simulated species pairs – parameter m (larger it is - more accurate the yielded distribution should be). In order to use the program copy all the code below (highlighted in red) and paste it in the R console. After that in the same console write the name of the function, put in it the parameters and push ENTER key. In the result you will see the histogram of probable gap sizes (in Ka) and the printed average, median gap sizes between ancestor and descendant species. Example: ancestor.descendant.gaps (n1=20, n2=30, p1=0.25, p2=0.3, m=1000, k=50, sd=20, gs=400) n1 – species' 1 duration in samples n2 – species' 2 duration in samples p1 – probability of sampling of species 1, found during the Gap Analysis (p=R) p2 – probability of sampling of species 2, found during the Gap Analysis (p=R) m – number of simulated species ancestor-descendant pairs k – average duration of a sampling interval (in Ka) sd – standard deviation of sampling interval (in Ka) gs – size of the gap which we want to know the probability of it (in Ka). ancestor.descendant.gaps<- function (n1, n2, p1, p2, m, k, sd, gs) { random.gaps<- function(n, m, p, k, sd) { a<-matrix(c(runif(n*m)),nrow=n, ncol=m) b<-ifelse(a <=p, 1, 0) c<-row(a, as.factor = FALSE) d<- b*c e<-ifelse(d<1, max(d), d) f<-apply(e, 2, min) f<-matrix (f, nrow=1) fkkk<- function (bb) { k1<-k sd1<-sd fff<-rnorm(bb,k1, sd1) fff1<- ifelse(fff <=0, 0, fff) zzz<-sum(fff1) } g <-apply(f, 2, fkkk) } set.seed(5) jj<- (random.gaps (n=n1, p=p1, m, k, sd))-rnorm(m, k, sd) j<- ifelse( jj<=0, 0, jj) set.seed(1115) s<- random.gaps (n=n2, p=p2, m, k, sd) w<-j+s hist(w, freq = TRUE, main="Distribution of gap sizes between ancestors and descendants", xlab="Gap size (Ka)", ylab="Frequency", col="green") ddd<- ifelse(w>= gs, 1, 0) probability.of.gap <-(sum(ddd))/(length(w)) small.p<-1/(length(w)) gap.prob<-ifelse(probability.of.gap==0, small.p , probability.of.gap) print("probability of chosen gap") print(gap.prob) summary(w) } ******************************************************************** Analytical method for homogenous probabilities This is analytical method that depends on the assumption that ancestor and descendant had the the same probability of preservation (therefore just one p). You can find probability of the gap with the size from st to ed using numerical integration of the equation (simplified case of eq 3.1 of Akkouchi, 2008): Example: analytical.homo.gap (0.15, 8, 800, 500, 4000 , 1, 1000) Parameters: p1 – probability of uccurence of species per sample k – average duration of a sampling interval (in Ka) l- length of sequence in Ka which you want to plot st – start of the interval of which probability you want to know ed – end of the interval of which probability you want to know sta1 – beginning of the interval where you searching mode of the distribution end1 – end of the interval where you searching mode of the distribution analytical.homo.gap <- function (p1, k, l, st, ed, sta1, end1 ) { lambda.1<- (-log(1-p1))/k time.Ka<-seq(1:l) homog.gap<- ((lambda.1)^2)* time.Ka*exp(-lambda.1* time.Ka) plot(time.Ka, homog.gap, type= "l", lwd = 5, xlab="Gap size (Ka)", ylab="Probability density") y <-function (time.Ka) { ((lambda.1)^2)* time.Ka*exp(-lambda.1* time.Ka) } print("Mode") mode<-optimize(y, interval=c(sta1, end1), maximum=TRUE) print (mode) print("probability") integrate(y, st, ed) } This approach assumes that probabilities of occurrence (p1, p2) are different for ancestors and descendants. All other parameters identical to the previous function. Studied equation (simplified case of eq 2.3 of Akkouchi, 2008): Example: analytical.hetero.gap (0.15, 0.16, 8, 1500, 300, 400 , 1, 1000) p1 – probability of sampling of species 1, found during the Gap Analysis (p=R) p2 – probability of sampling of species 2, found during the Gap Analysis (p=R) k – average duration of a sampling interval (in Ka) l- length of sequence in Ka which you want to plot st – start of the interval of which probability you want to know ed – end of the interval of which probability you want to know sta1 – beginning of the interval where you searching mode of the distribution end1 – end of the interval where you searching mode of the distribution analytical.hetero.gap<- function (p1, p2, k, l, st, ed, sta1, end1) { lambda.1<- (-log(1-p1))/k lambda.2<- (-log(1-p2))/k time.Ka<-seq(1:l) hetero.gap <-((lambda.1* lambda.2)/( lambda.2- lambda.1))*exp(-time.Ka *lambda.1)+((lambda.1* lambda.2)/( lambda.1-lambda.2))*exp(-time.Ka*lambda.2) plot(time.Ka, hetero.gap, type= "l", lwd = 5, xlab="Gap size (Ka)", ylab="Probability density") y <-function (time.Ka) { ((lambda.1* lambda.2)/( lambda.2- lambda.1))*exp(-time.Ka *lambda.1)+((lambda.1* lambda.2)/( lambda.1-lambda.2))*exp(-time.Ka*lambda.2) } print("Mode") mode<-optimize(y, interval=c(sta1, end1), maximum=TRUE) print (mode) print("probability") integrate(y, st, ed) }