function [sumY, varY] = mars_sum_func(y, sumfunc, wt)
% mars_sum_func - creates summary stats for region data
% FORMAT [sumY, varY] = mars_sum_func(y, sumfunc, wt)
  
if nargin < 2
  error('Need data matrix and summary function');
end
if nargin < 3
  wt = ones(1, size(y,2));
end

if any(size(y) == 1)
  sumY = y; varY = y*Inf;
  return
end

varY = ones(size(y, 1),1) * Inf;
switch sumfunc
 case 'mean'
  sumY = mean(y, 2);
  varY = var(y')';
 case 'wtmean'
  wmean = sum(y .* (ones(ny, 1) * wt), 2) ./sum(wt.^2);
  sumY = wmean;
  % Note that weighted variance is not adjusted for N-1
  varY = sum((y - (wmean * wt)).^2,2)/nwt;
 case 'median'
  sumY = median(y, 2);
 case 'eigen1'    
  % code taken from spm_regions.m l 230-247, with thanks
  % @(#)spm_regions.m	2.7 Karl Friston 00/10/04
  
% compute regional response in terms of first eigenvariate
%-----------------------------------------------------------------------
[m n]   = size(y);
if m > n
	[v s v] = svd(spm_atranspa(y));
	s       = diag(s);
	v       = v(:,1);
	u       = y*v/sqrt(s(1));
else
	[u s u] = svd(spm_atranspa(y'));
	s       = diag(s);
	u       = u(:,1);
	v       = y'*u/sqrt(s(1));
end
d       = sign(sum(v));
u       = u*d;
v       = v*d;
sumY       = u*sqrt(s(1)/n);

% end of paste from spm_regions
      
 otherwise
  error(['Do not recognize summary function ' sumfunc]);
end
