P-value correction

From Rosetta Code
P-value correction is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Given a list of 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 "FDR" After adjustment, the p-values will be significantly higher.

Task Description

Given one list of p-values, return the p-values correcting for multiple comparisons

   p = {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}

There are numerous implementations of how to do this, namely Benjamini-Hochberg, Benjamini-Yekutieli, Holm, Hochberg, Hommel, Bonferroni. Each of which has its own advantages and disadvantages.

C[edit]

Works with: C99
Translation of: R

This work is a translation of the R source code. In order to confirm that the new function is working correctly, each value is compared to R's output and a cumulative absolute error is returned.

The C function p_adjust is designed to work as similarly to the R function p.adjust as possible, and is able to do any one of the methods.

This program, for example, fdr.c, can be compiled by

gcc -o fdr fdr.c -Wall -pedantic -std=c11 -lm -O4

or

clang -o fdr fdr.c -Wall -pedantic -std=c11 -lm -O4.

Link with -lm

#include <stdio.h>//printf
#include <stdlib.h>//qsort
#include <math.h>//fabs
#include <stdbool.h>//bool data type
#include <strings.h>//strcasecmp
 
unsigned int *restrict seq_len(const size_t START, const size_t END) {
//named after R function of same name, but simpler function
size_t start = START;
size_t end = END;
if (START == END) {
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int));
if (sequence == NULL) {
printf("malloc failed at %s line %u\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
for (size_t i = 0; i < end; i++) {
sequence[i] = i+1;
}
return sequence;
}
if (START > END) {
end = START;
start = END;
}
const size_t LENGTH = end - start ;
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int));
if (sequence == NULL) {
printf("malloc failed at %s line %u\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
if (START < END) {
for (size_t index = 0; index <= LENGTH; index++) {
sequence[index] = start + index;
}
} else {
for (size_t index = 0; index <= LENGTH; index++) {
sequence[index] = end - index;
}
}
return sequence;
}
 
//modified from https://phoxis.org/2012/07/12/get-sorted-index-orderting-of-an-array/
 
double *restrict base_arr = NULL;
 
static int compar_increase (const void *restrict a, const void *restrict b) {
int aa = *((int *restrict ) a), bb = *((int *restrict) b);
if (base_arr[aa] < base_arr[bb]) {
return 1;
} else if (base_arr[aa] == base_arr[bb]) {
return 0;
} else {
return -1;
}
}
 
static int compar_decrease (const void *restrict a, const void *restrict b) {
int aa = *((int *restrict ) a), bb = *((int *restrict) b);
if (base_arr[aa] < base_arr[bb]) {
return -1;
} else if (base_arr[aa] == base_arr[bb]) {
return 0;
} else {
return 1;
}
}
 
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));
if (idx == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
base_arr = malloc(sizeof(double) * SIZE);
if (base_arr == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
for (unsigned int i = 0; i < SIZE; i++) {
base_arr[i] = ARRAY[i];
idx[i] = i;
}
if (DECREASING == false) {
qsort(idx, SIZE, sizeof(unsigned int), compar_decrease);
} else if (DECREASING == true) {
qsort(idx, SIZE, sizeof(unsigned int), compar_increase);
}
free(base_arr); base_arr = NULL;
return idx;
}
 
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
if (NO_OF_ARRAY_ELEMENTS < 1) {
puts("cummin function requires at least one element.\n");
printf("Failed at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
double *restrict output = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (output == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
double cumulative_min = ARRAY[0];
for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
if (ARRAY[i] < cumulative_min) {
cumulative_min = ARRAY[i];
}
output[i] = cumulative_min;
}
return output;
}
 
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
if (NO_OF_ARRAY_ELEMENTS < 1) {
puts("function requires at least one element.\n");
printf("Failed at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
double *restrict output = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (output == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
double cumulative_max = ARRAY[0];
for (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {
if (ARRAY[i] > cumulative_max) {
cumulative_max = ARRAY[i];
}
output[i] = cumulative_max;
}
return output;
}
 
double *restrict pminx(const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS, const double X) {
//named after the R function pmin
if (NO_OF_ARRAY_ELEMENTS < 1) {
puts("pmin requires at least one element.\n");
printf("Failed at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
double *restrict pmin_array = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (pmin_array == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
if (ARRAY[index] < X) {
pmin_array[index] = ARRAY[index];
} else {
pmin_array[index] = X;
}
}
return pmin_array;
}
 
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
printf("[1] %e", ARRAY[0]);
for (size_t i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
printf(" %.10f", ARRAY[i]);
if (((i+1) % 5) == 0) {
printf("\n[%zu]", i+1);
}
}
puts("\n");
}
 
/*void uint_say (const unsigned int *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
//for debugging
printf("%u", ARRAY[0]);
for (size_t i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) {
printf(",%u", ARRAY[i]);
}
puts("\n");
}*/

 
double *restrict uint2double (const unsigned int *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (doubleArray == NULL) {
printf("Failure to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
doubleArray[index] = (double)ARRAY[index];
}
return doubleArray;
}
 
double min2 (const double N1, const double N2) {
if (N1 < N2) {
return N1;
} else {
return N2;
}
}
 
double *restrict p_adjust (const double *restrict PVALUES, const size_t 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
if (NO_OF_ARRAY_ELEMENTS < 1) {
puts("p_adjust requires at least one element.\n");
printf("Failed at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
short int TYPE = -1;
if (strcasecmp(STRING, "BH") == 0) {
TYPE = 0;
} else if (strcasecmp(STRING, "fdr") == 0) {
TYPE = 0;
} else if (strcasecmp(STRING, "by") == 0) {
TYPE = 1;
} else if (strcasecmp(STRING, "Bonferroni") == 0) {
TYPE = 2;
} else if (strcasecmp(STRING, "hochberg") == 0) {
TYPE = 3;
} else if (strcasecmp(STRING, "holm") == 0) {
TYPE = 4;
} else if (strcasecmp(STRING, "hommel") == 0) {
TYPE = 5;
} else {
printf("%s doesn't match any accepted FDR methods.\n", STRING);
printf("Failed at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
if (TYPE == 2) {//Bonferroni method
double *restrict bonferroni = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (bonferroni == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS;
if (BONFERRONI >= 1.0) {
bonferroni[index] = 1.0;
} else if ((0.0 <= BONFERRONI) && (BONFERRONI < 1.0)) {
bonferroni[index] = BONFERRONI;
} else {
printf("%g is outside of the interval I planned.\n", BONFERRONI);
printf("Failure at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
}
return bonferroni;
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
} else if (TYPE == 4) {//Holm method
/*these values are computed separately from BH, BY, and Hochberg because they are
computed differently*/

unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, false);
//sorted in reverse of methods 0-3
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * PVALUES[o[index]];
// printf("cummax_input[%zu] = %e\n", index, cummax_input[index]);
}
free(o); o = NULL;
unsigned int *restrict ro = order(o2double, NO_OF_ARRAY_ELEMENTS, false);
free(o2double); o2double = NULL;
 
double *restrict cummax_output = cummax(cummax_input, NO_OF_ARRAY_ELEMENTS);
free(cummax_input); cummax_input = NULL;
 
double *restrict pmin = pminx(cummax_output, NO_OF_ARRAY_ELEMENTS, 1);
free(cummax_output); cummax_output = NULL;
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
qvalues[index] = pmin[ro[index]];
}
free(pmin); pmin = NULL;
free(ro); ro = NULL;
return qvalues;
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
} else if (TYPE == 5) {//Hommel method
//i <- seq_len(n)
//o <- order(p)
unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, false);//false is R's default
//p <- p[o]
double *restrict p = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (p == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
p[index] = PVALUES[o[index]];
}
//ro <- order(o)
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
free(o); o = NULL;
unsigned int *restrict ro = order(o2double, NO_OF_ARRAY_ELEMENTS, false);
free(o2double); o2double = NULL;
// puts("ro");
//q <- pa <- rep.int(min(n * p/i), n)
double *restrict q = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (q == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
double *restrict pa = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (pa == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0];
for (size_t index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (1+index);
if (TEMP < min) {
min = TEMP;
}
}
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
pa[index] = min;
q[index] = min;
}
// puts("q & pa");
// double_say(q, NO_OF_ARRAY_ELEMENTS);
/*for (j in (n - 1):2) {
ij <- seq_len(n - j + 1)
i2 <- (n - j + 2):n
q1 <- min(j * p[i2]/(2:j))
q[ij] <- pmin(j * p[ij], q1)
q[i2] <- q[n - j + 1]
pa <- pmax(pa, q)
}
*/

for (size_t j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) {
// printf("j = %zu\n", j);
unsigned int *restrict ij = seq_len(1,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_t 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 = j * p[i2[0]] / 2.0;
for (size_t i = 1; i < I2_LENGTH; i++) {//loop through 2:j
const double TEMP_Q1 = (double)j * p[i2[i]] / (2 + i);
if (TEMP_Q1 < q1) {
q1 = TEMP_Q1;
}
}
 
for (size_t i = 0; i < (NO_OF_ARRAY_ELEMENTS - j + 1); i++) {//q[ij] <- pmin(j * p[ij], q1)
q[ij[i]] = min2( j*p[ij[i]], q1);
}
free(ij); ij = NULL;
 
for (size_t 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_t i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q)
if (pa[i] < q[i]) {
pa[i] = q[i];
}
}
// printf("j = %zu, pa = \n", j);
// double_say(pa, N);
}//end j loop
free(p); p = NULL;
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
q[index] = pa[ro[index]];//Hommel q-values
}
//now free memory
free(ro); ro = NULL;
free(pa); pa = NULL;
return q;
}
//The methods are similarly computed and thus can be combined for clarity
unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, true);
if (o == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
double *restrict o_double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
if ((PVALUES[index] < 0) || (PVALUES[index] > 1)) {
printf("array[%u] = %lf, which is outside the interval [0,1]\n", index, PVALUES[index]);
printf("died at %s line %u\n", __FILE__, __LINE__);
exit(EXIT_FAILURE);
}
}
 
unsigned int *restrict ro = order(o_double, NO_OF_ARRAY_ELEMENTS, false);
if (ro == NULL) {
printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__);
perror("");
exit(EXIT_FAILURE);
}
free(o_double); o_double = NULL;
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (TYPE == 0) {//BH method
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (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_t index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) {
q += (double) 1.0/index;
}
for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (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_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
// pmin(1, cummin((n - i + 1L) * p[o]))[ro]
cummin_input[index] = (index + 1) * PVALUES[o[index]];
}
}
free(o); o = NULL;
double *restrict cummin_array = NULL;
cummin_array = cummin(cummin_input, NO_OF_ARRAY_ELEMENTS);
free(cummin_input); cummin_input = NULL;//I don't need this anymore
double *restrict pmin = pminx(cummin_array, NO_OF_ARRAY_ELEMENTS, 1);
free(cummin_array); cummin_array = NULL;
double *restrict q_array = malloc(NO_OF_ARRAY_ELEMENTS*sizeof(double));
for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) {
q_array[index] = pmin[ro[index]];
}
 
free(ro); ro = NULL;
free(pmin); pmin = NULL;
return q_array;
}
 
int main(void) {
const 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};//just the pvalues
const double CORRECT_ANSWERS[6][50] = {//each first index is type
{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-Hochberg
{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},//Benjamini & Yekutieli
{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},//Bonferroni
{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},//Hochberg
{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},//Holm
{ 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}//Hommel
};
//the following loop checks each type with R's answers
const char *restrict TYPES[] = {"bh", "by", "bonferroni", "hochberg", "holm", "hommel"};
for (unsigned short int type = 0; type <= 5; type++) {
double *restrict q = p_adjust(PVALUES, sizeof(PVALUES) / sizeof(*PVALUES), TYPES[type]);
double error = fabs(q[0] - CORRECT_ANSWERS[type][0]);
// printf("%e - %e = %g\n", q[0], CORRECT_ANSWERS[type][0], error);
// puts("p q");
// printf("%g\t%g\n", pvalues[0], q[0]);
for (unsigned int i = 1; i < sizeof(PVALUES) / sizeof(*PVALUES); i++) {
const double this_error = fabs(q[i] - CORRECT_ANSWERS[type][i]);
// printf("%e - %e = %g\n", q[i], CORRECT_ANSWERS[type][i], error);
error += this_error;
}
double_say(q, sizeof(PVALUES) / sizeof(*PVALUES));
free(q); q = NULL;
printf("\ntype %u = '%s' has cumulative error of %g\n", type, TYPES[type], error);
}
 
return 0;
}
 
Output:
[1] 6.126681e-01 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 cumulative error of 8.03053e-07
[1] 1.000000e+00 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 cumulative error of 3.64072e-07
[1] 1.000000e+00 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 cumulative error of 6.5e-08
[1] 9.991834e-01 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 cumulative error of 2.7375e-07
[1] 1.000000e+00 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 cumulative error of 2.8095e-07
[1] 9.991834e-01 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 cumulative error of 4.35302e-07

Kotlin[edit]

Translation of: C
// version 1.1.51
 
import java.util.Arrays
 
typealias IAE = IllegalArgumentException
 
fun seqLen(start: Int, end: Int) =
when {
start == end -> IntArray(end + 1) { it + 1 }
start < end -> IntArray(end - start + 1) { start + it }
else -> IntArray(start - end + 1) { start - it }
}
 
var baseArr: DoubleArray? = null
 
fun compareIncrease(a: Int, b: Int): Int = baseArr!![b].compareTo(baseArr!![a])
 
fun compareDecrease(a: Int, b: Int): Int = baseArr!![a].compareTo(baseArr!![b])
 
fun order(array: DoubleArray, decreasing: Boolean): IntArray {
val size = array.size
var idx = IntArray(size) { it }
baseArr = array.copyOf()
if (!decreasing) {
idx = Arrays.stream(idx)
.boxed()
.sorted { a, b -> compareDecrease(a, b) }
.mapToInt { it }
.toArray()
}
else {
idx = Arrays.stream(idx)
.boxed()
.sorted { a, b -> compareIncrease(a, b) }
.mapToInt { it }
.toArray()
}
baseArr = null
return idx
}
 
fun cummin(array: DoubleArray): DoubleArray {
val size = array.size
if (size < 1) throw IAE("cummin requires at least one element")
val output = DoubleArray(size)
var cumulativeMin = array[0]
for (i in 0 until size) {
if (array[i] < cumulativeMin) cumulativeMin = array[i]
output[i] = cumulativeMin
}
return output
}
 
fun cummax(array: DoubleArray): DoubleArray {
val size = array.size
if (size < 1) throw IAE("cummax requires at least one element")
val output = DoubleArray(size)
var cumulativeMax = array[0]
for (i in 0 until size) {
if (array[i] > cumulativeMax) cumulativeMax = array[i]
output[i] = cumulativeMax
}
return output
}
 
fun pminx(array: DoubleArray, x: Double): DoubleArray {
val size = array.size
if (size < 1) throw IAE("pmin requires at least one element")
return DoubleArray(size) { if (array[it] < x) array[it] else x }
}
 
fun doubleSay(array: DoubleArray) {
print("[ 1] %e".format(array[0]))
for (i in 1 until array.size) {
print(" %.10f".format(array[i]))
if ((i + 1) % 5 == 0) print("\n[%2d]".format(i + 1))
}
println()
}
 
fun intToDouble(array: IntArray) = DoubleArray(array.size) { array[it].toDouble() }
 
fun doubleArrayMin(array: DoubleArray) =
if (array.size < 1) throw IAE("pAdjust requires at least one element")
else array.min()!!
 
fun pAdjust(pvalues: DoubleArray, str: String): DoubleArray {
val size = pvalues.size
if (size < 1) throw IAE("pAdjust requires at least one element")
val type = when(str.toLowerCase()) {
"bh", "fdr" -> 0
"by" -> 1
"bonferroni" -> 2
"hochberg" -> 3
"holm" -> 4
"hommel" -> 5
else -> throw IAE("'$str' doesn't match any accepted FDR types")
}
if (type == 2) { // Bonferroni method
return DoubleArray(size) {
val b = pvalues[it] * size
when {
b >= 1 -> 1.0
0 <= b && b < 1 -> b
else -> throw RuntimeException("$b is outside [0, 1)")
}
}
}
else if (type == 4) { // Holm method
val o = order(pvalues, false)
val o2Double = intToDouble(o)
val cummaxInput = DoubleArray(size) { (size - it) * pvalues[o[it]] }
val ro = order(o2Double, false)
val cummaxOutput = cummax(cummaxInput)
val pmin = pminx(cummaxOutput, 1.0)
return DoubleArray(size) { pmin[ro[it]] }
}
else if (type == 5) { // Hommel method
val indices = seqLen(size, size)
val o = order(pvalues, false)
val p = DoubleArray(size) { pvalues[o[it]] }
val o2Double = intToDouble(o)
val ro = order(o2Double, false)
val q = DoubleArray(size)
val pa = DoubleArray(size)
val npi = DoubleArray(size) { p[it] * size / indices[it] }
val min = doubleArrayMin(npi)
q.fill(min)
pa.fill(min)
for (j in size - 1 downTo 2) {
val ij = seqLen(1, size - j + 1)
for (i in 0 until size - j + 1) ij[i]--
val i2Length = j - 1
val i2 = IntArray(i2Length) { size - j + 2 + it - 1 }
val pi2Length = i2Length
var q1 = j * p[i2[0]] / 2.0
for (i in 1 until pi2Length) {
val temp_q1 = p[i2[i]] * j / (2.0 + i)
if(temp_q1 < q1) q1 = temp_q1
}
for (i in 0 until size - j + 1) {
q[ij[i]] = minOf(p[ij[i]] * j, q1)
}
for (i in 0 until i2Length) q[i2[i]] = q[size - j]
for (i in 0 until size) if (pa[i] < q[i]) pa[i] = q[i]
}
for (index in 0 until size) q[index] = pa[ro[index]]
return q
}
val ni = DoubleArray(size)
val o = order(pvalues, true)
val oDouble = intToDouble(o)
for (index in 0 until size) {
if (pvalues[index] !in 0.0 .. 1.0) {
throw RuntimeException("array[$index] = ${pvalues[index]} is outside [0, 1]")
}
ni[index] = size.toDouble() / (size - index)
}
val ro = order(oDouble, false)
val cumminInput = DoubleArray(size)
if (type == 0) { // BH method
for (index in 0 until size) {
cumminInput[index] = ni[index] * pvalues[o[index]]
}
}
else if (type == 1) { // BY method
var q = 0.0
for (index in 1 until size + 1) q += 1.0 / index
for (index in 0 until size) {
cumminInput[index] = q * ni[index] * pvalues[o[index]]
}
}
else if (type == 3) { // Hochberg method
for (index in 0 until size) {
cumminInput[index] = (index + 1) * pvalues[o[index]]
}
}
val cumminArray = cummin(cumminInput)
val pmin = pminx(cumminArray, 1.0)
return DoubleArray(size) { pmin[ro[it]] }
}
 
fun main(args: Array<String>) {
val pvalues = doubleArrayOf(
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
)
 
val correctAnswers = listOf(
doubleArrayOf( // 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
),
doubleArrayOf( // 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
),
doubleArrayOf( // 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
),
doubleArrayOf( // 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
),
doubleArrayOf( // 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
),
doubleArrayOf( // 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
)
)
val types = listOf("bh", "by", "bonferroni", "hochberg", "holm", "hommel")
val f = "\ntype %d = '%s' has cumulative error of %g"
for (type in 0 until types.size) {
val q = pAdjust(pvalues, types[type])
var error = 0.0
for (i in 0 until pvalues.size) {
error += Math.abs(q[i] - correctAnswers[type][i])
}
doubleSay(q)
println(f.format(type, types[type], error))
}
}
Output:
[ 1] 6.126681e-01 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.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 cumulative error of 8.03053e-07
[ 1] 1.000000e+00 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 cumulative error of 3.64072e-07
[ 1] 1.000000e+00 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 cumulative error of 6.50000e-08
[ 1] 9.991834e-01 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 cumulative error of 2.73750e-07
[ 1] 1.000000e+00 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 cumulative error of 2.80950e-07
[ 1] 9.991834e-01 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 cumulative error of 4.35302e-07

Perl[edit]

Translation of: C
#!/usr/bin/env perl
 
use strict; use warnings; use Cwd;
my $TOP_DIRECTORY = getcwd();
local $SIG{__WARN__} = sub {#kill the program if there are any warnings
my $message = shift;
my $fail_filename = "$TOP_DIRECTORY/$0.FAIL";
open my $fh, '>', $fail_filename or die "Can't write $fail_filename: $!";
printf $fh ("$message @ %s\n", getcwd());
close $fh;
die "$message\n";
};#http://perlmaven.com/how-to-capture-and-save-warnings-in-perl
 
sub pmin {
my $array_ref = shift;
my $x = 1;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "cummin requires an array.\n";
die;
}
my @pmin_array;
my $n = scalar @$array_ref;
for (my $index = 0; $index < $n; $index++) {
if (@$array_ref[$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];
foreach my $p (@$array_ref) {
if ($p < $cumulative_min) {
$cumulative_min = $p;
}
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];
foreach my $p (@$array_ref) {
if ($p > $cumulative_max) {
$cumulative_max = $p;
}
push @cummax, $cumulative_max;
}
return @cummax;
}
 
sub order {#made to match R's "order"
my $array_ref = shift;
my $decreasing = 'false';
if (defined $_[0]) {
my $option = shift;
if ($option =~ m/true/i) {
$decreasing = 'true';
} elsif ($option =~ m/false/i) {
#do nothing, it's already set to false
} else {
print "2nd option should only be case-insensitive 'true' or 'false'";
die;
}
}
unless ((ref $array_ref) =~ m/ARRAY/) {
print "You should have entered an array.\n";
die;
}
my @array;
my $max_index = scalar @$array_ref-1;
if ($decreasing eq 'false') {
@array = sort { @$array_ref[$a] <=> @$array_ref[$b] } 0..$max_index;
} elsif ($decreasing eq 'true') {
@array = sort { @$array_ref[$b] <=> @$array_ref[$a] } 0..$max_index;
}
}
 
sub min {
my $array_ref = shift;
unless ((ref $array_ref) =~ m/ARRAY/) {
print "min requires an array.\n";
die;
}
my $min = @$array_ref[0];
foreach my $e (@$array_ref) {
if ($e < $min) {
$min = $e;
}
}
return $min;
}
 
sub min2 {
my $n1 = shift;
my $n2 = shift;
if ($n1 < $n2) {
return $n1;
} else {
return $n2;
}
}
 
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 = (
'bh' => 1,
'fdr' => 1,
'by' => 1,
'holm' => 1,
'hommel' => 1,
'bonferroni' => 1,
'hochberg' => 1
);
my $method_found = 'no';
foreach my $key (keys %methods) {
if ((uc $method) eq (uc $key)) {
$method = $key;
$method_found = 'yes';
last;
}
}
if ($method_found eq 'no') {
if ($method =~ m/benjamini-?\s*hochberg/i) {
$method = 'bh';
$method_found = 'yes';
} elsif ($method =~ m/benjamini-?\s*yekutieli/i) {
$method = 'by';
$method_found = 'yes';
}
}
if ($method_found eq 'no') {
print "No method could be determined from $method.\n";
die;
}
my $lp = scalar @$pvalues_ref;
my $n = $lp;
my @qvalues;
if ($method eq 'hochberg') {
my @o = order($pvalues_ref, 'TRUE');
my @cummin_input;
for (my $index = 0; $index < $n; $index++) {
$cummin_input[$index] = ($index+1)* @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o]
}
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') {
my @o = order($pvalues_ref, 'TRUE');
my @cummin_input;
for (my $index = 0; $index < $n; $index++) {
$cummin_input[$index] = ($n/($n-$index))* @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o]
}
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') {
my $q = 0.0;
my @o = order($pvalues_ref, 'TRUE');
my @ro = order(\@o);
for (my $index = 1; $index < ($n+1); $index++) {
$q += 1.0 / $index;
}
my @cummin_input;
for (my $index = 0; $index < $n; $index++) {
$cummin_input[$index] = $q * ($n/($n-$index)) * @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o]
}
my @cummin = cummin(\@cummin_input);
undef @cummin_input;
my @pmin = pmin(\@cummin);
undef @cummin;
@qvalues = @pmin[@ro];
} elsif ($method eq 'bonferroni') {
for (my $index = 0; $index < $n; $index++) {
my $q = @$pvalues_ref[$index]*$n;
if ((0 <= $q) && ($q < 1)) {
$qvalues[$index] = $q;
} elsif ($q >= 1) {
$qvalues[$index] = 1.0;
} else {
print "Failed to get Bonferroni adjusted p.";
die;
}
}
} elsif ($method eq 'holm') {
my @o = order($pvalues_ref);
my @cummax_input;
for (my $index = 0; $index < $n; $index++) {
$cummax_input[$index] = ($n - $index) * @$pvalues_ref[$o[$index]];
}
my @ro = order(\@o);
undef @o;
my @cummax = cummax(\@cummax_input);
undef @cummax_input;
my @pmin = pmin(\@cummax);
undef @cummax;
@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 @pa;
my @q;
my $min = $n*$p[0];
for (my $index = 0; $index < $n; $index++) {
my $temp = $n*$p[$index] / ($index + 1);
if ($temp < $min) {
$min = $temp;
}
}
for (my $index = 0; $index < $n; $index++) {
$pa[$index] = $min;#q <- pa <- rep.int(min(n * p/i), n)
$q[$index] = $min;#q <- pa <- rep.int(min(n * p/i), n)
}
for (my $j = ($n-1); $j >= 2; $j--) {
# 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;
for (my $i = 0; $i < $I2_LENGTH; $i++) {
$i2[$i] = $n-$j+2+$i-1;
#R's indices are 1-based, C's are 0-based, I added the -1
}
 
my $q1 = $j * $p[$i2[0]] / 2.0;
for (my $i = 1; $i < $I2_LENGTH; $i++) {#loop through 2:j
my $TEMP_Q1 = $j * $p[$i2[$i]] / (2 + $i);
if ($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]] = min2( $j*$p[$ij[$i]], $q1);
}
 
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
}
 
for (my $i = 0; $i < $n; $i++) {#pa <- pmax(pa, q)
if ($pa[$i] < $q[$i]) {
$pa[$i] = $q[$i];
}
}
# printf("j = %zu, pa = \n", j);
# double_say(pa, N);
}#end j loop
@qvalues = @pa[@ro];
} 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,
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);
 
my %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]);
 
 
foreach my $method ('Hochberg','Benjamini-Hochberg','Benjamini-Yekutieli', 'Bonferroni', 'Holm', 'Hommel') {
print "$method\n";
my @qvalues = p_adjust(\@pvalues, $method);
my $error = 0.0;
foreach my $q (0..$#qvalues) {
$error += abs($qvalues[$q] - $correct_answers{$method}[$q]);
}
printf("type $method has cumulative error of %g.\n", $error);
}
 
Output:
Hochberg
type Hochberg has cumulative error of 2.7375e-07.
Benjamini-Hochberg
type Benjamini-Hochberg has cumulative error of 8.03053e-07.
Benjamini-Yekutieli
type Benjamini-Yekutieli has cumulative error of 3.64072e-07.
Bonferroni
type Bonferroni has cumulative error of 6.5e-08.
Holm
type Holm has cumulative error of 2.8095e-07.
Hommel
type Hommel has cumulative error of 4.35302e-07.

Perl 6[edit]

Works with: Rakudo version 2017.10
########################### 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, $_
}
Output:
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

R[edit]

 
#p.adjust's source code is below
#> p.adjust
#function (p, method = p.adjust.methods, n = length(p))
#{
# method <- match.arg(method)
# if (method == "fdr")
# method <- "BH"
# nm <- names(p)
# p <- as.numeric(p)
# p0 <- setNames(p, nm)
# if (all(nna <- !is.na(p)))
# nna <- TRUE
# p <- p[nna]
# lp <- length(p)
# stopifnot(n >= lp)
# if (n <= 1)
# return(p0)
# if (n == 2 && method == "hommel")
# method <- "hochberg"
# p0[nna] <- switch(method, bonferroni = pmin(1, n * p), holm = {
# i <- seq_len(lp)
# o <- order(p)
# ro <- order(o)
# pmin(1, cummax((n - i + 1L) * p[o]))[ro]
# }, hommel = {
# if (n > lp) p <- c(p, rep.int(1, n - lp))
# i <- seq_len(n)
# o <- order(p)
# p <- p[o]
# ro <- order(o)
# q <- pa <- rep.int(min(n * p/i), n)
## for (j in (n - 1):2) {
# ij <- seq_len(n - j + 1)
# i2 <- (n - j + 2):n
# q1 <- min(j * p[i2]/(2:j))
# q[ij] <- pmin(j * p[ij], q1)
# q[i2] <- q[n - j + 1]
# pa <- pmax(pa, q)
# }
# pmax(pa, p)[if (lp < n) ro[1:lp] else ro]
# }, hochberg = {
# i <- lp:1L
# o <- order(p, decreasing = TRUE)
# ro <- order(o)
# pmin(1, cummin((n - i + 1L) * p[o]))[ro]
# }, BH = {
# i <- lp:1L
# o <- order(p, decreasing = TRUE)
# ro <- order(o)
# pmin(1, cummin(n/i * p[o]))[ro]
# }, BY = {
# i <- lp:1L
# o <- order(p, decreasing = TRUE)
# ro <- order(o)
# q <- sum(1L/(1L:n))
# pmin(1, cummin(q * n/i * p[o]))[ro]
# }, none = p)
# p0
#}
#<bytecode: 0x3a61b88>
#<environment: namespace:stats>
 
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,
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)
 
p.adjust(p, method = 'BH')
print("Benjamini-Hochberg")
writeLines("\n")
 
p.adjust(p, method = 'BY')
print("Benjamini & Yekutieli")
writeLines("\n")
 
p.adjust(p, method = 'bonferroni')
print("Bonferroni")
writeLines("\n")
 
p.adjust(p, method = 'hochberg')
print("Hochberg")
writeLines("\n");
 
p.adjust(p, method = 'hommel')
writeLines("Hommel\n")
 
Output:
 [1] 6.126681e-01 8.521710e-01 1.987205e-01 1.891595e-01 3.217789e-01
 [6] 9.301450e-01 4.870370e-01 9.301450e-01 6.049731e-01 6.826753e-01
[11] 6.482629e-01 7.253722e-01 5.280973e-01 8.769926e-01 4.705703e-01
[16] 9.241867e-01 6.049731e-01 7.856107e-01 4.887526e-01 1.136717e-01
[21] 4.991891e-01 8.769926e-01 9.991834e-01 3.217789e-01 9.301450e-01
[26] 2.304958e-01 5.832475e-01 3.899547e-02 8.521710e-01 1.476843e-01
[31] 1.683638e-02 2.562902e-03 3.516084e-02 6.250189e-02 3.636589e-03
[36] 2.562902e-03 2.946883e-02 6.166064e-03 3.899547e-02 2.688991e-03
[41] 4.502862e-04 1.252228e-05 7.881555e-02 3.142613e-02 4.846527e-03
[46] 2.562902e-03 4.846527e-03 1.101708e-03 7.252032e-02 2.205958e-02
[1] "Benjamini-Hochberg"


 [1] 1.000000e+00 1.000000e+00 8.940844e-01 8.510676e-01 1.000000e+00
 [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[16] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 5.114323e-01
[21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[26] 1.000000e+00 1.000000e+00 1.754486e-01 1.000000e+00 6.644618e-01
[31] 7.575031e-02 1.153102e-02 1.581959e-01 2.812089e-01 1.636176e-02
[36] 1.153102e-02 1.325863e-01 2.774239e-02 1.754486e-01 1.209832e-02
[41] 2.025930e-03 5.634031e-05 3.546073e-01 1.413926e-01 2.180552e-02
[46] 1.153102e-02 2.180552e-02 4.956812e-03 3.262838e-01 9.925057e-02
[1] "Benjamini & Yekutieli"


 [1] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
 [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[16] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[26] 1.000000e+00 1.000000e+00 7.019185e-01 1.000000e+00 1.000000e+00
[31] 2.020365e-01 1.516674e-02 5.625735e-01 1.000000e+00 2.909271e-02
[36] 1.537741e-02 4.125636e-01 6.782670e-02 6.803480e-01 1.882294e-02
[41] 9.005725e-04 1.252228e-05 1.000000e+00 4.713920e-01 4.395577e-02
[46] 1.088915e-02 4.846527e-02 3.305125e-03 1.000000e+00 2.867745e-01
[1] "Bonferroni"


 [1] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
 [6] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[11] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[16] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[21] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[26] 9.991834e-01 9.991834e-01 4.632662e-01 9.991834e-01 9.991834e-01
[31] 1.575885e-01 1.383967e-02 3.938014e-01 7.600230e-01 2.501973e-02
[36] 1.383967e-02 3.052971e-01 5.426136e-02 4.626366e-01 1.656419e-02
[41] 8.825610e-04 1.252228e-05 9.930759e-01 3.394022e-01 3.692284e-02
[46] 1.023581e-02 3.974152e-02 3.172920e-03 8.992520e-01 2.179486e-01
[1] "Hochberg"


 [1] 9.991834e-01 9.991834e-01 9.991834e-01 9.987624e-01 9.991834e-01
 [6] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[11] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[16] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.595180e-01
[21] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[26] 9.991834e-01 9.991834e-01 4.351895e-01 9.991834e-01 9.766522e-01
[31] 1.414256e-01 1.304340e-02 3.530937e-01 6.887709e-01 2.385602e-02
[36] 1.322457e-02 2.722920e-01 5.426136e-02 4.218158e-01 1.581127e-02
[41] 8.825610e-04 1.252228e-05 8.743649e-01 3.016908e-01 3.516461e-02
[46] 9.582456e-03 3.877222e-02 3.172920e-03 8.122276e-01 1.950067e-01
Hommel

zkl[edit]

Translation of: C
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),..
o,ro  := order(pvalues,True),order(o,False); # sort pvalues, sort indices
in  := psz.pump(List,'wrap(n){ n_i[n]*pvalues[o[n]] });
pmin  := cummin(in).apply((1.0).min); # (min(1,c[0]),min(1,c[1]),...)
ro.apply(pmin.get); # (pmin[ro[0]],pmin[ro[1]],...)
}
 
fcn by(pvalues){ // Benjamini & Yekutieli
psz,pszf := pvalues.len(), psz.toFloat();
o,ro  := order(pvalues,True),order(o,False); # sort pvalues, sort indices
n_i  := psz.pump(List,'wrap(n){ pszf/(psz - n) }); # N/(N-0),N/(N-1),..
q  := [1..psz].reduce(fcn(q,n){ q+=1.0/n },0.0);
in  := psz.pump(List,'wrap(n){ q * n_i[n] * pvalues[o[n]] });
cummin(in).apply((1.0).min) : ro.apply(_.get);
}
 
fcn hochberg(pvalues){
psz,pszf := pvalues.len(), psz.toFloat();
o,ro  := order(pvalues,True),order(o,False); # sort pvalues, sort indices
n_i  := psz.pump(List,'wrap(n){ pszf/(psz - n) }); # N/(N-0),N/(N-1),..
in  := psz.pump(List,'wrap(n){ pvalues[o[n]]*(n + 1) });
cummin(in).apply((1.0).min) : ro.apply(_.get);
}
 
fcn cummin(pvalues){ // R's cumulative minima --> list of mins
out,m := List.createLong(pvalues.len()), pvalues[0];
foreach pv in (pvalues){ out.append(m=m.min(pv)) }
out
}
fcn order(list,downUp){ // True==increasing, --> List(int) sorted indices
f:=(downUp) and fcn(a,b){ a[1]>b[1] } or fcn(a,b){ a[1]<b[1] };
[0..].zip(list).pump(List()).sort(f).pump(List,T("get",0))
}
 
fcn bonferroni(pvalues){ // -->List
sz,r := pvalues.len(),List();
foreach pv in (pvalues){
b:=pv*sz;
if(b>=1.0) r.append(1.0);
else if(0.0<=b<1.0) r.append(b);
else throw(Exception.ValueError(
"%g is outside of the interval I planned.".fmt(b)));
}
r
}
 
fcn hommel(pvalues){
psz,indices := pvalues.len(), [1..psz].walk(); // 1,2,3,4...
o,ro  := order(pvalues,False),order(o,False); # sort pvalues, sort indices
p  := o.apply('wrap(n){ pvalues[n] }).copy(); // pvalues[*o]
npi  := [1..].zip(p).apply('wrap([(n,p)]){ p*psz/n });
min  := (0.0).min(npi); // min value in npi
pa,q  := List.createLong(psz,min), pa.copy(); #(min,min,,,)
foreach j in ([psz - 1..2,-1]){
ij:=[0..psz - j].walk();
i2:=(j - 1).pump(List,'+(psz - j + 1));
q1:=(0.0).min((j-1).pump(List,'wrap(n){ p[i2[n]]*j/(2 + n) }));
foreach i in (psz - j + 1){ q[ij[i]] = q1.min(p[ij[i]]*j) }
foreach i in (j - 1){ q[i2[i]] = q[psz - j] }
foreach i in (psz){ pa[i] = pa[i].max(q[i]) }
}
psz.pump(List,'wrap(n){ pa[ro[n]] }); // Hommel q-values
}
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,
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);
 
bh(pvalues)  : format(_,"\nBenjamini-Hochberg");
by(pvalues)  : format(_,"\nBenjamini & Yekutieli");
bonferroni(pvalues) : format(_,"\nBonferroni");
hochberg(pvalues)  : format(_,"\nHochberg");
hommel(pvalues)  : format(_,"\nHommel");
 
fcn format(list,title){
print(title,":");
foreach n in ([1..list.len(),5]){
print("\n[%2d]:".fmt(n));
foreach x in (list[n-1,5]){ print(" %.6e".fmt(x)) }
}
println();
}
Output:
Benjamini-Hochberg:
[ 1]: 6.126681e-01 8.521710e-01 1.987205e-01 1.891595e-01 3.217789e-01
[ 6]: 9.301450e-01 4.870370e-01 9.301450e-01 6.049731e-01 6.826753e-01
[11]: 6.482629e-01 7.253722e-01 5.280973e-01 8.769926e-01 4.705703e-01
[16]: 9.241867e-01 6.049731e-01 7.856107e-01 4.887526e-01 1.136717e-01
[21]: 4.991891e-01 8.769926e-01 9.991834e-01 3.217789e-01 9.301450e-01
[26]: 2.304958e-01 5.832475e-01 3.899547e-02 8.521710e-01 1.476843e-01
[31]: 1.683638e-02 2.562902e-03 3.516084e-02 6.250189e-02 3.636589e-03
[36]: 2.562902e-03 2.946883e-02 6.166064e-03 3.899547e-02 2.688991e-03
[41]: 4.502862e-04 1.252228e-05 7.881555e-02 3.142613e-02 4.846527e-03
[46]: 2.562902e-03 4.846527e-03 1.101708e-03 7.252032e-02 2.205958e-02

Benjamini & Yekutieli:
[ 1]: 1.000000e+00 1.000000e+00 8.940844e-01 8.510676e-01 1.000000e+00
[ 6]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[11]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[16]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 5.114323e-01
[21]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[26]: 1.000000e+00 1.000000e+00 1.754486e-01 1.000000e+00 6.644618e-01
[31]: 7.575031e-02 1.153102e-02 1.581959e-01 2.812089e-01 1.636176e-02
[36]: 1.153102e-02 1.325863e-01 2.774239e-02 1.754486e-01 1.209832e-02
[41]: 2.025930e-03 5.634031e-05 3.546073e-01 1.413926e-01 2.180552e-02
[46]: 1.153102e-02 2.180552e-02 4.956812e-03 3.262838e-01 9.925057e-02

Bonferroni:
[ 1]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[ 6]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[11]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[16]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[21]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[26]: 1.000000e+00 1.000000e+00 7.019185e-01 1.000000e+00 1.000000e+00
[31]: 2.020365e-01 1.516674e-02 5.625735e-01 1.000000e+00 2.909271e-02
[36]: 1.537741e-02 4.125636e-01 6.782670e-02 6.803480e-01 1.882294e-02
[41]: 9.005725e-04 1.252228e-05 1.000000e+00 4.713920e-01 4.395577e-02
[46]: 1.088915e-02 4.846527e-02 3.305125e-03 1.000000e+00 2.867745e-01

Hochberg:
[ 1]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[ 6]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[11]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[16]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[21]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[26]: 9.991834e-01 9.991834e-01 4.632662e-01 9.991834e-01 9.991834e-01
[31]: 1.575885e-01 1.383967e-02 3.938014e-01 7.600230e-01 2.501973e-02
[36]: 1.383967e-02 3.052971e-01 5.426136e-02 4.626366e-01 1.656419e-02
[41]: 8.825610e-04 1.252228e-05 9.930759e-01 3.394022e-01 3.692284e-02
[46]: 1.023581e-02 3.974152e-02 3.172920e-03 8.992520e-01 2.179486e-01

Hommel:
[ 1]: 9.991834e-01 9.991834e-01 9.991834e-01 9.987624e-01 9.991834e-01
[ 6]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[11]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[16]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.595180e-01
[21]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01
[26]: 9.991834e-01 9.991834e-01 4.351895e-01 9.991834e-01 9.766522e-01
[31]: 1.414256e-01 1.304340e-02 3.530937e-01 6.887709e-01 2.385602e-02
[36]: 1.322457e-02 2.722920e-01 5.426136e-02 4.218158e-01 1.581127e-02
[41]: 8.825610e-04 1.252228e-05 8.743649e-01 3.016908e-01 3.516461e-02
[46]: 9.582456e-03 3.877222e-02 3.172920e-03 8.122276e-01 1.950067e-01