# HG changeset patch
# User jwe
# Date 1109827617 0
# Node ID 53f80b6d98d33aac8054b26bb5549c4f4ea851d8
# Parent  6758c11b5b9920daa8cd2565df955d20e31f6e87
[project @ 2005-03-03 05:25:23 by jwe]

diff -r 6758c11b5b99 -r 53f80b6d98d3 scripts/statistics/distributions/binomial_pdf.m
--- a/scripts/statistics/distributions/binomial_pdf.m	Thu Mar 03 05:18:04 2005 +0000
+++ b/scripts/statistics/distributions/binomial_pdf.m	Thu Mar 03 05:26:57 2005 +0000
@@ -33,31 +33,29 @@
     usage ("binomial_pdf (x, n, p)");
   endif
 
-  if (!isscalar (n) || !isscalar (p))
+  if (! isscalar (n) || ! isscalar (p))
     [retval, x, n, p] = common_size (x, n, p);
     if (retval > 0)
       error ("binomial_pdf: x, n and p must be of common size or scalar");
     endif
   endif
 
-  sz = size (x);
-  pdf = zeros (sz);
-
-  k = find (isnan (x) | !(n >= 0) | (n != round (n)) | !(p >= 0) | !(p <= 1));
-  if (any (k))
-    pdf(k) = NaN;
-  endif
+  k = ((x >= 0) & (x <= n)
+       & (x == round (x)) & (n == round (n))
+       & (p >= 0) & (p <= 1));
 
-  k = find ((x >= 0) & (x <= n) & (x == round (x))
-	    & (n == round (n)) & (p >= 0) & (p <= 1));
-  if (any (k))
-    if (isscalar (n) && isscalar (p))
-      pdf(k) = (bincoeff (n, x(k)) .* (p .^ x(k))
-		.* ((1 - p) .^ (n - x(k))));
-    else
-      pdf(k) = (bincoeff (n(k), x(k)) .* (p(k) .^ x(k))
-		.* ((1 - p(k)) .^ (n(k) - x(k))));
+  pdf = zeros (size (x));
+  pdf(! k) = NaN;
+  if (any (k(:)))
+    x = x(k);
+    if (! isscalar (n))
+      n = n(k);
     endif
+    if (! isscalar (p))
+      p = p(k);
+    endif
+    z = gammaln(n+1) - gammaln(x+1) - gammaln(n-x+1) + x.*log(p) + (n-x).*log(1-p);
+    pdf(k) = exp (z);
   endif
 
 endfunction