Factorial, tgamma & matrix building

Modeling
#1

Dear stan community, I am facing the following programming problem writing my own log probability mass function for the multivariate Fisher & Wallenius’ non central hypergeometric distributions.

Given two positive integers n and c>1, I would like to build a matrix of integers with c columns and
\frac{(n+c-1)!}{n!(c-1)!} rows.

To calculate the factorials i use the tgamma() function, so that \Gamma(n+1) = n!.
However tgamma returns real values and I cannot pass real values a s dimension of a matrix.
Any workaround?

Thanks in advance.

#2

Stan currently doesn’t have a function for converting reals to ints, but you can write one pretty simply:

int real_to_int(real x) {
    int i = 1;
    while(i < x) {
        i += 1;
    }
    return i;
}

Note that this is a pretty basic implementation and will always round up:

real_to_int(25.0) = 25
real_to_int(25.1) = 26

So keep that in mind if you end up working with reals that have decimals

2 Likes
Interpolation needs integer conversion
#3

Thanks @andrjohns. It works like a charm!

I thought, that another way would be me writing a function like

int factorial(n) {
... // recursive stuff
}

Do you see any drawback? Could the factorial function be faster than real_to_int in this case where n! grows really fast and so the loop could be time-consuming? Or am I saving time using tgamma(n+1) over my factorial()?
These of mine are just free thoughts: “premature optimization is the root of all evil”… ;-)

1 Like
#4

Yes, I’d say that implementing factorial (possibly with a loop) is a cleaner solution in this case. It might even be sensible for Stan to expose such a function. For maximum robustness (e.g. avoid overflow when one of n and c is large but the other small), you might want to split the computation into falling_factorial( n + c - 1, n + c - 1 - max(n, c - 1)) / factorial(min(n, c - 1)) as that will avoid some large intermediate values but maybe that’s even prematurer optimization :-D

1 Like
#5

@martinmodrak
Thanks for the reply:
here’s my quick&dirty hack for everyone in need for a “pack-animal-function” for int factorial(n).

/**
* Return the factorial of the non negative number n: the product of all 
* positive integers less than or equal to n, by using recursion.
*
* @param n number for wich the factorial is calculated.
*
* @return the factorial of n, $n!$.
*/    
int factorial(int n) {
   if (n >= 1)
       return n * factorial(n-1);
   else
       return 1;
}

I will implement all the lpmfs, then I will go back and try to do better with the factorial function and then share a decent function to all the stan community!

3 Likes