fdr_thresholds.m

来自「toolbox of BVQX, This is the access betw」· M 代码 · 共 79 行

function th = fdr_thresholds(plist, levels, lnmeth)
% fdr_thresholds  - compute probability thresholds for q(FDR) levels
%
% FORMAT:       th = fdr_thresholds(plist, levels [, lnmeth])
%
% Input fields:
%
%       plist       list of p-values to be thresholded
%       levels      q(FDR) levels
%       lnmeth      if given and true (boolean!), use
%                   c(V) = ln(V) + E  instead of  c(V) = 1 approximation
%
% Output fields:
%
%       th          p-thresholds for q(FDR)

% Version:  v0.6f
% Build:    7070509
% Date:     Jul-05 2007, 9:32 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   (~isa(plist, 'double') && ...
    ~isa(plist, 'single')) || ...
    numel(plist) ~= length(plist) || ...
    any(plist < 0) || ...
   (~isa(levels, 'double') && ...
    ~isa(levels, 'single')) || ...
    isempty(levels) || ...
    numel(levels) ~= length(levels) || ...
    any(levels < 0 | levels >= 1)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing first argument.' ...
    );
end

% sort plist
splist = sort(plist);
nlist = numel(splist);
levels = levels(:)';
th = zsz(levels);

% divide with FDR factor list for fast find access
fsplist = splist(:)' ./ ((1/nlist):(1/nlist):1);
if nargin > 2 && ...
    islogical(lnmeth) && ...
   ~isempty(lnmeth) && ...
    lnmeth(1)
    th = [th(:), th(:)];
    lnmeth = true;
else
    lnmeth = false;
end

% find thresholds
for thc = 1:numel(levels)
    thi = find(fsplist > levels(thc));
    if isempty(thi)
        th(thc) = splist(end);
    elseif thi(1) > 1
        th(thc) = splist(thi(1) - 1);
    else
        th(thc) = splist(1) / 2.001;
    end
    if lnmeth
        thi = find(fsplist > (levels(thc) / (log(nlist) + 0.5772)));
        if isempty(thi)
            th(thc, 2) = splist(end);
        elseif thi(1) > 1
            th(thc, 2) = splist(thi(1) - 1);
        else
            th(thc, 2) = splist(1) / 2.001;
        end
    end
end