more on gamma

[Michael Creel]

Thanks to a couple of you for suggestions. Two solutions to the gamma overflow 
problem are:

1) to write a wrapper function, as follows:

        function [f] = mc_gamma(z)

        test = z > 170;
        percentage = mean(test);
        if sum(test) > 0
                warning("mc_gamma: %f of the arguments were bounded", 
percentage);
        endif   
        z = 170*test + (1-test).*z;
        f = gamma(z);

        endfunction

This prevents overflows and gives you a warning when it occurs.


2) Run the input (x) to gamma(g(x)) through a nonlinear map that bounds it so 
that eps < g(x) < 170, e.g, 

        alpha = 165/(1 + exp(5 - psi(k+1,:))) + eps; # eps < alpha < 165

Now gamma(alpha) won't crash.

These two solutions together seem to be working for me, and were easy to 
implement.

Thanks, Michael



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------