Mercurial > forge

changeset 10180:93615b0e893f octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Relaxed condition for invalid rows in mnpdf and mnrnd so that larger numbers of categories are supported.
author asnelt
date Sun, 06 May 2012 14:32:55 +0000
parents 32ef62b6ff2a
children 7a304f7f9bd5
files main/statistics/NEWS main/statistics/inst/mnpdf.m main/statistics/inst/mnrnd.m
diffstat 3 files changed, 20 insertions(+), 11 deletions(-) [+]
line wrap: on
line diff
--- a/main/statistics/NEWS	Sat May 05 13:03:58 2012 +0000
+++ b/main/statistics/NEWS	Sun May 06 14:32:55 2012 +0000
@@ -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:
 -------------------------------------------------------------------
 
--- a/main/statistics/inst/mnpdf.m	Sat May 05 13:03:58 2012 +0000
+++ b/main/statistics/inst/mnpdf.m	Sun May 06 14:32:55 2012 +0000
@@ -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
 
--- a/main/statistics/inst/mnrnd.m	Sat May 05 13:03:58 2012 +0000
+++ b/main/statistics/inst/mnrnd.m	Sun May 06 14:32:55 2012 +0000
@@ -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