function [P,p,Em,En,EN] = spm_P_RF(c,k,Z,df,STAT,R,n)
% Returns the [un]corrected P value using unifed EC theory
% FORMAT [P p Em En EN] = spm_P_RF(c,k,Z,df,STAT,R,n)
%
% c - cluster number
% k - extent {RESELS}
% Z - height {minimum over n values}
% df - [df{interest} df{error}]
% STAT - Statisical feild
% 'Z' - Gaussian feild
% 'T' - T - feild
% 'X' - Chi squared feild
% 'F' - F - feild
% R - RESEL Count {defining search volume}
% n - number of component SPMs in conjunction
%
% P - corrected P value - P(n > kmax}
% p - uncorrected P value - P(n > k}
% Em - expected total number of maxima {m}
% En - expected total number of resels per cluster {n}
% EN - expected total number of voxels {N}
%
%___________________________________________________________________________
%
% spm_P_RF returns the probability of c or more clusters with more than
% k voxels in volume process of R RESELS thresholded at u. All p values
% can be considered special cases:
%
% spm_P_RF(1,0,Z,df,STAT,1,n) = uncorrected p value
% spm_P_RF(1,0,Z,df,STAT,R,n) = corrected p value {based on height Z)
% spm_P_RF(1,k,u,df,STAT,R,n) = corrected p value {based on extent k at u)
% spm_P_RF(c,k,u,df,STAT,R,n) = corrected p value {based on number c at k and u)
% spm_P_RF(c,0,u,df,STAT,R,n) = omnibus p value {based on number c at u)
%
% If n > 1 a conjunction probility over the n values of the statistic
% is returned
%
% Ref: Hasofer AM (1978) Upcrossings of random fields
% Suppl Adv Appl Prob 10:14-21
% Ref: Friston et al (1993) Comparing functional images: Assessing
% the spatial extent of activation foci
% Ref: Worsley KJ et al 1996, Hum Brain Mapp. 4:58-73
%___________________________________________________________________________
% @(#)spm_P_RF.m 1.4 Karl Friston 01/07/29
% get expectations
%===========================================================================
% get EC densities
%---------------------------------------------------------------------------
D = max(find(R));
R = R(1:D);
G = sqrt(pi)./gamma(([1:D])/2);
EC = spm_ECdensity(STAT,Z,df);
EC = EC([1:D]);
% corrected p value
%---------------------------------------------------------------------------
P = triu(toeplitz(EC'.*G))^n;
P = P(1,:)';
Em = (R./G)*P;
EN = P(1)*R(D);
En = G(D)*P(1)/(eps+P(D)); % i.e. En = EN/Em;
% get P{n > k}
%===========================================================================
% assume a Gaussian form for P{n > k} ~ exp(-beta*k^(2/D))
% Appropriate for SPM{Z} and high d.f. SPM{T}
%---------------------------------------------------------------------------
D = D - 1;
if ~k | ~D
p = 1;
elseif STAT == 'Z'
beta = (gamma(D/2 + 1)/En)^(2/D);
p = exp(-beta*(k^(2/D)));
elseif STAT == 'T'
beta = (gamma(D/2 + 1)/En)^(2/D);
p = exp(-beta*(k^(2/D)));
elseif STAT == 'X'
beta = (gamma(D/2 + 1)/En)^(2/D);
p = exp(-beta*(k^(2/D)));
elseif STAT == 'F'
beta = (gamma(D/2 + 1)/En)^(2/D);
p = exp(-beta*(k^(2/D)));
end
% Poisson clumping heuristic {for multiple clusters}
%===========================================================================
P = 1 - spm_Pcdf(c - 1,(Em + eps)*p);
% set P and p = [] for non-implemented cases
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if k > 0 & n > 1
P = []; p = [];
end
if k > 0 & (STAT == 'X' | STAT == 'F')
P = []; p = [];
end
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% spm_ECdensity
%===========================================================================
function [EC] = spm_ECdensity(STAT,t,df)
% Returns the EC density
%___________________________________________________________________________
%
% Reference : Worsley KJ et al 1996, Hum Brain Mapp. 4:58-73
%
%---------------------------------------------------------------------------
% EC densities (EC}
%---------------------------------------------------------------------------
t = t(:)';
if STAT == 'Z'
% Gaussian Field
%-------------------------------------------------------------------
a = 4*log(2);
b = exp(-t.^2/2);
EC(1,:) = 1 - spm_Ncdf(t);
EC(2,:) = a^(1/2)/(2*pi)*b;
EC(3,:) = a/((2*pi)^(3/2))*b.*t;
EC(4,:) = a^(3/2)/((2*pi)^2)*b.*(t.^2 - 1);
elseif STAT == 'T'
% T - Field
%-------------------------------------------------------------------
v = df(2);
a = 4*log(2);
b = exp(gammaln((v+1)/2) - gammaln(v/2));
c = (1+t.^2/v).^((1-v)/2);
EC(1,:) = 1 - spm_Tcdf(t,v);
EC(2,:) = a^(1/2)/(2*pi)*c;
EC(3,:) = a/((2*pi)^(3/2))*c.*t/((v/2)^(1/2))*b;
EC(4,:) = a^(3/2)/((2*pi)^2)*c.*((v-1)*(t.^2)/v - 1);
elseif STAT == 'X'
% X - Field
%-------------------------------------------------------------------
v = df(2);
a = (4*log(2))/(2*pi);
b = t.^(1/2*(v - 1)).*exp(-t/2-gammaln(v/2))/2^((v-2)/2);
EC(1,:) = 1 - spm_Xcdf(t,v);
EC(2,:) = a^(1/2)*b;
EC(3,:) = a*b.*(t-(v-1));
EC(4,:) = a^(3/2)*b.*(t.^2-(2*v-1)*t+(v-1)*(v-2));
elseif STAT == 'F'
% F Field
%-------------------------------------------------------------------
k = df(1);
v = df(2);
a = (4*log(2))/(2*pi);
b = gammaln(v/2) + gammaln(k/2);
EC(1,:) = 1 - spm_Fcdf(t,df);
EC(2,:) = a^(1/2)*exp(gammaln((v+k-1)/2)-b)*2^(1/2)...
*(k*t/v).^(1/2*(k-1)).*(1+k*t/v).^(-1/2*(v+k-2));
EC(3,:) = a*exp(gammaln((v+k-2)/2)-b)*(k*t/v).^(1/2*(k-2))...
.*(1+k*t/v).^(-1/2*(v+k-2)).*((v-1)*k*t/v-(k-1));
EC(4,:) = a^(3/2)*exp(gammaln((v+k-3)/2)-b)...
*2^(-1/2)*(k*t/v).^(1/2*(k-3)).*(1+k*t/v).^(-1/2*(v+k-2))...
.*((v-1)*(v-2)*(k*t/v).^2-(2*v*k-v-k-1)*(k*t/v)+(k-1)*(k-2));
end
