proba.R 3.25 KB
Newer Older
Eric Coissac committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
#' Logarithm of the sum of two exponentials
#' 
#' `log_add_exp` computes the logarithm of the sum
#' of the exponential of two value.
#' 
#' \deqn{log\_add\_exp(a,b) = \log(e^a + e^b)}{%
#'       log_add_exp(a,b) = log(exp(a) + exp(b))}
#' 
#' @param a first logarithmic values to add 
#' @param b second logarithmic values to add 
#' 
#' @return a numeric vector containing the logarithm of the sum
#' of the exponential of `a` and `b`.
#'
#' @author Eric Coissac
#' @export
log_add_exp = function(a, b) {
     a + log1p(exp(b - a))
}

#' Logarithm of the difference of two exponentials
#' 
#' `log_add_exp` computes the logarithm of the difference
#' of the exponential of two value.
#' 
#' \deqn{log\_add\_exp(a,b) = \log(e^a - e^b)}{%
#'       log_add_exp(a,b) = log(exp(a) - exp(b))}
#' 
#' @param a first logarithmic values to add 
#' @param b second logarithmic values to add 
#' 
#' @return a numeric vector containing the logarithm of the 
#' difference of the exponential of `a` and `b`.
#'
#' @author Eric Coissac
#' @export
log_diff_exp = function(a, b) {
    a + log(1 - exp(b - a))
}

#' Logarithm of the sum of exponentials
#' 
#' `log_sum_exp` computes the logarithm of the sum
#' of exponential of values.
#' 
#' \deqn{log\_add\_exp(x_1,x_2,\ldots) = \log(e^{x_1} + e^{x_2} + \ldots)}{%
#'       log_add_exp(x_1,x_2,\ldots) = log(exp(x_1) + exp(x_2) + \ldots)}
#' 
#' @param x logarithmic values to sum 
#' 
#' @return a numeric value containing the logarithm of the 
#' sum of the exponential of the values of `x`.
#'
#' @author Eric Coissac
#' @export
log_sum_exp = function(x) Reduce(logaddexp, x)

58
#' Density probability function of a multinomial distribution for nucleotides
Eric Coissac committed
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
#' 
#' Computes the probalility to observe a given number of `A`, `C`, `G`, `T`
#' when probalilities of each nucleotide are `pA`, `pC`, `pG`, `pT`.
#' 
#' @param A numeric; number of observed *A* at a given site.
#' @param C numeric; number of observed *C* at a given site.
#' @param G numeric; number of observed *G* at a given site.
#' @param T numeric; number of observed *T* at a given site.
#' @param pA numeric; probability to observed a *A* at a given site.
#' @param pC numeric; probability to observed a *C* at a given site.
#' @param pG numeric; probability to observed a *G* at a given site.
#' @param pT numeric; probability to observed a *T* at a given site.
#' @param log logical; if TRUE, probabilities p are given as log(p).
#' 
#' @return numeric; the density of probability of the log that probability
#'         of the observations given the probabilities of observation
#'         
#' @examples
77 78 79 80 81
#'   dmultinomialACGT(c(3,4),c(3,6),c(3,2),c(3,5),
#'                    0.25,0.25,0.25,0.25)
#'   dmultinomialACGT(c(3,4),c(3,6),c(3,2),c(3,5),
#'                    pA = 0.25, pC = 0.25, pG = 0.25, pT = 0.25,
#'                    log=TRUE)
Eric Coissac committed
82 83 84
#' 
#' @author Eric Coissac
#' @export
85
dmultinomialACGT = function(A,C,G,T,pA,pC,pG,pT,log=FALSE) {
Eric Coissac committed
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102

   lpcomb =(  log(pA) * A 
            + log(pC) * C 
            + log(pG) * G 
            + log(pT) * T)
    
   N = A + C + G + T
   
   tetranome = (lgamma(N+1) - lgamma(A+1) 
                            - lgamma(C+1) 
                            - lgamma(G+1) 
                            - lgamma(T+1))

   lp = lpcomb + tetranome
   
   if (log) lp else exp(lp)
}