# HG changeset patch
# User jwe
# Date 1081468365 0
# Node ID 265d566cc770d8b782a0cd45822065d3a2943e70
# Parent  499d2ca4698240fe49e481a5f82d35240df9c8fd
[project @ 2004-04-08 23:52:45 by jwe]

diff --git a/scripts/statistics/distributions/discrete_cdf.m b/scripts/statistics/distributions/discrete_cdf.m
--- a/scripts/statistics/distributions/discrete_cdf.m
+++ b/scripts/statistics/distributions/discrete_cdf.m
@@ -33,7 +33,7 @@
     usage ("discrete_cdf (x, v, p)");
   endif
 
-  [r, c] = size (x);
+  sz = size (x);
 
   if (! isvector (v))
     error ("discrete_cdf: v must be a vector");
@@ -43,16 +43,16 @@
     error ("discrete_cdf: p must be a nonzero, nonnegative vector");
   endif
 
-  n = r * c;
+  n = numel (x);
   m = length (v);
   x = reshape (x, n, 1);
   v = reshape (v, 1, m);
   p = reshape (p / sum (p), m, 1);
 
-  cdf = zeros (n, 1);
+  cdf = zeros (sz);
   k = find (isnan (x));
   if (any (k))
-    cdf (k) = NaN * ones (length (k), 1);
+    cdf (k) = NaN;
   endif
   k = find (!isnan (x));
   if (any (k))
@@ -60,6 +60,4 @@
     cdf (k) = ((x(k) * ones (1, m)) >= (ones (n, 1) * v)) * p;
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/discrete_inv.m b/scripts/statistics/distributions/discrete_inv.m
--- a/scripts/statistics/distributions/discrete_inv.m
+++ b/scripts/statistics/distributions/discrete_inv.m
@@ -33,7 +33,7 @@
     usage ("discrete_inv (x, v, p)");
   endif
 
-  [r, c] = size (x);
+  sz = size (x);
 
   if (! isvector (v))
     error ("discrete_inv: v must be a vector");
@@ -43,7 +43,7 @@
     error ("discrete_inv: p must be a nonzero, nonnegative vector");
   endif
 
-  n = r * c;
+  n = numel (x);
   x = reshape (x, 1, n);
   m = length (v);
   v = sort (v);
@@ -52,12 +52,12 @@
   ## Allow storage allocated for P to be reclaimed.
   p = [];
 
-  inv = NaN * ones (n, 1);
+  inv = NaN * ones (sz);
   if (any (k = find (x == 0)))
-    inv(k) = -Inf * ones (1, length (k));
+    inv(k) = -Inf;
   endif
   if (any (k = find (x == 1)))
-    inv(k) = v(m) * ones (1, length (k));
+    inv(k) = v(m) * ones (size (k));
   endif
 
   if (any (k = find ((x > 0) & (x < 1))))
@@ -68,15 +68,13 @@
     ##
     ## Vectorized code is:
     ##
-    ##   inv(k) = v(sum ((ones (m, 1) * x(k)) > (s * ones (1, n))) + 1);
+    ##     inv(k) = v(sum ((ones (m, 1) * x(k)) > (s * ones (1, n))) + 1);
 
     for q = 1:n
-      inv(q) = v(sum (x(q) > s) + 1);
+      inv(k(q)) = v(sum (x(k(q)) > s) + 1);
     endfor
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
 
 
diff --git a/scripts/statistics/distributions/discrete_pdf.m b/scripts/statistics/distributions/discrete_pdf.m
--- a/scripts/statistics/distributions/discrete_pdf.m
+++ b/scripts/statistics/distributions/discrete_pdf.m
@@ -33,7 +33,7 @@
     usage ("discrete_pdf (x, v, p)");
   endif
 
-  [r, c] = size (x);
+  sz = size (x);
 
   if (! isvector (v))
     error ("discrete_pdf: v must be a vector");
@@ -43,16 +43,16 @@
     error ("discrete_pdf: p must be a nonzero, nonnegative vector");
   endif
 
-  n = r * c;
+  n = numel (x);
   m = length (v);
   x = reshape (x, n, 1);
   v = reshape (v, 1, m);
   p = reshape (p / sum (p), m, 1);
 
-  pdf = zeros (n, 1);
+  pdf = zeros (sz);
   k = find (isnan (x));
   if (any (k))
-    pdf (k) = NaN * ones (length (k), 1);
+    pdf (k) = NaN;
   endif
   k = find (!isnan (x));
   if (any (k))
@@ -60,6 +60,4 @@
     pdf (k) = ((x(k) * ones (1, m)) == (ones (n, 1) * v)) * p;
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/discrete_rnd.m b/scripts/statistics/distributions/discrete_rnd.m
--- a/scripts/statistics/distributions/discrete_rnd.m
+++ b/scripts/statistics/distributions/discrete_rnd.m
@@ -19,24 +19,51 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} discrete_rnd (@var{n}, @var{v}, @var{p})
+## @deftypefnx {Function File} {} discrete_rnd (@var{v}, @var{p}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} discrete_rnd (@var{v}, @var{p}, @var{sz})
 ## Generate a row vector containing a random sample of size @var{n} from
 ## the univariate distribution which assumes the values in @var{v} with
-## probabilities @var{p}.
+## probabilities @var{p}. @var{n} must be a scalar.
 ##
-## Currently, @var{n} must be a scalar.
+## If @var{r} and @var{c} are given create a matrix with @var{r} rows and
+## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size
+## @var{sz}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Description: Random deviates from a discrete distribution
 
-function rnd = discrete_rnd (n, v, p)
+function rnd = discrete_rnd (v, p, r, c)
+
+  if (nargin == 4)
+    if (! (isscalar (r) && (r > 0) && (r == round (r))))
+      error ("discrete_rnd: r must be a positive integer");
+    endif
+    if (! (isscalar (c) && (c > 0) && (c == round (c))))
+      error ("discrete_rnd: c must be a positive integer");
+    endif
+    sz = [r, c];
+  elseif (nargin == 3)
+    ## A potential problem happens here if all args are scalar, as
+    ## we can't distiguish between the command syntax. Thankfully this
+    ## case doesn't make much sense. So we assume the first syntax
+    ## if the first arg is scalar
 
-  if (nargin != 3)
-    usage ("discrete_rnd (n, v, p)");
-  endif
-
-  if (! isscalar (n))
-    error ("discrete_rnd: n must be a scalar");
+    if (isscalar (v))
+      sz = [1, floor(v)];
+      v = p;
+      p = r;
+    else
+      if (isscalar (r) && (r > 0))
+	sz = [r, r];
+      elseif (isvector(r) && all (r > 0))
+	sz = r(:)';
+      else
+	error ("discrete_rnd: r must be a postive integer or vector");
+      endif
+    endif
+  else
+    usage ("discrete_rnd (n, v, p) | discrete_rnd (v, p, r, c)");
   endif
 
   if (! isvector (v))
@@ -47,10 +74,24 @@
     error ("discrete_rnd: p must be a nonzero, nonnegative vector");
   endif
 
+  n = prod (sz);
+  m = length (v);
   u = rand (1, n);
-  m = length (p);
   s = reshape (cumsum (p / sum (p)), m, 1);
 
+  ## The following loop is a space/time tradeoff in favor of space,
+  ## since the dataset may be large.
+  ##
+  ## Vectorized code is:
+  ##
   rnd = v (1 + sum ((s * ones (1, n)) <= ((ones (m, 1) * u))));
+  rnd = reshape (rnd, sz);
+  ##
+  ## Non-vectorized code is:
+  ##
+  ##  rnd = zeros (sz);
+  ##  for q=1:n
+  ##    rnd (q) = v (sum (s <= u (q)) + 1);
+  ##  endfor
 
 endfunction
diff --git a/scripts/statistics/distributions/empirical_rnd.m b/scripts/statistics/distributions/empirical_rnd.m
--- a/scripts/statistics/distributions/empirical_rnd.m
+++ b/scripts/statistics/distributions/empirical_rnd.m
@@ -19,19 +19,35 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} empirical_rnd (@var{n}, @var{data})
+## @deftypefnx {Function File} {} empirical_rnd (@var{data}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} empirical_rnd (@var{data}, @var{sz})
 ## Generate a bootstrap sample of size @var{n} from the empirical
 ## distribution obtained from the univariate sample @var{data}.
+##
+## If @var{r} and @var{c} are given create a matrix with @var{r} rows and
+## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size
+## @var{sz}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Description: Bootstrap samples from the empirical distribution
 
-function rnd = empirical_rnd (n, data)
+function rnd = empirical_rnd (data, r, c)
+
+  if (nargin == 2)
+    if (isscalar(data))
+      c = data;
+      data = r;
+      r = 1;
+    endif
+  elseif (nargin != 3)
+    usage ("empirical_rnd (n, data) | empirical_rnd (data, r, c)");
+  endif
 
   if (! isvector (data))
     error ("empirical_rnd: data must be a vector");
   endif
 
-  rnd = discrete_rnd (n, data, ones (size (data)) / length (data));
+  rnd = discrete_rnd (data, ones (size (data)) / length (data), r, c);
 
 endfunction
diff --git a/scripts/statistics/distributions/exponential_cdf.m b/scripts/statistics/distributions/exponential_cdf.m
--- a/scripts/statistics/distributions/exponential_cdf.m
+++ b/scripts/statistics/distributions/exponential_cdf.m
@@ -35,32 +35,40 @@
     usage ("exponential_cdf (x, lambda)");
   endif
 
-  [retval, x, l] = common_size (x, l);
-  if (retval > 0)
-    error ("exponential_cdf: x and lambda must be of common size or scalar");
+  if (!isscalar (x) && !isscalar(l))
+    [retval, x, l] = common_size (x, l);
+    if (retval > 0)
+      error ("exponential_cdf: x and lambda must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  l = reshape (l, 1, s);
-  cdf = zeros (1, s);
+  if (isscalar (x))
+    sz = size (l);
+  else
+    sz = size (x);
+  endif
+
+  cdf = zeros (sz);
 
   k = find (isnan (x) | !(l > 0));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (l > 0));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x > 0) & (x < Inf) & (l > 0));
   if (any (k))
-    cdf (k) = 1 - exp (- l(k) .* x(k));
+    if isscalar (l)
+      cdf (k) = 1 - exp (- l .* x(k));
+    elseif isscalar (x)
+      cdf (k) = 1 - exp (- l(k) .* x);
+    else
+      cdf (k) = 1 - exp (- l(k) .* x(k));
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/exponential_inv.m b/scripts/statistics/distributions/exponential_inv.m
--- a/scripts/statistics/distributions/exponential_inv.m
+++ b/scripts/statistics/distributions/exponential_inv.m
@@ -33,32 +33,40 @@
     usage ("exponential_inv (x, lambda)");
   endif
 
-  [retval, x, l] = common_size (x, l);
-  if (retval > 0)
-    error ("exponential_inv: x and lambda must be of common size or scalar");
+  if (!isscalar (x) && !isscalar(l))
+    [retval, x, l] = common_size (x, l);
+    if (retval > 0)
+      error ("exponential_inv: x and lambda must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  l = reshape (l, 1, s);
-  inv = zeros (1, s);
+  if (isscalar (x))
+    sz = size (l);
+  else
+    sz = size (x);
+  endif
+
+  inv = zeros (sz);
 
   k = find (!(l > 0) | (x < 0) | (x > 1) | isnan (x));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 1) & (l > 0));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (l > 0));
   if (any (k))
-    inv(k) = - log (1 - x(k)) ./ l(k);
+    if isscalar (l)
+      inv(k) = - log (1 - x(k)) ./ l;
+    elseif isscalar (x)
+      inv(k) = - log (1 - x) ./ l(k);
+    else
+      inv(k) = - log (1 - x(k)) ./ l(k);
+    endif
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/exponential_pdf.m b/scripts/statistics/distributions/exponential_pdf.m
--- a/scripts/statistics/distributions/exponential_pdf.m
+++ b/scripts/statistics/distributions/exponential_pdf.m
@@ -32,27 +32,34 @@
     usage ("exponential_pdf (x, lambda)");
   endif
 
-  [retval, x, l] = common_size (x, l);
-  if (retval > 0)
-    error ("exponential_pdf: x and lambda must be of common size or scalar");
+  if (!isscalar (x) && !isscalar(l))
+    [retval, x, l] = common_size (x, l);
+    if (retval > 0)
+      error ("exponential_pdf: x and lambda must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  l = reshape (l, 1, s);
-  pdf = zeros (1, s);
+  if (isscalar (x))
+    sz = size (l);
+  else
+    sz = size (x);
+  endif
+  pdf = zeros (sz);
 
   k = find (!(l > 0) | isnan (x));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find ((x > 0) & (x < Inf) & (l > 0));
   if (any (k))
-    pdf(k) = l(k) .* exp (- l(k) .* x(k));
+    if isscalar (l)
+      pdf(k) = l .* exp (- l .* x(k));
+    elseif isscalar (x)
+      pdf(k) = l(k) .* exp (- l(k) .* x);
+    else
+      pdf(k) = l(k) .* exp (- l(k) .* x(k));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/exponential_rnd.m b/scripts/statistics/distributions/exponential_rnd.m
--- a/scripts/statistics/distributions/exponential_rnd.m
+++ b/scripts/statistics/distributions/exponential_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} exponential_rnd (@var{lambda}, @var{r}, @var{c})
+## @deftypefn {Function File} {} exponential_rnd (@var{lambda}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the
 ## exponential distribution with parameter @var{lambda}, which must be a
-## scalar or of size @var{r} by @var{c}.
+## scalar or of size @var{r} by @var{c}. Or if @var{sz} is a vector, 
+## create a matrix of size @var{sz}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the size of @var{lambda}.
@@ -39,29 +41,48 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("exponential_rnd: c must be a positive integer");
     endif
-    [retval, l] = common_size (l, zeros (r, c));
-    if (retval > 0)
-      error ("exponential_rnd: lambda must be scalar or of size %d by %d",
-	     r, c);
+    sz = [r, c];
+
+    if (any (size (l) != 1) && 
+	((length (size (nl)) != length (sz)) || any (size (l) != sz)))
+      error ("exponential_rnd: lambda must be scalar or of size [r, c]");
     endif
-  elseif (nargin != 1)
+  elseif (nargin == 2)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("exponential_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (l) != 1) && 
+	((length (size (l)) != length (sz)) || any (size (l) != sz)))
+      error ("exponential_rnd: lambda must be scalar or of size sz");
+    endif
+  elseif (nargin == 1)
+    sz = size (l);
+  else
     usage ("exponential_rnd (lambda, r, c)");
   endif
 
-  [r, c] = size (l);
-  s = r * c;
-  l = reshape (l, 1, s);
-  rnd = zeros (1, s);
 
-  k = find (!(l > 0) | !(l < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
-  endif
-  k = find ((l > 0) & (l < Inf));
-  if (any (k))
-    rnd(k) = - log (1 - rand (1, length (k))) ./ l(k);
+  if (isscalar (l))
+    if ((l > 0) && (l < Inf))
+      rnd = - log (1 - rand (sz)) ./ l;
+    else
+      rnd = NaN * ones (sz);
+    endif
+  else
+    rnd = zeros (sz);
+    k = find (!(l > 0) | !(l < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+    k = find ((l > 0) & (l < Inf));
+    if (any (k))
+      rnd(k) = - log (1 - rand (size (k))) ./ l(k);
+    endif
   endif
 
-  rnd = reshape (rnd, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/f_cdf.m b/scripts/statistics/distributions/f_cdf.m
--- a/scripts/statistics/distributions/f_cdf.m
+++ b/scripts/statistics/distributions/f_cdf.m
@@ -33,33 +33,34 @@
     usage ("f_cdf (x, m, n)");
   endif
 
-  [retval, x, m, n] = common_size (x, m, n);
-  if (retval > 0)
-    error ("f_cdf: x, m and n must be of common size or scalar");
+  if (!isscalar (m) || !isscalar (n))
+    [retval, x, m, n] = common_size (x, m, n);
+    if (retval > 0)
+      error ("f_cdf: x, m and n 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);
-  n = reshape (n, 1, s);
-  cdf = zeros (1, s);
+  sz = size (x);
+  cdf = zeros (sz);
 
   k = find (!(m > 0) | !(n > 0) | isnan (x));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (m > 0) & (n > 0));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0));
   if (any (k))
-    cdf(k) = 1 - betainc (1 ./ (1 + m(k) .* x(k) ./ n(k)), n(k) / 2, m(k) / 2);
+    if (isscalar (m) && isscalar (n))
+      cdf(k) = 1 - betainc (1 ./ (1 + m .* x(k) ./ n), n / 2, m / 2);
+    else
+      cdf(k) = 1 - betainc (1 ./ (1 + m(k) .* x(k) ./ n(k)), n(k) / 2, 
+			    m(k) / 2);
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/f_inv.m b/scripts/statistics/distributions/f_inv.m
--- a/scripts/statistics/distributions/f_inv.m
+++ b/scripts/statistics/distributions/f_inv.m
@@ -33,34 +33,34 @@
     usage ("f_inv (x, m, n)");
   endif
 
-  [retval, x, m, n] = common_size (x, m, n);
-  if (retval > 0)
-    error ("f_inv: x, m and n must be of common size or scalar");
+  if (!isscalar (m) || !isscalar (n))
+    [retval, x, m, n] = common_size (x, m, n);
+    if (retval > 0)
+      error ("f_inv: x, m and n 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);
-  n = reshape (n, 1, s);
-  inv = zeros (1, s);
+  sz = size (x);
+  inv = zeros (sz);
 
   k = find ((x < 0) | (x > 1) | isnan (x) | !(m > 0) | !(n > 0));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 1) & (m > 0) & (n > 0));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (m > 0) & (n > 0));
   if (any (k))
-    inv(k) = ((1 ./ beta_inv (1 - x(k), n(k) / 2, m(k) / 2) - 1)
-	      .* n(k) ./ m(k));
+    if (isscalar (m) && isscalar (n))
+      inv(k) = ((1 ./ beta_inv (1 - x(k), n / 2, m / 2) - 1) .* n ./ m);
+    else
+      inv(k) = ((1 ./ beta_inv (1 - x(k), n(k) / 2, m(k) / 2) - 1)
+		.* n(k) ./ m(k));
+    endif
   endif
 
-  inv = reshape (inv, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/f_pdf.m b/scripts/statistics/distributions/f_pdf.m
--- a/scripts/statistics/distributions/f_pdf.m
+++ b/scripts/statistics/distributions/f_pdf.m
@@ -33,31 +33,34 @@
     usage ("f_pdf (x, m, n)");
   endif
 
-  [retval, x, m, n] = common_size (x, m, n);
-  if (retval > 0)
-    error ("f_pdf: x, m and n must be of common size or scalar");
+  if (!isscalar (m) || !isscalar (n))
+    [retval, x, m, n] = common_size (x, m, n);
+    if (retval > 0)
+      error ("f_pdf: x, m and n 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);
-  n = reshape (n, 1, s);
-  pdf = zeros (1, s);
+  sz = size (x);
+  pdf = zeros (sz);
 
   k = find (isnan (x) | !(m > 0) | !(n > 0));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0));
   if (any (k))
-    tmp = m(k) .* x(k) ./ n(k);
-    pdf(k) = (exp ((m(k) / 2 - 1) .* log (tmp)
-		   - ((m(k) + n(k)) / 2) .* log (1 + tmp))
-	      .* (m(k) ./ n(k)) ./ beta (m(k) / 2, n(k) / 2));
+    if (isscalar (m) && isscalar (n))
+      tmp = m / n * x(k);
+      pdf(k) = (exp ((m / 2 - 1) .* log (tmp)
+		     - ((m + n) / 2) .* log (1 + tmp))
+		.* (m / n) ./ beta (m / 2, n / 2));
+    else
+      tmp = m(k) .* x(k) ./ n(k);
+      pdf(k) = (exp ((m(k) / 2 - 1) .* log (tmp)
+		     - ((m(k) + n(k)) / 2) .* log (1 + tmp))
+		.* (m(k) ./ n(k)) ./ beta (m(k) / 2, n(k) / 2));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/f_rnd.m b/scripts/statistics/distributions/f_rnd.m
--- a/scripts/statistics/distributions/f_rnd.m
+++ b/scripts/statistics/distributions/f_rnd.m
@@ -19,9 +19,12 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} f_rnd (@var{m}, @var{n}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} f_rnd (@var{m}, @var{n}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the F
 ## distribution with @var{m} and @var{n} degrees of freedom.  Both
 ## @var{m} and @var{n} must be scalar or of size @var{r} by @var{c}.
+## If @var{sz} is a vector the random samples are in a matrix of 
+## size @var{sz}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{m} and @var{n}.
@@ -32,6 +35,16 @@
 
 function rnd = f_rnd (m, n, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(m) || !isscalar(n)) 
+      [retval, m, n] = common_size (m, n);
+      if (retval > 0)
+	error ("f_rnd: m and n must be of common size or scalar");
+      endif
+    endif
+  endif
+
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("f_rnd: r must be a positive integer");
@@ -39,37 +52,52 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("f_rnd: c must be a positive integer");
     endif
-    [retval, m, n] = common_size (m, n, zeros (r, c));
-    if (retval > 0)
-      error ("f_rnd: m and n must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+
+    if (any (size (m) != 1) && 
+	((length (size (m)) != length (sz)) || any (size (m) != sz)))
+      error ("f_rnd: m and n must be scalar or of size [r,c]");
+    endif
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("f_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (m) != 1) && 
+	((length (size (m)) != length (sz)) || any (size (m) != sz)))
+      error ("f_rnd: m and n must be scalar or of size sz");
     endif
   elseif (nargin == 2)
-    [retval, m, n] = common_size (m, n);
-    if (retval > 0)
-      error ("f_rnd: m and n must be of common size or scalar");
-    endif
+    sz = size(a);
   else
     usage ("f_rnd (m, n, r, c)");
   endif
 
-  [r, c] = size (m);
-  s = r * c;
-  m = reshape (m, 1, s);
-  n = reshape (n, 1, s);
-  rnd = zeros (1, s);
+
+  if (isscalar (m) && isscalar (n))
+    if ((m > 0) && (m < Inf) && (n > 0) && (n < Inf))
+      rnd =  f_inv (rand (sz), m, n);
+    else
+      rnd = NaN * ones (sz);
+    endif
+  else
+    rnd = zeros (sz);
 
-  k = find (!(m > 0) | !(m < Inf) |
-            !(n > 0) | !(n < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(m > 0) | !(m < Inf) |
+              !(n > 0) | !(n < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+
+    k = find ((m > 0) & (m < Inf) &
+              (n > 0) & (n < Inf));
+    if (any (k))
+      rnd(k) = f_inv (rand (size (k)), m(k), n(k));
+    endif
   endif
 
-  k = find ((m > 0) & (m < Inf) &
-            (n > 0) & (n < Inf));
-  if (any (k))
-    rnd(k) = f_inv (rand (1, length (k)), m(k), n(k));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/geometric_cdf.m b/scripts/statistics/distributions/geometric_cdf.m
--- a/scripts/statistics/distributions/geometric_cdf.m
+++ b/scripts/statistics/distributions/geometric_cdf.m
@@ -32,32 +32,34 @@
     usage ("geometric_cdf (x, p)");
   endif
 
-  [retval, x, p] = common_size (x, p);
-  if (retval > 0)
-    error ("geometric_cdf: x and p must be of common size or scalar");
+  if (!isscalar (x) && !isscalar (p))
+    [retval, x, p] = common_size (x, p);
+    if (retval > 0)
+      error ("geometric_cdf: x and p must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, 1, s);
-  p   = reshape (p, 1, s);
-  cdf = zeros (1, s);
+  cdf = zeros (size (x));
 
   k = find (isnan (x) | !(p >= 0) | !(p <= 1));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (p >= 0) & (p <= 1));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (p > 0) & (p <= 1));
   if (any (k))
-    cdf(k) = 1 - ((1 - p(k)) .^ (x(k) + 1));
+    if (isscalar (x))
+      cdf(k) = 1 - ((1 - p(k)) .^ (x + 1));
+    elseif (isscalar (p))
+      cdf(k) = 1 - ((1 - p) .^ (x(k) + 1));
+    else
+      cdf(k) = 1 - ((1 - p(k)) .^ (x(k) + 1));
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/geometric_inv.m b/scripts/statistics/distributions/geometric_inv.m
--- a/scripts/statistics/distributions/geometric_inv.m
+++ b/scripts/statistics/distributions/geometric_inv.m
@@ -32,33 +32,34 @@
     usage ("geometric_inv (x, p)");
   endif
 
-  [retval, x, p] = common_size (x, p);
-  if (retval > 0)
-    error ("geometric_inv: x and p must be of common size or scalar");
+  if (!isscalar (x) && !isscalar (p))
+    [retval, x, p] = common_size (x, p);
+    if (retval > 0)
+      error ("geometric_inv: x and p must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, 1, s);
-  p   = reshape (p, 1, s);
-  inv = zeros (1, s);
+  inv = zeros (size (x));
 
   k = find (!(x >= 0) | !(x <= 1) | !(p >= 0) | !(p <= 1));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 1) & (p >= 0) & (p <= 1));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (p > 0) & (p <= 1));
   if (any (k))
-    inv(k) = max (ceil (log (1 - x(k)) ./ log (1 - p(k))) - 1,
-		  zeros (1, length (k)));
+    if (isscalar (x))
+      inv(k) = max (ceil (log (1 - x) ./ log (1 - p(k))) - 1, 0);
+    elseif (isscalar (p))
+      inv(k) = max (ceil (log (1 - x(k)) / log (1 - p)) - 1, 0);
+    else
+      inv(k) = max (ceil (log (1 - x(k)) ./ log (1 - p(k))) - 1, 0);
+    endif
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/geometric_pdf.m b/scripts/statistics/distributions/geometric_pdf.m
--- a/scripts/statistics/distributions/geometric_pdf.m
+++ b/scripts/statistics/distributions/geometric_pdf.m
@@ -32,33 +32,35 @@
     usage ("geometric_pdf (x, p)");
   endif
 
-  [retval, x, p] = common_size (x, p);
-  if (retval > 0)
-    error ("geometric_pdf: x and p must be of common size or scalar");
+  if (!isscalar (x) && !isscalar (p))
+    [retval, x, p] = common_size (x, p);
+    if (retval > 0)
+      error ("geometric_pdf: x and p must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x   = reshape (x, 1, s);
-  p   = reshape (p, 1, s);
-  cdf = zeros (1, s);
+  pdf = zeros (size (x));
 
   k = find (isnan (x) | !(p >= 0) | !(p <= 1));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   ## Just for the fun of it ...
   k = find ((x == Inf) & (p == 0));
   if (any (k))
-    pdf(k) = ones (1, length (k));
+    pdf(k) = 1;
   endif
 
   k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (p > 0) & (p <= 1));
   if (any (k))
-    pdf(k) = p(k) .* ((1 - p(k)) .^ x(k));
+    if (isscalar (x))
+      pdf(k) = p(k) .* ((1 - p(k)) .^ x);
+    elseif (isscalar (p))
+      pdf(k) = p .* ((1 - p) .^ x(k));
+    else
+      pdf(k) = p(k) .* ((1 - p(k)) .^ x(k));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/geometric_rnd.m b/scripts/statistics/distributions/geometric_rnd.m
--- a/scripts/statistics/distributions/geometric_rnd.m
+++ b/scripts/statistics/distributions/geometric_rnd.m
@@ -19,12 +19,14 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} geometric_rnd (@var{p}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} geometric_rnd (@var{p}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the
 ## geometric distribution with parameter @var{p}, which 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{p}.
+## If @var{r} and @var{c} are given create a matrix with @var{r} rows and
+## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size
+## @var{sz}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
@@ -39,34 +41,59 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("geometric_rnd: c must be a positive integer");
     endif
-    [retval, p] = common_size (p, zeros (r, c));
-    if (retval > 0)
-      error ("geometric_rnd: p must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+
+    if (any (size (p) != 1) && ((length (size (p)) != length (sz)) ||
+				any (size (p) != sz)))
+      error ("geometric_rnd: p must be scalar or of size [r, c]");
     endif
+  elseif (nargin == 2)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("geometric_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (p) != 1) && ((length (size (p)) != length (sz)) ||
+				any (size (p) != sz)))
+      error ("geometric_rnd: n must be scalar or of size sz");
+    endif
+  elseif (nargin == 1)
+    sz = size(n);
   elseif (nargin != 1)
     usage ("geometric_rnd (p, r, c)");
   endif
 
-  [r, c] = size (p);
-  s = r * c;
-  p = reshape (p, 1, s);
-  rnd = zeros (1, s);
+
+  if (isscalar (p))
+    if (!(p >= 0) || !(p <= 1))
+      rnd = NaN * ones (sz);
+    elseif (p == 0)
+      rnd = Inf * ones (sz);
+    elseif ((p > 0) & (p < 1));
+      rnd = floor (log (rand (sz)) / log (1 - p));
+    else
+      rnd = zeros (sz);
+    endif
+  else
+    rnd = zeros (sz);
 
-  k = find (!(p >= 0) | !(p <= 1));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(p >= 0) | !(p <= 1));
+    if (any (k))
+      rnd(k) = NaN * ones (1, length (k));
+    endif
+
+    k = find (p == 0);
+    if (any (k))
+      rnd(k) = Inf * ones (1, length (k));
+    endif
+
+    k = find ((p > 0) & (p < 1));
+    if (any (k))
+      rnd(k) = floor (log (rand (size (k))) ./ log (1 - p(k)));
+    endif
   endif
 
-  k = find (p == 0);
-  if (any (k))
-    rnd(k) = Inf * ones (1, length (k));
-  endif
-
-  k = find ((p > 0) & (p < 1));
-  if (any (k))
-    rnd(k) = floor (log (rand (1, length (k))) ./ log (1 - p(k)));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/hypergeometric_cdf.m b/scripts/statistics/distributions/hypergeometric_cdf.m
--- a/scripts/statistics/distributions/hypergeometric_cdf.m
+++ b/scripts/statistics/distributions/hypergeometric_cdf.m
@@ -36,7 +36,11 @@
 function cdf = hypergeometric_cdf (x, m, t, n)
 
   if (nargin != 4)
-    usage ("hypergeometrix_cdf (x, m, t, n)");
+    usage ("hypergeometric_cdf (x, m, t, n)");
+  endif
+
+  if (!isscalar (m) || !isscalar (t) || !isscalar (n))
+    error ("hypergeometric_cdf: m, t and n must all be positive integers");
   endif
 
   if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) |
diff --git a/scripts/statistics/distributions/hypergeometric_inv.m b/scripts/statistics/distributions/hypergeometric_inv.m
--- a/scripts/statistics/distributions/hypergeometric_inv.m
+++ b/scripts/statistics/distributions/hypergeometric_inv.m
@@ -36,6 +36,10 @@
     usage ("hypergeometric_inv (x, m, t, n)");
   endif
 
+  if (!isscalar (m) || !isscalar (t) || !isscalar (n))
+    error ("hypergeometrix_inv: m, t and n must all be positive integers");
+  endif
+
   if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) |
       (t != round (t)) | (n != round (n)) | (m > t) | (n > t))
     inv = NaN * ones (size (x))
diff --git a/scripts/statistics/distributions/hypergeometric_pdf.m b/scripts/statistics/distributions/hypergeometric_pdf.m
--- a/scripts/statistics/distributions/hypergeometric_pdf.m
+++ b/scripts/statistics/distributions/hypergeometric_pdf.m
@@ -37,18 +37,15 @@
     usage ("hypergeometric_pdf (x, m, t, n)");
   endif
 
-  [retval, x, m, t, n] = common_size (x, m, t, n);
-  if (retval > 0)
-    error ("hypergeometric_pdf: x, m, t, and n must be of common size or scalar");
+  if (!isscalar (m) || !isscalar (t) || !isscalar (n))
+    [retval, x, m, t, n] = common_size (x, m, t, n);
+    if (retval > 0)
+      error ("hypergeometric_pdf: x, m, t, and n 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);
-  t = reshape (t, 1, s);
-  n = reshape (n, 1, s);
-  pdf = zeros * ones (1, s);
+  pdf = zeros (size (x));
+
   ## everything in i1 gives NaN
   i1 = ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) |
         (t != round (t)) | (n != round (n)) | (m > t) | (n > t));
@@ -56,12 +53,21 @@
   i2 = ((x != round (x)) | (x < 0) | (x > m) | (n < x) | (n-x > t-m));
   k = find (i1);
   if (any (k))
-    pdf (k) = NaN * ones (size (k));
+    if (isscalar (m) && isscalar (t) && isscalar (n))
+      pdf = NaN * ones ( size (x));
+    else
+      pdf (k) = NaN;
+    endif
   endif
   k = find (!i1 & !i2);
   if (any (k))
-    pdf (k) = (bincoeff (m(k), x(k)) .* bincoeff (t(k)-m(k), n(k)-x(k))
-               ./ bincoeff (t(k), n(k)));
+    if (isscalar (m) && isscalar (t) && isscalar (n))
+      pdf (k) = (bincoeff (m, x(k)) .* bincoeff (t-m, n-x(k))
+		 / bincoeff (t, n));
+    else
+      pdf (k) = (bincoeff (m(k), x(k)) .* bincoeff (t(k)-m(k), n(k)-x(k))
+		 ./ bincoeff (t(k), n(k)));
+    endif
   endif
 
 endfunction
diff --git a/scripts/statistics/distributions/hypergeometric_rnd.m b/scripts/statistics/distributions/hypergeometric_rnd.m
--- a/scripts/statistics/distributions/hypergeometric_rnd.m
+++ b/scripts/statistics/distributions/hypergeometric_rnd.m
@@ -19,25 +19,61 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} hypergeometric_rnd (@var{n_size}, @var{m}, @var{t}, @var{n})
-## Generate a row vector containing a random sample of size @var{n_size} from
-## the hypergeometric distribution with parameters @var{m}, @var{t}, and
-## @var{n}.
+## @deftypefnx {Function File} {} hypergeometric_rnd (@var{@var{m}, @var{t}, @var{n}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} hypergeometric_rnd (@var{@var{m}, @var{t}, @var{n}, @var{sz})
+## Generate a row vector containing a random sample of size @var{n_size}
+## from the hypergeometric distribution with parameters @var{m}, @var{t},
+## and @var{n}.
+##
+## If  @var{r} and @var{c} are given create a matrix with @var{r} rows and
+## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size
+## @var{sz}.
 ##
 ## The parameters @var{m}, @var{t}, and @var{n} must positive integers
 ## with @var{m} and @var{n} not greater than @var{t}.
 ## @end deftypefn
 
-function rnd = hypergeometric_rnd (N, m, t, n)
+## function rnd = hypergeometric_rnd (N, m, t, n)
+function rnd = hypergeometric_rnd (m, t, n, r, c)
 
-  if (nargin != 4)
-    usage ("hypergeometric_rnd (N, m, t, n)");
+  if (nargin == 5)
+    if (! (isscalar (r) && (r > 0) && (r == round (r))))
+      error ("hypergeometric_rnd: r must be a positive integer");
+    endif
+    if (! (isscalar (c) && (c > 0) && (c == round (c))))
+      error ("hypergeometric_rnd: c must be a positive integer");
+    endif
+    sz = [r, c];
+  elseif (nargin == 4)
+    ## A potential problem happens here if all args are scalar, as
+    ## we can distiguish between the command syntax. This is quite
+    ## ambigous! I assume that if the last arg is a vector then 
+    ## then third form is assumed. This means that you can't define
+    ## and r-by-r matrix with a single scalar!
+
+    if (isscalar (r))
+      sz = [1, floor(m)];
+      m = t;
+      t = n;
+      n = r;
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("hypergeometric_rnd: r must be a vector of positive integers");
+    endif
+  else
+    usage ("hypergeometric_rnd (N, m, t, n) | hypergeometric_rnd (m, t, n, r, c)");
+  endif
+
+  if (!isscalar (m) || !isscalar (t) || !isscalar (n))
+    error ("hypergeometric_cdf: m, t and n must all be positive integers");
   endif
 
   if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) |
       (t != round (t)) | (n != round (n)) | (m > t) | (n > t))
-    rnd = NaN * ones (1, N)
+    rnd = NaN * ones (sz)
   else
-    rnd = discrete_rnd (N, 0 : n, hypergeometric_pdf (0 : n, m, t, n));
+    rnd = discrete_rnd (0 : n, hypergeometric_pdf (0 : n, m, t, n), sz);
   endif
 
 endfunction
diff --git a/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m b/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m
--- a/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m
+++ b/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m
@@ -57,17 +57,20 @@
     endif
   endif
 
-  [nr, nc] = size (x);
-  if (min (nr, nc) == 0)
+  n = numel (x);
+  if (n == 0)
     error ("kolmogorov_smirnov_cdf: x must not be empty");
   endif
 
-  n   = nr * nc;
-  x   = reshape (x, 1, n);
-  cdf = zeros (1, n);
+  cdf = zeros (size (x));
+
   ind = find (x > 0);
   if (length (ind) > 0)
-    y   = x(ind);
+    if (size(ind,2) < size(ind,1))
+      y = x(ind.');
+    else
+      y   = x(ind);
+    endif
     K   = ceil (sqrt (- log (tol) / 2) / min (y));
     k   = (1:K)';
     A   = exp (- 2 * k.^2 * y.^2);
@@ -76,6 +79,4 @@
     cdf(ind) = 1 + 2 * sum (A);
   endif
 
-  cdf = reshape (cdf, nr, nc);
-
 endfunction
diff --git a/scripts/statistics/distributions/laplace_cdf.m b/scripts/statistics/distributions/laplace_cdf.m
--- a/scripts/statistics/distributions/laplace_cdf.m
+++ b/scripts/statistics/distributions/laplace_cdf.m
@@ -32,19 +32,16 @@
     usage ("laplace_cdf (x)");
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  cdf = zeros (1, s);
+  cdf = zeros (size (x));
 
   k = find (isnan (x));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find (x == Inf);
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x > -Inf) & (x < Inf));
@@ -52,6 +49,4 @@
     cdf(k) = (1 + sign (x(k)) .* (1 - exp (- abs (x(k))))) / 2;
   endif
 
-  cdf = reshape (cdf, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/laplace_inv.m b/scripts/statistics/distributions/laplace_inv.m
--- a/scripts/statistics/distributions/laplace_inv.m
+++ b/scripts/statistics/distributions/laplace_inv.m
@@ -32,19 +32,16 @@
     usage ("laplace_inv (x)");
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  inv = (-Inf) * ones (1, s);
+  inv = (-Inf) * ones (size (x));
 
   k = find (isnan (x) | (x < 0) | (x > 1));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find (x == 1);
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1));
@@ -53,6 +50,4 @@
 	      - (x(k) > 1/2) .* log (2 * (1 - x(k))));
   endif
 
-  inv = reshape (inv, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/laplace_pdf.m b/scripts/statistics/distributions/laplace_pdf.m
--- a/scripts/statistics/distributions/laplace_pdf.m
+++ b/scripts/statistics/distributions/laplace_pdf.m
@@ -32,14 +32,11 @@
     usage ("laplace_pdf (x)");
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  pdf = zeros (1, s);
+  pdf = zeros (size (x));
 
   k = find (isnan (x));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find ((x > -Inf) & (x < Inf));
@@ -47,6 +44,4 @@
     pdf(k) = exp (- abs (x(k))) / 2;
   endif
 
-  pdf = reshape (pdf, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/laplace_rnd.m b/scripts/statistics/distributions/laplace_rnd.m
--- a/scripts/statistics/distributions/laplace_rnd.m
+++ b/scripts/statistics/distributions/laplace_rnd.m
@@ -19,8 +19,10 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} laplace_rnd (@var{r}, @var{c})
+## @deftypefnx {Function File} {} laplace_rnd (@var{sz});
 ## Return an @var{r} by @var{c} matrix of random numbers from the
-## Laplace distribution.
+## Laplace distribution. Or is @var{sz} is a vector, create a matrix of
+## @var{sz}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
@@ -28,18 +30,27 @@
 
 function rnd = laplace_rnd (r, c)
 
-  if (nargin != 2)
+  if (nargin == 2)
+    if (! (isscalar (r) && (r > 0) && (r == round (r))))
+      error ("laplace_rnd: r must be a positive integer");
+    endif
+    if (! (isscalar (c) && (c > 0) && (c == round (c))))
+      error ("laplace_rnd: c must be a positive integer");
+    endif
+    sz = [r, c];
+  elseif (nargin == 1)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("laplace_rnd: r must be a postive integer or vector");
+    endif
+  else
     usage ("laplace_rnd (r, c)");
   endif
 
-  if (! (isscalar (r) && (r > 0) && (r == round (r))))
-    error ("laplace_rnd: r must be a positive integer");
-  endif
-  if (! (isscalar (c) && (c > 0) && (c == round (c))))
-    error ("laplace_rnd: c must be a positive integer");
-  endif
-
-  tmp = rand (r, c);
+  tmp = rand (sz);
   rnd = ((tmp < 1/2) .* log (2 * tmp)
          - (tmp > 1/2) .* log (2 * (1 - tmp)));
 
diff --git a/scripts/statistics/distributions/logistic_inv.m b/scripts/statistics/distributions/logistic_inv.m
--- a/scripts/statistics/distributions/logistic_inv.m
+++ b/scripts/statistics/distributions/logistic_inv.m
@@ -32,24 +32,21 @@
     usage ("logistic_inv (x)");
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  inv = zeros (1, s);
+  inv = zeros (size (x));
 
   k = find ((x < 0) | (x > 1) | isnan (x));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find (x == 0);
   if (any (k))
-    inv(k) = (-Inf) * ones (1, length (k));
+    inv(k) = -Inf;
   endif
 
   k = find (x == 1);
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1));
@@ -57,6 +54,4 @@
     inv (k) = - log (1 ./ x(k) - 1);
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/logistic_rnd.m b/scripts/statistics/distributions/logistic_rnd.m
--- a/scripts/statistics/distributions/logistic_rnd.m
+++ b/scripts/statistics/distributions/logistic_rnd.m
@@ -19,8 +19,10 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} logistic_rnd (@var{r}, @var{c})
+## @deftypefn {Function File} {} logistic_rnd (@var{sz})
 ## Return an @var{r} by @var{c} matrix of random numbers from the
-## logistic distribution.
+## logistic distribution. Or is @var{sz} is a vector, create a matrix of
+## @var{sz}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
@@ -28,17 +30,27 @@
 
 function rnd = logistic_rnd (r, c)
 
-  if (nargin != 2)
+
+  if (nargin == 2)
+    if (! (isscalar (r) && (r > 0) && (r == round (r))))
+      error ("logistic_rnd: r must be a positive integer");
+    endif
+    if (! (isscalar (c) && (c > 0) && (c == round (c))))
+      error ("logistic_rnd: c must be a positive integer");
+    endif
+    sz = [r, c];
+  elseif (nargin == 1)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("logistic_rnd: r must be a postive integer or vector");
+    endif
+  else
     usage ("logistic_rnd (r, c)");
   endif
 
-  if (! (isscalar (r) && (r > 0) && (r == round (r))))
-    error ("logistic_rnd: r must be a positive integer");
-  endif
-  if (! (isscalar (c) && (c > 0) && (c == round (c))))
-    error ("logistic_rnd: c must be a positive integer");
-  endif
-
-  rnd = - log (1 ./ rand (r, c) - 1);
+  rnd = - log (1 ./ rand (sz) - 1);
 
 endfunction
diff --git a/scripts/statistics/distributions/lognormal_cdf.m b/scripts/statistics/distributions/lognormal_cdf.m
--- a/scripts/statistics/distributions/lognormal_cdf.m
+++ b/scripts/statistics/distributions/lognormal_cdf.m
@@ -47,33 +47,32 @@
   ## cdf = normal_cdf (log (x), log (a), v);
   ## Hence ...
 
-  [retval, x, a, v] = common_size (x, a, v);
-  if (retval > 0)
-    error ("lognormal_cdf: x, a and v must be of common size or scalars");
+  if (!isscalar (a) || !isscalar (v))
+    [retval, x, a, v] = common_size (x, a, v);
+    if (retval > 0)
+      error ("lognormal_cdf: x, a and v 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);
-  v = reshape (v, 1, s);
-  cdf = zeros (1, s);
+  cdf = zeros (size (x));
 
   k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf));
   if (any (k))
-    cdf(k) = stdnormal_cdf ((log (x(k)) - log (a(k))) ./ sqrt (v(k)));
+    if (isscalar (a) && isscalar (v))
+      cdf(k) = stdnormal_cdf ((log (x(k)) - log (a)) / sqrt (v));
+    else
+      cdf(k) = stdnormal_cdf ((log (x(k)) - log (a(k))) ./ sqrt (v(k)));
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/lognormal_inv.m b/scripts/statistics/distributions/lognormal_inv.m
--- a/scripts/statistics/distributions/lognormal_inv.m
+++ b/scripts/statistics/distributions/lognormal_inv.m
@@ -47,34 +47,33 @@
   ## inv = exp (normal_inv (x, log (a), v));
   ## Hence ...
 
-  [retval, x, a, v] = common_size (x, a, v);
-  if (retval > 0)
-    error ("lognormal_inv: x, a and v must be of common size or scalars");
+  if (!isscalar (a) || !isscalar (v))
+    [retval, x, a, v] = common_size (x, a, v);
+    if (retval > 0)
+      error ("lognormal_inv: x, a and v 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);
-  v = reshape (v, 1, s);
-  inv = zeros (1, s);
+  inv = zeros (size (x));
 
   k = find (!(x >= 0) | !(x <= 1) | !(a > 0) | !(a < Inf)
 	    | !(v > 0) | !(v < Inf));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf));
   if (any (k))
-    inv(k) = a(k) .* exp (sqrt (v(k)) .* stdnormal_inv (x(k)));
+    if (isscalar (a) && isscalar (v))
+      inv(k) = a .* exp (sqrt (v) .* stdnormal_inv (x(k)));
+    else
+      inv(k) = a(k) .* exp (sqrt (v(k)) .* stdnormal_inv (x(k)));
+    endif
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/lognormal_pdf.m b/scripts/statistics/distributions/lognormal_pdf.m
--- a/scripts/statistics/distributions/lognormal_pdf.m
+++ b/scripts/statistics/distributions/lognormal_pdf.m
@@ -47,28 +47,27 @@
   ## pdf = (x > 0) ./ x .* normal_pdf (log (x), log (a), v);
   ## Hence ...
 
-  [retval, x, a, v] = common_size (x, a, v);
-  if (retval > 0)
-    error ("lognormal_pdf: x, a and v must be of common size or scalars");
+  if (!isscalar (a) || !isscalar (v))
+    [retval, x, a, v] = common_size (x, a, v);
+    if (retval > 0)
+      error ("lognormal_pdf: x, a and v 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);
-  v = reshape (v, 1, s);
-  pdf = zeros (1, s);
+  pdf = zeros (size (x));
 
   k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf));
   if (any (k))
-    pdf(k) = normal_pdf (log (x(k)), log (a(k)), v(k)) ./ x(k);
+    if (isscalar (a) && isscalar (v))
+      pdf(k) = normal_pdf (log (x(k)), log (a), v) ./ x(k);
+    else
+      pdf(k) = normal_pdf (log (x(k)), log (a(k)), v(k)) ./ x(k);
+  endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/lognormal_rnd.m b/scripts/statistics/distributions/lognormal_rnd.m
--- a/scripts/statistics/distributions/lognormal_rnd.m
+++ b/scripts/statistics/distributions/lognormal_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} lognormal_rnd (@var{a}, @var{v}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} lognormal_rnd (@var{a}, @var{v}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the
 ## lognormal distribution with parameters @var{a} and @var{v}. Both
 ## @var{a} and @var{v} must be scalar or of size @var{r} by @var{c}.
+## Or if @var{sz} is a vector, create a matrix of size @var{sz}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{a} and @var{v}.
@@ -32,6 +34,15 @@
 
 function rnd = lognormal_rnd (a, v, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(a) || !isscalar(v)) 
+      [retval, a, v] = common_size (a, v);
+      if (retval > 0)
+	error ("lognormal_rnd: a and v must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("lognormal_rnd: r must be a positive integer");
@@ -39,35 +50,51 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("lognormal_rnd: c must be a positive integer");
     endif
-    [retval, a, v] = common_size (a, v, zeros (r, c));
-    if (retval > 0)
-      error ("lognormal_rnd: a and v must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+
+    if (any (size (a) != 1) && 
+	((length (size (a)) != length (sz)) || any (size (a) != sz)))
+      error ("lognormal_rnd: a and b must be scalar or of size [r, c]");
+    endif
+
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("lognormal_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (a) != 1) && 
+	((length (size (a)) != length (sz)) || any (size (a) != sz)))
+      error ("lognormal_rnd: a and b must be scalar or of size sz");
     endif
   elseif (nargin == 2)
-    [retval, a, v] = common_size (a, v);
-    if (retval > 0)
-      error ("lognormal_rnd: a and v must be of common size or scalar");
-    endif
+    sz = size(a);
   else
     usage ("lognormal_rnd (a, v, r, c)");
   endif
 
-  [r, c] = size (a);
-  s = r * c;
-  a = reshape (a, 1, s);
-  v = reshape (v, 1, s);
-  rnd = zeros (1, s);
+  if (isscalar (a) && isscalar (v))
+    if  (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf))
+      rnd = NaN * ones (sz);
+    elseif find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf));
+      rnd = a * exp (sqrt (v) .* randn (sz));
+    else
+      rnd = zeros (sz);
+    endif
+  else
+    rnd = zeros (sz);
+    k = find (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf));
+    if (any (k))
+      rnd(k) = NaN * ones (1, length (k));
+    endif
 
-  k = find (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf));
+    if (any (k))
+      rnd(k) = a(k) .* exp (sqrt (v(k)) .* randn (1, length (k)));
+    endif
   endif
 
-  k = find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf));
-  if (any (k))
-    rnd(k) = a(k) .* exp (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/pascal_cdf.m b/scripts/statistics/distributions/pascal_cdf.m
--- a/scripts/statistics/distributions/pascal_cdf.m
+++ b/scripts/statistics/distributions/pascal_cdf.m
@@ -36,51 +36,58 @@
     usage ("pascal_cdf (x, n, p)");
   endif
 
-  [retval, x, n, p] = common_size (x, n, p);
-  if (retval > 0)
-    error ("pascal_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 ("pascal_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);
+  
+  cdf = zeros (size (x));
 
   k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n))
 	    | (p < 0) | (p > 1));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (n > 0) & (n < Inf) & (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 < Inf) & (x == round (x)) & (n > 0)
 	    & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1));
   if (any (k))
     ## Does anyone know a better way to do the summation?
-    m = zeros (1, length (k));
+    m = zeros (size (k));
     x = floor (x(k));
-    n = n(k);
-    p = p(k);
     y = cdf(k);
-    while (1)
-      l = find (m <= x);
-      if (any (l))
-        y(l) = y(l) + pascal_pdf (m(l), n(l), p(l));
-        m(l) = m(l) + 1;
-      else
-        break;
-      endif
-    endwhile
+    if (isscalar (n) && isscalar (p))
+      while (1)
+	l = find (m <= x);
+	if (any (l))
+          y(l) = y(l) + pascal_pdf (m(l), n, p);
+          m(l) = m(l) + 1;
+	else
+          break;
+	endif
+      endwhile
+    else
+      n = n(k);
+      p = p(k);
+      while (1)
+	l = find (m <= x);
+	if (any (l))
+          y(l) = y(l) + pascal_pdf (m(l), n(l), p(l));
+          m(l) = m(l) + 1;
+	else
+          break;
+	endif
+      endwhile
+    endif
     cdf(k) = y;
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/pascal_inv.m b/scripts/statistics/distributions/pascal_inv.m
--- a/scripts/statistics/distributions/pascal_inv.m
+++ b/scripts/statistics/distributions/pascal_inv.m
@@ -37,50 +37,58 @@
     usage ("pascal_inv (x, n, p)");
   endif
 
-  [retval, x, n, p] = common_size (x, n, p);
-  if (retval > 0)
-    error ("pascal_inv: 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 ("pascal_inv: 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);
-  inv = zeros (1, s);
+  inv = zeros (size (x));
 
   k = find (isnan (x) | (x < 0) | (x > 1) | (n < 1) | (n == Inf)
 	    | (n != round (n)) | (p < 0) | (p > 1));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 1) & (n > 0) & (n < Inf) & (n == round (n))
 	    & (p >= 0) & (p <= 1));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x >= 0) & (x < 1) & (n > 0) & (n < Inf)
 	    & (n == round (n)) & (p > 0) & (p <= 1));
   if (any (k))
+    m = zeros (size (k));
     x = x(k);
-    n = n(k);
-    p = p(k);
-    m = zeros (1, length (k));
-    s = p .^ n;
-    while (1)
-      l = find (s < x);
-      if (any (l))
-        m(l) = m(l) + 1;
-        s(l) = s(l) + pascal_pdf (m(l), n(l), p(l));
-      else
-        break;
-      endif
-    endwhile
+    if (isscalar (n) && isscalar (p))
+      s = p ^ n * ones (size(k));
+      while (1)
+	l = find (s < x);
+	if (any (l))
+          m(l) = m(l) + 1;
+          s(l) = s(l) + pascal_pdf (m(l), n, p);
+	else
+          break;
+	endif
+      endwhile
+    else
+      n = n(k);
+      p = p(k);
+      s = p .^ n;
+      while (1)
+	l = find (s < x);
+	if (any (l))
+          m(l) = m(l) + 1;
+          s(l) = s(l) + pascal_pdf (m(l), n(l), p(l));
+	else
+          break;
+	endif
+      endwhile
+    endif
     inv(k) = m;
   endif
 
-  inv = reshape (inv, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/pascal_pdf.m b/scripts/statistics/distributions/pascal_pdf.m
--- a/scripts/statistics/distributions/pascal_pdf.m
+++ b/scripts/statistics/distributions/pascal_pdf.m
@@ -37,37 +37,36 @@
     usage ("pascal_pdf (x, n, p)");
   endif
 
-  [retval, x, n, p] = common_size (x, n, p);
-  if (retval > 0)
-    error ("pascal_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 ("pascal_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);
+  pdf = zeros (size (x));
 
   k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n))
 	    | (p < 0) | (p > 1));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   ## Just for the fun of it ...
   k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n))
 	    & (p == 0));
   if (any (k))
-    pdf(k) = ones (1, length (k));
+    pdf(k) = 1;
   endif
 
   k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0)
 	    & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1));
   if (any (k))
-    pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k));
+    if (isscalar (n) && isscalar (p))
+      pdf(k) = bincoeff (-n, x(k)) .* (p ^ n) .* ((p - 1) .^ x(k));
+    else
+      pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/pascal_rnd.m b/scripts/statistics/distributions/pascal_rnd.m
--- a/scripts/statistics/distributions/pascal_rnd.m
+++ b/scripts/statistics/distributions/pascal_rnd.m
@@ -19,12 +19,14 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} pascal_rnd (@var{n}, @var{p}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} pascal_rnd (@var{n}, @var{p}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the Pascal
 ## (negative 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}. 
+## 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}.
+## the common size of @var{n} and @var{p}. Or if @var{sz} is a vector, 
+## create a matrix of size @var{sz}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
@@ -32,6 +34,15 @@
 
 function rnd = pascal_rnd (n, p, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(n) || !isscalar(p)) 
+      [retval, n, p] = common_size (n, p);
+      if (retval > 0)
+	error ("pascal_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 ("pascal_rnd: r must be a positive integer");
@@ -39,39 +50,71 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("pascal_rnd: c must be a positive integer");
     endif
-    [retval, n, p] = common_size (n, p, zeros (r, c));
-    if (retval > 0)
-      error ("pascal_rnd: n and p must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+
+    if (any (size (n) != 1) && 
+	((length (size (n)) != length (sz)) || any (size (n) != sz)))
+      error ("pascal_rnd: n and p must be scalar or of size [r, c]");
+    endif
+
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("pascal_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (n) != 1) && 
+	((length (size (n)) != length (sz)) || any (size (n) != sz)))
+      error ("pascal_rnd: n and p must be scalar or of size sz");
     endif
   elseif (nargin == 2)
-    [retval, n, p] = common_size (n, p);
-    if (retval > 0)
-      error ("pascal_rnd: n and p must be of common size or scalar");
-    endif
+    sz = size(n);
   else
     usage ("pascal_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 ((n < 1) || (n == Inf) || (n != round (n)) || (p <= 0) || (p > 1));
+      rnd = NaN * ones (sz)
+    elseif ((n > 0) && (n < Inf) && (n == round (n)) && 
+	    (p > 0) && (p <= 1))
+      L = prod (sz);
+      tmp = floor (log (rand (n, L)) / log (1 - p));
+      if (n == 1)
+	rnd = tmp;
+      else
+	ind = (1 : n)' * ones (1, L);
+	rnd = sum (tmp .* (ind <= n));
+      endif
+    else
+      rnd = zeros (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 < 1) || (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 = reshape (n, 1, prod (sz));
+      p = reshape (p, 1, prod (sz));
+      N = max (n(k));
+      L = length (k);
+      tmp = floor (log (rand (N, L)) ./ (ones (N, 1) * log (1 - p(k))));
+      if (N == 1)
+	rnd(k) = tmp;
+      else
+	ind = (1 : N)' * ones (1, L);
+	rnd(k) = sum (tmp .* (ind <= ones (N, 1) * n(k)));
+      endif
+    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 = floor (log (rand (N, L)) ./ (ones (N, 1) * log (1 - p(k))));
-    ind = (1 : N)' * ones (1, L);
-    rnd(k) = sum (tmp .* (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/poisson_cdf.m b/scripts/statistics/distributions/poisson_cdf.m
--- a/scripts/statistics/distributions/poisson_cdf.m
+++ b/scripts/statistics/distributions/poisson_cdf.m
@@ -33,32 +33,32 @@
     usage ("poisson_cdf (x, lambda)");
   endif
 
-  [retval, x, l] = common_size (x, l);
-  if (retval > 0)
-    error ("poisson_cdf: x and lambda must be of common size or scalar");
+  if (!isscalar (l))
+    [retval, x, l] = common_size (x, l);
+    if (retval > 0)
+      error ("poisson_cdf: x and lambda must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  l = reshape (l, 1, s);
-  cdf = zeros (1, s);
+  cdf = zeros (size (x));
 
   k = find (isnan (x) | !(l > 0));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (l > 0));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x >= 0) & (x < Inf) & (l > 0));
   if (any (k))
-    cdf(k) = 1 - gammainc (l(k), floor (x(k)) + 1);
+    if (isscalar (l))
+      cdf(k) = 1 - gammainc (l, floor (x(k)) + 1);
+    else
+      cdf(k) = 1 - gammainc (l(k), floor (x(k)) + 1);
+    endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/poisson_inv.m b/scripts/statistics/distributions/poisson_inv.m
--- a/scripts/statistics/distributions/poisson_inv.m
+++ b/scripts/statistics/distributions/poisson_inv.m
@@ -33,41 +33,45 @@
     usage ("poisson_inv (x, lambda)");
   endif
 
-  [retval, x, l] = common_size (x, l);
-  if (retval > 0)
-    error ("poisson_inv: x and lambda must be of common size or scalar");
+  if (!isscalar (l))
+    [retval, x, l] = common_size (x, l);
+    if (retval > 0)
+      error ("poisson_inv: x and lambda must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  l = reshape (l, 1, s);
-  inv = zeros (1, s);
+  inv = zeros (size (x));
 
   k = find ((x < 0) | (x > 1) | isnan (x) | !(l > 0));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 1) & (l > 0));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (l > 0));
   if (any (k))
-    cdf = exp (-l(k));
+    if (isscalar (l))
+      cdf = exp (-l) * ones (size (k));
+    else
+      cdf = exp (-l(k));
+    endif
     while (1)
       m = find (cdf < x(k));
       if (any (m))
         inv(k(m)) = inv(k(m)) + 1;
-        cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l(k(m)));
+	if (isscalar (l))
+          cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l);
+	else
+          cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l(k(m)));
+	endif
       else
         break;
       endif
     endwhile
   endif
 
-  inv = reshape (inv, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/poisson_pdf.m b/scripts/statistics/distributions/poisson_pdf.m
--- a/scripts/statistics/distributions/poisson_pdf.m
+++ b/scripts/statistics/distributions/poisson_pdf.m
@@ -32,27 +32,27 @@
     usage ("poisson_pdf (x, lambda)");
   endif
 
-  [retval, x, l] = common_size (x, l);
-  if (retval > 0)
-    error ("poisson_pdf: x and lambda must be of common size or scalar");
+  if (!isscalar (l))
+    [retval, x, l] = common_size (x, l);
+    if (retval > 0)
+      error ("poisson_pdf: x and lambda must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  l = reshape (l, 1, s);
-  pdf = zeros (1, s);
+  pdf = zeros (size (x));
 
   k = find (!(l > 0) | isnan (x));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   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));
+    if (isscalar (l))
+      pdf(k) = exp (x(k) .* log (l) - l - gammaln (x(k) + 1));
+    else
+      pdf(k) = exp (x(k) .* log (l(k)) - l(k) - gammaln (x(k) + 1));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/poisson_rnd.m b/scripts/statistics/distributions/poisson_rnd.m
--- a/scripts/statistics/distributions/poisson_rnd.m
+++ b/scripts/statistics/distributions/poisson_rnd.m
@@ -39,43 +39,76 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("poisson_rnd: c must be a positive integer");
     endif
-    [retval, l] = common_size (l, zeros (r, c));
-    if (retval > 0)
-      error ("poisson_rnd: lambda must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+
+    if (any (size (l) != 1) && 
+	((length (size (l)) != length (sz)) || any (size (l) != sz)))
+      error ("poisson_rnd: lambda must be scalar or of size [r, c]");
     endif
-  elseif (nargin != 1)
+  elseif (nargin == 2)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("poisson_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (l) != 1) && 
+	((length (size (l)) != length (sz)) || any (size (l) != sz)))
+      error ("poisson_rnd: lambda must be scalar or of size sz");
+    endif
+  elseif (nargin == 1)
+    sz = size (l);
+  else
     usage ("poisson_rnd (lambda, r, c)");
   endif
 
-  [r, c] = size (l);
-  s = r * c;
-  l = reshape (l, 1, s);
-  rnd = zeros (1, s);
+  if (isscalar (l))
+
+    if (!(l > 0) | !(l < Inf))
+      rnd = NaN * ones (sz);
+    elseif ((l > 0) & (l < Inf))
+      num = zeros (sz);
+      sum = - log (1 - rand (sz)) ./ l;
+      while (1)
+	ind = find (sum < 1);
+	if (any (ind))
+          sum(ind) = (sum(ind) - log (1 - rand (size (ind))) / l);
+          num(ind) = num(ind) + 1;
+	else
+          break;
+	endif
+      endwhile
+      rnd = num;
+    else
+      rnd = zeros (sz);
+    endif
+  else
+    rnd = zeros (sz);
 
-  k = find (!(l > 0) | !(l < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(l > 0) | !(l < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+
+    k = find ((l > 0) & (l < Inf));
+    if (any (k))
+      l = l(k);
+      num = zeros (size (k));
+      sum = - log (1 - rand (size (k))) ./ l;
+      while (1)
+	ind = find (sum < 1);
+	if (any (ind))
+          sum(ind) = (sum(ind)
+                      - log (1 - rand (1, length (ind))) ./ l(ind));
+          num(ind) = num(ind) + 1;
+	else
+          break;
+	endif
+      endwhile
+      rnd(k) = num;
+    endif
   endif
 
-  k = find ((l > 0) & (l < Inf));
-  if (any (k))
-    l = l(k);
-    len = length (k);
-    num = zeros (1, len);
-    sum = - log (1 - rand (1, len)) ./ l;
-    while (1)
-      ind = find (sum < 1);
-      if (any (ind))
-        sum(ind) = (sum(ind)
-                    - log (1 - rand (1, length (ind))) ./ l(ind));
-        num(ind) = num(ind) + 1;
-      else
-        break;
-      endif
-    endwhile
-    rnd(k) = num;
-  endif
-
-  rnd = reshape (rnd, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/t_cdf.m b/scripts/statistics/distributions/t_cdf.m
--- a/scripts/statistics/distributions/t_cdf.m
+++ b/scripts/statistics/distributions/t_cdf.m
@@ -33,36 +33,36 @@
     usage ("t_cdf (x, n)");
   endif
 
-  [retval, x, n] = common_size (x, n);
-  if (retval > 0)
-    error ("t_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 ("t_cdf: x and n 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);
-  cdf = zeros (1, s);
+  cdf = zeros (size (x));
 
   k = find (isnan (x) | !(n > 0));
   if (any (k))
-    cdf(k) = NaN * ones (1, length (k));
+    cdf(k) = NaN;
   endif
 
   k = find ((x == Inf) & (n > 0));
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
   k = find ((x > -Inf) & (x < Inf) & (n > 0));
   if (any (k))
-    cdf(k) = betainc (1 ./ (1 + x(k) .^ 2 ./ n(k)), n(k) / 2, 1 / 2) / 2;
+    if (isscalar (n))
+      cdf(k) = betainc (1 ./ (1 + x(k) .^ 2 ./ n), n / 2, 1 / 2) / 2;
+    else
+      cdf(k) = betainc (1 ./ (1 + x(k) .^ 2 ./ n(k)), n(k) / 2, 1 / 2) / 2;
+    endif
     ind = find (x(k) > 0);
     if (any (ind))
       cdf(k(ind)) = 1 - cdf(k(ind));
     endif
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/t_inv.m b/scripts/statistics/distributions/t_inv.m
--- a/scripts/statistics/distributions/t_inv.m
+++ b/scripts/statistics/distributions/t_inv.m
@@ -37,37 +37,41 @@
     usage ("t_inv (x, n)");
   endif
 
-  [retval, x, n] = common_size (x, n);
-  if (retval > 0)
-    error ("t_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 ("t_inv: x and n 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);
-  inv = zeros (1, s);
+  inv = zeros (size (x));
 
   k = find ((x < 0) | (x > 1) | isnan (x) | !(n > 0));
   if (any (k))
-    inv(k) = NaN * ones (1, length (k));
+    inv(k) = NaN;
   endif
 
   k = find ((x == 0) & (n > 0));
   if (any (k))
-    inv(k) = (-Inf) * ones (1, length (k));
+    inv(k) = -Inf;
   endif
 
   k = find ((x == 1) & (n > 0));
   if (any (k))
-    inv(k) = Inf * ones (1, length (k));
+    inv(k) = Inf;
   endif
 
   k = find ((x > 0) & (x < 1) & (n > 0) & (n < 10000));
   if (any (k))
-    inv(k) = (sign (x(k) - 1/2)
-	      .* sqrt (n(k) .* (1 ./ beta_inv (2*min (x(k), 1 - x(k)),
-					       n(k)/2, 1/2) - 1)));
+    if (isscalar (n))
+      inv(k) = (sign (x(k) - 1/2)
+		.* sqrt (n .* (1 ./ beta_inv (2*min (x(k), 1 - x(k)),
+						 n/2, 1/2) - 1)));
+    else
+      inv(k) = (sign (x(k) - 1/2)
+		.* sqrt (n(k) .* (1 ./ beta_inv (2*min (x(k), 1 - x(k)),
+						 n(k)/2, 1/2) - 1)));
+    endif
   endif
 
   ## For large n, use the quantiles of the standard normal
@@ -76,6 +80,4 @@
     inv(k) = stdnormal_inv (x(k));
   endif
 
-  inv = reshape (inv, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/t_pdf.m b/scripts/statistics/distributions/t_pdf.m
--- a/scripts/statistics/distributions/t_pdf.m
+++ b/scripts/statistics/distributions/t_pdf.m
@@ -33,28 +33,29 @@
     usage ("t_pdf (x, n)");
   endif
 
-  [retval, x, n] = common_size (x, n);
-  if (retval > 0)
-    error ("t_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 ("t_pdf: x and n 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);
-  pdf = zeros (1, s);
+  pdf = zeros (size (x));
 
   k = find (isnan (x) | !(n > 0) | !(n < Inf));
   if (any (k))
-    pdf(k) = NaN * ones (1, length (k));
+    pdf(k) = NaN;
   endif
 
   k = find (!isinf (x) & !isnan (x) & (n > 0) & (n < Inf));
   if (any (k))
-    pdf(k) = (exp (- (n(k) + 1) .* log (1 + x(k) .^ 2 ./ n(k))/2)
-	      ./ (sqrt (n(k)) .* beta (n(k)/2, 1/2)));
+    if (isscalar (n))
+      pdf(k) = (exp (- (n + 1) .* log (1 + x(k) .^ 2 ./ n)/2)
+		/ (sqrt (n) * beta (n/2, 1/2)));
+    else
+      pdf(k) = (exp (- (n(k) + 1) .* log (1 + x(k) .^ 2 ./ n(k))/2)
+		./ (sqrt (n(k)) .* beta (n(k)/2, 1/2)));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/t_rnd.m b/scripts/statistics/distributions/t_rnd.m
--- a/scripts/statistics/distributions/t_rnd.m
+++ b/scripts/statistics/distributions/t_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} t_rnd (@var{n}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} t_rnd (@var{n}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the t
 ## (Student) distribution with @var{n} degrees of freedom.  @var{n} must
-## be a scalar or of size @var{r} by @var{c}.
+## be a scalar or of size @var{r} by @var{c}. Or if @va{sz} is a
+## vector create a matrix of size @var{sz}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the size of @var{n}.
@@ -39,29 +41,51 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("t_rnd: c must be a positive integer");
     endif
-    [retval, n] = common_size (n, zeros (r, c));
-    if (retval > 0)
-      error ("t_rnd: n must be scalar or of size %d by %d", r, c);
+    sz = [r, c];
+
+    if (any (size (n) != 1) && 
+	((length (size (n)) != length (sz)) || any (size (n) != sz)))
+      error ("t_rnd: n must be scalar or of size sz");
     endif
-  elseif (nargin != 1)
+  elseif (nargin == 2)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("t_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (n) != 1) && 
+	((length (size (n)) != length (sz)) || any (size (n) != sz)))
+      error ("t_rnd: n must be scalar or of size sz");
+    endif
+  elseif (nargin == 1)
+    sz = size (n);
+  else
     usage ("t_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 (!(n > 0) || !(n < Inf))
+      rnd = NaN * ones (sz);
+    elseif ((n > 0) && (n < Inf))
+      rnd = t_inv (rand (sz), n);
+    else
+      rnd = zeros (size (n));
+    endif
+  else
+    rnd = zeros (size (n));
 
-  k = find (!(n > 0) | !(n < Inf));
-  if (any (k))
-    rnd(k) = NaN * ones (1, length (k));
+    k = find (!(n > 0) | !(n < Inf));
+    if (any (k))
+      rnd(k) = NaN;
+    endif
+
+    k = find ((n > 0) & (n < Inf));
+    if (any (k))
+      rnd(k) = t_inv (rand (size (k)), n(k));
+    endif
   endif
 
-  k = find ((n > 0) & (n < Inf));
-  if (any (k))
-    rnd(k) = t_inv (rand (1, length (k)), n(k));
-  endif
-
-  rnd = reshape (rnd, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/weibull_cdf.m b/scripts/statistics/distributions/weibull_cdf.m
--- a/scripts/statistics/distributions/weibull_cdf.m
+++ b/scripts/statistics/distributions/weibull_cdf.m
@@ -40,36 +40,34 @@
     usage ("weibull_cdf (x, alpha, sigma)");
   endif
 
-  [retval, x, shape, scale] = common_size (x, shape, scale);
-  if (retval > 0)
-    error ("weibull_cdf: x, alpha and sigma must be of common size or scalar");
+  if (!isscalar (shape) || !isscalar (scale))
+    [retval, x, shape, scale] = common_size (x, shape, scale);
+    if (retval > 0)
+      error ("weibull_cdf: x, alpha and sigma must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  shape = reshape (shape, 1, s);
-  scale = reshape (scale, 1, s);
-
-  cdf = NaN * ones (1, s);
+  cdf = NaN * ones (size (x));
 
   ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf));
 
   k = find ((x <= 0) & ok);
   if (any (k))
-    cdf(k) = zeros (1, length (k));
+    cdf(k) = 0;
   endif
 
   k = find ((x > 0) & (x < Inf) & ok);
   if (any (k))
-    cdf(k) = 1 - exp (- (x(k) ./ scale(k)) .^ shape(k));
+    if (isscalar (shape) && isscalar (scale))
+      cdf(k) = 1 - exp (- (x(k) / scale) .^ shape);
+    else
+      cdf(k) = 1 - exp (- (x(k) ./ scale(k)) .^ shape(k));
+    endif
   endif
 
   k = find ((x == Inf) & ok);
   if (any (k))
-    cdf(k) = ones (1, length (k));
+    cdf(k) = 1;
   endif
 
-  cdf = reshape (cdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/weibull_inv.m b/scripts/statistics/distributions/weibull_inv.m
--- a/scripts/statistics/distributions/weibull_inv.m
+++ b/scripts/statistics/distributions/weibull_inv.m
@@ -33,35 +33,34 @@
     usage ("weibull_inv (x, alpha, sigma)");
   endif
 
-  [retval, x, shape, scale] = common_size (x, shape, scale);
-  if (retval > 0)
-    error ("weibull_inv: x, alpha and sigma must be of common size or scalar");
+  if (!isscalar (shape) || !isscalar (scale))
+    [retval, x, shape, scale] = common_size (x, shape, scale);
+    if (retval > 0)
+      error ("weibull_inv: x, alpha and sigma must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  shape = reshape (shape, 1, s);
-  scale = reshape (scale, 1, s);
+  inv = NaN * ones (size (x));
 
-  inv = NaN * ones (1, s);
   ok = ((shape > 0) & (shape < 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) = scale(k) .* (- log (1 - x(k))) .^ (1 ./ shape(k));
+    if (isscalar (shape) && isscalar (scale))
+      inv(k) = scale * (- log (1 - x(k))) .^ (1 / shape);
+    else
+      inv(k) = scale(k) .* (- log (1 - x(k))) .^ (1 ./ shape(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/weibull_pdf.m b/scripts/statistics/distributions/weibull_pdf.m
--- a/scripts/statistics/distributions/weibull_pdf.m
+++ b/scripts/statistics/distributions/weibull_pdf.m
@@ -40,32 +40,32 @@
     usage ("weibull_pdf (x, alpha, sigma)");
   endif
 
-  [retval, x, shape, scale] = common_size (x, shape, scale);
-  if (retval > 0)
-    error ("weibull_pdf: x, alpha and sigma must be of common size or scalar");
+  if (!isscalar (shape) || !isscalar (scale))
+    [retval, x, shape, scale] = common_size (x, shape, scale);
+    if (retval > 0)
+      error ("weibull_pdf: x, alpha and sigma must be of common size or scalar");
+    endif
   endif
 
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  shape = reshape (shape, 1, s);
-  scale = reshape (scale, 1, s);
-
-  pdf = NaN * ones (1, s);
+  pdf = NaN * ones (size (x));
   ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf));
 
   k = find ((x > -Inf) & (x <= 0) & ok);
   if (any (k))
-    pdf(k) = zeros (1, length (k));
+    pdf(k) = 0;
   endif
 
   k = find ((x > 0) & (x < Inf) & ok);
   if (any (k))
-    pdf(k) = (shape(k) .* (scale(k) .^ -shape(k))
-              .* (x(k) .^ (shape(k) - 1))
-              .* exp(- (x(k) ./ scale(k)) .^ shape(k)));
+    if (isscalar (shape) && isscalar (scale))
+      pdf(k) = (shape .* (scale .^ -shape)
+		.* (x(k) .^ (shape - 1))
+		.* exp(- (x(k) / scale) .^ shape));
+    else
+      pdf(k) = (shape(k) .* (scale(k) .^ -shape(k))
+		.* (x(k) .^ (shape(k) - 1))
+		.* exp(- (x(k) ./ scale(k)) .^ shape(k)));
+    endif
   endif
 
-  pdf = reshape (pdf, r, c);
-
 endfunction
diff --git a/scripts/statistics/distributions/weibull_rnd.m b/scripts/statistics/distributions/weibull_rnd.m
--- a/scripts/statistics/distributions/weibull_rnd.m
+++ b/scripts/statistics/distributions/weibull_rnd.m
@@ -19,9 +19,11 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} weibull_rnd (@var{alpha}, @var{sigma}, @var{r}, @var{c})
+## @deftypefnx {Function File} {} weibull_rnd (@var{alpha}, @var{sigma}, @var{sz})
 ## Return an @var{r} by @var{c} matrix of random samples from the
 ## Weibull distribution with parameters @var{alpha} and @var{sigma}
-## which must be scalar or of size @var{r} by @var{c}.
+## which must be scalar or of size @var{r} by @var{c}. Or if @var{sz}
+## is a vector return a matrix of size @var{sz}.
 ##
 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
 ## the common size of @var{alpha} and @var{sigma}.
@@ -32,6 +34,15 @@
 
 function rnd = weibull_rnd (shape, scale, r, c)
 
+  if (nargin > 1)
+    if (!isscalar(shape) || !isscalar(scale)) 
+      [retval, shape, scale] = common_size (shape, scale);
+      if (retval > 0)
+	error ("weibull_rnd: shape and scale must be of common size or scalar");
+      endif
+    endif
+  endif
+
   if (nargin == 4)
     if (! (isscalar (r) && (r > 0) && (r == round (r))))
       error ("weibull_rnd: r must be a positive integer");
@@ -39,32 +50,46 @@
     if (! (isscalar (c) && (c > 0) && (c == round (c))))
       error ("weibull_rnd: c must be a positive integer");
     endif
-    [retval, shape, scale] = common_size (shape, scale, zeros (r, c));
-    if (retval > 0)
-      error ("weibull_rnd: alpha and sigma must be scalar or of size %d by %d",
-	     r, c);
+    sz = [r, c];
+
+    if (any (size (scale) != 1) && 
+	((length (size (scale)) != length (sz)) ||
+	 any (size (scale) != sz)))
+      error ("weilbull_rnd: shape and scale must be scalar or of size [r, c]");
+    endif
+  elseif (nargin == 3)
+    if (isscalar (r) && (r > 0))
+      sz = [r, r];
+    elseif (isvector(r) && all (r > 0))
+      sz = r(:)';
+    else
+      error ("weibull_rnd: r must be a postive integer or vector");
+    endif
+
+    if (any (size (scale) != 1) && 
+	((length (size (scale)) != length (sz)) ||
+	 any (size (scale) != sz)))
+      error ("weibull_rnd: shape and scale must be scalar or of size sz");
     endif
   elseif (nargin == 2)
-    [retval, shape, scale] = common_size (shape, scale);
-    if (retval > 0)
-      error ("weibull_rnd: alpha and sigma must be of common size or scalar");
-    endif
+    sz = size(shape);
   else
     usage ("weibull_rnd (alpha, sigma, r, c)");
   endif
 
-  [r, c] = size (shape);
-  s = r * c;
-  shape = reshape (shape, 1, s);
-  scale = reshape (scale, 1, s);
-
-  rnd = NaN * ones (1, s);
-  k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf));
-  if (any (k))
-    rnd(k) = (scale(k)
-              .* (- log (1 - rand (1, length (k)))) .^ (1 ./ shape(k)));
+  if (isscalar (shape) && isscalar (scale))
+    if ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf))
+      rnd = (scale * (- log (1 - rand (sz))) .^ (1 / shape));
+    else
+      rnd = NaN * ones (sz);
+    endif
+  else
+    rnd = NaN * ones (sz);
+    k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf));
+    if (any (k))
+      rnd(k) = (scale(k)
+		.* (- log (1 - rand (size (k)))) .^ (1 ./ shape(k)));
+    endif
   endif
 
-  rnd = reshape (rnd, r, c);
-
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/distributions/wiener_rnd.m b/scripts/statistics/distributions/wiener_rnd.m
--- a/scripts/statistics/distributions/wiener_rnd.m
+++ b/scripts/statistics/distributions/wiener_rnd.m
@@ -40,7 +40,11 @@
   elseif (nargin == 2)
     n = 1000;
   elseif (nargin > 3)
-    usage ("wiener_rnd (t, d,n)");
+    usage ("wiener_rnd (t, d, n)");
+  endif
+
+  if (!isscalar (t) || !isscalar (d) || !isscalar (n))
+    error ("wiener_rnd: t, d and n must all be positive integers");
   endif
 
   retval = randn (n * t, d);