--- a/main/statistics/INDEX	Mon Dec 30 15:20:20 2013 +0000
+++ b/main/statistics/INDEX	Mon Dec 30 18:37:39 2013 +0000
@@ -13,6 +13,7 @@
  geostat
  gevcdf gevfit gevfit_lmom gevinv gevlike gevpdf gevrnd gevstat 
  hygestat
+ iwishrnd
  jsucdf
  jsupdf
  lognstat
@@ -30,6 +31,7 @@
  unifstat
  vmpdf vmrnd
  wblstat
+ wishrnd
 Descriptive statistics
  nansum
  nanmax

--- a/main/statistics/NEWS	Mon Dec 30 15:20:20 2013 +0000
+++ b/main/statistics/NEWS	Mon Dec 30 18:37:39 2013 +0000
@@ -6,6 +6,10 @@
  ** Fixed second output of nanmax and nanmin.
 
  ** Corrected handle for outliers in boxplot.
+ 
+ ** The following functions are new:
+ 
+    wishrnd iwishrnd 
 
 Summary of important user-visible changes for statistics 1.2.2:
 -------------------------------------------------------------------
@@ -19,8 +23,7 @@
 
       pcares  pcacov  runstest  stepwisefit hist3
 
- ** dendogram now returns the leaf node numbers and order that the nodes were
-    displayed in.
+ ** dendogram now returns the leaf node numbers and order that the nodes were displayed in.
 
  ** New faster implementation of princomp.
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/statistics/inst/iwishrnd.m	Mon Dec 30 18:37:39 2013 +0000
@@ -0,0 +1,69 @@
+## Copyright (C) 2013 Nir Krakauer <nkrakauer@ccny.cuny.edu>
+##
+## This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License along with Octave; see the file COPYING.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {} [@var{W}[, @var{DI}]] = iwishrnd (@var{Psi}, @var{df}[, @var{DI}][, @var{n}=1])
+## Return a random matrix sampled from the inverse Wishart distribution with given parameters
+##
+## Inputs: the @var{p} x @var{p} positive definite matrix @var{Tau} and scalar degrees of freedom parameter @var{df} (and optionally the transposed Cholesky factor @var{DI} of @var{Sigma} = @code{inv(Tau)}).
+## @var{df} can be non-integer as long as @var{df} > @var{d}
+##
+## Output: a random @var{p} x @var{p}  matrix @var{W} from the inverse Wishart(@var{Tau}, @var{df}) distribution. (@code{inv(W)} is from the Wishart(@code{inv(Tau)}, @var{df}) distribution.) If @var{n} > 1, then @var{W} is @var{p} x @var{p} x @var{n} and holds @var{n} such random matrices. (Optionally, the transposed Cholesky factor @var{DI} of @var{Sigma} is also returned.)
+##
+## Averaged across many samples, the mean of @var{W} should approach @var{Tau} / (@var{df} - @var{p} - 1).
+##
+## Reference: Yu-Cheng Ku and Peter Bloomfield (2010), Generating Random Wishart Matrices with Fractional Degrees of Freedom in OX, http://www.gwu.edu/~forcpgm/YuChengKu-030510final-WishartYu-ChengKu.pdf
+## 
+## @seealso{wishrnd, iwishpdf}
+## @end deftypefn
+
+## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
+## Description: Random matrices from the inverse Wishart distribution
+
+function [W, DI] = iwishrnd(Tau, df, DI, n = 1)
+
+if (nargin < 2)
+  print_usage ();
+endif
+
+if nargin < 3 || isempty(DI)
+  try
+    D = chol(inv(Tau));
+  catch
+    error('Cholesky decomposition failed; Tau probably not positive definite')
+  end_try_catch
+  DI = D';
+else  
+  D = DI';  
+endif
+
+w = wishrnd([], df, D, n);
+
+if n > 1
+  p = size(D, 1);
+  W = nan(p, p, n);
+endif
+
+for i = 1:n
+  W(:, :, i) = inv(w(:, :, i));
+endfor
+
+endfunction
+
+
+
+%!assert(size (iwishrnd (1,2,1)), [1, 1]);
+%!assert(size (iwishrnd ([],2,1)), [1, 1]);
+%!assert(size (iwishrnd ([3 1; 1 3], 2.00001, [], 1)), [2, 2]);
+%!assert(size (iwishrnd (eye(2), 2, [], 3)), [2, 2, 3]);
+
+%% Test input validation
+%!error iwishrnd ()
+%!error iwishrnd (1)
+%!error iwishrnd ([-3 1; 1 3],1)
+%!error iwishrnd ([1; 1],1)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/statistics/inst/wishrnd.m	Mon Dec 30 18:37:39 2013 +0000
@@ -0,0 +1,88 @@
+## Copyright (C) 2013 Nir Krakauer <nkrakauer@ccny.cuny.edu>
+##
+## This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License along with Octave; see the file COPYING.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {} [@var{W}[, @var{D}]] = wishrnd (@var{Sigma}, @var{df}[, @var{D}][, @var{n}=1])
+## Return a random matrix sampled from the Wishart distribution with given parameters
+##
+## Inputs: the @var{p} x @var{p} positive definite matrix @var{Sigma} and scalar degrees of freedom parameter @var{df} (and optionally the Cholesky factor @var{D} of @var{Sigma}).
+## @var{df} can be non-integer as long as @var{df} > @var{p}
+##
+## Output: a random @var{p} x @var{p}  matrix @var{W} from the Wishart(@var{Sigma}, @var{df}) distribution. If @var{n} > 1, then @var{W} is @var{p} x @var{p} x @var{n} and holds @var{n} such random matrices. (Optionally, the Cholesky factor @var{D} of @var{Sigma} is also returned.)
+##
+## Averaged across many samples, the mean of @var{W} should approach @var{df}*@var{Sigma}, and the variance of each element @var{W}_{ij} should approach @var{df}*(@var{Sigma}_{ij}^2 + @var{Sigma}_{ii}*@var{Sigma}_{jj})
+##
+## Reference: Yu-Cheng Ku and Peter Bloomfield (2010), Generating Random Wishart Matrices with Fractional Degrees of Freedom in OX, http://www.gwu.edu/~forcpgm/YuChengKu-030510final-WishartYu-ChengKu.pdf
+## 
+## @seealso{iwishrnd, wishpdf}
+## @end deftypefn
+
+## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
+## Description: Compute the probability density function of the Wishart distribution
+
+function [W, D] = wishrnd(Sigma, df, D, n=1)
+
+if (nargin < 3)
+  print_usage ();
+endif
+
+if nargin < 3 || isempty(D)
+  try
+    D = chol(Sigma);
+  catch
+    error('Cholesky decomposition failed; Sigma probably not positive definite')
+  end_try_catch
+endif
+
+p = size(D, 1);
+
+if df < p
+  df = floor(df); #distribution not defined for small noninteger df
+  df_isint = 1;
+else 
+#check for integer degrees of freedom
+ df_isint = (df == floor(df));
+endif
+
+if ~df_isint
+  [ii, jj] = ind2sub([p, p], 1:(p*p));
+endif
+
+if n > 1
+  W = nan(p, p, n);
+endif
+
+for i = 1:n
+  if df_isint
+    Z = randn(df, p) * D;
+    W(:, :, i) = Z'*Z;
+  else
+    Z = diag(sqrt(chi2rnd(df - (0:(p-1))))); #fill diagonal
+    #note: chi2rnd(x) is equivalent to 2*randg(x/2), but the latter seems to offer no performance advantage
+    Z(ii > jj) = randn(p*(p-1)/2, 1); #fill lower triangle with normally distributed variates
+    Z = D * Z;
+    W(:, :, i) = Z*Z';
+  endif
+
+  
+endfor
+
+endfunction
+
+
+
+%!assert(size (wishrnd (1,2,1)), [1, 1]);
+%!assert(size (wishrnd ([],2,1)), [1, 1]);
+%!assert(size (wishrnd ([3 1; 1 3], 2.00001, [], 1)), [2, 2]);
+%!assert(size (wishrnd (eye(2), 2, [], 3)), [2, 2, 3]);
+
+%% Test input validation
+%!error wishrnd ()
+%!error wishrnd (1)
+%!error wishrnd ([-3 1; 1 3],1)
+%!error wishrnd ([1; 1],1)