P-value correction: Difference between revisions

m
(removed draft stauts)
m (→‎{{header|Wren}}: Minor tidy)
 
(16 intermediate revisions by 9 users not shown)
Line 1:
{{task|Probability and statistics}}
Given a list of [[wp:p-value|p-values]], adjust the p-values for multiple comparisons. This is done in order to control the false positive, or Type 1 error rate. This is also known as the "[[wp:False discovery rate|false discovery rate]]" (FDR). After adjustment, the p-values will be higher but still inside [0,1]. The adjusted p-values are sometimes called "q-values".
 
Given a list of [[wp:p-value|p-values]], adjust the p-values for multiple comparisons. This is done in order to control the false positive, or Type 1 error rate.
 
This is also known as the "[[wp:False discovery rate|false discovery rate]]" (FDR). After adjustment, the p-values will be higher but still inside [0,1].
 
The adjusted p-values are sometimes called "q-values".
 
 
;Task:
Line 14 ⟶ 21:
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}
 
 
There are several methods to do this, see:
Line 21 ⟶ 29:
* Yosef Hochberg, "A sharper Bonferroni procedure for multiple tests of significance", ''Biometrika'', Vol. 75, No. 4 (1988), pp 800–802, DOI:[https://doi.org/10.1093/biomet/75.4.800 10.1093/biomet/75.4.800] JSTOR:[https://www.jstor.org/stable/2336325 2336325]
* Gerhard Hommel, "A stagewise rejective multiple test procedure based on a modified Bonferroni test", ''Biometrika'', Vol. 75, No. 2 (1988), pp 383–386, DOI:[https://doi.org/10.1093/biomet/75.2.383 10.1093/biomet/75.2.383] JSTOR:[https://www.jstor.org/stable/2336190 2336190]
 
 
Each method has its own advantages and disadvantages.
<br><br>
 
=={{header|C}}==
Line 43 ⟶ 53:
 
Link with <code>-lm</code>
<langsyntaxhighlight Clang="c">#include <stdio.h>//printf
#include <stdlib.h>//qsort
#include <math.h>//fabs
#include <stdbool.h>//bool data type
#include <strings.h>//strcasecmp
#include <assert.h>//assert, necessary for random integer selection
 
unsigned int *restrict seq_len(const size_tunsigned int START, const size_tunsigned int END) {
//named after R function of same name, but simpler function
size_tunsigned start = (unsigned)START;
size_tunsigned end = (unsigned)END;
if (START == END) {
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int));
Line 60 ⟶ 71:
exit(EXIT_FAILURE);
}
for (size_tunsigned i = 0; i < end; i++) {
sequence[i] = i+1;
}
Line 66 ⟶ 77:
}
if (START > END) {
end = (unsigned)START;
start = (unsigned)END;
}
const size_tunsigned LENGTH = end - start ;
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int));
if (sequence == NULL) {
Line 77 ⟶ 88:
}
if (START < END) {
for (size_tunsigned index = 0; index <= LENGTH; index++) {
sequence[index] = start + index;
}
} else {
for (size_tunsigned index = 0; index <= LENGTH; index++) {
sequence[index] = end - index;
}
Line 114 ⟶ 125:
}
 
unsigned int *restrict order (const double *restrict ARRAY, const unsigned int SIZE, const bool DECREASING) {
//this has the same name as the same R function
unsigned int *restrict idx = malloc(SIZE * sizeof(unsigned int));
Line 141 ⟶ 152:
}
 
double *restrict cummin(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
//this takes the same name of the R function which it copies
//this requires a free() afterward where it is used
Line 165 ⟶ 176:
}
 
double *restrict cummax(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) {
//this takes the same name of the R function which it copies
//this requires a free() afterward where it is used
Line 180 ⟶ 191:
}
double cumulative_max = ARRAY[0];
for (size_tunsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
if (ARRAY[i] > cumulative_max) {
cumulative_max = ARRAY[i];
Line 189 ⟶ 200:
}
 
double *restrict pminx(const double *restrict ARRAY, const size_tunsigned int NO_OF_ARRAY_ELEMENTS, const double X) {
//named after the R function pmin
if (NO_OF_ARRAY_ELEMENTS < 1) {
Line 214 ⟶ 225:
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
printf("[1] %e", ARRAY[0]);
for (size_tunsigned int i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
printf(" %.10f", ARRAY[i]);
if (((i+1) % 5) == 0) {
printf("\n[%zuu]", i+1);
}
}
Line 232 ⟶ 243:
}*/
 
double *restrict uint2double (const unsigned int *restrict ARRAY, const size_tunsigned int NO_OF_ARRAY_ELEMENTS) {
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (doubleArray == NULL) {
Line 239 ⟶ 250:
exit(EXIT_FAILURE);
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
doubleArray[index] = (double)ARRAY[index];
}
Line 253 ⟶ 264:
}
 
double *restrict p_adjust (const double *restrict PVALUES, const size_tunsigned int NO_OF_ARRAY_ELEMENTS, const char *restrict STRING) {
//this function is a translation of R's p.adjust "BH" method
// i is always i[index] = NO_OF_ARRAY_ELEMENTS - index - 1
Line 262 ⟶ 273:
}
short int TYPE = -1;
if (strcasecmp(STRING, "BH") == 0NULL) {
TYPE = 0;
} else if (strcasecmp(STRING, "BH") == 0) {
TYPE = 0;
} else if (strcasecmp(STRING, "fdr") == 0) {
Line 290 ⟶ 303:
exit(EXIT_FAILURE);
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS;
if (BONFERRONI >= 1.0) {
Line 312 ⟶ 325:
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (size_tunsigned index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * (double)PVALUES[o[index]];
// printf("cummax_input[%zu] = %e\n", index, cummax_input[index]);
}
Line 326 ⟶ 339:
free(cummax_output); cummax_output = NULL;
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
qvalues[index] = pmin[ro[index]];
}
Line 345 ⟶ 358:
exit(EXIT_FAILURE);
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
p[index] = PVALUES[o[index]];
}
Line 368 ⟶ 381:
}
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0];
for (size_tunsigned index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (double)(1+index);
if (TEMP < min) {
min = TEMP;
}
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
pa[index] = min;
q[index] = min;
Line 389 ⟶ 402:
}
*/
for (size_tunsigned j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) {
// printf("j = %zu\n", j);
unsigned int *restrict ij = seq_len(10,NO_OF_ARRAY_ELEMENTS - j + 1);
for (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS - j + 1; i++) {
ij[i]--;//R's indices are 1-based, C's are 0-based
}
const size_t I2_LENGTH = j - 1;
unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int));
for (size_tunsigned i = 0; i < I2_LENGTH; i++) {
i2[i] = NO_OF_ARRAY_ELEMENTS-j+2+i-1;
//R's indices are 1-based, C's are 0-based, I added the -1
}
 
double q1 = (double)j * p[i2[0]] / 2.0;
for (size_tunsigned int i = 1; i < I2_LENGTH; i++) {//loop through 2:j
const double TEMP_Q1 = (double)j * p[i2[i]] / (double)(2 + i);
if (TEMP_Q1 < q1) {
q1 = TEMP_Q1;
Line 410 ⟶ 420:
}
 
for (size_tunsigned int i = 0; i < (NO_OF_ARRAY_ELEMENTS - j + 1); i++) {//q[ij] <- pmin(j * p[ij], q1)
q[ij[i]] = min2( (double)j*p[ij[i]], q1);
}
free(ij); ij = NULL;
 
for (size_tunsigned int i = 0; i < I2_LENGTH; i++) {//q[i2] <- q[n - j + 1]
q[i2[i]] = q[NO_OF_ARRAY_ELEMENTS - j];//subtract 1 because of starting index difference
}
free(i2); i2 = NULL;
 
for (size_tunsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q)
if (pa[i] < q[i]) {
pa[i] = q[i];
Line 429 ⟶ 439:
}//end j loop
free(p); p = NULL;
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
q[index] = pa[ro[index]];//Hommel q-values
}
Line 462 ⟶ 472:
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (TYPE == 0) {//BH method
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (double)(NO_OF_ARRAY_ELEMENTS - index);// n/i simplified
cummin_input[index] = NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o]
}
} else if (TYPE == 1) {//BY method
double q = 1.0;
for (size_tunsigned int index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) {
q += (double) 1.0/(double)index;
}
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (double)(NO_OF_ARRAY_ELEMENTS - index);// n/i simplified
cummin_input[index] = q * NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o]
}
} else if (TYPE == 3) {//Hochberg method
for (size_tunsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
// pmin(1, cummin((n - i + 1L) * p[o]))[ro]
cummin_input[index] = (double)(index + 1) * PVALUES[o[index]];
}
}
Line 496 ⟶ 506:
return q_array;
}
 
 
int main(void) {
Line 590 ⟶ 601:
return 0;
}
</syntaxhighlight>
</lang>
 
{{out}}
Line 681 ⟶ 692:
{{works with|C89}}
{{trans|Kotlin}}
To avoid licensing issues, this version is a translation of the Kotlin entry (Version 2) which is itself a partial translation of the Perl 6Raku entry. If using gcc, you need to link to the math library (-lm).
<langsyntaxhighlight Clang="c">#include <stdio.h>
#include <stdlib.h>
#include <math.h>
Line 876 ⟶ 887:
each_i(0, 8) adjusted(p_values, types[i]);
return 0;
}</langsyntaxhighlight>
 
{{output}}
Line 883 ⟶ 894:
</pre>
 
=={{header|C++ sharp|C#}}==
{{trans|Java}}
<syntaxhighlight lang="csharp">using System;
<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|C#|C sharp}}==
{{trans|Java}}
<lang csharp>using System;
using System.Collections.Generic;
using System.Linq;
Line 1,624 ⟶ 1,227:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 6.126681E-001 8.521710E-001 1.987205E-001 1.891595E-001 3.217789E-001
Line 1,703 ⟶ 1,306:
[50]
type 5 = 'hommel' has a cumulative error of 4.353024E-007</pre>
 
=={{header|C++}}==
{{trans|Java}}
<syntaxhighlight 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;
}</syntaxhighlight>
{{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}}==
{{trans|Kotlin}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight Dlang="d">import std.algorithm;
import std.conv;
import std.math;
Line 2,046 ⟶ 2,057:
writefln("\ntype %d = '%s' has a cumulative error of %g", type, types[type], error);
}
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286
Line 2,129 ⟶ 2,140:
=={{header|Go}}==
{{trans|Kotlin (Version 2)}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 2,432 ⟶ 2,443:
fmt.Println(s)
}
}</langsyntaxhighlight>
 
{{out}}
Line 2,535 ⟶ 2,546:
{{works with|Java|8}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight Javalang="java">import java.util.Arrays;
import java.util.Comparator;
 
Line 2,884 ⟶ 2,895:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286
Line 2,966 ⟶ 2,977:
 
=={{header|Julia}}==
<syntaxhighlight lang="julia">using MultipleTesting, IterTools, Printf
{{works with|Julia|0.6}}
 
<lang julia>using MultipleTesting
using IterTools
 
p = [4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
Line 2,993 ⟶ 3,001:
println("\n", corr)
printpvalues(adjust(p, corr))
end</langsyntaxhighlight>
 
{{out}}
Line 3,050 ⟶ 3,058:
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
 
<langsyntaxhighlight lang="scala">// version 1.1.51
 
import java.util.Arrays
Line 3,331 ⟶ 3,339:
println(f.format(type, types[type], error))
}
}</langsyntaxhighlight>
 
{{out}}
Line 3,416 ⟶ 3,424:
 
===Version 2 ===
{{trans|Perl 6Raku}}
 
To avoid licensing issues, this version follows the approach of the Perl 6Raku entry of which it is a partial translation. However, the correction routines themselves have been coded independently, common code factored out into separate functions (analogous to Perl 6Raku) and (apart from the Šidák method) agree with the Perl 6Raku results.
<langsyntaxhighlight lang="scala">// version 1.2.21
 
typealias DList = List<Double>
Line 3,564 ⟶ 3,572:
 
types.forEach { println(adjusted(pValues, it)) }
}</langsyntaxhighlight>
 
{{out}}
Same as Perl 6Raku entry except:
<pre>
....
Line 3,593 ⟶ 3,601:
Šidák
</pre>
 
=={{header|Nim}}==
{{trans|Kotlin (Version 2)}}
<syntaxhighlight lang="nim">import algorithm, math, sequtils, strformat, strutils, sugar
 
type
 
CorrectionType {.pure.} = enum
BenjaminiHochberg = "Benjamini-Hochberg"
BenjaminiYekutieli = "Benjamini-Yekutieli"
Bonferroni = "Bonferroni"
Hochberg = "Hochberg"
Holm = "Holm"
Hommel = "Hommel"
Šidák = "Šidák"
 
Direction {.pure.} = enum Up, Down
 
PValues = seq[float]
 
 
template newPValues(length: Natural): PValues =
## Create a PValues object of given length.
newSeq[float](length)
 
 
func ratchet(p: var PValues; dir: Direction) =
var m = p[0]
case dir
of Up:
for i in 1..p.high:
if p[i] > m: p[i] = m
m = p[i]
of Down:
for i in 1..p.high:
if p[i] < m: p[i] = m
m = p[i]
for i in 0..p.high:
if p[i] > 1: p[i] = 1
 
 
func schwartzian(p, mult: PValues; dir: Direction): PValues =
 
let length = p.len
let sortOrder = if dir == Up: Descending else: Ascending
let order1 = toSeq(p.pairs).sorted((x, y) => cmp(x.val, y.val), sortOrder).mapIt(it.key)
 
var pa = newPValues(length)
for i in 0..pa.high:
pa[i] = mult[i] * p[order1[i]]
 
ratchet(pa, dir)
 
let order2 = toSeq(order1.pairs).sortedByIt(it.val).mapIt(it.key)
for idx in order2:
result.add pa[idx]
 
 
proc adjust(p: PValues; ctype: CorrectionType): PValues =
let length = p.len
assert length > 0
let flength = length.toFloat
 
case ctype
 
of BenjaminiHochberg:
var mult = newPValues(length)
for i in 0..mult.high:
mult[i] = flength / (flength - i.toFloat)
return schwartzian(p, mult, Up)
 
of BenjaminiYekutieli:
var q = 0.0
for i in 1..length: q += 1 / i
var mult = newPValues(length)
for i in 0..mult.high:
mult[i] = (q * flength) / (flength - i.toFloat)
return schwartzian(p, mult, Up)
 
of Bonferroni:
result = newPValues(length)
for i in 0..result.high:
result[i] = min(p[i] * flength, 1)
return
 
of Hochberg:
var mult = newPValues(length)
for i in 0..mult.high:
mult[i] = i.toFloat + 1
return schwartzian(p, mult, Up)
 
of Holm:
var mult = newPValues(length)
for i in 0..mult.high:
mult[i] = flength - i.toFloat
return schwartzian(p, mult, Down)
 
of Hommel:
let order1 = toSeq(p.pairs).sortedByIt(it.val).mapIt(it.key)
let s = order1.mapIt(p[it])
var m = Inf
for i in 0..s.high:
m = min(m, s[i] * flength / (i + 1).toFloat)
var q, pa = repeat(m, length)
 
for j in countdown(length - 1, 2):
let lower = toSeq(0..length - j)
let upper = toSeq((length - j + 1)..<length)
var qmin = j.toFloat * s[upper[0]] / 2
for i in 1..upper.high:
let val = s[upper[i]] * j.toFloat / (i + 2).toFloat
if val < qmin: qmin = val
for idx in lower: q[idx] = min(s[idx] * j.toFloat, qmin)
for idx in upper: q[idx] = q[^j]
for i, val in q:
if pa[i] < val: pa[i] = val
 
let order2 = toSeq(order1.pairs).sortedByIt(it.val).mapIt(it.key)
return order2.mapIt(pa[it])
 
of Šidák:
result = newPValues(length)
for i in 0..result.high:
result[i] = 1 - (1 - p[i])^length
return
 
 
func pformat(p: PValues; cols = 5): string =
var lines: seq[string]
for i in countup(0, p.high, cols):
let fchunk = p[i..<(i + cols)]
var schunk = newSeq[string](fchunk.len)
for j in 0..<cols:
schunk[j] = fchunk[j].formatFloat(ffDecimal, 10)
lines.add &"[{i:2}] {schunk.join(\" \")}"
result = lines.join("\n")
 
 
func adjusted(p: PValues; ctype: CorrectionType): string =
doAssert p.len > 0 and min(p) >= 0 and max(p) <= 1, "p-values must be in range 0.0 to 1.0."
result = &"\n{ctype}\n{pformat(p.adjust(ctype))}"
 
when isMainModule:
 
const PVals = @[
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 ctype in CorrectionType:
echo adjusted(PVals, ctype)</syntaxhighlight>
 
{{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.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
 
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.0008825610 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.0008825610 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.0008825610 0.0000125223 0.8743649143 0.3016908480 0.0351646120
[45] 0.0095824564 0.0387722160 0.0031729200 0.8122276400 0.1950066600
 
Šidák
[ 0] 1.0000000000 1.0000000000 0.9946598274 0.9914285749 0.9999515274
[ 5] 1.0000000000 0.9999999688 1.0000000000 1.0000000000 1.0000000000
[10] 1.0000000000 1.0000000000 0.9999999995 1.0000000000 0.9999998801
[15] 1.0000000000 1.0000000000 1.0000000000 0.9999999855 0.9231179729
[20] 0.9999999956 1.0000000000 1.0000000000 0.9999317605 1.0000000000
[25] 0.9983109511 1.0000000000 0.5068253940 1.0000000000 0.9703301333
[30] 0.1832692440 0.0150545753 0.4320729669 0.6993672225 0.0286818157
[35] 0.0152621104 0.3391808707 0.0656206307 0.4959194266 0.0186503726
[40] 0.0009001752 0.0000125222 0.8142104886 0.3772612062 0.0430222116
[45] 0.0108312558 0.0473319661 0.0032997780 0.7705015898 0.2499384839</pre>
 
=={{header|Perl}}==
{{trans|C}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight lang="perl">#!/usr/bin/env perl
 
use strict;
use warnings FATAL => 'all';
use autodie ':all';
use List::Util 'min';
use feature 'say';
 
sub pmin {
my $array_refarray = shift;
my $x = 1;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @pmin_array;
my $n = scalar @$array_refarray;
for (my $index = 0; $index < $n; $index++) {
if$pmin_array[$index] = min(@$array_refarray[$index] <, $x) {;
$pmin_array[$index] = @$array_ref[$index];
} else {
$pmin_array[$index] = $x;
}
}
return @pmin_array;
}
 
sub cummin {
my $array_ref = shift;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @cummin;
my $cumulative_min = @$array_ref[0];
Line 3,635 ⟶ 3,878:
push @cummin, $cumulative_min;
}
return @cummin;
}
 
sub cummax {
my $array_ref = shift;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @cummax;
my $cumulative_max = @$array_ref[0];
Line 3,652 ⟶ 3,891:
push @cummax, $cumulative_max;
}
return @cummax;
}
 
Line 3,668 ⟶ 3,907:
die;
}
}
unless ((ref $array_ref) =~ m/ARRAY/) {
print "You should have entered an array.\n";
die;
}
my @array;
Line 3,680 ⟶ 3,915:
@array = sort { @$array_ref[$b] <=> @$array_ref[$a] } 0..$max_index;
}
@array
}
 
use List::Util 'min';
 
sub p_adjust {
my $pvalues_ref = shift;
unless ((ref $pvalues_ref) =~ m/ARRAY/) {
print "p_adjust requires an array.\n";
die;
}
my $method;
if (defined $_[0]) {
$method = shift;
} else {
$method = 'Holm';
}
my %methods = (
Line 3,710 ⟶ 3,941:
$method = $key;
$method_found = 'yes';
last;
}
}
Line 3,724 ⟶ 3,955:
if ($method_found eq 'no') {
print "No method could be determined from $method.\n";
die;
}
my $lp = scalar @$pvalues_ref;
Line 3,736 ⟶ 3,967:
}
my @cummin = cummin(\@cummin_input);
undef @cummin_input;
my @pmin = pmin(\@cummin);
undef @cummin;
my @ro = order(\@o);
undef @o;
@qvalues = @pmin[@ro];
} elsif ($method eq 'bh') {
Line 3,749 ⟶ 3,977:
}
my @ro = order(\@o);
undef @o;
my @cummin = cummin(\@cummin_input);
undef @cummin_input;
my @pmin = pmin(\@cummin);
undef @cummin;
@qvalues = @pmin[@ro];
} elsif ($method eq 'by') {
Line 3,766 ⟶ 3,991:
$cummin_input[$index] = $q * ($n/($n-$index)) * @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o]
}
# say join (',', @cummin_input);
# say '@cummin_input # of elements = ' . scalar @cummin_input;
my @cummin = cummin(\@cummin_input);
undef @cummin_input;
my @pmin = pmin(\@cummin);
undef @cummin;
@qvalues = @pmin[@ro];
} elsif ($method eq 'bonferroni') {
Line 3,779 ⟶ 4,005:
$qvalues[$index] = 1.0;
} else {
printsay "'Failed to get Bonferroni adjusted p."';
die;
}
Line 3,797 ⟶ 4,023:
@qvalues = @pmin[@ro];
} elsif ($method eq 'hommel') {
my @i = 1..$n;
my @o = order($pvalues_ref);
my @p = @$pvalues_ref[@o];
my @ro = order(\@o);
undef @o;
my (@q, @pa);
my @q;
my $min = $n*$p[0];
for (my $index = 0; $index < $n; $index++) {
my $temp = $n*$p[$index] / ($index + 1);
if ($tempmin <= min($min), {$temp);
$min = $temp;
}
}
for (my $index = 0; $index < $n; $index++) {
Line 3,816 ⟶ 4,038:
}
for (my $j = ($n-1); $j >= 2; $j--) {
my @ij = 0..($n - $j);#ij <- seq_len(n - j + 1)
# printf("j = %zu\n", j);
my @ij = 1..($n - $j + 1);#ij <- seq_len(n - j + 1)
for (my $i = 0; $i < $n - $j + 1; $i++) {
$ij[$i]--;#R's indices are 1-based, C's are 0-based
}
my $I2_LENGTH = $j - 1;
my @i2;
Line 3,831 ⟶ 4,049:
for (my $i = 1; $i < $I2_LENGTH; $i++) {#loop through 2:j
my $TEMP_Q1 = $j * $p[$i2[$i]] / (2 + $i);
if$q1 = min($TEMP_Q1 <, $q1) {;
$q1 = $TEMP_Q1;
}
}
 
for (my $i = 0; $i < ($n - $j + 1); $i++) {#q[ij] <- pmin(j * p[ij], q1)
$q[$ij[$i]] = min( $j*$p[$ij[$i]], $q1);
Line 3,841 ⟶ 4,056:
 
for (my $i = 0; $i < $I2_LENGTH; $i++) {#q[i2] <- q[n - j + 1]
$q[$i2[$i]] = $q[$n - $j];#subtract 1 because of starting index difference
}
 
Line 3,855 ⟶ 4,070:
} else {
print "$method doesn't fit my types.\n";
die;
}
return @qvalues;
}
my @pvalues = (4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
Line 3,944 ⟶ 4,159:
printf("type $method has cumulative error of %g.\n", $error);
}
</syntaxhighlight>
</lang>
 
{{out}}
Line 3,961 ⟶ 4,176:
type Hommel has cumulative error of 4.35302e-07.
</pre>
 
=={{header|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|Phix}}==
Translation of Kotlin (version 2), except for the Hommel part, which is translated from Go.<br>
<!--<syntaxhighlight lang="phix">(phixonline)-->
Note that sq_min(), extract(), and custom_sort() as used below require 0.8.0+
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<lang Phix>enum UP, DOWN
<span style="color: #008080;">enum</span> <span style="color: #000000;">UP</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">DOWN</span>
 
function ratchet(sequence p, integer direction)
<span style="color: #008080;">function</span> <span style="color: #000000;">ratchet</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">direction</span><span style="color: #0000FF;">)</span>
atom m = p[1]
<span style="color: #004080;">atom</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
for i=1 to length(p) do
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
if iff(direction=UP?p[i]>m:p[i]<m) then p[i] = m end if
<span style="color: #008080;">if</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">direction</span><span style="color: #0000FF;">=</span><span style="color: #000000;">UP</span><span style="color: #0000FF;">?</span><span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]></span><span style="color: #000000;">m</span><span style="color: #0000FF;">:</span><span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]<</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
m = p[i]
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return sq_min(p,1)
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sq_min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
function schwartzian(sequence p, mult, integer direction)
<span style="color: #008080;">function</span> <span style="color: #000000;">schwartzian</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">direction</span><span style="color: #0000FF;">)</span>
sequence order = custom_sort(p,tagset(length(p)))
<span style="color: #004080;">sequence</span> <span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">custom_sort</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)))</span>
if direction=UP then order = reverse(order) end if
<span style="color: #008080;">if</span> <span style="color: #000000;">direction</span><span style="color: #0000FF;">=</span><span style="color: #000000;">UP</span> <span style="color: #008080;">then</span> <span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">reverse</span><span style="color: #0000FF;">(</span><span style="color: #000000;">order</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
sequence pa = ratchet(sq_mul(mult,extract(p,order)), direction)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pa</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ratchet</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mult</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">order</span><span style="color: #0000FF;">)),</span> <span style="color: #000000;">direction</span><span style="color: #0000FF;">)</span>
return extract(pa,order,invert:=true)
<span style="color: #008080;">return</span> <span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">order</span><span style="color: #0000FF;">,</span><span style="color: #000000;">invert</span><span style="color: #0000FF;">:=</span><span style="color: #004600;">true</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
function adjust(sequence p, string method)
<span style="color: #008080;">function</span> <span style="color: #000000;">adjust</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">string</span> <span style="color: #000000;">method</span><span style="color: #0000FF;">)</span>
integer size = length(p)
<span style="color: #004080;">integer</span> <span style="color: #000000;">size</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
sequence mult = tagset(size)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">mult</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
switch method
<span style="color: #008080;">switch</span> <span style="color: #000000;">method</span>
 
case "Benjamini-Hochberg":
<span style="color: #008080;">case</span> <span style="color: #008000;">"Benjamini-Hochberg"</span><span style="color: #0000FF;">:</span>
mult = sq_div(size,sq_sub(size+1,mult))
<span style="color: #000000;">mult</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mult</span><span style="color: #0000FF;">))</span>
return schwartzian(p, mult, UP)
<span style="color: #008080;">return</span> <span style="color: #000000;">schwartzian</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">UP</span><span style="color: #0000FF;">)</span>
 
case "Benjamini-Yekutieli":
<span style="color: #008080;">case</span> <span style="color: #008000;">"Benjamini-Yekutieli"</span><span style="color: #0000FF;">:</span>
atom q = sum(sq_div(1,mult))
<span style="color: #004080;">atom</span> <span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mult</span><span style="color: #0000FF;">))</span>
mult = sq_div(q*size,sq_sub(size+1,mult))
<span style="color: #000000;">mult</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">*</span><span style="color: #000000;">size</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mult</span><span style="color: #0000FF;">))</span>
return schwartzian(p, mult, UP)
<span style="color: #008080;">return</span> <span style="color: #000000;">schwartzian</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">UP</span><span style="color: #0000FF;">)</span>
 
case "Bonferroni":
<span style="color: #008080;">case</span> <span style="color: #008000;">"Bonferroni"</span><span style="color: #0000FF;">:</span>
return sq_min(sq_mul(p,size),1)
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sq_min</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">size</span><span style="color: #0000FF;">),</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
 
case "Hochberg":
<span style="color: #008080;">case</span> <span style="color: #008000;">"Hochberg"</span><span style="color: #0000FF;">:</span>
return schwartzian(p, mult, UP)
<span style="color: #008080;">return</span> <span style="color: #000000;">schwartzian</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">UP</span><span style="color: #0000FF;">)</span>
case "Holm":
<span style="color: #008080;">case</span> <span style="color: #008000;">"Holm"</span><span style="color: #0000FF;">:</span>
mult = sq_sub(size+1,mult)
<span style="color: #000000;">mult</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mult</span><span style="color: #0000FF;">)</span>
return schwartzian(p, mult, DOWN)
<span style="color: #008080;">return</span> <span style="color: #000000;">schwartzian</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">DOWN</span><span style="color: #0000FF;">)</span>
case "Hommel":
<span style="color: #008080;">case</span> <span style="color: #008000;">"Hommel"</span><span style="color: #0000FF;">:</span>
sequence ivdx = repeat(0,size)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">ivdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
for i=1 to size do ivdx[i] = {p[i],i} end for
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">size</span> <span style="color: #008080;">do</span> <span style="color: #000000;">ivdx</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">i</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
ivdx = sort(ivdx)
<span style="color: #000000;">ivdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sort</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ivdx</span><span style="color: #0000FF;">)</span>
sequence s = vslice(ivdx,1),
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">vslice</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ivdx</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span>
m = sq_div(sq_mul(s,size),mult),
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #000000;">size</span><span style="color: #0000FF;">),</span><span style="color: #000000;">mult</span><span style="color: #0000FF;">),</span>
{q,pa} @= repeat(min(m),size),
<span style="color: #000000;">qh</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">),</span><span style="color: #000000;">size</span><span style="color: #0000FF;">),</span>
order = vslice(ivdx,2)
<span style="color: #000000;">pa</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">),</span><span style="color: #000000;">size</span><span style="color: #0000FF;">),</span>
for j=size-1 to 2 by -1 do
<span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">vslice</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ivdx</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
sequence lwr = tagset(size-j+1),
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">size</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">2</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
upr = sq_add(size-j+1,tagset(j-1))
<span style="color: #004080;">sequence</span> <span style="color: #000000;">lwr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">-</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span>
atom qmin = j*s[upr[1]]/2
<span style="color: #000000;">upr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">-</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span>
for i=2 to length(upr) do
<span style="color: #004080;">atom</span> <span style="color: #000000;">qmin</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">*</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">upr</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]]/</span><span style="color: #000000;">2</span>
qmin = min(s[upr[i]]*j/(i+1),qmin)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">upr</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">qmin</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">upr</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]]*</span><span style="color: #000000;">j</span><span style="color: #0000FF;">/(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span><span style="color: #000000;">qmin</span><span style="color: #0000FF;">)</span>
for i=1 to length(lwr) do
<span style="color: #008080;">end</span> q[lwr[i]]<span style="color: min(s[lwr[i]]*j, qmin)#008080;">for</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lwr</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">qh</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lwr</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lwr</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]]*</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">qmin</span><span style="color: #0000FF;">)</span>
for i=1 to length(upr) do
<span style="color: #008080;">end</span> q[upr[i]]<span style="color: q[size-j+1]#008080;">for</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">upr</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">qh</span><span style="color: #0000FF;">[</span><span style="color: #000000;">upr</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">qh</span><span style="color: #0000FF;">[</span><span style="color: #000000;">size</span><span style="color: #0000FF;">-</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
pa = sq_max(pa,q)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end for
<span style="color: #000000;">pa</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">qh</span><span style="color: #0000FF;">)</span>
return extract(pa,order,invert:=true)
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 
<span style="color: #008080;">return</span> <span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">order</span><span style="color: #0000FF;">,</span><span style="color: #000000;">invert</span><span style="color: #0000FF;">:=</span><span style="color: #004600;">true</span><span style="color: #0000FF;">)</span>
case "Sidak":
for i=1 to length(p) do
<span style="color: #008080;">case</span> <span style="color: #008000;">"Sidak"</span><span style="color: #0000FF;">:</span>
p[i] = 1 - power(1-p[i],size)
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
end for
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
return p
<span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> <span style="color: #0000FF;">-</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
else
<span style="color: #008080;">return</span> <span style="color: #000000;">p</span>
return {} -- (unknown method)
 
<span style="color: #008080;">else</span>
end switch
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{}</span> <span style="color: #000080;font-style:italic;">-- (unknown method)</span>
return p
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">switch</span>
 
<span style="color: #008080;">return</span> <span style="color: #000000;">p</span>
constant {types,correct_answers} = columnize({
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
{"Benjamini-Hochberg",
{6.126681e-01, 8.521710e-01, 1.987205e-01, 1.891595e-01, 3.217789e-01,
<span style="color: #008080;">constant</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">types</span><span style="color: #0000FF;">,</span><span style="color: #000000;">correct_answers</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">columnize</span><span style="color: #0000FF;">({</span>
9.301450e-01, 4.870370e-01, 9.301450e-01, 6.049731e-01, 6.826753e-01,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Benjamini-Hochberg"</span><span style="color: #0000FF;">,</span>
6.482629e-01, 7.253722e-01, 5.280973e-01, 8.769926e-01, 4.705703e-01,
<span style="color: #0000FF;">{</span><span style="color: #000000;">6.126681e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.521710e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.987205e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.891595e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.217789e-01</span><span style="color: #0000FF;">,</span>
9.241867e-01, 6.049731e-01, 7.856107e-01, 4.887526e-01, 1.136717e-01,
<span style="color: #000000;">9.301450e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.870370e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.301450e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.049731e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.826753e-01</span><span style="color: #0000FF;">,</span>
4.991891e-01, 8.769926e-01, 9.991834e-01, 3.217789e-01, 9.301450e-01,
<span style="color: #000000;">6.482629e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.253722e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.280973e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.769926e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.705703e-01</span><span style="color: #0000FF;">,</span>
2.304958e-01, 5.832475e-01, 3.899547e-02, 8.521710e-01, 1.476843e-01,
<span style="color: #000000;">9.241867e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.049731e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.856107e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.887526e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.136717e-01</span><span style="color: #0000FF;">,</span>
1.683638e-02, 2.562902e-03, 3.516084e-02, 6.250189e-02, 3.636589e-03,
<span style="color: #000000;">4.991891e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.769926e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.217789e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.301450e-01</span><span style="color: #0000FF;">,</span>
2.562902e-03, 2.946883e-02, 6.166064e-03, 3.899547e-02, 2.688991e-03,
<span style="color: #000000;">2.304958e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.832475e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.899547e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.521710e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.476843e-01</span><span style="color: #0000FF;">,</span>
4.502862e-04, 1.252228e-05, 7.881555e-02, 3.142613e-02, 4.846527e-03,
<span style="color: #000000;">1.683638e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.562902e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.516084e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.250189e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.636589e-03</span><span style="color: #0000FF;">,</span>
2.562902e-03, 4.846527e-03, 1.101708e-03, 7.252032e-02, 2.205958e-02}},
<span style="color: #000000;">2.562902e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.946883e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.166064e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.899547e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.688991e-03</span><span style="color: #0000FF;">,</span>
{"Benjamini-Yekutieli",
<span style="color: #000000;">4.502862e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.252228e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.881555e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.142613e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.846527e-03</span><span style="color: #0000FF;">,</span>
{1.000000e+00, 1.000000e+00, 8.940844e-01, 8.510676e-01, 1.000000e+00,
<span style="color: #000000;">2.562902e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.846527e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.101708e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.252032e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.205958e-02</span><span style="color: #0000FF;">}},</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Benjamini-Yekutieli"</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #0000FF;">{</span><span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.940844e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.510676e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 5.114323e-01,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.754486e-01, 1.000000e+00, 6.644618e-01,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.114323e-01</span><span style="color: #0000FF;">,</span>
7.575031e-02, 1.153102e-02, 1.581959e-01, 2.812089e-01, 1.636176e-02,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.153102e-02, 1.325863e-01, 2.774239e-02, 1.754486e-01, 1.209832e-02,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.754486e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.644618e-01</span><span style="color: #0000FF;">,</span>
2.025930e-03, 5.634031e-05, 3.546073e-01, 1.413926e-01, 2.180552e-02,
<span style="color: #000000;">7.575031e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.153102e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.581959e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.812089e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.636176e-02</span><span style="color: #0000FF;">,</span>
1.153102e-02, 2.180552e-02, 4.956812e-03, 3.262838e-01, 9.925057e-02}},
<span style="color: #000000;">1.153102e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.325863e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.774239e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.754486e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.209832e-02</span><span style="color: #0000FF;">,</span>
{"Bonferroni",
<span style="color: #000000;">2.025930e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.634031e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.546073e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.413926e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.180552e-02</span><span style="color: #0000FF;">,</span>
{1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.153102e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.180552e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.956812e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.262838e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.925057e-02</span><span style="color: #0000FF;">}},</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Bonferroni"</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #0000FF;">{</span><span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 7.019185e-01, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
2.020365e-01, 1.516674e-02, 5.625735e-01, 1.000000e+00, 2.909271e-02,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.537741e-02, 4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.019185e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01, 4.395577e-02,
<span style="color: #000000;">2.020365e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.516674e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.625735e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.909271e-02</span><span style="color: #0000FF;">,</span>
1.088915e-02, 4.846527e-02, 3.305125e-03, 1.000000e+00, 2.867745e-01}},
<span style="color: #000000;">1.537741e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.125636e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.782670e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.803480e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.882294e-02</span><span style="color: #0000FF;">,</span>
{"Hochberg",
<span style="color: #000000;">9.005725e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.252228e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.713920e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.395577e-02</span><span style="color: #0000FF;">,</span>
{9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #000000;">1.088915e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.846527e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.305125e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.867745e-01</span><span style="color: #0000FF;">}},</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Hochberg"</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #0000FF;">{</span><span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 4.632662e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
1.575885e-01, 1.383967e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
1.383967e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.632662e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02,
<span style="color: #000000;">1.575885e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.383967e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.938014e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.600230e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.501973e-02</span><span style="color: #0000FF;">,</span>
1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01}},
<span style="color: #000000;">1.383967e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.052971e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.426136e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.626366e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.656419e-02</span><span style="color: #0000FF;">,</span>
{"Holm",
<span style="color: #000000;">8.825610e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.252228e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.930759e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.394022e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.692284e-02</span><span style="color: #0000FF;">,</span>
{1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.023581e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.974152e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.172920e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.992520e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.179486e-01</span><span style="color: #0000FF;">}},</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Holm"</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #0000FF;">{</span><span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.000000e+00, 1.000000e+00, 4.632662e-01, 1.000000e+00, 1.000000e+00,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.575885e-01, 1.395341e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
1.395341e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02,
<span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.632662e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.000000e+00</span><span style="color: #0000FF;">,</span>
8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02,
<span style="color: #000000;">1.575885e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.395341e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.938014e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.600230e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.501973e-02</span><span style="color: #0000FF;">,</span>
1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01}},
<span style="color: #000000;">1.395341e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.052971e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.426136e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.626366e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.656419e-02</span><span style="color: #0000FF;">,</span>
{"Hommel",
<span style="color: #000000;">8.825610e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.252228e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.930759e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.394022e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.692284e-02</span><span style="color: #0000FF;">,</span>
{9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01, 9.991834e-01,
<span style="color: #000000;">1.023581e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.974152e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.172920e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.992520e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.179486e-01</span><span style="color: #0000FF;">}},</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Hommel"</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #0000FF;">{</span><span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.987624e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
9.991834e-01, 9.991834e-01, 4.351895e-01, 9.991834e-01, 9.766522e-01,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.595180e-01</span><span style="color: #0000FF;">,</span>
1.414256e-01, 1.304340e-02, 3.530937e-01, 6.887709e-01, 2.385602e-02,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span>
1.322457e-02, 2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02,
<span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.351895e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.766522e-01</span><span style="color: #0000FF;">,</span>
8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01, 3.516461e-02,
<span style="color: #000000;">1.414256e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.304340e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.530937e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.887709e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.385602e-02</span><span style="color: #0000FF;">,</span>
9.582456e-03, 3.877222e-02, 3.172920e-03, 8.122276e-01, 1.950067e-01}},
<span style="color: #000000;">1.322457e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.722920e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.426136e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.218158e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.581127e-02</span><span style="color: #0000FF;">,</span>
{"Sidak",
<span style="color: #000000;">8.825610e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.252228e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.743649e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.016908e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.516461e-02</span><span style="color: #0000FF;">,</span>
{1.0000000000, 1.0000000000, 0.9946598274, 0.9914285749, 0.9999515274,
<span style="color: #000000;">9.582456e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.877222e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.172920e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.122276e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.950067e-01</span><span style="color: #0000FF;">}},</span>
1.0000000000, 0.9999999688, 1.0000000000, 1.0000000000, 1.0000000000,
<span style="color: #0000FF;">{</span><span style="color: #008000;">"Sidak"</span><span style="color: #0000FF;">,</span>
1.0000000000, 1.0000000000, 0.9999999995, 1.0000000000, 0.9999998801,
<span style="color: #0000FF;">{</span><span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9946598274</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9914285749</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9999515274</span><span style="color: #0000FF;">,</span>
1.0000000000, 1.0000000000, 1.0000000000, 0.9999999855, 0.9231179729,
<span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9999999688</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span>
0.9999999956, 1.0000000000, 1.0000000000, 0.9999317605, 1.0000000000,
<span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9999999995</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9999998801</span><span style="color: #0000FF;">,</span>
0.9983109511, 1.0000000000, 0.5068253940, 1.0000000000, 0.9703301333,
<span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9999999855</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9231179729</span><span style="color: #0000FF;">,</span>
0.1832692440, 0.0150545753, 0.4320729669, 0.6993672225, 0.0286818157,
<span style="color: #000000;">0.9999999956</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9999317605</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span>
0.0152621104, 0.3391808707, 0.0656206307, 0.4959194266, 0.0186503726,
<span style="color: #000000;">0.9983109511</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.5068253940</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0000000000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.9703301333</span><span style="color: #0000FF;">,</span>
0.0009001752, 0.0000125222, 0.8142104886, 0.3772612062, 0.0430222116,
<span style="color: #000000;">0.1832692440</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0150545753</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.4320729669</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.6993672225</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0286818157</span><span style="color: #0000FF;">,</span>
0.0108312558, 0.0473319661, 0.0032997780, 0.7705015898, 0.2499384839}}})
<span style="color: #000000;">0.0152621104</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.3391808707</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0656206307</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.4959194266</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0186503726</span><span style="color: #0000FF;">,</span>
-- {"Unknown",{1}}})
<span style="color: #000000;">0.0009001752</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0000125222</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.8142104886</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.3772612062</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0430222116</span><span style="color: #0000FF;">,</span>
 
<span style="color: #000000;">0.0108312558</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0473319661</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.0032997780</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.7705015898</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.2499384839</span><span style="color: #0000FF;">}}})</span>
constant pValues = {4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
<span style="color: #000080;font-style:italic;">-- {"Unknown",{1<nowiki>}}</nowiki>})</span>
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,
<span style="color: #008080;">constant</span> <span style="color: #000000;">pValues</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">4.533744e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.296024e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.936026e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.079658e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.801962e-01</span><span style="color: #0000FF;">,</span>
8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02,
<span style="color: #000000;">8.752257e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.922222e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.115421e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.355806e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.324867e-01</span><span style="color: #0000FF;">,</span>
3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01,
<span style="color: #000000;">4.926798e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.802978e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.485442e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.883130e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.729308e-01</span><span style="color: #0000FF;">,</span>
1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02,
<span style="color: #000000;">8.502518e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4.268138e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.442008e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.030266e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.001555e-02</span><span style="color: #0000FF;">,</span>
4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04,
<span style="color: #000000;">3.194810e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.892933e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.991834e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.745691e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.037516e-01</span><span style="color: #0000FF;">,</span>
3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04,
<span style="color: #000000;">1.198578e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.966083e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.403837e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7.328671e-01</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.793476e-02</span><span style="color: #0000FF;">,</span>
1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04,
<span style="color: #000000;">4.040730e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.033349e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.125147e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.375072e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.818542e-04</span><span style="color: #0000FF;">,</span>
2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03}
<span style="color: #000000;">3.075482e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.251272e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.356534e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.360696e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.764588e-04</span><span style="color: #0000FF;">,</span>
 
<span style="color: #000000;">1.801145e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.504456e-07</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3.310253e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.427839e-03</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.791153e-04</span><span style="color: #0000FF;">,</span>
if length(pValues)=0 or min(pValues)<0 or max(pValues)>1 then
<span style="color: #000000;">2.177831e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9.693054e-04</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6.610250e-05</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.900813e-02</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5.735490e-03</span><span style="color: #0000FF;">}</span>
crash("p-values must be in range 0.0 to 1.0")
end if
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pValues</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pValues</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pValues</span><span style="color: #0000FF;">)></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
 
<span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"p-values must be in range 0.0 to 1.0"</span><span style="color: #0000FF;">)</span>
for i=1 to length(types) do
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
string ti = types[i]
sequence res = adjust(pValues,ti)
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">types</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
if res={} then
<span style="color: #004080;">string</span> <span style="color: #000000;">ti</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">types</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
printf(1,"\nSorry, do not know how to do %s correction.\n"&
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">adjust</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pValues</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ti</span><span style="color: #0000FF;">)</span>
"Perhaps you want one of these?:\n %s\n",
<span style="color: #008080;">if</span> <span style="color: #000000;">res</span><span style="color: #0000FF;">={}</span> <span style="color: #008080;">then</span>
{ti,join(types[1..$-1],"\n ")})
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\nSorry, do not know how to do %s correction.\n"</span><span style="color: #0000FF;">&</span>
exit
<span style="color: #008000;">"Perhaps you want one of these?:\n %s\n"</span><span style="color: #0000FF;">,</span>
end if
<span style="color: #0000FF;">{</span><span style="color: #000000;">ti</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">types</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #008000;">"\n "</span><span style="color: #0000FF;">)})</span>
-- printf(1,"%s\n",{ti})
<span style="color: #008080;">exit</span>
-- res = correct_answers[i] -- (for easier comparison only)
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
-- pp(res,{pp_FltFmt,"%13.10f",pp_IntFmt,"%13.10f",pp_Maxlen,75,pp_Pause,0})
<span style="color: #000080;font-style:italic;">-- printf(1,"%s\n",{ti})
atom error = sum(sq_abs(sq_sub(res,correct_answers[i])))
-- res = correct_answers[i] -- (for easier comparison only)
printf(1,"%s has cumulative error of %g\n", {ti,error})
-- pp(res,{pp_FltFmt,"%13.10f",pp_IntFmt,"%13.10f",pp_Maxlen,75,pp_Pause,0})</span>
end for</lang>
<span style="color: #004080;">atom</span> <span style="color: #000000;">error</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_abs</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">correct_answers</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s has cumulative error of %g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">ti</span><span style="color: #0000FF;">,</span><span style="color: #000000;">error</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
{{out}}
Matches Kotlin (etc) when some of those lines just above are uncommented.
Line 4,353 ⟶ 4,389:
{{trans|Perl}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight lang="python">from __future__ import division
import sys
 
Line 4,588 ⟶ 4,624:
error += abs(q[i] - correct_answers[key][i])
print '%s error = %g' % (key.upper(), error)
</syntaxhighlight>
</lang>
 
{{out}}
Line 4,603 ⟶ 4,639:
The '''p.adjust''' function is built-in, see [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html R manual].
 
<langsyntaxhighlight Rlang="rsplus">p <- c(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,
Line 4,631 ⟶ 4,667:
 
p.adjust(p, method = 'hommel')
writeLines("Hommel\n")</langsyntaxhighlight>
 
{{out}}
Line 4,698 ⟶ 4,734:
[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}}
 
<syntaxhighlight lang="raku" line>########################### 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, $_
}</syntaxhighlight>
 
{{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}}==
{{trans|Perl}}
<syntaxhighlight lang="ruby">def pmin(array)
x = 1
pmin_array = []
array.each_index do |i|
pmin_array[i] = [array[i], x].min
abort if pmin_array[i] > 1
end
pmin_array
end
 
def cummin(array)
cumulative_min = array[0]
arr_cummin = []
array.each do |p|
cumulative_min = [p, cumulative_min].min
arr_cummin.push(cumulative_min)
end
arr_cummin
end
 
def cummax(array)
cumulative_max = array[0]
arr_cummax = []
array.each do |p|
cumulative_max = [p, cumulative_max].max
arr_cummax.push(cumulative_max)
end
arr_cummax
end
 
# decreasing variable is optional
def order(array, decreasing = false)
if decreasing == false
array.sort.map { |n| array.index(n) }
else
array.sort.map { |n| array.index(n) }.reverse
end
end
 
def p_adjust(arr_pvalues, method = 'Holm')
lp = arr_pvalues.size
n = lp
if method.casecmp('hochberg').zero?
arr_o = order(arr_pvalues, true)
arr_cummin_input = []
(0..n).each do |index|
arr_cummin_input[index] = (index + 1) * arr_pvalues[arr_o[index].to_i]
end
arr_cummin = cummin(arr_cummin_input)
arr_pmin = pmin(arr_cummin)
arr_ro = order(arr_o)
return arr_pmin.values_at(*arr_ro)
elsif method.casecmp('bh').zero? || method.casecmp('benjamini-hochberg').zero?
arr_o = order(arr_pvalues, true)
arr_cummin_input = []
(0..(n - 1)).each do |i|
arr_cummin_input[i] = (n / (n - i).to_f) * arr_pvalues[arr_o[i]]
end
arr_ro = order(arr_o)
arr_cummin = cummin(arr_cummin_input)
arr_pmin = pmin(arr_cummin)
return arr_pmin.values_at(*arr_ro)
elsif method.casecmp('by').zero? || method.casecmp('benjamini-yekutieli').zero?
q = 0.0
arr_o = order(arr_pvalues, true)
arr_ro = order(arr_o)
(1..n).each do |index|
q += 1.0 / index
end
arr_cummin_input = []
(0..(n - 1)).each do |i|
arr_cummin_input[i] = q * (n / (n - i).to_f) * arr_pvalues[arr_o[i]]
end
arr_cummin = cummin(arr_cummin_input)
arr_pmin = pmin(arr_cummin)
return arr_pmin.values_at(*arr_ro)
elsif method.casecmp('bonferroni').zero?
arr_qvalues = []
(0..(n - 1)).each do |i|
q = arr_pvalues[i] * n
if (q >= 0) && (q < 1)
arr_qvalues[i] = q
elsif q >= 1
arr_qvalues[i] = 1.0
else
puts "Falied to get Bonferroni adjusted p for #{arr_pvalues[i]}"
end
end
return arr_qvalues
elsif method.casecmp('holm').zero?
o = order(arr_pvalues)
cummax_input = []
(0..(n - 1)).each do |index|
cummax_input[index] = (n - index) * arr_pvalues[o[index]]
end
ro = order(o)
arr_cummax = cummax(cummax_input)
arr_pmin = pmin(arr_cummax)
return arr_pmin.values_at(*ro)
elsif method.casecmp('hommel').zero?
o = order(arr_pvalues)
arr_p = arr_pvalues.values_at(*o)
ro = order(o)
q = []
pa = []
min = n * arr_p[0]
(0..(n - 1)).each do |index|
temp = n * arr_p[index] / (index + 1)
min = [min, temp].min
end
(0..(n - 1)).each do |index|
pa[index] = min
q[index] = min
end
j = n - 1
while j >= 2
ij = Array 0..(n - j)
i2_length = j - 1
i2 = []
(0..(i2_length - 1)).each do |i|
i2[i] = n - j + 2 + i - 1 # R's indices are 1-based, C's are 0-based
end
q1 = j * arr_p[i2[0]] / 2.0
(1..(i2_length - 1)).each do |i|
temp_q1 = j * arr_p[i2[i]] / (2 + i)
q1 = [temp_q1, q1].min
end
(0..(n - j)).each do |i|
tmp = j * arr_p[ij[i]]
q[ij[i]] = [tmp, q1].min
end
(0..(i2_length - 1)).each do |i|
q[i2[i]] = q[n - j]
end
(0..(n - 1)).each do |i|
pa[i] = q[i] if pa[i] < q[i]
end
j -= 1
end
return pa.values_at(*ro)
else
puts "#{method} isn't accepted."
abort
end
end
 
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]
 
correct_answers = {
'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]
}
 
# correct_answers.each do |method, answers|
methods = ['Benjamini-Yekutieli', 'Benjamini-Hochberg', 'Hochberg',
'Bonferroni', 'Holm', 'Hommel']
methods.each do |method|
puts method
error = 0.0
arr_q = p_adjust(pvalues, method)
arr_q.each_index do |p|
error += (correct_answers[method][p] - arr_q[p])
end
puts "total error for #{method} = #{error}"
end
</syntaxhighlight>
{{out}}
<pre>Benjamini-Yekutieli
total error for Benjamini-Yekutieli = -1.7373780825929845e-07
Benjamini-Hochberg
total error for Benjamini-Hochberg = -1.4736877299143964e-08
Hochberg
total error for Hochberg = -1.2354999978398105e-07
Bonferroni
total error for Bonferroni = 4.49999999152751e-08
Holm
total error for Holm = -1.163499997815258e-07
Hommel
total error for Hommel = 1.1483094955369324e-07
</pre>
 
=={{header|Rust}}==
<syntaxhighlight lang="rust">
use std::iter;
 
#[rustfmt::skip]
const PVALUES:[f64;50] = [
4.533_744e-01, 7.296_024e-01, 9.936_026e-02, 9.079_658e-02, 1.801_962e-01,
8.752_257e-01, 2.922_222e-01, 9.115_421e-01, 4.355_806e-01, 5.324_867e-01,
4.926_798e-01, 5.802_978e-01, 3.485_442e-01, 7.883_130e-01, 2.729_308e-01,
8.502_518e-01, 4.268_138e-01, 6.442_008e-01, 3.030_266e-01, 5.001_555e-02,
3.194_810e-01, 7.892_933e-01, 9.991_834e-01, 1.745_691e-01, 9.037_516e-01,
1.198_578e-01, 3.966_083e-01, 1.403_837e-02, 7.328_671e-01, 6.793_476e-02,
4.040_730e-03, 3.033_349e-04, 1.125_147e-02, 2.375_072e-02, 5.818_542e-04,
3.075_482e-04, 8.251_272e-03, 1.356_534e-03, 1.360_696e-02, 3.764_588e-04,
1.801_145e-05, 2.504_456e-07, 3.310_253e-02, 9.427_839e-03, 8.791_153e-04,
2.177_831e-04, 9.693_054e-04, 6.610_250e-05, 2.900_813e-02, 5.735_490e-03
];
 
#[derive(Debug)]
enum CorrectionType {
BenjaminiHochberg,
BenjaminiYekutieli,
Bonferroni,
Hochberg,
Holm,
Hommel,
Sidak,
}
 
enum SortDirection {
Increasing,
Decreasing,
}
 
/// orders **input** vector by value and multiplies with **multiplier** vector
/// Finally returns the multiplied values in the original order of **input**
fn ordered_multiply(input: &[f64], multiplier: &[f64], direction: &SortDirection) -> Vec<f64> {
let order_by_value = match direction {
SortDirection::Increasing => {
|a: &(f64, usize), b: &(f64, usize)| b.0.partial_cmp(&a.0).unwrap()
}
SortDirection::Decreasing => {
|a: &(f64, usize), b: &(f64, usize)| a.0.partial_cmp(&b.0).unwrap()
}
};
 
let cmp_minmax = match direction {
SortDirection::Increasing => |a: f64, b: f64| a.gt(&b),
SortDirection::Decreasing => |a: f64, b: f64| a.lt(&b),
};
 
// add original order index
let mut input_indexed = input
.iter()
.enumerate()
.map(|(idx, &p_value)| (p_value, idx))
.collect::<Vec<_>>();
 
// order by value desc/asc
input_indexed.sort_unstable_by(order_by_value);
 
// do the multiplication in place, clamp it at 1.0,
// keep the original index in place
for i in 0..input_indexed.len() {
input_indexed[i] = (
f64::min(1.0, input_indexed[i].0 * multiplier[i]),
input_indexed[i].1,
);
}
 
// make vector strictly monotonous increasing/decreasing in place
for i in 1..input_indexed.len() {
if cmp_minmax(input_indexed[i].0, input_indexed[i - 1].0) {
input_indexed[i] = (input_indexed[i - 1].0, input_indexed[i].1);
}
}
 
// re-sort back to original order
input_indexed.sort_unstable_by(|a: &(f64, usize), b: &(f64, usize)| a.1.cmp(&b.1));
 
// remove ordering index
let (resorted, _): (Vec<_>, Vec<_>) = input_indexed.iter().cloned().unzip();
resorted
}
 
#[allow(clippy::cast_precision_loss)]
fn hommel(input: &[f64]) -> Vec<f64> {
// using algorith described:
// http://stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/MultipleComparision/Writght92.pdf
 
// add original order index
let mut input_indexed = input
.iter()
.enumerate()
.map(|(idx, &p_value)| (p_value, idx))
.collect::<Vec<_>>();
 
// order by value asc
input_indexed
.sort_unstable_by(|a: &(f64, usize), b: &(f64, usize)| a.0.partial_cmp(&b.0).unwrap());
 
let (p_values, order): (Vec<_>, Vec<_>) = input_indexed.iter().cloned().unzip();
 
let n = input.len();
 
// initial minimal n*p/i values
// get the smalles of these values
let min_result = (0..n)
.map(|i| ((p_values[i] * n as f64) / (i + 1) as f64))
.fold(1. / 0. /* -inf */, f64::min);
 
// // initialize result vector with minimal values
let mut result = iter::repeat(min_result).take(n).collect::<Vec<_>>();
 
for m in (2..n).rev() {
let cmin: f64;
let m_as_float = m as f64;
let mut a = p_values.clone();
// println!("\nn: {}", m);
{
// split p-values into two group
let (_, second) = p_values.split_at(n - m + 1);
 
// calculate minumum of m*p/i for this second group
cmin = second
.iter()
.zip(2..=m)
.map(|(p, i)| (m_as_float * p) / i as f64)
.fold(1. / 0. /* inf */, f64::min);
}
 
// replace p values if p<cmin in the second group
((n - m + 1)..n).for_each(|i| a[i] = a[i].max(cmin));
 
// replace p values if min(cmin, m*p) > p
(0..=(n - m)).for_each(|i| a[i] = a[i].max(f64::min(cmin, m_as_float * p_values[i])));
 
// store in the result vector if any adjusted p is higher than the current one
(0..n).for_each(|i| result[i] = result[i].max(a[i]));
}
 
// re-sort into the original order
let mut result = result
.into_iter()
.zip(order.into_iter())
.map(|(p, idx)| (p, idx))
.collect::<Vec<_>>();
result.sort_unstable_by(|a: &(f64, usize), b: &(f64, usize)| a.1.cmp(&b.1));
let (result, _): (Vec<_>, Vec<_>) = result.iter().cloned().unzip();
result
}
#[allow(clippy::cast_precision_loss)]
fn p_value_correction(p_values: &[f64], ctype: &CorrectionType) -> Vec<f64> {
let p_vec = p_values.to_vec();
if p_values.is_empty() {
return p_vec;
}
 
let fsize = p_values.len() as f64;
 
match ctype {
CorrectionType::BenjaminiHochberg => {
let multiplier = (0..p_values.len())
.map(|index| fsize / (fsize - index as f64))
.collect::<Vec<_>>();
 
ordered_multiply(&p_vec, &multiplier, &SortDirection::Increasing)
}
CorrectionType::BenjaminiYekutieli => {
let q: f64 = (1..=p_values.len()).map(|index| 1. / index as f64).sum();
let multiplier = (0..p_values.len())
.map(|index| q * fsize / (fsize - index as f64))
.collect::<Vec<_>>();
 
ordered_multiply(&p_vec, &multiplier, &SortDirection::Increasing)
}
CorrectionType::Bonferroni => p_vec
.iter()
.map(|p| f64::min(p * fsize, 1.0))
.collect::<Vec<_>>(),
CorrectionType::Hochberg => {
let multiplier = (0..p_values.len())
.map(|index| 1. + index as f64)
.collect::<Vec<_>>();
ordered_multiply(&p_vec, &multiplier, &SortDirection::Increasing)
}
CorrectionType::Holm => {
let multiplier = (0..p_values.len())
.map(|index| fsize - index as f64)
.collect::<Vec<_>>();
 
ordered_multiply(&p_vec, &multiplier, &SortDirection::Decreasing)
}
CorrectionType::Sidak => p_vec
.iter()
.map(|x| 1. - (1. - x).powf(fsize))
.collect::<Vec<_>>(),
CorrectionType::Hommel => hommel(&p_vec),
}
}
 
// prints array into a nice table, max 5 floats/row
fn array_to_string(a: &[f64]) -> String {
a.chunks(5)
.enumerate()
.map(|(index, e)| {
format!(
"[{:>2}]: {}",
index * 5,
e.iter()
.map(|x| format!("{:>1.10}", x))
.collect::<Vec<_>>()
.join(", ")
)
})
.collect::<Vec<_>>()
.join("\n")
}
fn main() {
let ctypes = [
CorrectionType::BenjaminiHochberg,
CorrectionType::BenjaminiYekutieli,
CorrectionType::Bonferroni,
CorrectionType::Hochberg,
CorrectionType::Holm,
CorrectionType::Sidak,
CorrectionType::Hommel,
];
 
for ctype in &ctypes {
println!("\n{:?}:", ctype);
println!("{}", array_to_string(&p_value_correction(&PVALUES, ctype)));
}
}
</syntaxhighlight>
{{out}}
<pre style="height:60ex;overflow:scroll;">
BenjaminiHochberg:
[ 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.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
 
BenjaminiYekutieli:
[ 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.0008825610, 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.0008825610, 0.0000125223, 0.9930759000, 0.3394022040, 0.0369228426
[45]: 0.0102358057, 0.0397415214, 0.0031729200, 0.8992520300, 0.2179486200
 
Sidak:
[ 0]: 1.0000000000, 1.0000000000, 0.9946598274, 0.9914285749, 0.9999515274
[ 5]: 1.0000000000, 0.9999999688, 1.0000000000, 1.0000000000, 1.0000000000
[10]: 1.0000000000, 1.0000000000, 0.9999999995, 1.0000000000, 0.9999998801
[15]: 1.0000000000, 1.0000000000, 1.0000000000, 0.9999999855, 0.9231179729
[20]: 0.9999999956, 1.0000000000, 1.0000000000, 0.9999317605, 1.0000000000
[25]: 0.9983109511, 1.0000000000, 0.5068253940, 1.0000000000, 0.9703301333
[30]: 0.1832692440, 0.0150545753, 0.4320729669, 0.6993672225, 0.0286818157
[35]: 0.0152621104, 0.3391808707, 0.0656206307, 0.4959194266, 0.0186503726
[40]: 0.0009001752, 0.0000125222, 0.8142104886, 0.3772612062, 0.0430222116
[45]: 0.0108312558, 0.0473319661, 0.0032997780, 0.7705015898, 0.2499384839
 
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.0008825610, 0.0000125223, 0.8743649143, 0.3016908480, 0.0351646120
[45]: 0.0095824564, 0.0387722160, 0.0031729200, 0.8122276400, 0.1950066600
</pre>
 
=={{header|SAS}}==
 
<langsyntaxhighlight lang="sas">data pvalues;
input raw_p @@;
cards;
Line 4,718 ⟶ 5,542:
 
proc multtest pdata=pvalues bon sid hom hoc holm;
run;</langsyntaxhighlight>
 
'''output'''
Line 4,795 ⟶ 5,619:
First, install the package with:
 
<syntaxhighlight lang ="stata">ssc install qqvalue</langsyntaxhighlight>
 
Given a dataset containing the p-values in a variable, the qqvalue command generates another variable with the adjusted p-values. Here is an example showing the result with all implemented methods:
 
<langsyntaxhighlight lang="stata">clear
 
#delimit ;
Line 4,821 ⟶ 5,645:
}
 
list</langsyntaxhighlight>
 
'''output'''
Line 4,888 ⟶ 5,712:
50. | .00573549 .2867745 .24993848 .21794862 .19633763 .21794862 .02205958 .09925057 |
+-----------------------------------------------------------------------------------------------+</pre>
 
=={{header|Wren}}==
{{trans|Kotlin (version 2)}}
{{libheader|Wren-dynamic}}
{{libheader|Wren-fmt}}
{{libheader|Wren-seq}}
{{libheader|Wren-math}}
{{libheader|Wren-sort}}
<syntaxhighlight lang="wren">import "./dynamic" for Enum
import "./fmt" for Fmt
import "./seq" for Lst
import "./math" for Nums
import "./sort" for Sort
 
var Direction = Enum.create("Direction", ["UP", "DOWN"])
 
// test also for 'Unknown' correction type
var types = [
"Benjamini-Hochberg", "Benjamini-Yekutieli", "Bonferroni", "Hochberg",
"Holm", "Hommel", "Šidák", "Unknown"
]
 
var pFormat = Fn.new { |p, cols|
var i = -cols
var fmt = "$1.10f"
return Lst.chunks(p, cols).map { |chunk|
i = i + cols
return Fmt.swrite("[$2d $s", i, chunk.map { |v| Fmt.swrite(fmt, v) }.join(" "))
}.join("\n")
}
 
var check = Fn.new { |p|
if (p.count == 0 || Nums.min(p) < 0 || Nums.max(p) > 1) {
Fiber.abort("p-values must be in range 0 to 1")
}
return p
}
 
var ratchet = Fn.new { |p, dir|
var pp = p.toList
var m = pp[0]
if (dir == Direction.UP) {
for (i in 1...pp.count) {
if (pp[i] > m) pp[i] = m
m = pp[i]
}
} else {
for (i in 1...pp.count) {
if (pp[i] < m) pp[i] = m
m = pp[i]
}
}
return pp.map { |v| (v < 1) ? v : 1 }.toList
}
 
var schwartzian = Fn.new { |p, mult, dir|
var size = p.count
var pwi = List.filled(size, null)
for (i in 0...size) pwi[i] = [i, p[i]]
var cmp = (dir == Direction.UP) ? Fn.new { |a, b| (b[1] - a[1]).sign } :
Fn.new { |a, b| (a[1] - b[1]).sign }
var order = Sort.merge(pwi, cmp).map { |e| e[0] }.toList
var pa = List.filled(size, 0)
for (i in 0...size) pa[i] = mult[i] * p[order[i]]
pa = ratchet.call(pa, dir)
var owi = List.filled(order.count, null)
for (i in 0...order.count) owi[i] = [i, order[i]]
cmp = Fn.new { |a, b| (a[1] - b[1]).sign }
var order2 = Sort.merge(owi, cmp).map { |e| e[0] }.toList
var res = List.filled(size, 0)
for (i in 0...size) res[i] = pa[order2[i]]
return res
}
 
var adjust = Fn.new { |p, type|
var size = p.count
if (size == 0) Fiber.abort("List cannot be empty.")
if (type == "Benjamini-Hochberg") {
var mult = List.filled(size, 0)
for (i in 0...size) mult[i] = size / (size - i)
return schwartzian.call(p, mult, Direction.UP)
 
} else if (type == "Benjamini-Yekutieli") {
var q = (1..size).reduce { |acc, i| acc + 1/i }
var mult = List.filled(size, 0)
for (i in 0...size) mult[i] = q * size / (size - i)
return schwartzian.call(p, mult, Direction.UP)
 
} else if (type == "Bonferroni") {
return p.map { |v| (v * size).min(1) }.toList
 
} else if (type == "Hochberg") {
var mult = List.filled(size, 0)
for (i in 0...size) mult[i] = i + 1
return schwartzian.call(p, mult, Direction.UP)
 
} else if (type == "Holm") {
var mult = List.filled(size, 0)
for (i in 0...size) mult[i] = size - i
return schwartzian.call(p, mult, Direction.DOWN)
 
} else if (type == "Hommel") {
var pwi = List.filled(size, null)
for (i in 0...size) pwi[i] = [i, p[i]]
var cmp = Fn.new { |a, b| (a[1] - b[1]).sign }
var order = Sort.merge(pwi, cmp).map { |e| e[0] }.toList
var s = List.filled(size, 0)
for (i in 0...size) s[i] = p[order[i]]
var m = List.filled(size, 0)
for (i in 0...size) m[i] = s[i] * size / (i + 1)
var min = Nums.min(m)
var q = List.filled(size, min)
var pa = List.filled(size, min)
for (j in size-1..2) {
var lower = List.filled(size - j + 1, 0) // lower indices
for (i in 0...lower.count) lower[i] = i
var upper = List.filled(j - 1, 0) // upper indices
for (i in 0...upper.count) upper[i] = size - j + 1 + i
var qmin = j * s[upper[0]] / 2
for (i in 1...upper.count) {
var temp = s[upper[i]] * j / (2 + i)
if (temp < qmin) qmin = temp
}
for (i in 0...lower.count) {
q[lower[i]] = qmin.min(s[lower[i]] * j)
}
for (i in 0...upper.count) q[upper[i]] = q[size - j]
for (i in 0...size) if (pa[i] < q[i]) pa[i] = q[i]
}
var owi = List.filled(order.count, null)
for (i in 0...order.count) owi[i] = [i, order[i]]
var order2 = Sort.merge(owi, cmp).map { |e| e[0] }.toList
var res = List.filled(size, 0)
for (i in 0...size) res[i] = pa[order2[i]]
return res
 
} else if (type == "Šidák") {
return p.map { |v| 1 - (1 - v).pow(size) }.toList
 
} else {
System.print("\nSorry, do not know how to do '%(type)' correction.\n" +
"Perhaps you want one of these?:\n" +
types[0...-1].map { |t| " %(t)" }.join("\n")
)
Fiber.suspend()
}
}
 
var adjusted = Fn.new { |p, type| "\n%(type)\n%(pFormat.call(adjust.call(check.call(p), type), 5))" }
 
var 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
]
types.each { |type| System.print(adjusted.call(pValues, type)) }</syntaxhighlight>
 
{{out}}
<pre>
Same as Kotlin (version 2) entry.
</pre>
 
=={{header|zkl}}==
Line 4,893 ⟶ 5,885:
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
 
<langsyntaxhighlight lang="zkl">fcn bh(pvalues){ // Benjamini-Hochberg
psz,pszf := pvalues.len(), psz.toFloat();
n_i := psz.pump(List,'wrap(n){ pszf/(psz - n) }); # N/(N-0),N/(N-1),..
Line 4,957 ⟶ 5,949:
}
psz.pump(List,'wrap(n){ pa[ro[n]] }); // Hommel q-values
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">pvalues:=T(
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,
Line 4,983 ⟶ 5,975:
}
println();
}</langsyntaxhighlight>
{{out}}
<pre style="height:45ex">
9,476

edits