P-value correction

From Rosetta Code
Revision as of 00:56, 5 November 2017 by Thundergnat (talk | contribs) (→‎{{header|Perl 6}}: Simplify and comment Hommel routine)
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

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 as possible, and is able to do any one of the methods.

Link with -lm <lang C>

  1. include <stdio.h>//printf
  2. include <stdlib.h>//qsort
  3. include <math.h>//fabs
  4. include <stdbool.h>//bool data type
  5. 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 double_array_min(const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { 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); } double min = ARRAY[0]; for (size_t index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) { if (ARRAY[index] < min) { min = ARRAY[index]; } } return min; }

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 types.\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) unsigned int *restrict indices = seq_len(NO_OF_ARRAY_ELEMENTS, NO_OF_ARRAY_ELEMENTS); //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 *restrict npi = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (npi == 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++) { npi[index] = (double)NO_OF_ARRAY_ELEMENTS * p[index] / indices[index]; } free(indices); indices = NULL; const double min = double_array_min(npi, NO_OF_ARRAY_ELEMENTS); free(npi); npi = NULL; 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 }

const size_t PI2_LENGTH = I2_LENGTH;

double q1 = j * p[i2[0]] / 2.0; for (size_t i = 1; i < PI2_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 double *restrict n_i = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);//R's n/i if (n_i == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } 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); } n_i[index] = (double)NO_OF_ARRAY_ELEMENTS / (NO_OF_ARRAY_ELEMENTS - index);// n/i simplified }

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++) { cummin_input[index] = n_i[index] * PVALUES[o[index]];//PVALUES[o[index]] is p[o] } } else if (TYPE == 1) {//BY method double q = 0.0; for (size_t index = 1; index < (1+NO_OF_ARRAY_ELEMENTS); index++) { q += (double) 1.0/index; } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { cummin_input[index] = q * n_i[index] * 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; free(n_i); n_i = 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; } </lang>

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

Translation of: C

<lang scala>// 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))
   }

}</lang>

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 6

Works with: Rakudo version 2017.10

<lang perl6>########################### Helper subs ########################### sub correct (@p, $type) { 

   say "\n$type";
   pretty-print adjust( @p, $type );

}

sub pretty-print ( @p ) { 

   my $cols = 5;
   my $i = -$cols;
   my $fmt = "%6e ";
   printf "[%2d]  {$fmt xx $cols}\n", $i+=5, $_ for @p.rotor($cols);

}

                          1. The various p-value correction routines #############
  1. Factor out common code; perform a Schwartzian transform

sub schwartzian ( @p, &transform ) { 

   @p.kv.map({[$^v, $^k]}).sort(-*.[0]).values.map(&transform).sort(*.[1]).values»[0]

}

multi adjust( @p, 'Benjamini-Hochberg' ) { 

   schwartzian @p, { [$_.[0] * @p / (@p - $++) min 1, $_.[1]] }

}

multi adjust( @p, 'Benjamini-Yekutieli' ) { 

   schwartzian @p, { state $q=(1..@p).map({1/$_}).sum; [($q*$_.[0]*@p/(@p-$++)) min 1, $_.[1]] }

}

multi adjust( @p, 'Hochberg' ) { 

   schwartzian @p, { state $m = @p.max; [$_.[0] * (1+$++) min $m, $_.[1]] }

}

multi adjust( @p, 'Holm' ) { 

   schwartzian @p, { [$_.[0] * (@p - $_.[1]) min 1, $_.[1]] }

}

  1. Order unimportant for Bonferroni & Sidak corrections

multi adjust( @p, 'Bonferroni' ) { 

   @p.map: { $_ * @p min 1 }

}

multi adjust( @p, 'Sidak' ) { 

   @p.map: { 1 - (1 - $_) }

}

  1. Hommel correction can't be easily reduced to a one pass transform

multi adjust( @p, 'Hommel' ) { 

   my @s = @p.kv.map( { [$^v, $^k] } ).sort( *.[0] ).values; # sorted
   my \z = +@p; # elements
   my @q;       # scratch array
   my @pa = @s»[0].kv.map( {$^v * z / ($^k+1)} ).min xx z; # p adjusted
   for @p - 1, {$_-1} ... 1 -> $i {
       my @i1  = (^z)[^(z - $i + 1)     ]; # lower indices
       my @i2  = (^z)[  z - $i + 1 .. * ]; # upper indices
       my $q   = @s[@i2]»[0].kv.map( { $^v * $i / (2 + $^k) } ).min;
       @q[@i1] = @s[@i1]»[0].map: { [min] $_ * $i, $q, @s[*-1][0] };
       @pa     = (^z).map: { [max] @pa[$_], @q[$_] }
   }
   @pa[@s[*;1].kv.map( {[$^v, $^k]} ).sort( *.[0] ).values»[1]]

}

                                                      1. 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

correct @p-values, 'Benjamini-Hochberg'; correct @p-values, 'Benjamini-Yekutieli'; correct @p-values, 'Bonferroni'; correct @p-values, 'Hochberg'; correct @p-values, 'Holm'; correct @p-values, 'Hommel'; correct @p-values, 'Sidak';</lang>

Output:
Benjamini-Hochberg
[ 0]  6.126681e-01  8.685743e-01  1.987205e-01  1.891595e-01  3.217789e-01 
[ 5]  9.310912e-01  4.870370e-01  9.301450e-01  6.049731e-01  6.826753e-01 
[10]  6.482629e-01  7.253723e-01  5.280973e-01  8.958102e-01  4.705703e-01 
[15]  9.241867e-01  6.097340e-01  7.856107e-01  4.887526e-01  1.136717e-01 
[20]  4.991891e-01  8.769926e-01  9.991834e-01  3.232761e-01  9.414079e-01 
[25]  2.304958e-01  5.832475e-01  3.899547e-02  8.521710e-01  1.476843e-01 
[30]  1.683638e-02  3.033349e-03  3.516084e-02  6.250189e-02  3.636589e-03 
[35]  2.562902e-03  2.946883e-02  6.166064e-03  4.002047e-02  2.688991e-03 
[40]  4.502863e-04  1.252228e-05  7.881555e-02  3.142613e-02  4.883974e-03 
[45]  2.722289e-03  4.846527e-03  1.101708e-03  7.252033e-02  2.205958e-02 

Benjamini-Yekutieli
[ 0]  1.000000e+00  1.000000e+00  8.940844e-01  8.510676e-01  1.000000e+00 
[ 5]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[10]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[15]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  5.114323e-01 
[20]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[25]  1.000000e+00  1.000000e+00  1.754486e-01  1.000000e+00  6.644618e-01 
[30]  7.575031e-02  1.364766e-02  1.581959e-01  2.812089e-01  1.636176e-02 
[35]  1.153102e-02  1.325863e-01  2.774239e-02  1.800603e-01  1.209832e-02 
[40]  2.025930e-03  5.634031e-05  3.546073e-01  1.413926e-01  2.197400e-02 
[45]  1.224814e-02  2.180552e-02  4.956812e-03  3.262838e-01  9.925057e-02 

Bonferroni
[ 0]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[ 5]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[10]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[15]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[20]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[25]  1.000000e+00  1.000000e+00  7.019185e-01  1.000000e+00  1.000000e+00 
[30]  2.020365e-01  1.516675e-02  5.625735e-01  1.000000e+00  2.909271e-02 
[35]  1.537741e-02  4.125636e-01  6.782670e-02  6.803480e-01  1.882294e-02 
[40]  9.005725e-04  1.252228e-05  1.000000e+00  4.713920e-01  4.395577e-02 
[45]  1.088916e-02  4.846527e-02  3.305125e-03  1.000000e+00  2.867745e-01 

Hochberg
[ 0]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[ 5]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[10]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[15]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[20]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[25]  9.991834e-01  9.991834e-01  4.632662e-01  9.991834e-01  9.991834e-01 
[30]  1.575885e-01  1.395341e-02  3.938015e-01  7.600230e-01  2.501973e-02 
[35]  1.383967e-02  3.052971e-01  5.426136e-02  4.626366e-01  1.656419e-02 
[40]  8.825611e-04  1.252228e-05  9.930759e-01  3.394022e-01  3.692284e-02 
[45]  1.023581e-02  3.974152e-02  3.172920e-03  8.992520e-01  2.179486e-01 

Holm
[ 0]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[ 5]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[10]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[15]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[20]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 
[25]  1.000000e+00  1.000000e+00  3.228825e-01  1.000000e+00  1.000000e+00 
[30]  8.081460e-02  5.763363e-03  2.025265e-01  4.037622e-01  9.309667e-03 
[35]  4.613223e-03  1.155178e-01  1.763494e-02  1.632835e-01  4.141047e-03 
[40]  1.801145e-04  2.254010e-06  2.648202e-01  6.599487e-02  5.274692e-03 
[45]  1.088916e-03  3.877222e-03  1.983075e-04  5.801626e-02  5.735490e-03 

Hommel
[ 0]  9.991834e-01  9.991834e-01  9.991834e-01  9.987624e-01  9.991834e-01 
[ 5]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[10]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[15]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.595180e-01 
[20]  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01  9.991834e-01 
[25]  9.991834e-01  9.991834e-01  4.351895e-01  9.991834e-01  9.766523e-01 
[30]  1.414256e-01  1.304340e-02  3.530937e-01  6.887709e-01  2.385602e-02 
[35]  1.322457e-02  2.722920e-01  5.426136e-02  4.218158e-01  1.581127e-02 
[40]  8.825611e-04  1.252228e-05  8.743649e-01  3.016908e-01  3.516461e-02 
[45]  9.582456e-03  3.877222e-02  3.172920e-03  8.122276e-01  1.950067e-01 

Sidak
[ 0]  4.533744e-01  7.296024e-01  9.936026e-02  9.079658e-02  1.801962e-01 
[ 5]  8.752257e-01  2.922222e-01  9.115421e-01  4.355806e-01  5.324867e-01 
[10]  4.926798e-01  5.802978e-01  3.485442e-01  7.883130e-01  2.729308e-01 
[15]  8.502518e-01  4.268138e-01  6.442008e-01  3.030266e-01  5.001555e-02 
[20]  3.194810e-01  7.892933e-01  9.991834e-01  1.745691e-01  9.037516e-01 
[25]  1.198578e-01  3.966083e-01  1.403837e-02  7.328671e-01  6.793476e-02 
[30]  4.040730e-03  3.033349e-04  1.125147e-02  2.375072e-02  5.818542e-04 
[35]  3.075482e-04  8.251272e-03  1.356534e-03  1.360696e-02  3.764588e-04 
[40]  1.801145e-05  2.504456e-07  3.310253e-02  9.427839e-03  8.791153e-04 
[45]  2.177831e-04  9.693054e-04  6.610250e-05  2.900813e-02  5.735490e-03

R

<lang R>

  1. p.adjust's source code is below
  2. > p.adjust
  3. function (p, method = p.adjust.methods, n = length(p))
  4. {
  5. method <- match.arg(method)
  6. if (method == "fdr")
  7. method <- "BH"
  8. nm <- names(p)
  9. p <- as.numeric(p)
  10. p0 <- setNames(p, nm)
  11. if (all(nna <- !is.na(p)))
  12. nna <- TRUE
  13. p <- p[nna]
  14. lp <- length(p)
  15. stopifnot(n >= lp)
#   if (n <= 1) 
#       return(p0)
  1. if (n == 2 && method == "hommel")
  2. method <- "hochberg"
  3. p0[nna] <- switch(method, bonferroni = pmin(1, n * p), holm = {
  4. i <- seq_len(lp)
#       o <- order(p)
  1. ro <- order(o)
  2. pmin(1, cummax((n - i + 1L) * p[o]))[ro]
  3. }, hommel = {
#       if (n > lp) p <- c(p, rep.int(1, n - lp))
  1. i <- seq_len(n)
#       o <- order(p)
  1. p <- p[o]
  2. ro <- order(o)
  3. q <- pa <- rep.int(min(n * p/i), n)
    1. for (j in (n - 1):2) {
#           ij <- seq_len(n - j + 1)
  1. i2 <- (n - j + 2):n
  2. q1 <- min(j * p[i2]/(2:j))
  3. q[ij] <- pmin(j * p[ij], q1)
  4. q[i2] <- q[n - j + 1]
  5. pa <- pmax(pa, q)
  6. }
  7. pmax(pa, p)[if (lp < n) ro[1:lp] else ro]
  8. }, hochberg = {
#       i <- lp:1L
  1. o <- order(p, decreasing = TRUE)
  2. ro <- order(o)
  3. pmin(1, cummin((n - i + 1L) * p[o]))[ro]
  4. }, BH = {
  5. i <- lp:1L
  6. o <- order(p, decreasing = TRUE)
  7. ro <- order(o)
  8. pmin(1, cummin(n/i * p[o]))[ro]
  9. }, BY = {
  10. i <- lp:1L
  11. o <- order(p, decreasing = TRUE)
  12. ro <- order(o)
  13. q <- sum(1L/(1L:n))
#       pmin(1, cummin(q * n/i * p[o]))[ro]
  1. }, none = p)
  2. p0
  3. }
  4. <bytecode: 0x3a61b88>
  5. <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") </lang>

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

Translation of: C

<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),..
  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

}</lang> <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,
  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();

}</lang>

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
Retrieved from "https://rosettacode.org/wiki/P-value_correction?oldid=278814"