function [pID,pN,pA] = fdr(p,q)
% FORMAT pt = FDR(p,q)
% 
% p   - vector of p-values
% q   - False Discovery Rate level
%
% pID - p-value threshold based on independence or positive dependence
% pN  - Nonparametric p-value threshold
%pA - adaptive p-value threshold


p = sort(p(:));
V = length(p);
I = (1:V)';

cVID = 1;
cVN = sum(1./(1:V));

pID = p(max(find(p<=I/V*q/cVID)));

pN = p(max(find(p<=I/V*q/cVN)));

cA = (V - max(find(p<=I/V*q/cVID)))/(V*(1-q));

if cA<1 & cA>0
    pA = p(max(find(p<=I/V*q/cA)));
else
    pA = pID;
end