lgamma->gammaln

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

lgamma->gammaln

Paul Kienzle-2

Paul Kienzle <[hidden email]>

        * miscellaneous/bincoeff.m: Replace lgamma with gammaln.
        * specfun/beta.m: Likewise.
        * special-matrix/invhilb.m: Likewise (but it is only in a comment).
        * statistics/distributions/gamma_pdf.m: Likewise.
        * statistics/distributions/poisson_pdf.m: Likewise.

Index: miscellaneous/bincoeff.m
===================================================================
RCS file: /cvs/octave/scripts/miscellaneous/bincoeff.m,v
retrieving revision 1.5
diff -c -r1.5 bincoeff.m
*** miscellaneous/bincoeff.m 2002/05/01 06:48:35 1.5
--- miscellaneous/bincoeff.m 2002/08/09 17:29:49
***************
*** 83,96 ****
 
    ind = find ((k > 0) & ((n == real (round (n))) & (n < 0)));
    if any (ind)
!     b(ind) = (-1) .^ k(ind) .* exp (lgamma (abs (n(ind)) + k(ind)) ...
!         - lgamma (k(ind) + 1) - lgamma (abs (n(ind))));
    endif
 
    ind = find ((k > 0) & ((n != real (round (n))) | (n >= k)));
    if (length (ind) > 0)
!     b(ind) = exp (lgamma (n(ind) + 1) - lgamma (k(ind) + 1) ...
!         - lgamma (n(ind) - k(ind) + 1));
    endif
 
    ## clean up rounding errors
--- 83,96 ----
 
    ind = find ((k > 0) & ((n == real (round (n))) & (n < 0)));
    if any (ind)
!     b(ind) = (-1) .^ k(ind) .* exp (gammaln (abs (n(ind)) + k(ind)) ...
!         - gammaln (k(ind) + 1) - gammaln (abs (n(ind))));
    endif
 
    ind = find ((k > 0) & ((n != real (round (n))) | (n >= k)));
    if (length (ind) > 0)
!     b(ind) = exp (gammaln (n(ind) + 1) - gammaln (k(ind) + 1) ...
!         - gammaln (n(ind) - k(ind) + 1));
    endif
 
    ## clean up rounding errors
Index: specfun/beta.m
===================================================================
RCS file: /cvs/octave/scripts/specfun/beta.m,v
retrieving revision 1.15
diff -c -r1.15 beta.m
*** specfun/beta.m 2000/01/13 08:40:27 1.15
--- specfun/beta.m 2002/08/09 17:29:49
***************
*** 45,50 ****
      usage ("beta (a, b)");
    endif
 
!   retval = exp (lgamma (a) + lgamma (b) - lgamma (a+b));
 
  endfunction
--- 45,50 ----
      usage ("beta (a, b)");
    endif
 
!   retval = exp (gammaln (a) + gammaln (b) - gammaln (a+b));
 
  endfunction
Index: special-matrix/invhilb.m
===================================================================
RCS file: /cvs/octave/scripts/special-matrix/invhilb.m,v
retrieving revision 1.16
diff -c -r1.16 invhilb.m
*** special-matrix/invhilb.m 2002/04/04 23:03:15 1.16
--- special-matrix/invhilb.m 2002/08/09 17:29:49
***************
*** 85,91 ****
      ## machine number, the result is also exact.  Otherwise we calculate
      ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)).
      ##
!     ## The Octave bincoeff routine uses transcendental functions (lgamma
      ## and exp) rather than multiplications, for the sake of speed.  
      ## However, it rounds the answer to the nearest integer, which
      ## justifies the claim about exactness made above.
--- 85,91 ----
      ## machine number, the result is also exact.  Otherwise we calculate
      ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)).
      ##
!     ## The Octave bincoeff routine uses transcendental functions (gammaln
      ## and exp) rather than multiplications, for the sake of speed.  
      ## However, it rounds the answer to the nearest integer, which
      ## justifies the claim about exactness made above.
Index: statistics/distributions/gamma_pdf.m
===================================================================
RCS file: /cvs/octave/scripts/statistics/distributions/gamma_pdf.m,v
retrieving revision 1.5
diff -c -r1.5 gamma_pdf.m
*** statistics/distributions/gamma_pdf.m 2002/06/27 14:14:08 1.5
--- statistics/distributions/gamma_pdf.m 2002/08/09 17:29:49
***************
*** 59,65 ****
    k = find ((x > 0) & (a > 1) & (b > 0));
    if (any (k))
      pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
!  - b(k) .* x(k) - lgamma (a(k)));
    endif
 
    pdf = reshape (pdf, r, c);
--- 59,65 ----
    k = find ((x > 0) & (a > 1) & (b > 0));
    if (any (k))
      pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
!  - b(k) .* x(k) - gammaln (a(k)));
    endif
 
    pdf = reshape (pdf, r, c);
Index: statistics/distributions/poisson_pdf.m
===================================================================
RCS file: /cvs/octave/scripts/statistics/distributions/poisson_pdf.m,v
retrieving revision 1.4
diff -c -r1.4 poisson_pdf.m
*** statistics/distributions/poisson_pdf.m 2002/05/01 06:48:36 1.4
--- statistics/distributions/poisson_pdf.m 2002/08/09 17:29:49
***************
*** 50,56 ****
 
    k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0));
    if (any (k))
!     pdf(k) = exp (x(k) .* log (l(k)) - l(k) - lgamma (x(k) + 1));
    endif
 
    pdf = reshape (pdf, r, c);
--- 50,56 ----
 
    k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0));
    if (any (k))
!     pdf(k) = exp (x(k) .* log (l(k)) - l(k) - gammaln (x(k) + 1));
    endif
 
    pdf = reshape (pdf, r, c);


Reply | Threaded
Open this post in threaded view
|

lgamma->gammaln

John W. Eaton-6
On  9-Aug-2002, Paul Kienzle <[hidden email]> wrote:

| Paul Kienzle <[hidden email]>
|
| * miscellaneous/bincoeff.m: Replace lgamma with gammaln.
| * specfun/beta.m: Likewise.
| * special-matrix/invhilb.m: Likewise (but it is only in a comment).
| * statistics/distributions/gamma_pdf.m: Likewise.
| * statistics/distributions/poisson_pdf.m: Likewise.

I applied this patch.

Thanks,

jwe