Fix #64

* fail on weird values in VBA_orth
* simplifies nan suite

core/diagnostics/VBA_getDiagnostics.m

@@ -166,7 +166,7 @@
     dy(iS).dy = res(:);
     dy(iS).R = VBA_spm_autocorr(res);
     dy(iS).m = VBA_nanmean(dy(iS).dy);
-    dy(iS).v = VBA_nanvar(dy(iS).dy);
+    dy(iS).v = var(dy(iS).dy,'omitnan');
     [dy(iS).ny,dy(iS).nx] = hist(dy(iS).dy,10);
     dy(iS).ny = dy(iS).ny./sum(dy(iS).ny);
     d = diff(dy(iS).nx);

demos/3_behavioural/demo_ToMgames.m

@@ -115,7 +115,7 @@
 if noisy
     na = 2;
     mP2 = VBA_nanmean(Perf2,3);
-    sP2 = VBA_nanstd(Perf2,[],3);
+    sP2 = std(Perf2,0,3,'omitnan');
 else
     na = 1;
 end

thrid-party/spm/VBA_spm_autocorr.m

@@ -15,7 +15,7 @@
 R = zeros(n,t*2);
 % standardize y
 my = VBA_nanmean(y,2);
-sy = VBA_nanstd(y,[],2);
+sy = std(y,0,2,'omitnan');
 sy(sy==0) = 1; % correct for identically constant time series
 y = y - repmat(my,1,t);
 y = diag(1./sy)*y;

utils/VBA_nanmean.m

@@ -1,50 +1,11 @@
-function y = VBA_nanmean(x,dim)
-% FORMAT: Y = VBA_NANMEAN(X,DIM)
-% 
-%    Average or mean value ignoring NaNs
+function y = VBA_nanmean(varargin)
+% // VBA toolbox //////////////////////////////////////////////////////////
 %
-%    This function enhances the functionality of NANMEAN as distributed in
-%    the MATLAB Statistics Toolbox and is meant as a replacement (hence the
-%    identical name).  
+% y = VBA_nanmean(varargin)
+% mean ignoring NaNs
 %
-%    NANMEAN(X,DIM) calculates the mean along any dimension of the N-D
-%    array X ignoring NaNs.  If DIM is omitted NANMEAN averages along the
-%    first non-singleton dimension of X.
+% see mean() for the arguments
 %
-%    Similar replacements exist for NANSTD, NANMEDIAN, NANMIN, NANMAX, and
-%    NANSUM which are all part of the NaN-suite.
-%
-%    See also MEAN
-
-% -------------------------------------------------------------------------
-%    author:      Jan Gl?scher
-%    affiliation: Neuroimage Nord, University of Hamburg, Germany
-%    email:       [email protected]
-%    
-%    $Revision: 1.1 $ $Date: 2004/07/15 22:42:13 $
-
-if isempty(x)
-	y = NaN;
-	return
-end
-
-if nargin < 2
-	dim = min(find(size(x)~=1));
-	if isempty(dim)
-		dim = 1;
-	end
-end
-
-% Replace NaNs with zeros.
-nans = isnan(x);
-x(isnan(x)) = 0; 
-
-% denominator
-count = size(x,dim) - sum(nans,dim);
-
-% Protect against a  all NaNs in one dimension
-i = find(count==0);
-count(i) = ones(size(i));
+% /////////////////////////////////////////////////////////////////////////
 
-y = sum(x,dim)./count;
-y(i) = i + NaN;
+y = mean(varargin{:},'omitnan');

utils/VBA_orth.m

@@ -1,4 +1,4 @@
-function A = VBA_orth(A,norm)
+function A = VBA_orth (A, norm)
 % serial orthogonalisation of matrix A
 % function A = VBA_orth(A,norm)
 % This function considers each column serially, projecting it in the null
@@ -8,17 +8,26 @@
 %   - norm: a flag for normalization of A's columns
 % OUT:
 %   - A: the serially orthogonalized matrix A
-try,norm;catch,norm=0;end
+
+if nargin < 2
+    norm = false;
+end
+
+if VBA_isWeird(A) 
+    error('***VBA_orth: your matrix should not include any weird values (Nans, Inf, etc).');
+end
+
 [n,p] = size(A);
 I = eye(n);
-if norm && std(A(:,1))~=0
-    A(:,1) = (A(:,1)-mean(A(:,1)))./std(A(:,1));
+
+if norm
+    A(:,1) = VBA_zscore (A(:, 1));
 end
 for i=2:p
-    X = A(:,1:i-1);
-    nX = I - X*pinv(X'*X)*X';
-    A(:,i) = nX*A(:,i);
-    if norm && std(A(:,i))~=0
-        A(:,i) = (A(:,i)-mean(A(:,i)))./std(A(:,i));
+    X = A(:, 1 : i - 1);
+    nX = I - X * pinv(X' * X) * X';
+    A(:,i) = nX * A(:,i);
+    if norm
+        A(:,i) = VBA_zscore (A(:,i));
     end
 end

utils/VBA_zscore.m

@@ -0,0 +1,31 @@
+function [z, mu, sigma] = VBA_zscore (x, w, dim, nanflag)
+% // VBA toolbox //////////////////////////////////////////////////////////
+%
+% [z, mu, sigma] = VBA_zscore (x, flag, dim, nanflag)
+% standard score, aka z-score
+%
+% Compute the stardardized zscore of x, ie z = (x - mean(x)) / std(x)
+%
+% /////////////////////////////////////////////////////////////////////////
+
+if nargin < 2 || isempty (w)
+	w = 0;
+end
+
+if nargin < 3 || isempty (dim)
+	dim = find (size(x) > 1, 1);
+	if isempty(dim)
+		dim = 1;
+	end
+end
+
+if nargin < 4
+    nanflag = 'includenan';
+end
+
+mu = mean (x, dim, nanflag);
+sigma = std (x, w, dim, nanflag);
+sigma(sigma == 0) = 1;
+
+z = bsxfun(@minus, x, mu);
+z = bsxfun(@rdivide, z, sigma);