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')