mismatch <- function(matrix_pol,populations,THETA){
npops <- length(populations)
popnames <- paste("pop",1:npops)
if(npops!=length(THETA)){warning("You have to choose a theta value for every popualtion ")}
RE <- matrix(NaN,3,npops)
rownames(RE) <- c("MDSD","MDG1","MDG2")
colnames(RE) <- popnames
for(xx in 1:npops){
if(length(populations[[xx]])==0){next;}
thetaT <- THETA[xx]
popmatrix <- matrix_pol[populations[[xx]],,drop=FALSE]
poplength <- length(populations[[xx]])
if(poplength>1){seqpairs <- combn(poplength,2)}else{seqpairs<-rbind(1,1)}# only one sequence
ncomb <- choose(poplength,2)
res <- apply(seqpairs,2,function(x){
seq1 <- popmatrix[x[1],]
seq2 <- popmatrix[x[2],]
id <- which(!is.na(seq1) & !is.na(seq2) & seq1!=seq2)
freq1 <- length(id)
moment <- freq1 - thetaT
freq2 <- moment*moment
freq3 <- moment*moment*moment
freq4 <- moment*moment*moment*moment
return(c(freq2,freq3,freq4))
})
rownames(res) <- c("freq2","freq3","freq4")
E <- rowSums(res)
s2 <- E["freq2"]
g1 <- E["freq3"]
g2 <- E["freq4"]
value <- s2/(ncomb-1)
if(!is.na(value)){
s <- sqrt(s2/(ncomb-1))
RE["MDSD",xx] <- s
}else{s <- FALSE}
if(s){
RE["MDG1",xx] <- g1*ncomb/((ncomb-2)*((ncomb-1)*s*s*s))
RE["MDG2",xx] <- g2* (ncomb-1) * ncomb /((ncomb-3)*(ncomb-2)*(ncomb-1)*s*s*s*s)-(3*(ncomb-1)*(ncomb-1))/((ncomb-3)*(ncomb-2))
}
}# END for over all POPS
return(RE)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.