[Octave-cvsupdate] SF.net SVN: octave:[10366] trunk/octave-forge/main/statistics

[<as...@us...>]

Revision: 10366
          http://octave.svn.sourceforge.net/octave/?rev=10366&view=rev
Author:   asnelt
Date:     2012-05-06 14:32:55 +0000 (Sun, 06 May 2012)
Log Message:
-----------
Relaxed condition for invalid rows in mnpdf and mnrnd so that larger numbers of
categories are supported.

Modified Paths:
--------------
    trunk/octave-forge/main/statistics/NEWS
    trunk/octave-forge/main/statistics/inst/mnpdf.m
    trunk/octave-forge/main/statistics/inst/mnrnd.m

Modified: trunk/octave-forge/main/statistics/NEWS
===================================================================
--- trunk/octave-forge/main/statistics/NEWS	2012-05-05 13:03:58 UTC (rev 10365)
+++ trunk/octave-forge/main/statistics/NEWS	2012-05-06 14:32:55 UTC (rev 10366)
@@ -1,3 +1,9 @@
+Summary of important user-visible changes for statistics 1.1.3:
+-------------------------------------------------------------------
+
+ ** The functions mnpdf and mnrnd are now also usable for greater numbers
+    of categories for which the rows do not exactly sum to 1.
+
 Summary of important user-visible changes for statistics 1.1.2:
 -------------------------------------------------------------------
 

Modified: trunk/octave-forge/main/statistics/inst/mnpdf.m
===================================================================
--- trunk/octave-forge/main/statistics/inst/mnpdf.m	2012-05-05 13:03:58 UTC (rev 10365)
+++ trunk/octave-forge/main/statistics/inst/mnpdf.m	2012-05-06 14:32:55 UTC (rev 10366)
@@ -115,7 +115,9 @@
   t = x .* log (p);
   t(x == 0) = 0;
   y = exp (gammaln (n+1) - sum (gammaln (x+1), 2) + sum (t, 2));
-  y(sum (p, 2) != 1) = NaN;
+  # Set invalid rows to NaN
+  k = (abs (sum (p, 2) - 1) > 1e-6);
+  y(k) = NaN;
 
 endfunction
 

Modified: trunk/octave-forge/main/statistics/inst/mnrnd.m
===================================================================
--- trunk/octave-forge/main/statistics/inst/mnrnd.m	2012-05-05 13:03:58 UTC (rev 10365)
+++ trunk/octave-forge/main/statistics/inst/mnrnd.m	2012-05-06 14:32:55 UTC (rev 10366)
@@ -135,23 +135,24 @@
 
   # Upper bounds of categories
   ub = cumsum (p, 2);
+  # Make sure that the greatest upper bound is 1
+  gub = ub(:, end);
+  ub(:, end) = 1;
   # Lower bounds of categories
   lb = [zeros(sz(1), 1) ub(:, 1:(end-1))];
 
   # Draw multinomial samples
   x = zeros (sz);
   for i = 1:sz(1)
-    # Check if the probabilities sum to 1
-    if (ub(i, end) == 1)
-      # Draw uniform random numbers
-      r = repmat (rand (n(i), 1), 1, sz(2));
-      # Compare the random numbers of r to the cumulated probabilities of p and
-      # count the number of samples for each category
-      x(i, :) =  sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1);
-    else
-      x(i, :) = NaN;
-    endif
+    # Draw uniform random numbers
+    r = repmat (rand (n(i), 1), 1, sz(2));
+    # Compare the random numbers of r to the cumulated probabilities of p and
+    # count the number of samples for each category
+    x(i, :) =  sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1);
   endfor
+  # Set invalid rows to NaN
+  k = (abs (gub - 1) > 1e-6);
+  x(k, :) = NaN;
 
 endfunction
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.