302 lines
7.8 KiB
Plaintext
302 lines
7.8 KiB
Plaintext
#include "string.bi"
|
|
|
|
#define MIN(a, b) iif((a) < (b), (a), (b))
|
|
|
|
Enum Direction
|
|
UP
|
|
DOWN
|
|
End Enum
|
|
|
|
Type Sequence
|
|
As Integer length
|
|
As Double values(Any)
|
|
End Type
|
|
|
|
' Constants for correction types
|
|
Dim As String correctionTypes(7) = { _
|
|
"Benjamini-Hochberg", "Benjamini-Yekutieli", "Bonferroni", "Hochberg", _
|
|
"Holm", "Hommel", "Sidak", "Unknown" }
|
|
|
|
' Test p-values
|
|
Dim Shared As Double test_values(49) = { _
|
|
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 }
|
|
|
|
Function minimum(p As Sequence) As Double
|
|
Dim m As Double = p.values(0)
|
|
For i As Integer = 1 To p.length - 1
|
|
If p.values(i) < m Then m = p.values(i)
|
|
Next
|
|
Return m
|
|
End Function
|
|
|
|
Function maximum(p As Sequence) As Double
|
|
Dim m As Double = p.values(0)
|
|
For i As Integer = 1 To p.length - 1
|
|
If p.values(i) > m Then m = p.values(i)
|
|
Next
|
|
Return m
|
|
End Function
|
|
|
|
Sub ratchet(p As Sequence, direcc As Integer)
|
|
Dim m As Double = p.values(0)
|
|
Dim i As Integer
|
|
|
|
If direcc = UP Then
|
|
For i = 1 To p.length - 1
|
|
' Corrected logic: only update if greater than minimum
|
|
If p.values(i) > m Then p.values(i) = m
|
|
m = p.values(i)
|
|
Next
|
|
Else
|
|
For i = 1 To p.length - 1
|
|
If p.values(i) < m Then p.values(i) = m
|
|
m = p.values(i)
|
|
Next
|
|
End If
|
|
|
|
' Cap at 1.0
|
|
For i = 0 To p.length - 1
|
|
If p.values(i) > 1.0 Then p.values(i) = 1.0
|
|
Next
|
|
End Sub
|
|
|
|
Function schwartzian(p As Sequence, mult As Sequence, direcc As Integer) As Sequence
|
|
Dim As Sequence result
|
|
result.length = p.length
|
|
Redim result.values(p.length-1)
|
|
Dim As Integer i, j
|
|
|
|
' Sort with indices
|
|
Dim As Integer indices(p.length-1)
|
|
For i = 0 To p.length-1
|
|
indices(i) = i
|
|
Next
|
|
|
|
' Sort based on direccection
|
|
For i = 0 To p.length-2
|
|
For j = 0 To p.length-2-i
|
|
Dim As Boolean cond
|
|
If direcc = UP Then
|
|
cond = p.values(indices(j)) < p.values(indices(j+1))
|
|
Else
|
|
cond = p.values(indices(j)) > p.values(indices(j+1))
|
|
End If
|
|
If cond Then Swap indices(j), indices(j+1)
|
|
Next
|
|
Next
|
|
|
|
' Apply multipliers
|
|
For i = 0 To p.length-1
|
|
result.values(i) = p.values(indices(i)) * mult.values(i)
|
|
Next
|
|
|
|
ratchet(result, direcc)
|
|
|
|
' Restore original order
|
|
Dim As Sequence final
|
|
final.length = p.length
|
|
Redim final.values(p.length-1)
|
|
|
|
For i = 0 To p.length-1
|
|
final.values(indices(i)) = result.values(i)
|
|
Next
|
|
|
|
Return final
|
|
End Function
|
|
|
|
Function adjust(p As Sequence, method As String) As Sequence
|
|
Dim As Sequence result
|
|
result.length = p.length
|
|
Redim result.values(p.length-1)
|
|
|
|
Dim As Sequence mult
|
|
mult.length = p.length
|
|
Redim mult.values(p.length-1)
|
|
|
|
Dim As Integer i
|
|
Dim As Double tmp
|
|
|
|
Select Case method
|
|
Case "Benjamini-Hochberg"
|
|
For i = 0 To p.length-1
|
|
mult.values(i) = p.length / (p.length - i)
|
|
Next
|
|
Return schwartzian(p, mult, UP)
|
|
|
|
Case "Benjamini-Yekutieli"
|
|
tmp = 0
|
|
For i = 1 To p.length
|
|
tmp += 1.0 / i
|
|
Next
|
|
|
|
For i = 0 To p.length-1
|
|
mult.values(i) = tmp * p.length / (p.length - i)
|
|
Next
|
|
Return schwartzian(p, mult, UP)
|
|
|
|
Case "Bonferroni"
|
|
For i = 0 To p.length-1
|
|
result.values(i) = Min(p.values(i) * p.length, 1.0)
|
|
Next
|
|
Return result
|
|
|
|
Case "Hochberg"
|
|
For i = 0 To p.length-1
|
|
mult.values(i) = i + 1
|
|
Next
|
|
Return schwartzian(p, mult, UP)
|
|
|
|
Case "Holm"
|
|
For i = 0 To p.length-1
|
|
mult.values(i) = p.length - i
|
|
Next
|
|
Return schwartzian(p, mult, DOWN)
|
|
|
|
Case "Hommel"
|
|
Dim As Integer order(p.length-1), j
|
|
Dim As Double s(p.length-1)
|
|
|
|
' Sort and get order
|
|
For i = 0 To p.length-1
|
|
order(i) = i
|
|
Next
|
|
For i = 0 To p.length-2
|
|
For j = 0 To p.length-2-i
|
|
If p.values(order(j)) > p.values(order(j+1)) Then Swap order(j), order(j+1)
|
|
Next
|
|
Next
|
|
|
|
' Get sorted values
|
|
For i = 0 To p.length-1
|
|
s(i) = p.values(order(i))
|
|
Next
|
|
|
|
' Calculate initial minimum
|
|
Dim As Double m(p.length-1)
|
|
For i = 0 To p.length-1
|
|
m(i) = s(i) * p.length / (i + 1)
|
|
Next
|
|
|
|
Dim As Double min_val = m(0)
|
|
For i = 1 To p.length-1
|
|
If m(i) < min_val Then min_val = m(i)
|
|
Next
|
|
|
|
Dim As Double pa(p.length-1), q(p.length-1)
|
|
For i = 0 To p.length-1
|
|
pa(i) = min_val
|
|
q(i) = min_val
|
|
Next
|
|
|
|
' Hommel algorithm
|
|
For j = p.length-1 To 2 Step -1
|
|
Dim As Integer lower_count = p.length - j + 1
|
|
Dim As Integer upper_count = j - 1
|
|
|
|
' Initialize qmin with first upper value
|
|
Dim As Double qmin = j * s(p.length - j + 1) / 2.0
|
|
|
|
' Check remaining upper values
|
|
For i = 1 To upper_count-1
|
|
Dim As Double tmp = s(p.length - j + 1 + i) * j / (2.0 + i)
|
|
qmin = MIN(qmin, tmp)
|
|
Next
|
|
|
|
' Update lower values
|
|
For i = 0 To lower_count-1
|
|
q(i) = MIN(s(i) * j, qmin)
|
|
Next
|
|
|
|
' Update upper values
|
|
For i = lower_count To p.length-1
|
|
q(i) = q(p.length - j)
|
|
Next
|
|
|
|
' Update pa with maximum between current and new values
|
|
For i = 0 To p.length-1
|
|
If q(i) > pa(i) Then pa(i) = q(i)
|
|
Next
|
|
Next
|
|
|
|
' Map back to original order
|
|
For i = 0 To p.length-1
|
|
result.values(order(i)) = MIN(pa(i), 1.0)
|
|
Next
|
|
Return result
|
|
|
|
Case "Šidák", "Sidak"
|
|
For i = 0 To p.length-1
|
|
result.values(i) = 1.0 - (1.0 - p.values(i)) ^ p.length
|
|
Next
|
|
Return result
|
|
|
|
Case Else
|
|
result.length = 0
|
|
End Select
|
|
|
|
Return result
|
|
End Function
|
|
|
|
Function formatOutput(values As Sequence) As String
|
|
Dim As String result = ""
|
|
Dim As Integer i, j
|
|
|
|
For i = 0 To values.length -1 Step 5
|
|
result &= "[" & Format(i, "00") & "] "
|
|
For j = 0 To 4
|
|
If i + j < values.length Then
|
|
result &= Format(values.values(i+j), "0.0000000000") & " "
|
|
End If
|
|
Next
|
|
result &= !"\n"
|
|
Next
|
|
Return result
|
|
End Function
|
|
|
|
'Main program
|
|
Dim As Sequence p
|
|
p.length = 50
|
|
Redim p.values(49)
|
|
Dim As Integer i
|
|
|
|
' Initialize p-values
|
|
For i = 0 To 49
|
|
p.values(i) = test_values(i)
|
|
Next
|
|
|
|
' Check p-values
|
|
If p.length = 0 Orelse minimum(p) < 0 Orelse maximum(p) > 1 Then
|
|
Print "p-values must be in range 0.0 to 1.0"
|
|
End
|
|
End If
|
|
|
|
' Apply each correction method
|
|
For i = 0 To 7
|
|
Dim As String method = correctionTypes(i)
|
|
|
|
Print method
|
|
|
|
Dim As Sequence result = adjust(p, method)
|
|
|
|
If result.length > 0 Then
|
|
Print formatOutput(result)
|
|
Else
|
|
Print "Sorry, do not know how to do '" & method & "' correction."
|
|
Print "Perhaps you want one of these?:"
|
|
For j As Integer = 0 To 6
|
|
Print " " & correctionTypes(j)
|
|
Next
|
|
End If
|
|
Next
|
|
|
|
Sleep
|