example_synth.m

Optimal Spatial Regelarisation of Autocorrelation Estimates in fMRI Analysis,very useful,if you like

% Example analysis on the synthetical data sets described in Gautama and
% Van Hulle (2004).

clear all

ar_order = 1;		% order of AR-model used for correction
alpha = .05;		% Statistical threshold (uncorrected level)
batch_mode = 2;

switch (batch_mode)
  case 1
	gamma_act = 1;
	n_act = 500;
	add_ac = 0;
	ar_order = 1;
	sx = 40;
	sy = 50;

  case 2
	gamma_act = 1;
	n_act = 500;
	add_ac = .3;
	ar_order = 1;
	sx = 40;
	sy = 50;

  case 3
	gamma_act = 0;
	n_act = 500;
	add_ac = .3;
	ar_order = 1;
	sx = 40;
	sy = 50;

end


		%%%%%%%%%%%%%%%%%
		% Generate Data %
		%%%%%%%%%%%%%%%%%
[XX C INDact X] = generate_synth (gamma_act, n_act, add_ac, ar_order, sx, sy);
R = [1 0 0; 0 1 0];
[st sx sy] = size(XX);
df1 = size(R,1);
df2 = st-size(R,2);
f_thres = finv(1-alpha,df1, df2);	% Threshold for uncorrected F-testing


		%%%%%%%%%%
		% F-test %
		%%%%%%%%%%
[Funbr,Fols,Cunbr,Cunb] = analyse_f (XX, X, R, ar_order);


		%%%%%%%%%%%%%%%%%%
		% The Art Corner %
		%%%%%%%%%%%%%%%%%%
m1 = max(max(max(cat(3,C(:,:,2),Cunb(:,:,2),Cunbr(:,:,2)))));
subplot (2,3,1)
pcolor (C(:,:,2))
shading flat
axis image
caxis ([0 m1]);
title ('Desired AC(1)');

subplot (2,3,2)
pcolor (Cunb(:,:,2))
shading flat
axis image
caxis ([0 m1]);
title ('Unb AC(1)');

subplot (2,3,3)
pcolor (Cunbr(:,:,2))
shading flat
axis image
caxis ([0 m1]);
title ('UnbR AC(1)');

subplot (2,3,4)
pcolor (Fols)
shading flat
axis image
title ('F OLS');
caxis ([0 f_thres]);
colorbar

subplot (2,3,6)
pcolor (Funbr)
shading flat
axis image
title ('F UnbR');
caxis ([0 f_thres]);
colorbar


		%%%%%%%%%%%%%%
		% Evaluation %
		%%%%%%%%%%%%%%
INDnact = (1:(sx*sy))';
INDnact(INDact) = [];

% MSE-AC
mse_ac_unbr = mean(mean((C(:,:,2)-Cunbr(:,:,2)).^2));
mse_ac_unb = mean(mean((C(:,:,2)-Cunb(:,:,2)).^2));

% MSE-F
Fdes = finv((1:length(INDnact))'./(length(INDnact)+1), df1, df2);
mse_f_unbr = mean((Fdes(:)-sort(Funbr(INDnact))).^2);
mse_f_ols = mean((Fdes(:)-sort(Fols(INDnact))).^2);

% False positives
fp_unbr = mean(Funbr(INDnact)>f_thres);
fp_ols = mean(Fols(INDnact)>f_thres);

% True positives
if (~isempty(INDact))
	tp_unbr = mean(Funbr(INDact)>f_thres);
	tp_ols = mean(Fols(INDact)>f_thres);
else
	tp_unbr = NaN;
	tp_ols = NaN;
end

% Evaluation measures, see paper
fprintf ('\n*** Evaluation ***\n\n');
fprintf ('MSE-AC\t%.4f(Unb)\t%.4f(UnbR)\n', mse_ac_unb, mse_ac_unbr);
fprintf ('MSE-F\t%.4f(OLS)\t%.4f(UnbR)\n', mse_f_ols, mse_f_unbr);
fprintf ('FP\t%.4f(OLS)\t%.4f(UnbR)\n', fp_ols, fp_unbr);
fprintf ('TP\t%.4f(OLS)\t%.4f(UnbR)\n', tp_ols, tp_unbr);

% Unb	unbiased AC estimate
% UnbR	unbiased, regularised AC estimate
% OLS	ordinary least-squares regression: uncorrected