P-value correction: Difference between revisions
Content deleted Content added
→Version 1: new version does not show compiler warnings about ignored qualifiers or possible conversion errors. Eliminated reindexing step with Hommel method, so code is shorter and should be slightly faster with that method |
Thundergnat (talk | contribs) Rename Perl 6 -> Raku, alphabetize, minor clean-up |
||
Line 894:
</pre>
=={{header|C
{{trans|Java}}
<lang csharp>using System;
Line 1,714 ⟶ 1,306:
[50]
type 5 = 'hommel' has a cumulative error of 4.353024E-007</pre>
=={{header|C++}}==
{{trans|Java}}
<lang cpp>#include <algorithm>
#include <functional>
#include <iostream>
#include <numeric>
#include <vector>
std::vector<int> seqLen(int start, int end) {
std::vector<int> result;
if (start == end) {
result.resize(end + 1);
std::iota(result.begin(), result.end(), 1);
} else if (start < end) {
result.resize(end - start + 1);
std::iota(result.begin(), result.end(), start);
} else {
result.resize(start - end + 1);
std::iota(result.rbegin(), result.rend(), end);
}
return result;
}
std::vector<int> order(const std::vector<double>& arr, bool decreasing) {
std::vector<int> idx(arr.size());
std::iota(idx.begin(), idx.end(), 0);
std::function<bool(int, int)> cmp;
if (decreasing) {
cmp = [&arr](int a, int b) { return arr[b] < arr[a]; };
} else {
cmp = [&arr](int a, int b) { return arr[a] < arr[b]; };
}
std::sort(idx.begin(), idx.end(), cmp);
return idx;
}
std::vector<double> cummin(const std::vector<double>& arr) {
if (arr.empty()) throw std::runtime_error("cummin requries at least one element");
std::vector<double> output(arr.size());
double cumulativeMin = arr[0];
std::transform(arr.cbegin(), arr.cend(), output.begin(), [&cumulativeMin](double a) {
if (a < cumulativeMin) cumulativeMin = a;
return cumulativeMin;
});
return output;
}
std::vector<double> cummax(const std::vector<double>& arr) {
if (arr.empty()) throw std::runtime_error("cummax requries at least one element");
std::vector<double> output(arr.size());
double cumulativeMax = arr[0];
std::transform(arr.cbegin(), arr.cend(), output.begin(), [&cumulativeMax](double a) {
if (cumulativeMax < a) cumulativeMax = a;
return cumulativeMax;
});
return output;
}
std::vector<double> pminx(const std::vector<double>& arr, double x) {
if (arr.empty()) throw std::runtime_error("pmin requries at least one element");
std::vector<double> result(arr.size());
std::transform(arr.cbegin(), arr.cend(), result.begin(), [&x](double a) {
if (a < x) return a;
return x;
});
return result;
}
void doubleSay(const std::vector<double>& arr) {
printf("[ 1] %.10f", arr[0]);
for (size_t i = 1; i < arr.size(); ++i) {
printf(" %.10f", arr[i]);
if ((i + 1) % 5 == 0) printf("\n[%2d]", i + 1);
}
}
std::vector<double> pAdjust(const std::vector<double>& pvalues, const std::string& str) {
if (pvalues.empty()) throw std::runtime_error("pAdjust requires at least one element");
size_t size = pvalues.size();
int type;
if ("bh" == str || "fdr" == str) {
type = 0;
} else if ("by" == str) {
type = 1;
} else if ("bonferroni" == str) {
type = 2;
} else if ("hochberg" == str) {
type = 3;
} else if ("holm" == str) {
type = 4;
} else if ("hommel" == str) {
type = 5;
} else {
throw std::runtime_error(str + " doesn't match any accepted FDR types");
}
// Bonferroni method
if (2 == type) {
std::vector<double> result(size);
for (size_t i = 0; i < size; ++i) {
double b = pvalues[i] * size;
if (b >= 1) {
result[i] = 1;
} else if (0 <= b && b < 1) {
result[i] = b;
} else {
throw std::runtime_error("a value is outside [0, 1)");
}
}
return result;
}
// Holm method
else if (4 == type) {
auto o = order(pvalues, false);
std::vector<double> o2Double(o.begin(), o.end());
std::vector<double> cummaxInput(size);
for (size_t i = 0; i < size; ++i) {
cummaxInput[i] = (size - i) * pvalues[o[i]];
}
auto ro = order(o2Double, false);
auto cummaxOutput = cummax(cummaxInput);
auto pmin = pminx(cummaxOutput, 1.0);
std::vector<double> result(size);
std::transform(ro.cbegin(), ro.cend(), result.begin(), [&pmin](int a) { return pmin[a]; });
return result;
}
// Hommel
else if (5 == type) {
auto indices = seqLen(size, size);
auto o = order(pvalues, false);
std::vector<double> p(size);
std::transform(o.cbegin(), o.cend(), p.begin(), [&pvalues](int a) { return pvalues[a]; });
std::vector<double> o2Double(o.begin(), o.end());
auto ro = order(o2Double, false);
std::vector<double> q(size);
std::vector<double> pa(size);
std::vector<double> npi(size);
for (size_t i = 0; i < size; ++i) {
npi[i] = p[i] * size / indices[i];
}
double min = *std::min_element(npi.begin(), npi.end());
std::fill(q.begin(), q.end(), min);
std::fill(pa.begin(), pa.end(), min);
for (int j = size; j >= 2; --j) {
auto ij = seqLen(1, size - j + 1);
std::transform(ij.cbegin(), ij.cend(), ij.begin(), [](int a) { return a - 1; });
int i2Length = j - 1;
std::vector<int> i2(i2Length);
for (int i = 0; i < i2Length; ++i) {
i2[i] = size - j + 2 + i - 1;
}
double q1 = j * p[i2[0]] / 2.0;
for (int i = 1; i < i2Length; ++i) {
double temp_q1 = p[i2[i]] * j / (2.0 + i);
if (temp_q1 < q1) q1 = temp_q1;
}
for (size_t i = 0; i < size - j + 1; ++i) {
q[ij[i]] = std::min(p[ij[i]] * j, q1);
}
for (int i = 0; i < i2Length; ++i) {
q[i2[i]] = q[size - j];
}
for (size_t i = 0; i < size; ++i) {
if (pa[i] < q[i]) {
pa[i] = q[i];
}
}
}
std::transform(ro.cbegin(), ro.cend(), q.begin(), [&pa](int a) { return pa[a]; });
return q;
}
std::vector<double> ni(size);
std::vector<int> o = order(pvalues, true);
std::vector<double> od(o.begin(), o.end());
for (size_t i = 0; i < size; ++i) {
if (pvalues[i] < 0 || pvalues[i]>1) {
throw std::runtime_error("a value is outside [0, 1]");
}
ni[i] = (double)size / (size - i);
}
auto ro = order(od, false);
std::vector<double> cumminInput(size);
if (0 == type) { // BH method
for (size_t i = 0; i < size; ++i) {
cumminInput[i] = ni[i] * pvalues[o[i]];
}
} else if (1 == type) { // BY method
double q = 0;
for (size_t i = 1; i < size + 1; ++i) {
q += 1.0 / i;
}
for (size_t i = 0; i < size; ++i) {
cumminInput[i] = q * ni[i] * pvalues[o[i]];
}
} else if (3 == type) { // Hochberg method
for (size_t i = 0; i < size; ++i) {
cumminInput[i] = (i + 1) * pvalues[o[i]];
}
}
auto cumminArray = cummin(cumminInput);
auto pmin = pminx(cumminArray, 1.0);
std::vector<double> result(size);
for (size_t i = 0; i < size; ++i) {
result[i] = pmin[ro[i]];
}
return result;
}
int main() {
using namespace std;
vector<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
};
vector<vector<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
}
};
vector<string> types{ "bh", "by", "bonferroni", "hochberg", "holm", "hommel" };
for (size_t type = 0; type < types.size(); ++type) {
auto q = pAdjust(pvalues, types[type]);
double error = 0.0;
for (size_t i = 0; i < pvalues.size(); ++i) {
error += abs(q[i] - correctAnswers[type][i]);
}
doubleSay(q);
printf("\ntype = %d = '%s' has a cumulative error of %g\n\n\n", type, types[type].c_str(), error);
}
return 0;
}</lang>
{{out}}
<pre>[ 1] 0.6126681081 0.8521710465 0.1987205200 0.1891595417 0.3217789286
[ 5] 0.9301450000 0.4870370000 0.9301450000 0.6049730556 0.6826752564
[10] 0.6482628947 0.7253722500 0.5280972727 0.8769925556 0.4705703448
[15] 0.9241867391 0.6049730556 0.7856107317 0.4887525806 0.1136717045
[20] 0.4991890625 0.8769925556 0.9991834000 0.3217789286 0.9301450000
[25] 0.2304957692 0.5832475000 0.0389954722 0.8521710465 0.1476842609
[30] 0.0168363750 0.0025629017 0.0351608437 0.0625018947 0.0036365888
[35] 0.0025629017 0.0294688286 0.0061660636 0.0389954722 0.0026889914
[40] 0.0004502862 0.0000125223 0.0788155476 0.0314261300 0.0048465270
[45] 0.0025629017 0.0048465270 0.0011017083 0.0725203250 0.0220595769
[50]
type = 0 = 'bh' has a cumulative error of 8.03053e-07
[ 1] 1.0000000000 1.0000000000 0.8940844244 0.8510676197 1.0000000000
[ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 0.5114323399
[20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[25] 1.0000000000 1.0000000000 0.1754486368 1.0000000000 0.6644618149
[30] 0.0757503083 0.0115310209 0.1581958559 0.2812088585 0.0163617595
[35] 0.0115310209 0.1325863108 0.0277423864 0.1754486368 0.0120983246
[40] 0.0020259303 0.0000563403 0.3546073326 0.1413926119 0.0218055202
[45] 0.0115310209 0.0218055202 0.0049568120 0.3262838334 0.0992505663
[50]
type = 1 = 'by' has a cumulative error of 3.64072e-07
[ 1] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[25] 1.0000000000 1.0000000000 0.7019185000 1.0000000000 1.0000000000
[30] 0.2020365000 0.0151667450 0.5625735000 1.0000000000 0.0290927100
[35] 0.0153774100 0.4125636000 0.0678267000 0.6803480000 0.0188229400
[40] 0.0009005725 0.0000125223 1.0000000000 0.4713919500 0.0439557650
[45] 0.0108891550 0.0484652700 0.0033051250 1.0000000000 0.2867745000
[50]
type = 2 = 'bonferroni' has a cumulative error of 6.5e-08
[ 1] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[25] 0.9991834000 0.9991834000 0.4632662100 0.9991834000 0.9991834000
[30] 0.1575884700 0.0138396690 0.3938014500 0.7600230400 0.0250197306
[35] 0.0138396690 0.3052970640 0.0542613600 0.4626366400 0.0165641872
[40] 0.0008825610 0.0000125223 0.9930759000 0.3394022040 0.0369228426
[45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200
[50]
type = 3 = 'hochberg' has a cumulative error of 2.7375e-07
[ 1] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[25] 1.0000000000 1.0000000000 0.4632662100 1.0000000000 1.0000000000
[30] 0.1575884700 0.0139534054 0.3938014500 0.7600230400 0.0250197306
[35] 0.0139534054 0.3052970640 0.0542613600 0.4626366400 0.0165641872
[40] 0.0008825610 0.0000125223 0.9930759000 0.3394022040 0.0369228426
[45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200
[50]
type = 4 = 'holm' has a cumulative error of 2.8095e-07
[ 1] 0.9991834000 0.9991834000 0.9991834000 0.9987623800 0.9991834000
[ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9595180000
[20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[25] 0.9991834000 0.9991834000 0.4351894700 0.9991834000 0.9766522500
[30] 0.1414255500 0.0130434007 0.3530936533 0.6887708800 0.0238560222
[35] 0.0132245726 0.2722919760 0.0542613600 0.4218157600 0.0158112696
[40] 0.0008825610 0.0000125223 0.8743649143 0.3016908480 0.0351646120
[45] 0.0095824564 0.0387722160 0.0031729200 0.8122276400 0.1950066600
[50]
type = 5 = 'hommel' has a cumulative error of 4.35302e-07</pre>
=={{header|D}}==
Line 3,932:
type Hommel has cumulative error of 4.35302e-07.
</pre>
=={{header|Phix}}==
Line 4,669 ⟶ 4,486:
[46] 9.582456e-03 3.877222e-02 3.172920e-03 8.122276e-01 1.950067e-01
Hommel</pre>
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2019.03.1}}
<lang perl6>########################### Helper subs ###########################
sub adjusted (@p, $type) { "\n$type\n" ~ format adjust( check(@p), $type ) }
sub format ( @p, $cols = 5 ) {
my $i = -$cols;
my $fmt = "%1.10f";
join "\n", @p.rotor($cols, :partial).map:
{ sprintf "[%2d] { join ' ', $fmt xx $_ }", $i+=$cols, $_ };
}
sub check ( @p ) { die 'p-values must be in range 0.0 to 1.0' if @p.min < 0 or 1 < @p.max; @p }
multi ratchet ( 'up', @p ) { my $m; @p[$_] min= $m, $m = @p[$_] for ^@p; @p }
multi ratchet ( 'dn', @p ) { my $m; @p[$_] max= $m, $m = @p[$_] for ^@p .reverse; @p }
sub schwartzian ( @p, &transform, :$ratchet ) {
my @pa = @p.map( {[$_, $++]} ).sort( -*.[0] ).map: { [transform(.[0]), .[1]] };
@pa[*;0] = ratchet($ratchet, @pa»[0]);
@pa.sort( *.[1] )»[0]
}
############# The various p-value correction routines #############
multi adjust( @p, 'Benjamini-Hochberg' ) {
@p.&schwartzian: * * @p / (@p - $++) min 1, :ratchet('up')
}
multi adjust( @p, 'Benjamini-Yekutieli' ) {
my \r = ^@p .map( { 1 / ++$ } ).sum;
@p.&schwartzian: * * r * @p / (@p - $++) min 1, :ratchet('up')
}
multi adjust( @p, 'Hochberg' ) {
my \m = @p.max;
@p.&schwartzian: * * ++$ min m, :ratchet('up')
}
multi adjust( @p, 'Holm' ) {
@p.&schwartzian: * * ++$ min 1, :ratchet('dn')
}
multi adjust( @p, 'Šidák' ) {
@p.&schwartzian: 1 - (1 - *) ** ++$, :ratchet('dn')
}
multi adjust( @p, 'Bonferroni' ) {
@p.map: * * @p min 1
}
# Hommel correction can't be easily reduced to a one pass transform
multi adjust( @p, 'Hommel' ) {
my @s = @p.map( {[$_, $++]} ).sort: *.[0] ; # sorted
my \z = +@p; # array si(z)e
my @pa = @s»[0].map( * * z / ++$ ).min xx z; # p adjusted
my @q; # scratch array
for (1 ..^ z).reverse -> $i {
my @L = 0 .. z - $i; # lower indices
my @U = z - $i ^..^ z; # upper indices
my $q = @s[@U]»[0].map( { $_ * $i / (2 + $++) } ).min;
@q[@L] = @s[@L]»[0].map: { min $_ * $i, $q, @s[*-1][0] };
@pa = ^z .map: { max @pa[$_], @q[$_] }
}
@pa[@s[*;1].map( {[$_, $++]} ).sort( *.[0] )»[1]]
}
multi adjust ( @p, $unknown ) {
note "\nSorry, do not know how to do $unknown correction.\n" ~
"Perhaps you want one of these?:\n" ~
<Benjamini-Hochberg Benjamini-Yekutieli Bonferroni Hochberg
Holm Hommel Šidák>.join("\n");
exit
}
########################### The task ###########################
my @p-values =
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
;
for < Benjamini-Hochberg Benjamini-Yekutieli Bonferroni Hochberg Holm Hommel Šidák >
{
say adjusted @p-values, $_
}</lang>
{{out}}
<pre style="height:60ex;overflow:scroll;">Benjamini-Hochberg
[ 0] 0.6126681081 0.8521710465 0.1987205200 0.1891595417 0.3217789286
[ 5] 0.9301450000 0.4870370000 0.9301450000 0.6049730556 0.6826752564
[10] 0.6482628947 0.7253722500 0.5280972727 0.8769925556 0.4705703448
[15] 0.9241867391 0.6049730556 0.7856107317 0.4887525806 0.1136717045
[20] 0.4991890625 0.8769925556 0.9991834000 0.3217789286 0.9301450000
[25] 0.2304957692 0.5832475000 0.0389954722 0.8521710465 0.1476842609
[30] 0.0168363750 0.0025629017 0.0351608438 0.0625018947 0.0036365888
[35] 0.0025629017 0.0294688286 0.0061660636 0.0389954722 0.0026889914
[40] 0.0004502863 0.0000125223 0.0788155476 0.0314261300 0.0048465270
[45] 0.0025629017 0.0048465270 0.0011017083 0.0725203250 0.0220595769
Benjamini-Yekutieli
[ 0] 1.0000000000 1.0000000000 0.8940844244 0.8510676197 1.0000000000
[ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 0.5114323399
[20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[25] 1.0000000000 1.0000000000 0.1754486368 1.0000000000 0.6644618149
[30] 0.0757503083 0.0115310209 0.1581958559 0.2812088585 0.0163617595
[35] 0.0115310209 0.1325863108 0.0277423864 0.1754486368 0.0120983246
[40] 0.0020259303 0.0000563403 0.3546073326 0.1413926119 0.0218055202
[45] 0.0115310209 0.0218055202 0.0049568120 0.3262838334 0.0992505663
Bonferroni
[ 0] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[25] 1.0000000000 1.0000000000 0.7019185000 1.0000000000 1.0000000000
[30] 0.2020365000 0.0151667450 0.5625735000 1.0000000000 0.0290927100
[35] 0.0153774100 0.4125636000 0.0678267000 0.6803480000 0.0188229400
[40] 0.0009005725 0.0000125223 1.0000000000 0.4713919500 0.0439557650
[45] 0.0108891550 0.0484652700 0.0033051250 1.0000000000 0.2867745000
Hochberg
[ 0] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[25] 0.9991834000 0.9991834000 0.4632662100 0.9991834000 0.9991834000
[30] 0.1575884700 0.0138396690 0.3938014500 0.7600230400 0.0250197306
[35] 0.0138396690 0.3052970640 0.0542613600 0.4626366400 0.0165641872
[40] 0.0008825611 0.0000125223 0.9930759000 0.3394022040 0.0369228426
[45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200
Holm
[ 0] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000
[25] 1.0000000000 1.0000000000 0.4632662100 1.0000000000 1.0000000000
[30] 0.1575884700 0.0139534054 0.3938014500 0.7600230400 0.0250197306
[35] 0.0139534054 0.3052970640 0.0542613600 0.4626366400 0.0165641872
[40] 0.0008825611 0.0000125223 0.9930759000 0.3394022040 0.0369228426
[45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200
Hommel
[ 0] 0.9991834000 0.9991834000 0.9991834000 0.9987623800 0.9991834000
[ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9595180000
[20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000
[25] 0.9991834000 0.9991834000 0.4351894700 0.9991834000 0.9766522500
[30] 0.1414255500 0.0130434007 0.3530936533 0.6887708800 0.0238560222
[35] 0.0132245726 0.2722919760 0.0542613600 0.4218157600 0.0158112696
[40] 0.0008825611 0.0000125223 0.8743649143 0.3016908480 0.0351646120
[45] 0.0095824564 0.0387722160 0.0031729200 0.8122276400 0.1950066600
Šidák
[ 0] 0.9998642526 0.9999922727 0.9341844137 0.9234670175 0.9899922294
[ 5] 0.9999922727 0.9992955735 0.9999922727 0.9998642526 0.9998909746
[10] 0.9998642526 0.9999288207 0.9995533892 0.9999922727 0.9990991210
[15] 0.9999922727 0.9998642526 0.9999674876 0.9992955735 0.7741716825
[20] 0.9993332472 0.9999922727 0.9999922727 0.9899922294 0.9999922727
[25] 0.9589019598 0.9998137104 0.3728369461 0.9999922727 0.8605248833
[30] 0.1460714182 0.0138585952 0.3270159382 0.5366136349 0.0247164330
[35] 0.0138585952 0.2640282766 0.0528503728 0.3723753774 0.0164308228
[40] 0.0008821796 0.0000125222 0.6357389664 0.2889497995 0.0362651575
[45] 0.0101847015 0.0389807074 0.0031679962 0.5985019850 0.1963376344</pre>
=={{header|Ruby}}==
Line 4,951 ⟶ 4,952:
total error for Hommel = 1.1483094955369324e-07
</pre>
=={{header|SAS}}==
|