# HG changeset patch
# User jwe
# Date 1081305486 0
# Node ID 4b0f3b055331379dc6d3d5a1d7826e22cdb074cc
# Parent  66b3cce2bf37f8980c3e2e9a70118d5f3c1073be
[project @ 2004-04-07 02:37:05 by jwe]

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,5 +1,39 @@
 2004-04-06  David Bateman  <dbateman@free.fr>
 
+	* general/common_size.m, miscellaneous/bincoeff.m,
+	statistics/distributions/beta_cdf.m,
+	statistics/distributions/beta_inv.m,
+	statistics/distributions/beta_pdf.m,
+	statistics/distributions/beta_rnd.m,
+	statistics/distributions/binomial_cdf.m,
+	statistics/distributions/binomial_inv.m,
+	statistics/distributions/binomial_pdf.m,
+	statistics/distributions/binomial_rnd.m,
+	statistics/distributions/cauchy_cdf.m,
+	statistics/distributions/cauchy_inv.m,
+	statistics/distributions/cauchy_pdf.m,
+	statistics/distributions/cauchy_rnd.m,
+	statistics/distributions/chisquare_cdf.m,
+	statistics/distributions/chisquare_inv.m,
+	statistics/distributions/chisquare_pdf.m,
+	statistics/distributions/chisquare_rnd.m,
+	statistics/distributions/gamma_cdf.m,
+	statistics/distributions/gamma_inv.m,
+	statistics/distributions/gamma_pdf.m,
+	statistics/distributions/gamma_rnd.m,
+	statistics/distributions/normal_cdf.m,
+	statistics/distributions/normal_inv.m,
+	statistics/distributions/normal_pdf.m,
+	statistics/distributions/normal_rnd.m,
+	statistics/distributions/stdnormal_cdf.m,
+	statistics/distributions/stdnormal_pdf.m,
+	statistics/distributions/stdnormal_rnd.m,
+	statistics/distributions/uniform_cdf.m,
+	statistics/distributions/uniform_inv.m,
+	statistics/distributions/uniform_pdf.m,
+	statistics/distributions/uniform_rnd.m:
+	Allow the inputs to be N-d arrays.
+
 	* statistics/base/var.m: Update for N-d arrays.  Allow dimension arg.
 	* statistics/base/median.m: Likewise.
 
diff --git a/scripts/general/common_size.m b/scripts/general/common_size.m
--- a/scripts/general/common_size.m
+++ b/scripts/general/common_size.m
@@ -49,8 +49,20 @@
     error ("common_size: only makes sense if nargin >= 2");
   endif
 
+  len = 2;
   for i = 1 : nargin
-    s(i,:) = size (varargin{i});
+    sz =  size (varargin{i});
+    if (length (sz) < len)
+      s(i,:) = [sz, ones(1,len - length(sz))];
+    else
+      if (length (sz) > len)
+	if (i > 1)
+	  s = [s, ones(size(s,1), length(sz) - len)];
+	endif
+	len = length (sz);
+      endif
+      s(i,:) = sz;
+    endif
   endfor
 
   m = max (s);
diff --git a/scripts/miscellaneous/bincoeff.m b/scripts/miscellaneous/bincoeff.m
--- a/scripts/miscellaneous/bincoeff.m
+++ b/scripts/miscellaneous/bincoeff.m
@@ -65,20 +65,17 @@
     error ("bincoeff: n and k must be of common size or scalars");
   endif
 
-  [r, c] = size (n);
-  s = r * c;
-  n   = reshape (n, s, 1);
-  k   = reshape (k, s, 1);
-  b   = zeros (s, 1);
+  sz = size (n);
+  b   = zeros (sz);
 
   ind = find (! (k >= 0) | (k != real (round (k))) | isnan (n));
   if (any (ind))
-    b(ind) = NaN * ones (length (ind), 1);
+    b(ind) = NaN;
   endif
 
   ind = find (k == 0);
   if (any (ind))
-    b(ind) = ones (length (ind), 1);
+    b(ind) = 1;
   endif
 
   ind = find ((k > 0) & ((n == real (round (n))) & (n < 0)));
@@ -99,7 +96,5 @@
     b(ind) = round (b(ind));
   endif
 
-  b = reshape (b, r, c);
-
 endfunction
 
diff --git a/scripts/statistics/distributions/beta_cdf.m b/scripts/statistics/distributions/beta_cdf.m
--- a/scripts/statistics/distributions/beta_cdf.m
+++ b/scripts/statistics/distributions/beta_cdf.m
@@ -33,33 +33,33 @@
     usage ("beta_cdf (a, b, x)");
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("beta_cdf: x, a and b must be of common size or scalar");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("beta_cdf: x, a and b must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, s, 1);
-  a   = reshape (a, s, 1);
-  b   = reshape (b, s, 1);
-  cdf = zeros (s, 1);
+  sz = size(x);
+  cdf = zeros (sz);
 
   k = find (!(a > 0) | !(b > 0) | isnan (x));
   if (any (k))
-    cdf (k) = NaN * ones (length (k), 1);
+    cdf (k) = NaN;
   endif
 
   k = find ((x >= 1) & (a > 0) & (b > 0));
   if (any (k))
-    cdf (k) = ones (length (k), 1);
+    cdf (k) = 1;
   endif
 
   k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
   if (any (k))
-    cdf (k) = betainc (x(k), a(k), b(k));
+    if (isscalar (a) && isscalar(b))
+      cdf (k) = betainc (x(k), a, b);
+    else
+      cdf (k) = betainc (x(k), a(k), b(k));
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/beta_inv.m b/scripts/statistics/distributions/beta_inv.m
--- a/scripts/statistics/distributions/beta_inv.m
+++ b/scripts/statistics/distributions/beta_inv.m
@@ -33,34 +33,36 @@
     usage ("beta_inv (x, a, b)");
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("beta_inv: x, a and b must be of common size or scalars");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("beta_inv: x, a and b must be of common size or scalars");
+    endif
   endif
-
-  [r, c] = size (x);
-  s   = r * c;
-  x   = reshape (x, s, 1);
-  a   = reshape (a, s, 1);
-  b   = reshape (b, s, 1);
-  inv = zeros (s, 1);
+  
+  sz = size (x);
+  inv = zeros (sz);
 
   k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x));
   if (any (k))
-    inv (k) = NaN * ones (length (k), 1);
+    inv (k) = NaN;
   endif
 
   k = find ((x == 1) & (a > 0) & (b > 0));
   if (any (k))
-    inv (k) = ones (length (k), 1);
+    inv (k) = 1;
   endif
 
   k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
   if (any (k))
-    a = a (k);
-    b = b (k);
+    if (!isscalar(a) || !isscalar(b))
+      a = a (k);
+      b = b (k);
+      y = a ./ (a + b);
+    else
+      y = a / (a + b) * ones (size (k));
+    endif
     x = x (k);
-    y = a ./ (a + b);
     l = find (y < eps);
     if (any (l))
       y(l) = sqrt (eps) * ones (length (l), 1);
@@ -92,6 +94,4 @@
     inv (k) = y_new;
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/beta_pdf.m b/scripts/statistics/distributions/beta_pdf.m
--- a/scripts/statistics/distributions/beta_pdf.m
+++ b/scripts/statistics/distributions/beta_pdf.m
@@ -31,30 +31,31 @@
   if (nargin != 3)
     usage ("beta_pdf (a, b, x)");
   endif
-
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("beta_pdf: x, a and b must be of common size or scalar");
+  
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("beta_pdf: x, a and b must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, s, 1);
-  a   = reshape (a, s, 1);
-  b   = reshape (b, s, 1);
-  pdf = zeros (s, 1);
+  sz = size (x);
+  pdf = zeros (sz);
 
   k = find (!(a > 0) | !(b > 0) | isnan (x));
   if (any (k))
-    pdf (k) = NaN * ones (length (k), 1);
+    pdf (k) = NaN;
   endif
 
   k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
   if (any (k))
-    pdf(k) = exp ((a(k) - 1) .* log (x(k))
-		  + (b(k) - 1) .* log (1 - x(k))) ./ beta (a(k), b(k));
+    if (isscalar(a) && isscalar(b))
+      pdf(k) = exp ((a - 1) .* log (x(k))
+		    + (b - 1) .* log (1 - x(k))) ./ beta (a, b);
+    else
+      pdf(k) = exp ((a(k) - 1) .* log (x(k))
+		    + (b(k) - 1) .* log (1 - x(k))) ./ beta (a(k), b(k));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/beta_rnd.m b/scripts/statistics/distributions/beta_rnd.m
--- a/scripts/statistics/distributions/beta_rnd.m
+++ b/scripts/statistics/distributions/beta_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} beta_rnd (@var{a}, @var{b}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the Beta
-## distribution with parameters @var{a} and @var{b}.  Both @var{a} and
-## @var{b} must be scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} beta_rnd (@var{a}, @var{b}, @var{sz})
+## Return an @var{r} by @var{c} or @code{size (@var{sz})} matrix of 
+## random samples from the Beta distribution with parameters @var{a} and
+## @var{b}.  Both @var{a} and @var{b} must be scalar or of size @var{r}
+##  by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{a} and @var{b}.
@@ -32,6 +34,15 @@
 
 function rnd = beta_rnd (a, b, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(a) || !isscalar(b)) 
+      [retval, a, b] = common_size (a, b);
+      if (retval > 0)
+	error ("beta_rnd: a and b must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("beta_rnd: r must be a positive integer");
@@ -39,35 +50,39 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("beta_rnd: c must be a positive integer");
     endif
-    [retval, a, b] = common_size (a, b, zeros (r, c));
-    if (retval > 0)
-      error ("beta_rnd: a and b must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("beta_rnd: r must be a postive integer or vector");
     endif
   elseif (nargin == 2)
-    [retval, a, b] = common_size (a, b);
-    if (retval > 0)
-      error ("beta_rnd: a and b must be of common size or scalar");
-    endif
+    sz = size(a);
   else
     usage ("beta_rnd (a, b, r, c)");
   endif
 
-  [r, c] = size (a);
-  s = r * c;
-  a = reshape (a, 1, s);
-  b = reshape (b, 1, s);
-  rnd = zeros (1, s);
+  if (isscalar(a) && isscalar(b))
+    if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
+      rnd = NaN * ones (sz);
+    else
+      rnd = beta_inv (rand(sz), a, b);
+    endif
+  else
+    rnd = zeros (sz);
 
-  k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
+    if (any (k))
+      rnd(k) = NaN * ones (size (k));
+    endif
+
+    k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
+    if (any (k))
+      rnd(k) = beta_inv (rand (size (k)), a(k), b(k));
+    endif
   endif
 
-  k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
-  if (any (k))
-    rnd(k) = beta_inv (rand (1, length (k)), a(k), b(k));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/binomial_cdf.m b/scripts/statistics/distributions/binomial_cdf.m
--- a/scripts/statistics/distributions/binomial_cdf.m
+++ b/scripts/statistics/distributions/binomial_cdf.m
@@ -32,37 +32,37 @@
     usage ("binomial_cdf (x, n, p)");
   endif
 
-  [retval, x, n, p] = common_size (x, n, p);
-  if (retval > 0)
-    error ("binomial_cdf: x, n and p must be of common size or scalar");
+  if (!isscalar (n) || !isscalar (p))
+    [retval, x, n, p] = common_size (x, n, p);
+    if (retval > 0)
+      error ("binomial_cdf: x, n and p must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, 1, s);
-  n   = reshape (n, 1, s);
-  p   = reshape (p, 1, s);
-  cdf = zeros (1, s);
+  sz = size (x);
+  cdf = zeros (sz);
 
   k = find (isnan (x) | !(n >= 0) | (n != round (n))
 	    | !(p >= 0) | !(p <= 1));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x >= n) & (n >= 0) & (n == round (n))
 	    & (p >= 0) & (p <= 1));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x >= 0) & (x < n) & (n == round (n))
 	    & (p >= 0) & (p <= 1));
   if (any (k))
     tmp = floor (x(k));
-    cdf(k) = 1 - betainc (p(k), tmp + 1, n(k) - tmp);
+    if (isscalar (n) && isscalar (p))
+      cdf(k) = 1 - betainc (p, tmp + 1, n - tmp);
+    else
+      cdf(k) = 1 - betainc (p(k), tmp + 1, n(k) - tmp);
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/binomial_inv.m b/scripts/statistics/distributions/binomial_inv.m
--- a/scripts/statistics/distributions/binomial_inv.m
+++ b/scripts/statistics/distributions/binomial_inv.m
@@ -32,39 +32,48 @@
     usage ("binomial_inv (x, n, p)");
   endif
 
-  [retval, x, n, p] = common_size (x, n, p);
-  if (retval > 0)
-    error ("binomial_inv: x, n and p must be of common size or scalars");
+  if (!isscalar (n) || !isscalar (p))
+    [retval, x, n, p] = common_size (x, n, p);
+    if (retval > 0)
+      error ("binomial_inv: x, n and p must be of common size or scalars");
+    endif
   endif
-
-  [r, c] = size (x);
-  s   = r * c;
-  x   = reshape (x, 1, s);
-  n   = reshape (n, 1, s);
-  p   = reshape (p, 1, s);
-  inv = zeros (1, s);
+  
+  sz = size (x);
+  inv = zeros (sz);
 
   k = find (!(x >= 0) | !(x <= 1) | !(n >= 0) | (n != round (n))
 	    | !(p >= 0) | !(p <= 1));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x >= 0) & (x <= 1) & (n >= 0) & (n == round (n))
 	    & (p >= 0) & (p <= 1));
   if (any (k))
-    cdf = binomial_pdf (0, n(k), p(k));
-    while (any (inv(k) < n(k)))
-      m = find (cdf < x(k));
-      if (any (m))
-        inv(k(m)) = inv(k(m)) + 1;
-        cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n(k(m)), p(k(m)));
-      else
-        break;
-      endif
-    endwhile
+    if (isscalar (n) && isscalar (p))
+      cdf = binomial_pdf (0, n, p) * ones (size(k));
+      while (any (inv(k) < n))
+	m = find (cdf < x(k));
+	if (any (m))
+          inv(k(m)) = inv(k(m)) + 1;
+          cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n, p);
+	else
+          break;
+	endif
+      endwhile
+    else 
+      cdf = binomial_pdf (0, n(k), p(k));
+      while (any (inv(k) < n(k)))
+	m = find (cdf < x(k));
+	if (any (m))
+          inv(k(m)) = inv(k(m)) + 1;
+          cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n(k(m)), p(k(m)));
+	else
+          break;
+	endif
+      endwhile
+    endif
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/binomial_pdf.m b/scripts/statistics/distributions/binomial_pdf.m
--- a/scripts/statistics/distributions/binomial_pdf.m
+++ b/scripts/statistics/distributions/binomial_pdf.m
@@ -33,30 +33,31 @@
     usage ("binomial_pdf (x, n, p)");
   endif
 
-  [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");
+  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
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, 1, s);
-  n   = reshape (n, 1, s);
-  p   = reshape (p, 1, s);
-  cdf = zeros (1, s);
+  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 * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find ((x >= 0) & (x <= n) & (x == round (x))
 	    & (n == round (n)) & (p >= 0) & (p <= 1));
   if (any (k))
-    pdf(k) = (bincoeff (n(k), x(k)) .* (p(k) .^ x(k))
-	      .* ((1 - p(k)) .^ (n(k) - x(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))));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/binomial_rnd.m b/scripts/statistics/distributions/binomial_rnd.m
--- a/scripts/statistics/distributions/binomial_rnd.m
+++ b/scripts/statistics/distributions/binomial_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} binomial_rnd (@var{n}, @var{p}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the
-## binomial distribution with parameters @var{n} and @var{p}.  Both
-## @var{n} and @var{p} must be scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} binomial_rnd (@var{n}, @var{p}, @var{sz})
+## Return an @var{r} by @var{c}  or a @code{size (@var{sz})} matrix of 
+## random samples from the binomial distribution with parameters @var{n}
+## and @var{p}.  Both @var{n} and @var{p} must be scalar or of size
+## @var{r} by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{n} and @var{p}.
@@ -32,6 +34,15 @@
 
 function rnd = binomial_rnd (n, p, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(n) || !isscalar(p)) 
+      [retval, n, p] = common_size (n, p);
+      if (retval > 0)
+	error ("binomial_rnd: n and p must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("binomial_rnd: r must be a positive integer");
@@ -39,41 +50,51 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("binomial_rnd: c must be a positive integer");
     endif
-    [retval, n, p] = common_size (n, p, zeros (r, c));
-    if (retval > 0)
-      error ("binomial_rnd: n and p must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("binomial_rnd: r must be a postive integer or vector");
     endif
   elseif (nargin == 2)
-    [retval, n, p] = common_size (n, p);
-    if (retval > 0)
-      error ("binomial_rnd: n and p must be of common size or scalar");
-    endif
+    sz = size(n);
   else
     usage ("binomial_rnd (n, p, r, c)");
   endif
 
-  [r, c] = size (n);
-  s = r * c;
-  n = reshape (n, 1, s);
-  p = reshape (p, 1, s);
-  rnd = zeros (1, s);
+  if (isscalar (n) && isscalar (p))
+    if (find (!(n > 0) | !(n < Inf) | !(n == round (n)) |
+              !(p >= 0) | !(p <= 1)))
+      rnd = NaN * ones (sz);
+    else
+      nel = prod (sz);
+      tmp = rand (n, nel);
+      ind = (1 : n)' * ones (1, nel);
+      rnd = sum ((tmp < ones (n, nel) * p) &
+                    (ind <= ones (n, nel) * n));
+      rnd = reshape(rnd, sz);
+    endif
+  else
+    rnd = zeros (sz);
 
-  k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) |
-            !(p <= 0) | !(p >= 1));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) |
+              !(p >= 0) | !(p <= 1));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+
+    k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1));
+    if (any (k))
+      N = max (n(k));
+      L = length (k);
+      tmp = rand (N, L);
+      ind = (1 : N)' * ones (1, L);
+      rnd(k) = sum ((tmp < ones (N, 1) * p(k)(:)') &
+                    (ind <= ones (N, 1) * n(k)(:)'));
+    endif
   endif
 
-  k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1));
-  if (any (k))
-    N = max (n(k));
-    L = length (k);
-    tmp = rand (N, L);
-    ind = (1 : N)' * ones (1, L);
-    rnd(k) = sum ((tmp < ones (N, 1) * p(k)) &
-                  (ind <= ones (N, 1) * n(k)));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/cauchy_cdf.m b/scripts/statistics/distributions/cauchy_cdf.m
--- a/scripts/statistics/distributions/cauchy_cdf.m
+++ b/scripts/statistics/distributions/cauchy_cdf.m
@@ -39,24 +39,24 @@
     scale = 1;
   endif
 
-  [retval, x, location, scale] = common_size (x, location, scale);
-  if (retval > 0)
-    error ("cauchy_cdf: x, lambda and sigma must be of common size or scalar");
+  if (!isscalar (location) || !isscalar (scale)) 
+    [retval, x, location, scale] = common_size (x, location, scale);
+    if (retval > 0)
+      error ("cauchy_cdf: x, lambda and sigma must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  location = reshape (location, 1, s);
-  scale = reshape (scale, 1, s);
-  cdf = NaN * ones (1, s);
+  sz = size (x);
+  cdf = NaN * ones (sz);
 
   k = find ((x > -Inf) & (x < Inf) & (location > -Inf) &
             (location < Inf) & (scale > 0) & (scale < Inf));
   if (any (k))
-    cdf(k) = 0.5 + atan ((x(k) - location(k)) ./ scale(k)) / pi;
+    if (isscalar (location) && isscalar (scale)) 
+      cdf(k) = 0.5 + atan ((x(k) - location) ./ scale) / pi;
+    else
+      cdf(k) = 0.5 + atan ((x(k) - location(k)) ./ scale(k)) / pi;
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/cauchy_inv.m b/scripts/statistics/distributions/cauchy_inv.m
--- a/scripts/statistics/distributions/cauchy_inv.m
+++ b/scripts/statistics/distributions/cauchy_inv.m
@@ -39,37 +39,36 @@
     scale = 1;
   endif
 
-  [retval, x, location, scale] = common_size (x, location, scale);
-  if (retval > 0)
-    error ("cauchy_inv: x, lambda and sigma must be of common size or scalar");
+  if (!isscalar (location) || !isscalar (scale)) 
+    [retval, x, location, scale] = common_size (x, location, scale);
+    if (retval > 0)
+      error ("cauchy_inv: x, lambda and sigma must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  location = reshape (location, 1, s);
-  scale = reshape (scale, 1, s);
-
-  inv = NaN * ones (1, s);
+  sz = size (x);
+  inv = NaN * ones (sz);
 
   ok = ((location > -Inf) & (location < Inf) &
        (scale > 0) & (scale < Inf));
 
   k = find ((x == 0) & ok);
   if (any (k))
-    inv(k) = -Inf * ones (1, length (k));
+    inv(k) = -Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & ok);
   if (any (k))
-    inv(k) = location(k) - scale(k) .* cot (pi * x(k));
+    if (isscalar (location) && isscalar (scale)) 
+      inv(k) = location - scale .* cot (pi * x(k));
+    else
+      inv(k) = location(k) - scale(k) .* cot (pi * x(k));
+    endif
   endif
 
   k = find ((x == 1) & ok);
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/cauchy_pdf.m b/scripts/statistics/distributions/cauchy_pdf.m
--- a/scripts/statistics/distributions/cauchy_pdf.m
+++ b/scripts/statistics/distributions/cauchy_pdf.m
@@ -39,26 +39,26 @@
     scale = 1;
   endif
 
-  [retval, x, location, scale] = common_size (x, location, scale);
-  if (retval > 0)
-    error ("cauchy_pdf: x, lambda and sigma must be of common size or scalar");
+  if (!isscalar (location) || !isscalar (scale)) 
+    [retval, x, location, scale] = common_size (x, location, scale);
+    if (retval > 0)
+      error ("cauchy_pdf: x, lambda and sigma must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  location = reshape (location, 1, s);
-  scale = reshape (scale, 1, s);
-
-  pdf = NaN * ones (1, s);
+  sz = size (x);
+  pdf = NaN * ones (sz);
 
   k = find ((x > -Inf) & (x < Inf) & (location > -Inf) &
             (location < Inf) & (scale > 0) & (scale < Inf));
   if (any (k))
-    pdf(k) = ((1 ./ (1 + ((x(k) - location(k)) ./ scale(k)) .^ 2))
-	      / pi ./ scale(k));
+    if (isscalar (location) && isscalar (scale)) 
+      pdf(k) = ((1 ./ (1 + ((x(k) - location) ./ scale) .^ 2))
+		/ pi ./ scale);
+    else
+      pdf(k) = ((1 ./ (1 + ((x(k) - location(k)) ./ scale(k)) .^ 2))
+		/ pi ./ scale(k));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/cauchy_rnd.m b/scripts/statistics/distributions/cauchy_rnd.m
--- a/scripts/statistics/distributions/cauchy_rnd.m
+++ b/scripts/statistics/distributions/cauchy_rnd.m
@@ -19,9 +19,10 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} cauchy_rnd (@var{lambda}, @var{sigma}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the Cauchy
-## distribution with parameters @var{lambda} and @var{sigma} which must
-## both be scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} cauchy_rnd (@var{lambda}, @var{sigma}, @var{sz})
+## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of 
+## random samples from the Cauchy distribution with parameters @var{lambda} 
+## and @var{sigma} which must both be scalar or of size @var{r} by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{lambda} and @var{sigma}.
@@ -32,6 +33,15 @@
 
 function rnd = cauchy_rnd (l, scale, r, c)
 
+  if (nargin > 1)
+    if (!isscalar (l) || !isscalar (scale)) 
+      [retval, l, scale] = common_size (l, scale);
+      if (retval > 0)
+	error ("cauchy_rnd: lambda and sigma must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("cauchy_rnd: r must be a positive integer");
@@ -39,32 +49,33 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("cauchy_rnd: c must be a positive integer");
     endif
-    [retval, l, scale] = common_size (l, scale, zeros (r, c));
-    if (retval > 0)
-      error ("cauchy_rnd: lambda and sigma must be scalar or of size %d by %d",
-	     r, c); 
+    sz = [r, c];
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("cauchy_rnd: r must be a postive integer or vector");
     endif
   elseif (nargin == 2)
-    [retval, l, scale] = common_size (l, scale);
-    if (retval > 0)
-      error ("cauchy_rnd: lambda and sigma must be of common size or scalar");
-    endif
-    [r, c] = size (l);
+    sz = size(l);
   else
     usage ("cauchy_rnd (lambda, sigma, r, c)");
   endif
 
-  s = r * c;
-  l = reshape (l, 1, s);
-  scale = reshape (scale, 1, s);
-
-  rnd = NaN * ones (1, s);
-
-  k = find ((l > -Inf) & (l < Inf) & (scale > 0) & (scale < Inf));
-  if (any (k))
-    rnd(k) = l(k) - cot (pi * rand (1, length (k))) .* scale(k);
+  if (isscalar (l) && isscalar (scale)) 
+    if (find (!(l > -Inf) | !(l < Inf) | !(scale > 0) | !(scale < Inf)))
+      rnd = NaN * ones (sz);
+    else
+      rnd = l - cot (pi * rand (sz)) .* scale;
+    endif
+  else
+    rnd = NaN * ones (sz);
+    k = find ((l > -Inf) & (l < Inf) & (scale > 0) & (scale < Inf));
+    if (any (k))
+      rnd(k) = l(k)(:) - cot (pi * rand (size (k))) .* scale(k)(:);
+    endif
   endif
 
-  rnd = reshape (rnd, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/chisquare_cdf.m b/scripts/statistics/distributions/chisquare_cdf.m
--- a/scripts/statistics/distributions/chisquare_cdf.m
+++ b/scripts/statistics/distributions/chisquare_cdf.m
@@ -33,9 +33,11 @@
     usage ("chisquare_cdf (x, n)");
   endif
 
-  [retval, x, n] = common_size (x, n);
-  if (retval > 0)
-    error ("chisquare_cdf: x and n must be of common size or scalar");
+  if (!isscalar (n))
+    [retval, x, n] = common_size (x, n);
+    if (retval > 0)
+      error ("chisquare_cdf: x and n must be of common size or scalar");
+    endif
   endif
 
   cdf = gamma_cdf (x, n / 2, 1 / 2);
diff --git a/scripts/statistics/distributions/chisquare_inv.m b/scripts/statistics/distributions/chisquare_inv.m
--- a/scripts/statistics/distributions/chisquare_inv.m
+++ b/scripts/statistics/distributions/chisquare_inv.m
@@ -33,9 +33,11 @@
     usage ("chisquare_inv (x, n)");
   endif
 
-  [retval, x, n] = common_size (x, n);
-  if (retval > 0)
-    error ("chisquare_inv: x and n must be of common size or scalar");
+  if (!isscalar (n))
+    [retval, x, n] = common_size (x, n);
+    if (retval > 0)
+      error ("chisquare_inv: x and n must be of common size or scalar");
+    endif
   endif
 
   inv = gamma_inv (x, n / 2, 1 / 2);
diff --git a/scripts/statistics/distributions/chisquare_pdf.m b/scripts/statistics/distributions/chisquare_pdf.m
--- a/scripts/statistics/distributions/chisquare_pdf.m
+++ b/scripts/statistics/distributions/chisquare_pdf.m
@@ -33,9 +33,11 @@
     usage ("chisquare_pdf (x, n)");
   endif
 
-  [retval, x, n] = common_size (x, n);
-  if (retval > 0)
-    error ("chisquare_pdf: x and n must be of common size or scalar");
+  if (!isscalar (n))
+    [retval, x, n] = common_size (x, n);
+    if (retval > 0)
+      error ("chisquare_pdf: x and n must be of common size or scalar");
+    endif
   endif
 
   pdf = gamma_pdf (x, n / 2, 1 / 2);
diff --git a/scripts/statistics/distributions/chisquare_rnd.m b/scripts/statistics/distributions/chisquare_rnd.m
--- a/scripts/statistics/distributions/chisquare_rnd.m
+++ b/scripts/statistics/distributions/chisquare_rnd.m
@@ -19,9 +19,10 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} chisquare_rnd (@var{n}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the
-## chisquare distribution with @var{n} degrees of freedom.  @var{n} must
-## be a scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} chisquare_rnd (@var{n}, @var{sz})
+## Return an @var{r} by @var{c}  or a @code{size (@var{sz})} matrix of 
+## random samples from the chisquare distribution with @var{n} degrees 
+## of freedom.  @var{n} must be a scalar or of size @var{r} by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the size of @var{n}.
@@ -39,29 +40,43 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("chisquare_rnd: c must be a positive integer");
     endif
-    [retval, n] = common_size (n, zeros (r, c));
-    if (retval > 0)
-      error ("chisquare_rnd: n must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+  elseif (nargin == 2)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("chisquare_rnd: r must be a postive integer or vector");
     endif
-  elseif (nargin != 1)
+  elseif (nargin == 1)
+    sz = size(n);
+  else
     usage ("chisquare_rnd (n, r, c)");
   endif
 
-  [r, c] = size (n);
-  s = r * c;
-  n = reshape (n, 1, s);
-  rnd = zeros (1, s);
+  if (isscalar (n))
+     if (find (!(n > 0) | !(n < Inf)))
+       rnd = NaN * ones (sz);
+     else
+       rnd =  chisquare_inv (rand (sz), n);
+     endif
+  else
+    [retval, n, dummy] = common_size (n, ones (sz));
+    if (retval > 0)
+      error ("chisquare_rnd: a and b must be of common size or scalar");
+    endif
 
-  k = find (!(n > 0) | !(n < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    rnd = zeros (sz);
+    k = find (!(n > 0) | !(n < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+
+    k = find ((n > 0) & (n < Inf));
+    if (any (k))
+      rnd(k) = chisquare_inv (rand (size (k)), n(k));
+    endif
   endif
 
-  k = find ((n > 0) & (n < Inf));
-  if (any (k))
-    rnd(k) = chisquare_inv (rand (1, length (k)), n(k));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/gamma_cdf.m b/scripts/statistics/distributions/gamma_cdf.m
--- a/scripts/statistics/distributions/gamma_cdf.m
+++ b/scripts/statistics/distributions/gamma_cdf.m
@@ -33,28 +33,28 @@
     usage ("gamma_cdf (x, a, b)");
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("gamma_cdf: x, a and b must be of common size or scalars");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("gamma_cdf: x, a and b must be of common size or scalars");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, s, 1);
-  a   = reshape (a, s, 1);
-  b   = reshape (b, s, 1);
-  cdf = zeros (s, 1);
+  sz = size (x);
+  cdf = zeros (sz);
 
   k = find (!(a > 0) | !(b > 0) | isnan (x));
   if (any (k))
-    cdf (k) = NaN * ones (length (k), 1);
+    cdf (k) = NaN;
   endif
 
   k = find ((x > 0) & (a > 0) & (b > 0));
   if (any (k))
-    cdf (k) = gammainc (b(k) .* x(k), a(k));
+    if (isscalar (a) && isscalar(b))
+      cdf (k) = gammainc (b * x(k), a);
+    else
+      cdf (k) = gammainc (b(k) .* x(k), a(k));
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/gamma_inv.m b/scripts/statistics/distributions/gamma_inv.m
--- a/scripts/statistics/distributions/gamma_inv.m
+++ b/scripts/statistics/distributions/gamma_inv.m
@@ -33,34 +33,36 @@
     usage ("gamma_inv (x, a, b)");
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("gamma_inv: x, a and b must be of common size or scalars");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("gamma_inv: x, a and b must be of common size or scalars");
+    endif
   endif
 
-  [r, c] = size (x);
-  s   = r * c;
-  x   = reshape (x, s, 1);
-  a   = reshape (a, s, 1);
-  b   = reshape (b, s, 1);
-  inv = zeros (s, 1);
+  sz = size (x);
+  inv = zeros (sz);
 
   k = find ((x < 0) | (x > 1) | isnan (x) | !(a > 0) | !(b > 0));
   if (any (k))
-    inv (k) = NaN * ones (length (k), 1);
+    inv (k) = NaN;
   endif
 
   k = find ((x == 1) & (a > 0) & (b > 0));
   if (any (k))
-    inv (k) = Inf * ones (length (k), 1);
+    inv (k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
   if (any (k))
-    a = a (k);
-    b = b (k);
+    if (!isscalar(a) || !isscalar(b))
+      a = a (k);
+      b = b (k);
+      y = a ./ b;
+    else
+      y = a / b * ones (size (k));
+    endif
     x = x (k);
-    y = a ./ b;
     l = find (x < eps);
     if (any (l))
       y(l) = sqrt (eps) * ones (length (l), 1);
@@ -84,6 +86,4 @@
     inv (k) = y_new;
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/gamma_pdf.m b/scripts/statistics/distributions/gamma_pdf.m
--- a/scripts/statistics/distributions/gamma_pdf.m
+++ b/scripts/statistics/distributions/gamma_pdf.m
@@ -33,35 +33,41 @@
     usage ("gamma_pdf (x, a, b)");
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("gamma_pdf: x, a and b must be of common size or scalars");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("gamma_pdf: x, a and b must be of common size or scalars");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, s, 1);
-  a   = reshape (a, s, 1);
-  b   = reshape (b, s, 1);
-  pdf = zeros (s, 1);
+  sz = size(x);
+  pdf = zeros (sz);
 
   k = find (!(a > 0) | !(b > 0) | isnan (x));
   if (any (k))
-    pdf (k) = NaN * ones (length (k), 1);
+    pdf (k) = NaN;
   endif
 
   k = find ((x > 0) & (a > 0) & (a <= 1) & (b > 0));
   if (any (k))
-    pdf(k) = ((b(k) .^ a(k)) .* (x(k) .^ (a(k) - 1))
-	      .* exp(-b(k) .* x(k)) ./ gamma (a(k)));
+    if (isscalar(a) && isscalar(b))
+      pdf(k) = ((b .^ a) .* (x(k) .^ (a - 1))
+		.* exp(-b .* x(k)) ./ gamma (a));
+    else
+      pdf(k) = ((b(k) .^ a(k)) .* (x(k) .^ (a(k) - 1))
+		.* exp(-b(k) .* x(k)) ./ gamma (a(k)));
+    endif
   endif
 
   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)));
+    if (isscalar(a) && isscalar(b))
+      pdf(k) = exp (a .* log (b) + (a-1) .* log (x(k))
+		    - b .* x(k) - gammaln (a));
+    else
+      pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
+		    - b(k) .* x(k) - gammaln (a(k)));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/gamma_rnd.m b/scripts/statistics/distributions/gamma_rnd.m
--- a/scripts/statistics/distributions/gamma_rnd.m
+++ b/scripts/statistics/distributions/gamma_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} gamma_rnd (@var{a}, @var{b}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the Gamma
-## distribution with parameters @var{a} and @var{b}.  Both @var{a} and
-## @var{b} must be scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} gamma_rnd (@var{a}, @var{b}, @var{sz})
+## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of 
+## random samples from the Gamma distribution with parameters @var{a}
+## and @var{b}.  Both @var{a} and @var{b} must be scalar or of size 
+## @var{r} by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{a} and @var{b}.
@@ -32,6 +34,15 @@
 
 function rnd = gamma_rnd (a, b, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(a) || !isscalar(b)) 
+      [retval, a, b] = common_size (a, b);
+      if (retval > 0)
+	error ("gamma_rnd: a and b must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("gamma_rnd: r must be a positive integer");
@@ -39,34 +50,38 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("gamma_rnd: c must be a positive integer");
     endif
-    [retval, a, b] = common_size (a, b, zeros (r, c));
-    if (retval > 0)
-      error ("gamma_rnd: a and b must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("uniform_rnd: r must be a postive integer or vector");
     endif
   elseif (nargin == 2)
-    [retval, a, b] = common_size (a, b);
-    if (retval > 0)
-      error ("gamma_rnd: a and b must be of common size or scalar");
-    endif
+    sz = size(a);
   else
     usage ("gamma_rnd (a, b, r, c)");
   endif
 
-  [r, c] = size (a);
-  s = r * c;
-  a = reshape (a, 1, s);
-  b = reshape (b, 1, s);
-  rnd = zeros (1, s);
+  rnd = zeros (sz);
 
-  k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
-  endif
-  k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
-  if (any (k))
-    rnd(k) = gamma_inv (rand (1, length (k)), a(k), b(k));
+  if (isscalar (a) && isscalar(b))
+    if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
+      rnd = NaN * ones (sz);
+    else
+      rnd =  gamma_inv (rand (sz), a, b);
+    endif
+  else 
+    k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+    k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
+    if (any (k))
+      rnd(k) = gamma_inv (rand (size (k)), a(k), b(k));
+    endif
   endif
 
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/normal_cdf.m b/scripts/statistics/distributions/normal_cdf.m
--- a/scripts/statistics/distributions/normal_cdf.m
+++ b/scripts/statistics/distributions/normal_cdf.m
@@ -40,30 +40,34 @@
     v = 1;
   endif
 
-  [retval, x, m, v] = common_size (x, m, v);
-  if (retval > 0)
-    error ("normal_cdf: x, m and v must be of common size or scalars");
+  if (!isscalar (m) || !isscalar(v))
+    [retval, x, m, v] = common_size (x, m, v);
+    if (retval > 0)
+      error ("normal_cdf: x, m and v must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  m = reshape (m, 1, s);
-  v = reshape (v, 1, s);
-  cdf = zeros (1, s);
+  sz = size (x);
+  cdf = zeros (sz);
 
-  k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf));
-  if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
-  endif
+  if (isscalar (m) && isscalar(v))
+    if (find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)))
+      cdf = NaN * ones (sz);
+    else
+      cdf =  stdnormal_cdf ((x - m) ./ sqrt (v));
+    endif
+  else
+    k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf));
+    if (any (k))
+      cdf(k) = NaN;
+    endif
 
-  k = find (!isinf (m) & !isnan (m) & (v >= 0) & (v < Inf));
-  if (any (k))
-    cdf(k) = stdnormal_cdf ((x(k) - m(k)) ./ sqrt (v(k)));
+    k = find (!isinf (m) & !isnan (m) & (v >= 0) & (v < Inf));
+    if (any (k))
+      cdf(k) = stdnormal_cdf ((x(k) - m(k)) ./ sqrt (v(k)));
+    endif
   endif
 
   cdf((v == 0) & (x == m)) = 0.5;
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/normal_inv.m b/scripts/statistics/distributions/normal_inv.m
--- a/scripts/statistics/distributions/normal_inv.m
+++ b/scripts/statistics/distributions/normal_inv.m
@@ -31,7 +31,7 @@
 
 function inv = normal_inv (x, m, v)
 
-  if (! ((nargin == 1) || (nargin == 3)))
+  if (nargin != 1 && nargin != 3)
     usage ("normal_inv (x, m, v)");
   endif
 
@@ -40,26 +40,32 @@
     v = 1;
   endif
 
-  [retval, x, m, v] = common_size (x, m, v);
-  if (retval > 0)
-    error ("normal_inv: x, m and v must be of common size or scalars");
+  if (!isscalar (m) || !isscalar(v))
+    [retval, x, m, v] = common_size (x, m, v);
+    if (retval > 0)
+      error ("normal_inv: x, m and v must be of common size or scalars");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  m = reshape (m, 1, s);
-  v = reshape (v, 1, s);
-  inv = zeros (1, s);
+  sz = size (x);
+  inv = zeros (sz);
 
-  k = find (isinf (m) | isnan (m) | !(v > 0) | !(v < Inf));
-  if (any (k))
-    inv(k) = NaN * ones (1, length (k));
-  endif
+  if (isscalar (m) && isscalar(v))
+    if (find (isinf (m) | isnan (m) | !(v > 0) | !(v < Inf)))
+      inv = NaN * ones (sz);
+    else
+      inv =  m + sqrt (v) .* stdnormal_inv (x);
+    endif
+  else
+    k = find (isinf (m) | isnan (m) | !(v > 0) | !(v < Inf));
+    if (any (k))
+      inv(k) = NaN;
+    endif
 
-  k = find (!isinf (m) & !isnan (m) & (v > 0) & (v < Inf));
-  if (any (k))
-    inv(k) = m(k) + sqrt (v(k)) .* stdnormal_inv (x(k));
+    k = find (!isinf (m) & !isnan (m) & (v > 0) & (v < Inf));
+    if (any (k))
+      inv(k) = m(k) + sqrt (v(k)) .* stdnormal_inv (x(k));
+    endif
   endif
 
   k = find ((v == 0) & (x > 0) & (x < 1));
@@ -70,6 +76,4 @@
   inv((v == 0) & (x == 0)) = -Inf;
   inv((v == 0) & (x == 1)) = Inf;
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/normal_pdf.m b/scripts/statistics/distributions/normal_pdf.m
--- a/scripts/statistics/distributions/normal_pdf.m
+++ b/scripts/statistics/distributions/normal_pdf.m
@@ -31,7 +31,7 @@
 
 function pdf = normal_pdf (x, m, v)
 
-  if (! ((nargin == 1) || (nargin == 3)))
+  if (nargin != 1 && nargin != 3)
     usage ("normal_pdf (x, m, v)");
   endif
 
@@ -40,31 +40,35 @@
     v = 1;
   endif
 
-  [retval, x, m, v] = common_size (x, m, v);
-  if (retval > 0)
-    error ("normal_pdf: x, m and v must be of common size or scalars");
+  if (!isscalar (m) || !isscalar(v))
+    [retval, x, m, v] = common_size (x, m, v);
+    if (retval > 0)
+      error ("normal_pdf: x, m and v must be of common size or scalars");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  m = reshape (m, 1, s);
-  v = reshape (v, 1, s);
-  pdf = zeros (1, s);
+  sz = size (x);
+  pdf = zeros (sz);
 
-  k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf));
-  if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
-  endif
+  if (isscalar (m) && isscalar(v))
+    if (find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)))
+      pdf = NaN * ones (sz);
+    else
+      pdf = stdnormal_pdf ((x - m) ./ sqrt (v)) ./ sqrt (v);
+    endif
+  else
+    k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf));
+    if (any (k))
+      pdf(k) = NaN;
+    endif
 
-  k = find (!isinf (m) & !isnan (m) & (v >= 0) & (v < Inf));
-  if (any (k))
-    pdf(k) = stdnormal_pdf ((x(k) - m(k)) ./ sqrt (v(k))) ./ sqrt (v(k));
+    k = find (!isinf (m) & !isnan (m) & (v >= 0) & (v < Inf));
+    if (any (k))
+      pdf(k) = stdnormal_pdf ((x(k) - m(k)) ./ sqrt (v(k))) ./ sqrt (v(k));
+    endif
   endif
 
   pdf((v == 0) & (x == m)) = Inf;
   pdf((v == 0) & ((x < m) | (x > m))) = 0;
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/normal_rnd.m b/scripts/statistics/distributions/normal_rnd.m
--- a/scripts/statistics/distributions/normal_rnd.m
+++ b/scripts/statistics/distributions/normal_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} normal_rnd (@var{m}, @var{v}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the
-## normal distribution with parameters @var{m} and @var{v}.  Both
-## @var{m} and @var{v} must be scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} normal_rnd (@var{m}, @var{v}, @var{sz})
+## Return an @var{r} by @var{c}  or @code{size (@var{sz})} matrix of
+## random samples from the normal distribution with parameters @var{m} 
+## and @var{v}.  Both @var{m} and @var{v} must be scalar or of size 
+## @var{r} by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{m} and @var{v}.
@@ -32,6 +34,15 @@
 
 function rnd = normal_rnd (m, v, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(m) || !isscalar(v)) 
+      [retval, m, v] = common_size (m, v);
+      if (retval > 0)
+	error ("normal_rnd: m and v must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("normal_rnd: r must be a positive integer");
@@ -39,35 +50,33 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("normal_rnd: c must be a positive integer");
     endif
-    [retval, m, v] = common_size (m, v, zeros (r, c));
-    if (retval > 0)
-      error ("normal_rnd: m and v must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("normal_rnd: r must be a postive integer or vector");
     endif
   elseif (nargin == 2)
-    [retval, m, v] = common_size (m, v);
-    if (retval > 0)
-      error ("normal_rnd: m and v must be of common size or scalar");
-    endif
+    sz = size(m);
   else
     usage ("normal_rnd (m, v, r, c)");
   endif
 
-  [r, c] = size (m);
-  s = r * c;
-  m = reshape (m, 1, s);
-  v = reshape (v, 1, s);
-  rnd = zeros (1, s);
-
-  k = find (isnan (m) | isinf (m) | !(v > 0) | !(v < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+  if (isscalar (m) && isscalar (v))
+    if (find (isnan (m) | isinf (m) | !(v > 0) | !(v < Inf)))
+      rnd = NaN * ones (sz);
+    else
+      rnd =  m + sqrt (v) .* randn (sz);
+    endif
+  else
+    rnd = m + sqrt (v) .* randn (sz);
+    k = find (isnan (m) | isinf (m) | !(v > 0) | !(v < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
   endif
 
-  k = find ((m > -Inf) & (m < Inf) & (v > 0) & (v < Inf));
-  if (any (k))
-    rnd(k) = m(k) + sqrt (v(k)) .* randn (1, length (k));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/stdnormal_cdf.m b/scripts/statistics/distributions/stdnormal_cdf.m
--- a/scripts/statistics/distributions/stdnormal_cdf.m
+++ b/scripts/statistics/distributions/stdnormal_cdf.m
@@ -32,12 +32,12 @@
     usage ("stdnormal_cdf (x)");
   endif
 
-  [r_x, c_x] = size (x);
-  if (r_x * c_x == 0)
+  sz = size (x);
+  if (numel(x) == 0)
     error ("stdnormal_cdf: x must not be empty");
   endif
 
-  cdf = (ones (r_x, c_x) + erf (x / sqrt (2))) / 2;
+  cdf = (ones (sz) + erf (x / sqrt (2))) / 2;
 
 endfunction
 
diff --git a/scripts/statistics/distributions/stdnormal_pdf.m b/scripts/statistics/distributions/stdnormal_pdf.m
--- a/scripts/statistics/distributions/stdnormal_pdf.m
+++ b/scripts/statistics/distributions/stdnormal_pdf.m
@@ -32,14 +32,12 @@
     usage ("stdnormal_pdf (x)");
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  pdf = zeros (1, s);
+  sz = size(x);
+  pdf = zeros (sz);
 
   k = find (isnan (x));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find (!isinf (x));
@@ -47,6 +45,4 @@
     pdf (k) = (2 * pi)^(- 1/2) * exp (- x(k) .^ 2 / 2);
   endif
 
-  pdf = reshape (pdf, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/stdnormal_rnd.m b/scripts/statistics/distributions/stdnormal_rnd.m
--- a/scripts/statistics/distributions/stdnormal_rnd.m
+++ b/scripts/statistics/distributions/stdnormal_rnd.m
@@ -19,8 +19,9 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} stdnormal_rnd (@var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random numbers from the
-## standard normal distribution.
+## @deftypefnx {Function File} {} stdnormal_rnd (@var{sz})
+## Return an @var{r} by @var{c} or @code{size (@var{sz})} matrix of 
+## random numbers from the standard normal distribution.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
@@ -28,17 +29,28 @@
 
 function rnd = stdnormal_rnd (r, c)
 
-  if (nargin != 2)
+  if (nargin != 1 && nargin != 2)
     usage ("stdnormal_rnd (r, c)");
   endif
 
-  if (! (isscalar (r) && (r > 0) && (r == round (r))))
-    error ("stdnormal_rnd: r must be a positive integer");
-  endif
-  if (! (isscalar (c) && (c > 0) && (c == round (c))))
-    error ("stdnormal_rnd: c must be a positive integer");
+  if (nargin == 2)
+    if (! (isscalar (r) && (r > 0) && (r == round (r))))
+      error ("stdnormal_rnd: r must be a positive integer");
+    endif
+    if (! (isscalar (c) && (c > 0) && (c == round (c))))
+      error ("stdnormal_rnd: c must be a positive integer");
+    endif
+    sz = [r, c];
+  else
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("stdnormal_rnd: r must be a postive integer or vector");
+    endif
   endif
 
-  rnd = randn (r, c);
+  rnd = randn (sz);
 
 endfunction
diff --git a/scripts/statistics/distributions/uniform_cdf.m b/scripts/statistics/distributions/uniform_cdf.m
--- a/scripts/statistics/distributions/uniform_cdf.m
+++ b/scripts/statistics/distributions/uniform_cdf.m
@@ -30,7 +30,7 @@
 
 function cdf = uniform_cdf (x, a, b)
 
-  if (! (nargin == 1 || nargin == 3))
+  if (nargin != 1 && nargin != 3)
     usage ("uniform_cdf (x, a, b)");
   endif
 
@@ -39,33 +39,33 @@
     b = 1;
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("uniform_cdf: x, a and b must be of common size or scalar");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("uniform_cdf: x, a and b must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  a = reshape (a, 1, s);
-  b = reshape (b, 1, s);
-  cdf = zeros (1, s);
+  sz = size (x);
+  cdf = zeros (sz);
 
   k = find (isnan (x) | !(a < b));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x >= b) & (a < b));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
-
+  
   k = find ((x > a) & (x < b));
   if (any (k))
-    cdf(k) = (x(k) < b(k)) .* (x(k) - a(k)) ./ (b(k) - a(k));
+    if (isscalar (a) && isscalar(b))
+      cdf(k) = (x(k) < b) .* (x(k) - a) ./ (b - a);
+    else
+      cdf(k) = (x(k) < b(k)) .* (x(k) - a(k)) ./ (b(k) - a(k));
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/uniform_inv.m b/scripts/statistics/distributions/uniform_inv.m
--- a/scripts/statistics/distributions/uniform_inv.m
+++ b/scripts/statistics/distributions/uniform_inv.m
@@ -30,7 +30,7 @@
 
 function inv = uniform_inv (x, a, b)
 
-  if (! (nargin == 1 || nargin == 3))
+  if (nargin != 1 && nargin != 3)
     usage ("uniform_inv (x, a, b)");
   endif
 
@@ -39,28 +39,28 @@
     b = 1;
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("uniform_inv: x, a and b must be of common size or scalars");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("uniform_cdf: x, a and b must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s   = r * c;
-  x   = reshape (x, 1, s);
-  a   = reshape (a, 1, s);
-  b   = reshape (b, 1, s);
-  inv = zeros (1, s);
+  sz = size (x);
+  inv = zeros (sz);
 
   k = find ((x < 0) | (x > 1) | isnan (x) | !(a < b));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x >= 0) & (x <= 1) & (a < b));
   if (any (k))
-    inv(k) = a(k) + x(k) .* (b(k) - a(k));
+    if (isscalar (a) && isscalar(b))
+      inv(k) = a + x(k) .* (b - a);
+    else
+      inv(k) = a(k) + x(k) .* (b(k) - a(k));
+    endif
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/uniform_pdf.m b/scripts/statistics/distributions/uniform_pdf.m
--- a/scripts/statistics/distributions/uniform_pdf.m
+++ b/scripts/statistics/distributions/uniform_pdf.m
@@ -30,7 +30,7 @@
 
 function pdf = uniform_pdf (x, a, b)
 
-  if (! (nargin == 1 || nargin == 3))
+  if (nargin != 1 && nargin != 3)
     usage ("uniform_pdf (x, a, b)");
   endif
 
@@ -39,28 +39,28 @@
     b = 1;
   endif
 
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("uniform_pdf: x, a and b must be of common size or scalars");
+  if (!isscalar (a) || !isscalar(b))
+    [retval, x, a, b] = common_size (x, a, b);
+    if (retval > 0)
+      error ("uniform_pdf: x, a and b must be of common size or scalars");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  a = reshape (a, 1, s);
-  b = reshape (b, 1, s);
-  pdf = zeros (1, s);
+  sz = size (x);
+  pdf = zeros (sz);
 
   k = find (isnan (x) | !(a < b));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find ((x > a) & (x < b));
   if (any (k))
-    pdf(k) = 1 ./ (b(k) - a(k));
+    if (isscalar (a) && isscalar(b))
+      pdf(k) = 1 ./ (b - a);
+    else
+      pdf(k) = 1 ./ (b(k) - a(k));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/uniform_rnd.m b/scripts/statistics/distributions/uniform_rnd.m
--- a/scripts/statistics/distributions/uniform_rnd.m
+++ b/scripts/statistics/distributions/uniform_rnd.m
@@ -19,9 +19,10 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} uniform_rnd (@var{a}, @var{b}, @var{r}, @var{c})
-## Return an @var{r} by @var{c} matrix of random samples from the
-## uniform distribution on [@var{a}, @var{b}].  Both @var{a} and @var{b}
-## must be scalar or of size @var{r} by @var{c}.
+## @deftypefnx {Function File} {} uniform_rnd (@var{a}, @var{b}, @var{sz})
+## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of 
+## random samples from the uniform distribution on [@var{a}, @var{b}]. 
+## Both @var{a} and @var{b} must be scalar or of size @var{r} by @var{c}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{a} and @var{b}.
@@ -32,6 +33,15 @@
 
 function rnd = uniform_rnd (a, b, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(a) || !isscalar(b)) 
+      [retval, a, b] = common_size (a, b);
+      if (retval > 0)
+	error ("uniform_rnd: a and b must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("uniform_rnd: r must be a positive integer");
@@ -39,35 +49,34 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("uniform_rnd: c must be a positive integer");
     endif
-    [retval, a, b] = common_size (a, b, zeros (r, c));
-    if (retval > 0)
-      error ("uniform_rnd: a and b must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("uniform_rnd: r must be a postive integer or vector");
     endif
   elseif (nargin == 2)
-    [retval, a, b] = common_size (a, b);
-    if (retval > 0)
-      error ("uniform_rnd: a and b must be of common size or scalar");
-    endif
+    sz = size(a);
   else
     usage ("uniform_rnd (a, b, r, c)");
   endif
 
-  [r, c] = size (a);
-  s = r * c;
-  a = reshape (a, 1, s);
-  b = reshape (b, 1, s);
-  rnd = zeros (1, s);
+  if (isscalar(a) && isscalar(b))
+    if (find (!(-Inf < a) | !(a < b) | !(b < Inf)))
+      rnd = NaN * ones(sz);
+    else
+      rnd =  a + (b - a) .* rand (sz);
+    endif
+  else
+    rnd =  a + (b - a) .* rand (sz);
 
-  k = find (!(-Inf < a) | !(a < b) | !(b < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(-Inf < a) | !(a < b) | !(b < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
   endif
 
-  k = find ((-Inf < a) & (a < b) & (b < Inf));
-  if (any (k))
-    rnd(k) = a(k) + (b(k) - a(k)) .* rand (1, length (k));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/test/octave.test/stats/median-3.m b/test/octave.test/stats/median-3.m
--- a/test/octave.test/stats/median-3.m
+++ b/test/octave.test/stats/median-3.m
@@ -1,1 +1,1 @@
-median (1, 2)
+median (1, 2, 3)