(phixonline)-->
with javascript_semantics
enum UP, DOWN
function ratchet(sequence p, integer direction)
atom m = p[1]
for i=1 to length(p) do
if iff(direction=UP?p[i]>m:p[i]<m) then p[i] = m end if
m = p[i]
end for
return sq_min(p,1)
end function
function schwartzian(sequence p, mult, integer direction)
sequence order = custom_sort(p,tagset(length(p)))
if direction=UP then order = reverse(order) end if
sequence pa = ratchet(sq_mul(mult,extract(p,order)), direction)
return extract(pa,order,invert:=true)
end function
function adjust(sequence p, string method)
integer size = length(p)
sequence mult = tagset(size)
switch method
case "Benjamini-Hochberg":
mult = sq_div(size,sq_sub(size+1,mult))
return schwartzian(p, mult, UP)
case "Benjamini-Yekutieli":
atom q = sum(sq_div(1,mult))
mult = sq_div(q*size,sq_sub(size+1,mult))
return schwartzian(p, mult, UP)
case "Bonferroni":
return sq_min(sq_mul(p,size),1)
case "Hochberg":
return schwartzian(p, mult, UP)
case "Holm":
mult = sq_sub(size+1,mult)
return schwartzian(p, mult, DOWN)
case "Hommel":
sequence ivdx = repeat(0,size)
for i=1 to size do ivdx[i] = {p[i],i} end for
ivdx = sort(ivdx)
sequence s = vslice(ivdx,1),
m = sq_div(sq_mul(s,size),mult),
qh = repeat(min(m),size),
pa = repeat(min(m),size),
order = vslice(ivdx,2)
for j=size-1 to 2 by -1 do
sequence lwr = tagset(size-j+1),
upr = sq_add(size-j+1,tagset(j-1))
atom qmin = j*s[upr[1]]/2
for i=2 to length(upr) do
qmin = min(s[upr[i]]*j/(i+1),qmin)
end for
for i=1 to length(lwr) do
qh[lwr[i]] = min(s[lwr[i]]*j, qmin)
end for
for i=1 to length(upr) do
qh[upr[i]] = qh[size-j+1]
end for
pa = sq_max(pa,qh)
end for
return extract(pa,order,invert:=true)
case "Sidak":
p = deep_copy(p)
for i=1 to length(p) do
p[i] = 1 - power(1-p[i],size)
end for
return p
else
return {} -- (unknown method)
end switch
return p
end function
constant {types,correct_answers} = columnize({
{"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}},
{"Sidak",
{1.0000000000, 1.0000000000, 0.9946598274, 0.9914285749, 0.9999515274,
1.0000000000, 0.9999999688, 1.0000000000, 1.0000000000, 1.0000000000,
1.0000000000, 1.0000000000, 0.9999999995, 1.0000000000, 0.9999998801,
1.0000000000, 1.0000000000, 1.0000000000, 0.9999999855, 0.9231179729,
0.9999999956, 1.0000000000, 1.0000000000, 0.9999317605, 1.0000000000,
0.9983109511, 1.0000000000, 0.5068253940, 1.0000000000, 0.9703301333,
0.1832692440, 0.0150545753, 0.4320729669, 0.6993672225, 0.0286818157,
0.0152621104, 0.3391808707, 0.0656206307, 0.4959194266, 0.0186503726,
0.0009001752, 0.0000125222, 0.8142104886, 0.3772612062, 0.0430222116,
0.0108312558, 0.0473319661, 0.0032997780, 0.7705015898, 0.2499384839}}})
-- {"Unknown",{1}}})
constant 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}
if length(pValues)=0 or min(pValues)<0 or max(pValues)>1 then
crash("p-values must be in range 0.0 to 1.0")
end if
for i=1 to length(types) do
string ti = types[i]
sequence res = adjust(pValues,ti)
if res={} then
printf(1,"\nSorry, do not know how to do %s correction.\n"&
"Perhaps you want one of these?:\n %s\n",
{ti,join(types[1..$-1],"\n ")})
exit
end if
-- printf(1,"%s\n",{ti})
-- res = correct_answers[i] -- (for easier comparison only)
-- pp(res,{pp_FltFmt,"%13.10f",pp_IntFmt,"%13.10f",pp_Maxlen,75,pp_Pause,0})
atom error = sum(sq_abs(sq_sub(res,correct_answers[i])))
printf(1,"%s has cumulative error of %g\n", {ti,error})
end for