341 lines
13 KiB
D
341 lines
13 KiB
D
import std.algorithm;
|
|
import std.conv;
|
|
import std.math;
|
|
import std.stdio;
|
|
import std.string;
|
|
|
|
int[] seqLen(int start, int end) {
|
|
int[] result;
|
|
if (start == end) {
|
|
result.length = end+1;
|
|
for (int i; i<result.length; i++) {
|
|
result[i] = i+1;
|
|
}
|
|
} else if (start < end) {
|
|
result.length = end - start + 1;
|
|
for (int i; i<result.length; i++) {
|
|
result[i] = start+i;
|
|
}
|
|
} else {
|
|
result.length = start - end + 1;
|
|
for (int i; i<result.length; i++) {
|
|
result[i] = start-i;
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
int[] order(double[] array, bool decreasing) {
|
|
int size = array.length;
|
|
int[] idx;
|
|
idx.length = size;
|
|
double[] baseArr;
|
|
baseArr.length = size;
|
|
for (int i; i<size; i++) {
|
|
baseArr[i] = array[i];
|
|
idx[i] = i;
|
|
}
|
|
if (!decreasing) {
|
|
alias comp = (a,b) => baseArr[a] < baseArr[b];
|
|
idx.sort!comp;
|
|
} else {
|
|
alias comp = (a,b) => baseArr[b] < baseArr[a];
|
|
idx.sort!comp;
|
|
}
|
|
return idx;
|
|
}
|
|
|
|
double[] cummin(double[] array) {
|
|
int size = array.length;
|
|
if (size < 1) throw new Exception("cummin requires at least one element");
|
|
double[] output;
|
|
output.length = size;
|
|
auto cumulativeMin = array[0];
|
|
foreach (i; 0..size) {
|
|
if (array[i] < cumulativeMin) cumulativeMin = array[i];
|
|
output[i] = cumulativeMin;
|
|
}
|
|
return output;
|
|
}
|
|
|
|
double[] cummax(double[] array) {
|
|
auto size = array.length;
|
|
if (size < 1) throw new Exception("cummax requires at least one element");
|
|
double[] output;
|
|
output.length = size;
|
|
auto cumulativeMax = array[0];
|
|
foreach (i; 0..size) {
|
|
if (array[i] > cumulativeMax) cumulativeMax = array[i];
|
|
output[i] = cumulativeMax;
|
|
}
|
|
return output;
|
|
}
|
|
|
|
double[] pminx(double[] array, double x) {
|
|
auto size = array.length;
|
|
if (size < 1) throw new Exception("pmin requires at least one element");
|
|
double[] result;
|
|
result.length = size;
|
|
foreach (i; 0..size) {
|
|
if (array[i] < x) {
|
|
result[i] = array[i];
|
|
} else {
|
|
result[i] = x;
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
void doubleSay(double[] array) {
|
|
writef("[ 1] %e", array[0]);
|
|
foreach (i; 1..array.length) {
|
|
writef(" %.10f", array[i]);
|
|
if ((i+1) % 5 == 0) writef("\n[%2d]", i+1);
|
|
}
|
|
writeln;
|
|
}
|
|
|
|
auto toArray(T,U)(U[] array) {
|
|
T[] result;
|
|
result.length = array.length;
|
|
foreach(i; 0..array.length) {
|
|
result[i] = to!T(array[i]);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
double[] pAdjust(double[] pvalues, string str) {
|
|
auto size = pvalues.length;
|
|
if (size < 1) throw new Exception("pAdjust requires at least one element");
|
|
int type = str.toLower.predSwitch!"a==b"(
|
|
"bh", 0,
|
|
"fdr", 0,
|
|
"by", 1,
|
|
"bonferroni", 2,
|
|
"hochberg", 3,
|
|
"holm", 4,
|
|
"hommel", 5,
|
|
{ throw new Exception(text("'",str,"' doesn't match any accepted FDR types")); }()
|
|
);
|
|
if (type == 2) { // Bonferroni method
|
|
double[] result;
|
|
result.length = size;
|
|
foreach (i; 0..size) {
|
|
auto b = pvalues[i] * size;
|
|
if (b >= 1) {
|
|
result[i] = 1;
|
|
} else if (0 <= b && b < 1) {
|
|
result[i] = b;
|
|
} else {
|
|
throw new Exception(text(b," is outside [0, 1)"));
|
|
}
|
|
}
|
|
return result;
|
|
} else if (type == 4) { // Holm method
|
|
auto o = order(pvalues, false);
|
|
auto o2Double = toArray!(double,int)(o);
|
|
double[] cummaxInput;
|
|
cummaxInput.length = size;
|
|
foreach (i; 0..size) {
|
|
cummaxInput[i] = (size-i) * pvalues[o[i]];
|
|
}
|
|
auto ro = order(o2Double, false);
|
|
auto cummaxOutput = cummax(cummaxInput);
|
|
auto pmin = pminx(cummaxOutput, 1.0);
|
|
double[] result;
|
|
result.length = size;
|
|
foreach (i; 0..size) {
|
|
result[i] = pmin[ro[i]];
|
|
}
|
|
return result;
|
|
} else if (type == 5) {
|
|
auto indices = seqLen(size, size);
|
|
auto o = order(pvalues, false);
|
|
double[] p;
|
|
p.length = size;
|
|
foreach (i; 0..size) {
|
|
p[i] = pvalues[o[i]];
|
|
}
|
|
auto o2Double = toArray!double(o);
|
|
auto ro = order(o2Double, false);
|
|
double[] q;
|
|
q.length = size;
|
|
double[] pa;
|
|
pa.length = size;
|
|
double[] npi;
|
|
npi.length = size;
|
|
foreach (i; 0..size) {
|
|
npi[i] = p[i] * size / indices[i];
|
|
}
|
|
auto min_ = reduce!min(npi);
|
|
q[] = min_;
|
|
pa[] = min_;
|
|
foreach_reverse (j; 2..size) {
|
|
auto ij = seqLen(1, size - j + 1);
|
|
foreach (i; 0..size-j+1) {
|
|
ij[i]--;
|
|
}
|
|
auto i2Length = j-1;
|
|
int[] i2;
|
|
i2.length = i2Length;
|
|
foreach(i; 0..i2Length) {
|
|
i2[i] = size-j+2+i-1;
|
|
}
|
|
auto pi2Length = i2Length;
|
|
double q1 = j*p[i2[0]] / 2.0;
|
|
foreach (i; 1..pi2Length) {
|
|
auto temp_q1 = p[i2[i]] * j / (2.0 + i);
|
|
if (temp_q1 < q1) q1 = temp_q1;
|
|
}
|
|
foreach (i; 0..size-j+1) {
|
|
q[ij[i]] = min(p[ij[i]] * j, q1);
|
|
}
|
|
foreach(i; 0..i2Length) {
|
|
q[i2[i]] = q[size-j];
|
|
}
|
|
foreach(i; 0..size) if (pa[i] < q[i]) pa[i] = q[i];
|
|
}
|
|
foreach (index; 0..size) {
|
|
q[index] = pa[ro[index]];
|
|
}
|
|
return q;
|
|
}
|
|
|
|
double[] ni;
|
|
ni.length = size;
|
|
auto o = order(pvalues, true);
|
|
auto oDouble = toArray!double(o);
|
|
foreach (index; 0..size) {
|
|
if (pvalues[index] < 0 || pvalues[index] > 1) {
|
|
throw new Exception(text("array[", index, "] = ", pvalues[index], " is outside [0, 1]"));
|
|
}
|
|
ni[index] = cast(double) size / (size - index);
|
|
}
|
|
auto ro = order(oDouble, false);
|
|
double[] cumminInput;
|
|
cumminInput.length = size;
|
|
if (type == 0) { // BH method
|
|
foreach (index; 0..size) {
|
|
cumminInput[index] = ni[index] * pvalues[o[index]];
|
|
}
|
|
} else if (type == 1) { // BY method
|
|
double q = 0;
|
|
foreach (index; 1..size+1) q += 1.0 / index;
|
|
foreach (index; 0..size) {
|
|
cumminInput[index] = q * ni[index] * pvalues[o[index]];
|
|
}
|
|
} else if (type == 3) { // Hochberg method
|
|
foreach (index; 0..size) {
|
|
cumminInput[index] = (index + 1) * pvalues[o[index]];
|
|
}
|
|
}
|
|
auto cumminArray =cummin(cumminInput);
|
|
auto pmin = pminx(cumminArray, 1.0);
|
|
double[] result;
|
|
result.length = size;
|
|
foreach (i; 0..size) {
|
|
result[i] = pmin[ro[i]];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
void main() {
|
|
double[] pvalues = [
|
|
4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
|
|
8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01,
|
|
4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01,
|
|
8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02,
|
|
3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01,
|
|
1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02,
|
|
4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04,
|
|
3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04,
|
|
1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04,
|
|
2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03
|
|
];
|
|
|
|
double[][] correctAnswers = [
|
|
[ // Benjamini-Hochberg
|
|
6.126681e-01, 8.521710e-01, 1.987205e-01, 1.891595e-01, 3.217789e-01,
|
|
9.301450e-01, 4.870370e-01, 9.301450e-01, 6.049731e-01, 6.826753e-01,
|
|
6.482629e-01, 7.253722e-01, 5.280973e-01, 8.769926e-01, 4.705703e-01,
|
|
9.241867e-01, 6.049731e-01, 7.856107e-01, 4.887526e-01, 1.136717e-01,
|
|
4.991891e-01, 8.769926e-01, 9.991834e-01, 3.217789e-01, 9.301450e-01,
|
|
2.304958e-01, 5.832475e-01, 3.899547e-02, 8.521710e-01, 1.476843e-01,
|
|
1.683638e-02, 2.562902e-03, 3.516084e-02, 6.250189e-02, 3.636589e-03,
|
|
2.562902e-03, 2.946883e-02, 6.166064e-03, 3.899547e-02, 2.688991e-03,
|
|
4.502862e-04, 1.252228e-05, 7.881555e-02, 3.142613e-02, 4.846527e-03,
|
|
2.562902e-03, 4.846527e-03, 1.101708e-03, 7.252032e-02, 2.205958e-02
|
|
],
|
|
[ // Benjamini & Yekutieli
|
|
1.000000e+00, 1.000000e+00, 8.940844e-01, 8.510676e-01, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 5.114323e-01,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.754486e-01, 1.000000e+00, 6.644618e-01,
|
|
7.575031e-02, 1.153102e-02, 1.581959e-01, 2.812089e-01, 1.636176e-02,
|
|
1.153102e-02, 1.325863e-01, 2.774239e-02, 1.754486e-01, 1.209832e-02,
|
|
2.025930e-03, 5.634031e-05, 3.546073e-01, 1.413926e-01, 2.180552e-02,
|
|
1.153102e-02, 2.180552e-02, 4.956812e-03, 3.262838e-01, 9.925057e-02
|
|
],
|
|
[ // Bonferroni
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 7.019185e-01, 1.000000e+00, 1.000000e+00,
|
|
2.020365e-01, 1.516674e-02, 5.625735e-01, 1.000000e+00, 2.909271e-02,
|
|
1.537741e-02, 4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02,
|
|
9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01, 4.395577e-02,
|
|
1.088915e-02, 4.846527e-02, 3.305125e-03, 1.000000e+00, 2.867745e-01
|
|
],
|
|
[ // Hochberg
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 4.632662e-01, 9.991834e-01, 9.991834e-01,
|
|
1.575885e-01, 1.383967e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02,
|
|
1.383967e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
|
|
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02,
|
|
1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01
|
|
],
|
|
[ // Holm
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
|
|
1.000000e+00, 1.000000e+00, 4.632662e-01, 1.000000e+00, 1.000000e+00,
|
|
1.575885e-01, 1.395341e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02,
|
|
1.395341e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
|
|
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02,
|
|
1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01
|
|
],
|
|
[ // Hommel
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01,
|
|
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
|
|
9.991834e-01, 9.991834e-01, 4.351895e-01, 9.991834e-01, 9.766522e-01,
|
|
1.414256e-01, 1.304340e-02, 3.530937e-01, 6.887709e-01, 2.385602e-02,
|
|
1.322457e-02, 2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02,
|
|
8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01, 3.516461e-02,
|
|
9.582456e-03, 3.877222e-02, 3.172920e-03, 8.122276e-01, 1.950067e-01
|
|
]
|
|
];
|
|
auto types = ["bh", "by", "bonferroni", "hochberg", "holm", "hommel"];
|
|
foreach (type; 0..types.length) {
|
|
auto q = pAdjust(pvalues, types[type]);
|
|
double error = 0.0;
|
|
foreach (i; 0..pvalues.length) {
|
|
error += abs(q[i] - correctAnswers[type][i]);
|
|
}
|
|
doubleSay(q);
|
|
writefln("\ntype %d = '%s' has a cumulative error of %g", type, types[type], error);
|
|
}
|
|
}
|