pr_sum_func

 creates summary stats for region data
 FORMAT [sumY, varY] = pr_sum_func(y, sumfunc, wt)
 
 $Id$

0001 function [sumY, varY] = pr_sum_func(y, sumfunc, wt)
0002 % creates summary stats for region data
0003 % FORMAT [sumY, varY] = pr_sum_func(y, sumfunc, wt)
0004 %
0005 % $Id$
0006   
0007 if nargin < 2
0008   error('Need data matrix and summary function');
0009 end
0010 [m n]   = size(y);
0011 if any([m n] == 0)
0012   error('Data vector is empty');
0013 end
0014 if nargin < 3
0015   wt = [];
0016 end
0017 if isempty(wt)
0018   wt = ones(n,1);
0019 end
0020 if size(wt,1)==1
0021   wt = wt';
0022 end
0023 if n == 1 % only 1 sample in ROI
0024   sumY = y; varY = y*Inf;
0025   return
0026 end
0027 
0028 varY = ones(m,1) * Inf;
0029 switch sumfunc
0030  case 'mean'
0031   sumY = mean(y, 2);
0032   ssy  = sum(y.^2, 2);
0033   varY = (ssy - n*sumY.^2)/(n-1);
0034  case 'wtmean'
0035   % Formulae from GNU scientific library
0036   % $\hat\mu = (\sum w_i x_i) / (\sum w_i)$
0037   % \hat\sigma^2 = ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2)))
0038   %                \sum w_i (x_i - \hat\mu)^2$
0039   swt = sum(wt);
0040   sumY = y*wt/swt;
0041   varY = (swt/(swt.^2 - wt'*wt)) * ((y - (sumY * ones(1,n))).^2 * wt);
0042 
0043   % original formula
0044   %nwt = sum(wt ~= 0);
0045   %varY = (y - (sumY * ones(1,n))).^2 * wt / ((nwt-1) * swt / nwt);
0046  case 'median'
0047   sumY = median(y, 2);
0048  case 'eigen1'    
0049   % code taken from spm_regions.m l 230-247, with thanks
0050   % @(#)spm_regions.m    2.7 Karl Friston 00/10/04
0051   
0052 % compute regional response in terms of first eigenvariate
0053 %-----------------------------------------------------------------------
0054 if m > n
0055     [v s v] = svd(spm_atranspa(y));
0056     s       = diag(s);
0057     v       = v(:,1);
0058     u       = y*v/sqrt(s(1));
0059 else
0060     [u s u] = svd(spm_atranspa(y'));
0061     s       = diag(s);
0062     u       = u(:,1);
0063     v       = y'*u/sqrt(s(1));
0064 end
0065 d       = sign(sum(v));
0066 u       = u*d;
0067 v       = v*d;
0068 sumY       = u*sqrt(s(1)/n);
0069 
0070 % end of paste from spm_regions
0071       
0072  otherwise
0073   error(['Do not recognize summary function ' sumfunc]);
0074 end