--- a/extra/NaN/inst/zscore.m	Fri Sep 12 17:46:13 2014 +0000
+++ b/extra/NaN/inst/zscore.m	Fri Sep 12 17:51:52 2014 +0000
@@ -1,15 +1,15 @@
-function [i,v,m] = zscore(i,DIM)
+function [i,m,s] = zscore(i,OPT, DIM, W)
 % ZSCORE removes the mean and normalizes data 
 % to a variance of 1. Can be used for pre-whitening of data, too. 
 %
-% [z,r,m] = zscore(x,DIM)
+% [z,mu, sigma] = zscore(x [,OPT [, DIM])
 %   z   z-score of x along dimension DIM
-%   r   is the inverse of the standard deviation
-%   m   is the mean of x
+%   sigma is the inverse of the standard deviation
+%   mu   is the mean of x
 %
 % The data x can be reconstucted with 
-%     x = z*diag(1./r) + repmat(m,size(z)./size(m))  
-%     z = x*diag(r) - repmat(m.*v,size(z)./size(m))  
+%     x = z*diag(sigma) + repmat(m, size(z)./size(m))  
+%     z = x*diag(1./sigma) - repmat(m.*v, size(z)./size(m))  
 %
 % DIM	dimension
 %	1: STATS of columns
@@ -44,9 +44,12 @@
 if any(size(i)==0); return; end;
 
 if nargin<2
+        OPT=[]; 
+end
+if nargin<3
         DIM=[]; 
 end
-if nargin<3
+if nargin<4
         W = []; 
 end
 if isempty(DIM), 
@@ -56,9 +59,11 @@
 
 
 % pre-whitening
-m = mean(i,DIM,W);
-i = i-repmat(m,size(i)./size(m));  % remove mean
-v = 1./sqrt(mean(i.^2,DIM,W));
-i = i.*repmat(v,size(i)./size(v)); % scale to var=1
+[S,N,SSQ] = sumskipnan(i, DIM, W);
+m = S./N; 
+i = i-repmat(m, size(i)./size(m));  % remove mean
+s = std (i, OPT, DIM, W);
+s(s==0)=1;
+i = i ./ repmat(s,size(i)./size(s)); % scale to var=1