[adegenet-commits] r425 - pkg/R

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2009-08-27 20:40:07 +0200 (Thu, 27 Aug 2009)
New Revision: 425

Modified:
   pkg/R/sequences.R
Log:
added function to compute distance based on transition proba.
Have to stick to -logs(proba): raw proba won't work because of computer numeric limitations.


Modified: pkg/R/sequences.R
===================================================================
--- pkg/R/sequences.R	2009-06-28 13:17:47 UTC (rev 424)
+++ pkg/R/sequences.R	2009-08-27 18:40:07 UTC (rev 425)
@@ -68,3 +68,60 @@
 
     return(res)
 } # end DNAbin2genind
+
+
+
+
+
+
+
+###############
+## transiProb
+###############
+##
+## distance based on transition prob from one sequence to another
+## time is taken into account
+## output: matrix with term proba(rowIdx to colIdx)
+##
+transiProb <- function(x, mu, dates){
+    ## MISC CHECKS ##
+    if(!inherits(x,"DNAbin")) stop("x is not a DNAbin object")
+    if(!require(ape)) stop("The package ape is required.")
+
+
+    ## COMPUTATIONS ##
+
+    ## get numbers of differing nucleotides between sequences
+    seq.length <- ncol(as.matrix(x))
+    D <- as.matrix(dist.dna(x, model="raw")) * seq.length
+
+    ## compute matrix T
+    if(inherits(dates,"POSIXct")){ # dates in POSIXct format
+        temp <- outer(dates, dates, difftime, unit="days")
+        T <- -matrix(as.numeric(temp),ncol=length(dates))
+    } else { # dates are numeric
+        T <- -outer(dates, dates, "-")
+    }
+
+    ## spot negative times
+    toSetToInf <- T < 1e-15
+
+    ## compute term1 and term2
+    mu <- mu/365
+    term1 <- exp(-mu * T) + (1 - exp(-mu * T))/4
+    term1[toSetToInf] <- 0
+
+    term2 <- (1 - exp(-mu * T))/4
+    term2[toSetToInf] <- 0
+
+    ## compute res
+    ## res <- term1^D * term2^(seq.length-D)
+    res <- log(term1)*D + log(term2)*(seq.length-D)
+    res <- -res
+
+    res[toSetToInf] <- 1e15
+    diag(res) <- 0
+
+    ## RETURN RESULT ##
+    return(res)
+} # end transiProb