spm_P_RF.m 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112
  1. function [P,p,Ec,Ek] = spm_P_RF(c,k,Z,df,STAT,R,n)
  2. % Returns the [un]corrected P value using unifed EC theory
  3. % FORMAT [P p Ec Ek] = spm_P_RF(c,k,z,df,STAT,R,n)
  4. %
  5. % c - cluster number
  6. % k - extent {RESELS}
  7. % z - height {minimum over n values}
  8. % df - [df{interest} df{error}]
  9. % STAT - Statistical field
  10. % 'Z' - Gaussian field
  11. % 'T' - T - field
  12. % 'X' - Chi squared field
  13. % 'F' - F - field
  14. % R - RESEL Count {defining search volume}
  15. % n - number of component SPMs in conjunction
  16. %
  17. % P - corrected P value - P(C >= c | K >= k}
  18. % p - uncorrected P value
  19. % Ec - expected number of clusters (maxima)
  20. % Ek - expected number of resels per cluster
  21. %
  22. %__________________________________________________________________________
  23. %
  24. % spm_P_RF returns the probability of c or more clusters with more than
  25. % k resels in volume process of R RESELS thresholded at u. All p values
  26. % can be considered special cases:
  27. %
  28. % spm_P_RF(1,0,z,df,STAT,1,n) = uncorrected p value
  29. % spm_P_RF(1,0,z,df,STAT,R,n) = corrected p value {based on height z)
  30. % spm_P_RF(1,k,u,df,STAT,R,n) = corrected p value {based on extent k at u)
  31. % spm_P_RF(c,k,u,df,STAT,R,n) = corrected p value {based on number c at k and u)
  32. % spm_P_RF(c,0,u,df,STAT,R,n) = omnibus p value {based on number c at u)
  33. %
  34. % If n > 1 a conjunction probility over the n values of the statistic
  35. % is returned.
  36. %__________________________________________________________________________
  37. %
  38. % References:
  39. %
  40. % [1] Hasofer AM (1978) Upcrossings of random fields
  41. % Suppl Adv Appl Prob 10:14-21
  42. % [2] Friston KJ et al (1994) Assessing the Significance of Focal Activations
  43. % Using Their Spatial Extent
  44. % Human Brain Mapping 1:210-220
  45. % [3] Worsley KJ et al (1996) A Unified Statistical Approach for Determining
  46. % Significant Signals in Images of Cerebral Activation
  47. % Human Brain Mapping 4:58-73
  48. %__________________________________________________________________________
  49. % Copyright (C) 1999-2013 Wellcome Trust Centre for Neuroimaging
  50. % Karl Friston
  51. % $Id: spm_P_RF.m 5770 2013-11-27 20:12:29Z karl $
  52. %-Get expectations
  53. %==========================================================================
  54. % get EC densities
  55. %--------------------------------------------------------------------------
  56. D = find(R,1,'last');
  57. R = R(1:D);
  58. G = sqrt(pi)./gamma((1:D)/2);
  59. EC = spm_ECdensity(STAT,Z,df);
  60. EC = max(EC(1:D),eps);
  61. % corrected p value
  62. %--------------------------------------------------------------------------
  63. P = triu(toeplitz(EC'.*G))^n;
  64. P = P(1,:);
  65. EM = (R./G).*P; % <maxima> over D dimensions
  66. Ec = sum(EM); % <maxima>
  67. EN = P(1)*R(D); % <resels>
  68. Ek = EN/EM(D); % Ek = EN/EM(D);
  69. %-Get P{n > k}
  70. %==========================================================================
  71. % assume a Gaussian form for P{n > k} ~ exp(-beta*k^(2/D))
  72. % Appropriate for SPM{Z} and high d.f. SPM{T}
  73. %--------------------------------------------------------------------------
  74. D = D - 1;
  75. if ~k || ~D
  76. p = 1;
  77. elseif STAT == 'Z'
  78. beta = (gamma(D/2 + 1)/Ek)^(2/D);
  79. p = exp(-beta*(k^(2/D)));
  80. elseif STAT == 'T'
  81. beta = (gamma(D/2 + 1)/Ek)^(2/D);
  82. p = exp(-beta*(k^(2/D)));
  83. elseif STAT == 'X'
  84. beta = (gamma(D/2 + 1)/Ek)^(2/D);
  85. p = exp(-beta*(k^(2/D)));
  86. elseif STAT == 'F'
  87. beta = (gamma(D/2 + 1)/Ek)^(2/D);
  88. p = exp(-beta*(k^(2/D)));
  89. end
  90. %-Poisson clumping heuristic {for multiple clusters}
  91. %==========================================================================
  92. P = 1 - spm_Pcdf(c - 1,(Ec + eps)*p);