[R] strange behavior of lchoose in combinatorics problem

Gonçalo Ferraz gferraz29 at gmail.com
Sat Jun 25 16:13:54 CEST 2016


Hi,

I am working on interactions between animals, studying whether animal 1 is attracted to animal 2 (or vice-versa). I looked for the two animals in N=2178 sampling occasions, finding animal 1 a total of N1=165 times, and animal 2 a total of N2=331 times. In J=97 occasions, I saw both animals at the same time. 

The more frequently I see the two animals in the same sampling occasion, the more I will believe that they are attracted to each other. Therefore, I want to calculate the probability of finding J<=97 when the two animals are moving around independently of each other. The higher this probability, the stronger the attraction.

Following Veech (Journal of Biogeography 2014, 41: 1029-1035) I compute the log probability of obtaining a number n of encounters between animals as ‘lpn’ in the function lveech:

lveech <- function(N,N1,N2,n) {
    a <- lchoose(N,n)
    b <- lchoose(N-n,N2-n)
    c <- lchoose(N-N2,N1-n)
    d <- lchoose(N,N2)
    e <- lchoose(N,N1)
    lpn <- (a+b+c)-(d+e)
    return(lpn)
}

I have tested this function with shorter, intuitive numbers, and it works as expected. I use log probabilities to stay away from intractable numbers.

Next, I use function lveech to obtain a vector ‘lpvec’ that gives me all the log probabilities of getting n=0,1,2,…,97 simultaneous observations of the two animals:

lpvec <- rep(NA,J+1)
for(i in 1:(J+1)) lpvec[i] <- lveech(N,N1,N2,nseq[i])

PROBLEM: the sum of the probabilities in lpvec should be <=1, but it is not. The sum is something on the order of 1.48e-13. That happens whether I sum exp(lpvec) or use a function for summing probabilities in log-prob space. In most applications of lveech, the sum of the probabilities is <=1, as it should be, but occasionally I get this problem. Why should this be? Does anyone know of any issues with lchoose that could explain this?

I want to solve this because if I round up the sum to 1 I will not be able to quantify the attraction between animals and compare it with the attractions between other pairs of animals in my data.

Thank you so much for reading the long question.

G.








I have two binary vectors, v1 and v2, both with length N=2178. The sum of vector v1 is N1=165 and the sum of vector v2 is N2=331. I need to compute the probability that 


More information about the R-help mailing list