P-value correction: Difference between revisions
Content deleted Content added
→{{header|Perl}}: removed unused variable and unnecessary checks |
→Version 1: new version does not show compiler warnings about ignored qualifiers or possible conversion errors. Eliminated reindexing step with Hommel method, so code is shorter and should be slightly faster with that method |
||
Line 58:
#include <stdbool.h>//bool data type
#include <strings.h>//strcasecmp
#include <assert.h>//assert, necessary for random integer selection
unsigned int *
//named after R function of same name, but simpler function
if (START == END) {
unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int));
Line 70 ⟶ 71:
exit(EXIT_FAILURE);
}
for (
sequence[i] = i+1;
}
Line 76 ⟶ 77:
}
if (START > END) {
end = (unsigned)START;
start = (unsigned)END;
}
const
unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int));
if (sequence == NULL) {
Line 87 ⟶ 88:
}
if (START < END) {
for (
sequence[index] = start + index;
}
} else {
for (
sequence[index] = end - index;
}
Line 124 ⟶ 125:
}
unsigned int *
//this has the same name as the same R function
unsigned int *restrict idx = malloc(SIZE * sizeof(unsigned int));
Line 151 ⟶ 152:
}
double *
//this takes the same name of the R function which it copies
//this requires a free() afterward where it is used
Line 175 ⟶ 176:
}
double *
//this takes the same name of the R function which it copies
//this requires a free() afterward where it is used
Line 190 ⟶ 191:
}
double cumulative_max = ARRAY[0];
for (
if (ARRAY[i] > cumulative_max) {
cumulative_max = ARRAY[i];
Line 199 ⟶ 200:
}
double *
//named after the R function pmin
if (NO_OF_ARRAY_ELEMENTS < 1) {
Line 224 ⟶ 225:
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) {
printf("[1] %e", ARRAY[0]);
for (
printf(" %.10f", ARRAY[i]);
if (((i+1) % 5) == 0) {
printf("\n[%
}
}
Line 242 ⟶ 243:
}*/
double *
double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (doubleArray == NULL) {
Line 249 ⟶ 250:
exit(EXIT_FAILURE);
}
for (
doubleArray[index] = (double)ARRAY[index];
}
Line 263 ⟶ 264:
}
double *
//this function is a translation of R's p.adjust "BH" method
// i is always i[index] = NO_OF_ARRAY_ELEMENTS - index - 1
Line 272 ⟶ 273:
}
short int TYPE = -1;
if
TYPE = 0;
} else if (strcasecmp(STRING, "BH") == 0) {
TYPE = 0;
} else if (strcasecmp(STRING, "fdr") == 0) {
Line 300 ⟶ 303:
exit(EXIT_FAILURE);
}
for (
const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS;
if (BONFERRONI >= 1.0) {
Line 322 ⟶ 325:
double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS);
double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (
cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * (double)PVALUES[o[index]];
// printf("cummax_input[%zu] = %e\n", index, cummax_input[index]);
}
Line 336 ⟶ 339:
free(cummax_output); cummax_output = NULL;
double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
for (
qvalues[index] = pmin[ro[index]];
}
Line 355 ⟶ 358:
exit(EXIT_FAILURE);
}
for (
p[index] = PVALUES[o[index]];
}
Line 378 ⟶ 381:
}
double min = (double)NO_OF_ARRAY_ELEMENTS * p[0];
for (
const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (double)(1+index);
if (TEMP < min) {
min = TEMP;
}
}
for (
pa[index] = min;
q[index] = min;
Line 399 ⟶ 402:
}
*/
for (
// printf("j = %zu\n", j);
unsigned int *restrict ij = seq_len(
const size_t I2_LENGTH = j - 1;
unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int));
for (
i2[i] = NO_OF_ARRAY_ELEMENTS-j+2+i-1;
//R's indices are 1-based, C's are 0-based, I added the -1
}
double q1 = (double)j * p[i2[0]] / 2.0;
for (
const double TEMP_Q1 = (double)j * p[i2[i]] / (double)(2 + i);
if (TEMP_Q1 < q1) {
q1 = TEMP_Q1;
Line 420:
}
for (
q[ij[i]] = min2( (double)j*p[ij[i]], q1);
}
free(ij); ij = NULL;
for (
q[i2[i]] = q[NO_OF_ARRAY_ELEMENTS - j];//subtract 1 because of starting index difference
}
free(i2); i2 = NULL;
for (
if (pa[i] < q[i]) {
pa[i] = q[i];
Line 439:
}//end j loop
free(p); p = NULL;
for (
q[index] = pa[ro[index]];//Hommel q-values
}
Line 472:
double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);
if (TYPE == 0) {//BH method
for (
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (double)(NO_OF_ARRAY_ELEMENTS - index);// n/i simplified
cummin_input[index] = NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o]
}
} else if (TYPE == 1) {//BY method
double q = 1.0;
for (
q +=
}
for (
const double NI = (double)NO_OF_ARRAY_ELEMENTS / (double)(NO_OF_ARRAY_ELEMENTS - index);// n/i simplified
cummin_input[index] = q * NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o]
}
} else if (TYPE == 3) {//Hochberg method
for (
// pmin(1, cummin((n - i + 1L) * p[o]))[ro]
cummin_input[index] = (double)(index + 1) * PVALUES[o[index]];
}
}
Line 506:
return q_array;
}
int main(void) {
|