[Octave-cvsupdate] octave-forge/main/statistics nanmedian.m,1.2,1.3 nanstd.m,1.4,1.5

[David B. <ad...@us...>]

Update of /cvsroot/octave/octave-forge/main/statistics
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv8591

Modified Files:
	nanmedian.m nanstd.m 
Log Message:
NDArray update for nanstd.m and nanmedian.m

Index: nanstd.m
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/statistics/nanstd.m,v
retrieving revision 1.4
retrieving revision 1.5
diff -u -d -r1.4 -r1.5
--- nanstd.m	2 Apr 2004 13:21:19 -0000	1.4
+++ nanstd.m	6 Apr 2004 09:49:43 -0000	1.5
@@ -25,29 +25,31 @@
 ##		this provides the square root of the second moment around the mean
 ##
 ## See also: nanmin, nanmax, nansum, nanmedian, nanmean
-function v = nanstd (X, opt, dim)
+function v = nanstd (X, opt, varargin)
   if nargin < 1
     usage ("v = nanstd(X [, opt [, dim]])");
   else
     if nargin < 3
       dim = min(find(size(X)>1));
       if isempty(dim), dim=1; endif;
+    else
+      dim = varargin{1};
     endif
     if ((nargin < 2) || isempty(opt))
       opt = 0;
     endif
 
     ## determine the number of non-missing points in each data set
-    n = sum (!isnan(X), dim);
+    n = sum (!isnan(X), varargin{:});
     
     ## replace missing data with zero and compute the mean
     X(isnan(X)) = 0;
-    meanX = sum (X, dim) ./ n;
+    meanX = sum (X, varargin{:}) ./ n;
     
     ## subtract the mean from the data and compute the sum squared
     sz = ones(1,length(size(X)));
     sz(dim) = size(X,dim);
-    v = sumsq (X - repmat(meanX,sz), dim);
+    v = sumsq (X - repmat(meanX,sz), varargin{:});
     
     ## because the missing data was set to zero each missing data
     ## point will contribute (-meanX)^2 to sumsq, so remove these

Index: nanmedian.m
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/statistics/nanmedian.m,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -d -r1.2 -r1.3
--- nanmedian.m	12 Sep 2003 14:22:42 -0000	1.2
+++ nanmedian.m	6 Apr 2004 09:49:43 -0000	1.3
@@ -20,18 +20,19 @@
 ## [Is this behaviour compatible?]
 ##
 ## See also: nanmin, nanmax, nansum, nanmean
-function v = nanmedian (X, dim)
+function v = nanmedian (X, varargin)
   if nargin < 1 || nargin > 2
-    usage ("v = nanmean(X [, dim])");
+    usage ("v = nanmedian(X [, dim])");
+  endif
+  if nargin < 2
+    dim = min(find(size(X)>1));
+    if isempty(dim), dim=1; endif;
   else
-    if nargin == 1
-      if size(X,1) == 1
-	dim = 2; 
-      else
-        dim = 1;
-      endif
-    endif
-    if (dim == 2) X = X.'; endif
+    dim = varargin{:};
+  endif
+
+  sz = size (X);
+  if (prod (sz) > 1)
     try dfi = do_fortran_indexing;
     catch dfi = 0;
     end
@@ -43,32 +44,35 @@
       warn_fortran_indexing = 0;
       ## Find lengths of datasets after excluding NaNs; valid datasets
       ## are those that are not empty after you remove all the NaNs
-      n = size(X,1) - sum (isnan(X));
-      valid = find(n!=0);
+      n = sz(dim) - sum (isnan(X),varargin{:});
 
-      ## Extract all non-empty datasets and sort, replacing NaN with Inf
-      ## so that the invalid elements go toward the ends of the columns
-      X (isnan(X)) = Inf;
-      X = sort ( X (:, valid) );
+      ## When n is equal to zero, force it to one, so that median
+      ## picks up a NaN value below
+      n (n==0) = 1;
 
-      ## Determine the offset for each remaining column in single index mode
-      colidx = (0:size(X,2)-1)*size(X,1);
+      ## Sort the datasets, with the NaN going to the end of the data
+      X = sort (X, varargin{:});
 
-      ## Assume the median for all datasets will be NaNs
-      v = NaN*ones(size(n));
+      ## Determine the offset for each column in single index mode
+      colidx = reshape((0:(prod(sz) / sz(dim) - 1)), size(n)); 
+      colidx = floor(colidx / prod(sz(1:dim-1))) * prod(sz(1:dim)) + ...
+	  mod(colidx,prod(sz(1:dim-1)));
+      stride = prod(sz(1:dim-1));
 
       ## Average the two central values of the sorted list to compute
       ## the median, but only do so for valid rows.  If the dataset
       ## is odd length, the single central value will be used twice.
       ## E.g., 
-      ##   for n==5, ceil(2.5+0.4) is 3 and floor(2.5+0.6) is also 3
-      ##   for n==6, ceil(3.0+0.4) is 4 and floor(3.0+0.6) is 3
-      v(valid) = ( X (colidx + floor(n(valid)./2+0.6)) ...
-		 + X (colidx + ceil(n(valid)./2+0.4)) ) ./ 2;
+      ##   for n==5, ceil(2.5+0.5) is 3 and floor(2.5+0.5) is also 3
+      ##   for n==6, ceil(3.0+0.5) is 4 and floor(3.0+0.5) is 3
+      ## correction made for stride of data "stride*ceil(2.5-0.5)+1"
+      v = (X(colidx + stride*ceil(n./2-0.5) + 1)  + ...
+	   X(colidx + stride*floor(n./2-0.5) + 1)) ./ 2;
     unwind_protect_cleanup
       do_fortran_indexing = dfi;
       warn_fortran_indexing = wfi;
     end_unwind_protect
-    if (dim == 2) v = v.'; endif
+  else
+    error ("nanmedian: invalid matrix argument");
   endif
 endfunction