1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
|
unit fdr;
interface
uses define_types;
procedure EstimateFDR2(lnTests: integer; var Ps: SingleP; var lFDR05, lFDR01,lnegFDR05, lnegFDR01: double);
procedure qsort(lower, upper : integer; var Data:SingleP);
implementation
procedure qsort(lower, upper : integer; var Data:SingleP);
//40ms - very recursive...
var
left, right : integer;
pivot,lswap: single;
begin
pivot:=Data^[(lower+upper) div 2];
left:=lower;
right:=upper;
while left<=right do begin
while Data^[left] < pivot do left:=left+1; { Parting for left }
while Data^[right] > pivot do right:=right-1;{ Parting for right}
if left<=right then begin { Validate the change }
lswap := Data^[left];
Data^[left] := Data^[right];
Data^[right] := lswap;
left:=left+1;
right:=right-1;
end; //validate
end;//while left <=right
if right>lower then qsort(lower,right,Data); { Sort the LEFT part }
if upper>left then qsort(left ,upper,data); { Sort the RIGHT part }
end;
procedure EstimateFDR2(lnTests: integer; var Ps: SingleP; var lFDR05, lFDR01,lnegFDR05, lnegFDR01: double);
var
lInc: integer;
lrPs,Qs: SingleP;
begin
//rank Pvalues
//ShaQuickSort(lnTests,Singlep0(Ps[1]));
qSort(1,lnTests,Ps);
//qcksrt(1,lnTests,Ps);
GetMem(Qs,lnTests*sizeof(single));
//next findcrit FDR05
for lInc := 1 to lnTests do
Qs^[lInc] := (0.05*lInc)/lnTests;
lFDR05 := 0;
for lInc := 1 to lnTests do
if Ps^[lInc] <= Qs^[lInc] then
lFDR05 := Ps^[lInc];
//next findcrit FDR01
for lInc := 1 to lnTests do
Qs^[lInc] := (0.01*lInc)/lnTests;
lFDR01 := 0;
for lInc := 1 to lnTests do
if Ps^[lInc] <= Qs^[lInc] then
lFDR01 := Ps^[lInc];
//reverse
GetMem(lrPs,lnTests*sizeof(single));
for lInc := 1 to lnTests do
lrPs^[lInc] := 1- Ps^[lnTests-lInc+1];
for lInc := 1 to lnTests do
Qs^[lInc] := (0.05*lInc)/lnTests;
lnegFDR05 := 0;
for lInc := 1 to lnTests do
if lrPs^[lInc] <= Qs^[lInc] then
lnegFDR05 := lrPs^[lInc];
//next findcrit FDR01
for lInc := 1 to lnTests do
Qs^[lInc] := (0.01*lInc)/lnTests;
lnegFDR01 := 0;
for lInc := 1 to lnTests do
if lrPs^[lInc] <= Qs^[lInc] then
lnegFDR01 := lrPs^[lInc];
FreeMem(lrPs);
Freemem(Qs);
end;
end.
|
