proba.R 3.25 KB
 Eric Coissac committed Oct 19, 2018 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) Eric Coissac committed Oct 20, 2018 58 #' Density probability function of a multinomial distribution for nucleotides Eric Coissac committed Oct 19, 2018 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 Eric Coissac committed Oct 20, 2018 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 Oct 19, 2018 82 83 84 #' #' @author Eric Coissac #' @export Eric Coissac committed Oct 20, 2018 85 dmultinomialACGT = function(A,C,G,T,pA,pC,pG,pT,log=FALSE) { Eric Coissac committed Oct 19, 2018 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) }