function [P] = spm_P_FDR(Z,df,STAT,n,Ps)
% Returns the corrected FDR P value 
% FORMAT [P] = spm_P_FDR(Z,df,STAT,n,Ps)
%
% Z     - height (minimum of n statistics)
% df    - [df{interest} df{error}]
% STAT  - Statisical field
%		'Z' - Gaussian field
%		'T' - T - field
%		'X' - Chi squared field
%		'F' - F - field
%		'P' - P-values
% n     - number of component SPMs in conjunction
% Vs    - Vector of sorted (ascending) p-values in search volume
%
% P     - corrected FDR   P value
%
%___________________________________________________________________________
%
% The Benjamini & Hochberch (1995) False Discovery Rate (FDR) procedure
% finds a threshold u such that the expected FDR is at most q.  spm_P_FDR
% returns the smallest q such that Z>u. 
%
% If abs(n) > 1, a P-value for a minimum of n values of the statistic 
% is returned.  If n>0, the P-value assesses the conjunction null
% hypothesis of one or more of the null hypotheses being true.  If n<0,
% then the P-value assess the global null of all nulls being true. 
%
% FDR Background
%
% For a given threshold on a statistic image, the False Discovery Rate
% is the proportion of suprathreshold voxels which are false positives.
% Recall that the thresholding of each voxel consists of a hypothesis
% test, where the null hypothesis is rejected if the statistic is larger
% than threshold.  In this terminology, the FDR is the proportion of
% rejected tests where the null hypothesis is actually true.
%
% A FDR proceedure produces a threshold that controls the expected FDR
% at or below q.  The FDR adjusted p-value for a voxel is the smallest q
% such that the voxel would be suprathreshold.
%
% In comparison, a traditional multiple comparisons proceedure
% (e.g. Bonferroni or random field methods) controls Familywise Error
% rate (FWER) at or below alpha.  FWER is the *chance* of one or more
% false positives anywhere (whereas FDR is a *proportion* of false
% positives).  A FWER adjusted p-value for a voxel is the smallest alpha
% such that the voxel would be suprathreshold. 
%
% If there is truely no signal in the image anywhere, then a FDR
% proceedure controls FWER, just as Bonferroni and random field methods
% do. (Precisely, controlling E(FDR) yields weak control of FWE).  If
% there is some signal in the image, a FDR method should be more powerful
% than a traditional method.
%
% For careful definition of FDR-adjusted p-values (and distinction between
% corrected and adjusted p-values) see Yekutieli & Benjamini (1999).
%
%
% References
%
% Benjamini & Hochberg (1995), "Controlling the False Discovery Rate: A
% Practical and Powerful Approach to Multiple Testing". J Royal Stat Soc,
% Ser B.  57:289-300.
%
% Benjamini & Yekutieli (2001), "The Control of the false discovery rate
% in multiple testing under dependency". To appear, Annals of Statistics.
% Available at http://www.math.tau.ac.il/~benja 
%
% Yekutieli & Benjamini (1999). "Resampling-based false discovery rate
% controlling multiple test procedures for correlated test
% statistics".  J of Statistical Planning and Inference, 82:171-196.
%___________________________________________________________________________
% @(#)spm_P_FDR.m	2.4 Thomas Nichols 03/01/31
% $Id: spm_P_FDR.m,v 1.2 2004/07/06 19:55:55 nichols Exp $ UM Biostat

if n>0
  Cnj_n = 1;       % Inf on Conjunction Null
else
  Cnj_n = abs(n);  % Inf on Global Null
end

% Set Benjamini & Yeuketeli cV for independence/PosRegDep case
%-----------------------------------------------------------------------
cV   = 1; 
S    = length(Ps);

% Calculate p value of Z
%-----------------------------------------------------------------------
if      STAT == 'Z'
  PZ = (1 - spm_Ncdf(Z)).^Cnj_n;
elseif  STAT == 'T'
  PZ = (1 - spm_Tcdf(Z,df(2))).^Cnj_n;
elseif  STAT == 'X'
  PZ = (1 - spm_Xcdf(Z,df(2))).^Cnj_n;
elseif  STAT == 'F'
  PZ = (1-spm_Fcdf(Z,df)).^Cnj_n;
elseif  STAT == 'P'
  PZ = Z;
end

%-Calculate FDR p values
%-----------------------------------------------------------------------
% If Z is a value in the statistic image, then the adjusted p-value
% defined in Yekutieli & Benjamini (1999) (eqn 3) is obtained.  If Z
% isn't a value in the image, then the adjusted p-value for the next
% smallest statistic value (next largest uncorrected p) is returned.

%-"Corrected" p-values
%-----------------------------------------------------------------------
Qs    = Ps*S./(1:S)'*cV;    

% -"Adjusted" p-values
%-----------------------------------------------------------------------
P     = zeros(size(Z));
for i = 1:length(Z)

	% Find PZ(i) in Ps, or smallest Ps(j) s.t. Ps(j) >= PZ(i)
	%---------------------------------------------------------------
	I = min(find(Ps>=PZ(i)));  

	% -"Adjusted" p-values
	%---------------------------------------------------------------
	if isempty(I)
		P(i) = 1;
	else
		P(i) = min(Qs(I:S));
	end
end