close all; clear; clc
Ns = 10; % number of scans (fMRI) or subjects (VBM)
Nv = 10*10*10; % number of voxels (e.g. for 10*10*10 VOI)
randn('state', 0);
tol = 1e-12; % max absolute error tolerance

f = figure;
set(f,'Position', [60, 60, 800, 600]);
p = 1; % (subplot index)
for mn = [0 4 8] % different grand means over voxels and scans/subjects
    for vr = [5 3 1] % different variabilities of voxel data
        % Simulate random scans/subjects
        y = mn + randn(Ns, 1);
        % Simulate noisy data for VOI
        y = repmat(y, 1, Nv) + vr*randn(Ns, Nv);

        % Mean of VOI (mean over voxels, so one value per scan/subject)
        m = mean(y, 2);

        % Scaled largest eigenvariate (voxels treated as repeats like
        % subject covariance matrix, so one value per scan/subject)
        [V E] = eig(y*y');
        [mx ix] = max(diag(E));
        Vev = V(:, ix);
        Vev = Vev * sign(sum(y'*Vev)+eps); % (ensure in same dir as mean)
        % Scaling is by largest eigenvalue of y*y' divided by Nv
        Vsc = Vev * sqrt(mx / Nv);
        
        % Compare with version modified from spm_regions
        [u s u] = svd(spm_atranspa(y'));
        s       = diag(s);
        u       = u(:, 1);
        u       = u * sign(sum(y'*u));
        Y       = u * sqrt(s(1)/Nv);
        if max(abs(Y-Vsc)) > tol; keyboard; end
        
        subplot(3, 3, p)
        plot(1:Ns,m,'rs', 1:Ns,Y,'bx')
        title(sprintf('mn %g, vr %2.2g', mn, vr))
        p = p + 1;
    end
end
legend('mean', 'ev', 'Location','SouthEast')