[adegenet-commits] r429 - pkg/R

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

Author: jombart
Date: 2009-09-01 12:32:28 +0200 (Tue, 01 Sep 2009)
New Revision: 429

Modified:
   pkg/R/sequences.R
Log:
New Petery version for transition model.


Modified: pkg/R/sequences.R
===================================================================
--- pkg/R/sequences.R	2009-08-31 21:59:26 UTC (rev 428)
+++ pkg/R/sequences.R	2009-09-01 10:32:28 UTC (rev 429)
@@ -83,19 +83,23 @@
 ## time is taken into account
 ## output: matrix with term proba(rowIdx to colIdx)
 ##
-transiProb <- function(x, mu, dates){
+transiProb <- function(x, mu, dates, result=c("prob","dist")){
     ## MISC CHECKS ##
     if(!inherits(x,"DNAbin")) stop("x is not a DNAbin object")
     if(!require(ape)) stop("The package ape is required.")
+    result <- match.arg(result)
 
-
     ## COMPUTATIONS ##
 
     ## get numbers of differing nucleotides between sequences
     seq.length <- ncol(as.matrix(x))
     D <- as.matrix(dist.dna(x, model="raw")) * seq.length
+    if(sum(D-round(D)) > 1e-10){ # make sure we've got integers there
+        warning("Number of nucleotides are not all integers")
+    }
+    D <- round(D)
 
-    ## compute matrix T
+    ## compute matrix T (time between sequences)
     if(inherits(dates,"POSIXct")){ # dates in POSIXct format
         temp <- outer(dates, dates, difftime, unit="days")
         T <- -matrix(as.numeric(temp),ncol=length(dates))
@@ -104,24 +108,22 @@
     }
 
     ## spot negative times
-    toSetToInf <- T < 1e-15
+    toSetToNull <- T < 1e-15
 
-    ## compute term1 and term2
-    mu <- mu/365
-    term1 <- exp(-mu * T) + (1 - exp(-mu * T))/4
-    term1[toSetToInf] <- 0
+    ## compute proba(no change @ a site) term
+    mu <- mu/365 # express mu per day
+    p1 <- exp(-T*mu) + (1-exp(-T*mu))/4
+    res <- dbinom(D, size=seq.length, prob=(1-p1))
 
-    term2 <- (1 - exp(-mu * T))/4
-    term2[toSetToInf] <- 0
+    ## PROCESS/RETURN RESULT
+    if(result=="prob"){ # return probabilities
+        res[toSetToNull] <- 0
+        diag(res) <- 1
+    } else { # return d = -log(proba)
+        res <- -log(res)
+        res[toSetToNull] <- 1e15
+        diag(res) <- 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