Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftelschow
ftelschow committed Jun 25, 2018
1 parent 8718dc5 commit 95cc2d5
Show file tree
Hide file tree
Showing 21 changed files with 2,728 additions and 0 deletions.
96 changes: 96 additions & 0 deletions Adler_peakFDR.m
View file
@@ -0,0 +1,96 @@
function [Ps, I_crit, u, Loc] = Adler_peakFDR( Z, D, v, q, Loc )
%__________________________________________________________________________
% Computes p-values, the critical threshold and detected peaks in smooth Gaussian
% fields using the overshoot distribution in Theorem 3.6.1 Adler(1981) for
% stationary fields. Note that it has been recently established in CS(2015)
% that it this approximation is also valid for non-stationary fields.
% Inference is done by the BH procedure.
%
% Input:
% Z - Field over R^D or if 'Loc' is provided this is a vector containing
% the values of the local maxima!
% D - dimension of the domain of the field Z , only necessary input, if
% 'Loc' is provided
% v - pre-threshold for peaks
% q - FDR control threshold, number between 0 and 1
% Loc - If this is provided these are the locations of the local maxima
% and it is assumed that Z is a vector of local maxima instead of
% a field! (optional)
%
% Output:
% Ps - vector containing the p-values of the local maxima of Zs which
% are larger than v
% I_crit - critical index for significant peaks, all i <= I_crit are
% detections
% u - height of smallest significant peak after FDR control
% Loc - index location of the sorted p-values (use ind2sub on this vector
% to obtain the coordinates of the peaks)
%__________________________________________________________________________
% References:
% Adler, R. "The Geometry of Random Fields.", (1981), Vol. 62, SIAM.
% Cheng, D., Schwartzman, A., (2015), Distribution of the height of local maxima of Gaussian
% random fields, Extremes 18 (2), 213-240.
%__________________________________________________________________________
% Author: Fabian Telschow (ftelschow@ucsd.edu)
% Last changes: 06/19/18
%__________________________________________________________________________
%
% Start of function
%
if nargin == 5
% Check Input in case Loc is provided
if length(size(Z))> 2 || (all(size(Z))==1 && D==2)
error('Z, must be a vector containing the local maxima values, if Loc is provided.')
end
elseif nargin == 4
sZ = size(Z);

% Find local maxima above threshold
Maxima = zeros(sZ);
Imax = imregionalmax(Z);
% Cut off peaks at the boundary
switch D
case 2
Maxima( 2:(end-1), 2:(end-1)) = Imax( 2:(end-1), 2:(end-1) );
case 3
Maxima( 2:(end-1), 2:(end-1), 2:(end-1)) = Imax( 2:(end-1), ...
2:(end-1), ...
2:(end-1));
end

Z = Z( logical(Maxima) );
I = find(Z >= v);
Z = Z(I);

% save location of the maxima
Loc = find(Maxima);
Loc = Loc(I);
else
error('Please, input the appropriate number of input arguments.')
end

% Initialize p-value vector
Ps = zeros(1, length(Z));

% Compute p-values using Adler 2010
for k = 1:length(Z)
Ps(k) = (Z(k)/v)^(D-1) .* exp(-Z(k)^2/2) ./ exp(-v^2/2);
end

% sort p-values and report Locations, P values and critical threshold
[Ps, I] = sort(Ps);
Loc = Loc(I);

%%%% Which peaks are the discoveries?
% Find FDR threshold
S = length(Ps);
Fi = (1:S)/S*q;
% maximal index declaring detections after FDR
I_crit = find( Ps<= Fi, 1,'last' );
u = Z(I_crit);
if isempty(I_crit)
I_crit = 0;
end
if isempty(u)
u = NaN;
end
95 changes: 95 additions & 0 deletions CS_peakFDR.m
View file
@@ -0,0 +1,95 @@
function [Ps, I_crit, u, Loc] = CS_peakFDR( Z, D, v, q, kappa, Loc )
%__________________________________________________________________________
% Computes p-values, the critical threshold and detected peaks in smooth Gaussian
% fields using the STEM algorithm as in Cheng Schwartzman (2017)
%
% Input:
% Z - Field over R^D or if 'Loc' is provided this is a vector containing
% the values of the local maxima!
% D - dimension of the domain
% v - pre-threshold for peaks
% q - FDR control threshold (number between 0 and 1)
% kappa - estimated or true value if known of kappa for Z.
% Loc - If this is provided these are the locations of the local maxima
% and it is assumed that Z is a vector of local maxima instead of a
% field! (optional)
%
% Output:
% Ps - vector containing the p-values of the local maxima which are
% larger than v
% I_crit - critical index for significant peaks, all i <= I_crit are
% detections
% u - height of smallest significant peak after FDR control
% Loc - index location of the sorted p-values (use ind2sub on this vector
% to obtain the coordinates of the peaks)
%__________________________________________________________________________
% References:
% Cheng, Dan, and Armin Schwartzman. "Multiple testing of local maxima for
% detection of peaks in random fields." The Annals of Statistics 45.2
% (2017): 529-556.
%__________________________________________________________________________
% Author: Fabian Telschow (ftelschow@ucsd.edu)
% Last changes: 06/19/2018
%__________________________________________________________________________
%
% Start of function
%
if nargin == 6
% Check Input in case Loc is provided
if length(size(Z))> 2 || (all(size(Z))==1 && D==2)
error('Z, must be a vector containing the local maxima values, if Loc is provided.')
end
elseif nargin == 5
% Find local maxima above threshold
Maxima = zeros(size(Z));
Imax = imregionalmax(Z);
% Cut off peaks at the boundary
switch D
case 2
Maxima( 2:(end-1), 2:(end-1)) = Imax( 2:(end-1), 2:(end-1) );
case 3
Maxima( 2:(end-1), 2:(end-1), 2:(end-1)) = Imax( 2:(end-1), ...
2:(end-1), ...
2:(end-1));
end
Z = Z( logical(Maxima) );
I = find(Z >= v);
Z = Z(I);

% save location of the maxima
Loc = find(Maxima);
Loc = Loc(I);
else
error('Please, input the appropriate number of input arguments.')
end


% Define the peak height distribution as given in Proposition 6 CS(2017)
peakHeightDistr = @(x, D, kappa) integral( peakHeightDensity(D, kappa), x,...
Inf, 'ArrayValued',true);

% compute the p values of the peaks using Proposition 6 CS(2017)
Fkappa_v = peakHeightDistr(v, D, kappa);
Ps = zeros(1, length(I));

for k = 1:length(I)
Ps(k) = peakHeightDistr(Z(k), D, kappa) / Fkappa_v;
end

% sort p-values and report Locations, P values and critical threshold
[Ps, I] = sort(Ps);
Loc = Loc(I);

%%%% Which peaks are the discoveries?
% Find FDR threshold
S = length(Ps);
Fi = (1:S)/S*q;
% maximal index declaring detections after FDR
I_crit = find( Ps <= Fi, 1,'last' );
u = Z(I_crit);
if isempty(I_crit)
I_crit = 0;
end
if isempty(u)
u = NaN;
end
112 changes: 112 additions & 0 deletions Chumbley_peakFDR.m
View file
@@ -0,0 +1,112 @@
function [Ps, I_crit, u, Loc] = Chumbley_peakFDR( Z, D, v, q, STAT, df, Loc )
%__________________________________________________________________________
% Computes p-values, the critical threshold and detected peaks in smooth Gaussian
% fields using the method in Chumbley (2010) and applies BH procedure for
% inference.
%
% Input:
% Z - Field over R^D or if 'Loc' is provided this is a vector containing
% the values of the local maxima!
% D - dimension of the domain
% v - pre-threshold for peaks
% q - FDR control threshold (number between 0 and 1)
% STAT - either 'Z' or 'T' specifying whether the GKF for Gaussian or
% T-fields is used.
% df - degrees of freedom of the 'T'-field, put arbitrary number if 'Z'
% Loc - If this is provided these are the locations of the local maxima
% and it is assumed that Z is a vector of local maxima instead of a
% field! (optional)
%
% Output:
% Ps - vector containing the p-values of the local maxima which are
% larger than v
% I_crit - critical index for significant peaks, all i <= I_crit are
% detections
% u - height of smallest significant peak after FDR control
% Loc - index location of the sorted p-values (use ind2sub on this vector
% to obtain the coordinates of the peaks)
%__________________________________________________________________________
% References:
% Chumbley, J., et al. "Topological FDR for neuroimaging." Neuroimage 49.4
% (2010): 3057-3064.
%__________________________________________________________________________
% Author: Fabian Telschow (ftelschow@ucsd.edu)
% Last changes: 06/19/2018
%__________________________________________________________________________
%
% Start of function
%
if nargin == 7
% Check Input in case Loc is provided
if length(size(Z))> 2 || (all(size(Z))==1 && D==2)
error('Z, must be a vector containing the local maxima values, if Loc is provided.')
end
elseif nargin == 6
% Find local maxima above threshold
Maxima = zeros(size(Z));
Imax = imregionalmax(Z);
% Cut off peaks at the boundary
switch D
case 2
Maxima( 2:(end-1), 2:(end-1)) = Imax( 2:(end-1), 2:(end-1) );
case 3
Maxima( 2:(end-1), 2:(end-1), 2:(end-1)) = Imax( 2:(end-1), ...
2:(end-1), ...
2:(end-1));
end
Z = Z( logical(Maxima) );
I = find(Z >= v);
Z = Z(I);

% save location of the maxima
Loc = find(Maxima);
Loc = Loc(I);
else
error('Please, input the appropriate number of input arguments.')
end

% Initialize p-value vector
Ps = zeros(1, length(Z));

% Compute p-svalues using Chumbley 2010
switch STAT,
case 'Z',
if(D == 2),
H = @(x) x;
else
H = @(x) x.^2 - 1;
end
EECratio = @(u, v) ( H(u).*exp(-u.^2/2) ) ./ ( H(v).*exp(-v.^2/2) );
for k = 1:length(Z)
Ps(k) = EECratio(Z(k),v);
end
case 'T',
if(D == 2),
H = @(x) x;
else
H = @(x) (df-1)*x.^2/df - 1;
end
EECratio = @(u, v) ( H(u).*(1+u.^2/df)^(-(df-1)/2) ) ./ ...
( H(v).*(1+v.^2/df)^(-(df-1)/2) );
for k = 1:length(Z)
Ps(k) = EECratio(Z(k),v);
end
end

% sort p-values and report Locations, P values and critical threshold
[Ps, I] = sort(Ps);
Loc = Loc(I);

%%%% Which peaks are the discoveries?
% Find FDR threshold
S = length(Ps);
Fi = (1:S)/S*q;
% maximal index declaring detections after FDR
I_crit = find( Ps<= Fi, 1,'last' );
u = Z(I_crit);
if isempty(I_crit)
I_crit = 0;
end
if isempty(u)
u = NaN;
end
56 changes: 56 additions & 0 deletions PvalueTable_heightDistr.m
View file
@@ -0,0 +1,56 @@
function pValueTable = PvalueTable_heightDistr(D, kappa, du, umin, umax, show)
%__________________________________________________________________________
% Computes a p-value table for the height distribution using the formulas
% derived in Cheng Schwartzman (2015/2018) and is merely used to significantly
% speed up the computation time of CS_PeakFDR.m by using it as an input into
% PvalueTable_heightDistr.m, since CS_PeakFDR.m is slow for large numbers of peaks.
%
% Input:
% D - dimension of the domain
% kappa - value of kappa
% du - stepsize for integration
% umin - lower bound of integral
% umax - upper bound of integral
% show - boolean for plotting figures
%
% Output:
% pValueTable - first column are the heights and second column are the
% corresponding p-values
%
%__________________________________________________________________________
% References:
% Cheng, D., Schwartzman, A., 2015. On the explicit height distribution and expected number
% of local maxima of isotropic Gaussian random fields. arXiv preprint arXiv:1503.01328.
% Cheng, D., Schwartzman, A., 2018. Expected number and height distribution of critical
% points of smooth isotropic Gaussian random fields. Bernoulli 24 (4B), 3422-3446.
%__________________________________________________________________________
% Depends on,
% - peakHeightDensity.m
%__________________________________________________________________________
% Author: Fabian Telschow (ftelschow@ucsd.edu)
% Last changes: 06/19/2018
%__________________________________________________________________________
%
% Start of function
%
if nargin < 6
% estimate variance of the field, if unknown
show = 0;
end

density = peakHeightDensity(D, kappa);

u = umin:du:umax;
d = density(u);

Int = cumtrapz(d*du);

% plot the density to check whether the min and max are large enough
if show
figure(1)
plot(u,density(u))
figure(2)
plot(u,1-Int)
end

pValueTable = [ u' 1-Int ];
Binary file added README.pdf
View file
Binary file not shown.

0 comments on commit 95cc2d5

Please sign in to comment.