Berechnung der Wilcoxon Wahrscheinlichkeiten

3 Minuten zum Lesen

Aufzählung der Kombinationen

Die Häufigkeiten der Wilcoxon-Rangsummen sind Grundlage für statistische Tests, beispielsweise zur Lage des Medians (Vorzeichen-Test, Rangsummen-Test). Ein Rang ist hierbei die Position einer Zahl in einer geordneten Liste. Für die Wilcoxon-Rangsummen sind die Positionen positiver Zahlen (aus Differenzen) relevant.

Eine brute force-Möglichkeit, um zu den Häufigkeiten der Rangsummen zu gelangen, ist die Auflistung aller möglichen Listen. Die erste Möglichkeit ist, dass keine positive Zahl in der Liste enthalten ist. Wir betrachten eine Liste der Länge 5.

(-, -, -, -, -)

Die entsprechende Rangsumme hat den Wert 0.

Die Auflistung aller möglichen Listen entspricht allen binären Darstellungen der Zahlen zwischen 0 und 2 hoch minus 1.

In der folgenden Funktion wird eine Zahl hochgezählt, die zugeordnete Rangsumme bestimmt und mit rscount die Häufigkeiten festgestellt. Code-Beispiel Julia:

function hwilcoxoncount(n::Int)
    combination = UInt(0)
    ncomb = 2^n
    if ncomb > typemax(combination)
        throw(ErrorException("n is too big"))
    end
    rscount = zeros(Int, Int(n * (n + 1)/2 + 1))
    function is1(c::Char)::Bool
        c == '1'
    end
    while combination < ncomb
        val = sum(findall(is1, string(combination, base=2, pad=n)))
        rscount[val+1] += 1
        combination += 1
    end
    rscount
end

Ergebnis

Da die Zahl der möglichen Kombination sehr schnell wächst, muss man schon für je nach Hardware einige Sekunden auf das Ergebnis warten.

Übergang von zu

Für eine effizientere Berechnung betrachtet man den Übergang von zu . Welche neuen Möglichkeiten ergeben sich für eine Summe durch die zusätzliche Stelle? Über die Möglichkeiten wird nach der Zahl der beteiligten Summanden Buch geführt. Für jede Summe wird einer Zahl von Summanden die Zahl der Kombinationen zugeordnet, welche die Summe ergeben. Für ergibt sich:

Rangsumme Zuordnungen
0 Dict(0=>1.0)
1 Dict(1=>1.0)
2 Dict(1=>1.0)
3 Dict(2=>1.0,1=>1.0)
4 Dict(2=>1.0,1=>1.0)
5 Dict(2=>2.0)
6 Dict(2=>1.0,3=>1.0)
7 Dict(2=>1.0,3=>1.0)
8 Dict(3=>1.0)
9 Dict(3=>1.0)
10 Dict(4=>1.0)

Lesebeispiel: Für die Zahl sieben gibt es für zwei Summanden eine Kombination und eine weitere Kombination für drei Summanden.

Beispiel: Übergang von zu

Die Zahl der Möglichkeiten zur Bildung von Rangsummen zwischen 0 und 4 ändert sich nicht.

Für die Rangsumme 5 bestehen für die Kombinationen:

(+, -, -, +)
(-, +, +, -)

Diese bestehen auch weiterhin als:

(+, -, -, +, -)
(-, +, +, -, -)

Es kommt folgende Möglichkeit hinzu:

(-, -, -, -, +)

Dabei entspricht die Kombination (-, -, -, -) den Möglichkeiten für den Wert 0 für . Aus der Kombination für den Wert 0 mit 0 beteiligten Summanden wird eine Kombination mit dem Wert 5 mit einem beteiligten Summanden.

Aus der Zuordnung von der Zahl der Summanden zur Zahl der Möglichkeiten für und

2 => 2

wird für und

1 => 1, 2 => 2

Für die Rangsumme sind neben den bestehenden Kombinationen für , die schon 6 ergeben, die Kombinationen für die Rangsumme 1 zu betrachten, die zusammen mit einem + an der neuen fünften Stelle den Zielwert 6 ergeben. Es gibt nur eine Kombination für die Rangsumme 1, die einen Summanden hat. Daraus wird nun eine Kombination mit 2 Summanden.

Aus der Zuordnung für und

2 => 1, 3 => 1

wird für und

2 => 2, 3 => 1

Code-Beispiel Julia:

function cwilcoxonr(n::Int)
    @assert(n>0)
    if n == 1
        return res = [DictIntFloat(0 => 1),DictIntFloat(1 => 1)]
    end
    res = cwilcoxonr(n-1)
    len = length(res)
    lennext = Int(n * (n+1) / 2 + 1)
    resnext = Vector{DictIntFloat}(undef, lennext)
    function incrkey(p::eltype(DictIntFloat))
        typeof(p)(p[1] + 1, p[2])
    end
    resnext[1:n] = res[1:n]
    for i = n+1:len
        dadd = DictIntFloat(incrkey(p) for p = pairs(res[i-n]))
        resnext[i] = add(res[i], dadd)
    end
    for i = len+1:lennext
        resnext[i] = DictIntFloat(incrkey(p) for p = pairs(res[i-n]))
    end
    return resnext
end

Ausgabe für mit den entsprechenden Werten für die Summe:

0	Dict(0=>1.0)
1	Dict(1=>1.0)
2	Dict(1=>1.0)
3	Dict(2=>1.0,1=>1.0)
4	Dict(2=>1.0,1=>1.0)
5	Dict(2=>2.0)
6	Dict(2=>1.0,3=>1.0)
7	Dict(2=>1.0,3=>1.0)
8	Dict(3=>1.0)
9	Dict(3=>1.0)
10	Dict(4=>1.0)

Ausgabe für

0	Dict(0=>1.0)
1	Dict(1=>1.0)
2	Dict(1=>1.0)
3	Dict(2=>1.0,1=>1.0)
4	Dict(2=>1.0,1=>1.0)
5	Dict(2=>2.0,1=>1.0)
6	Dict(2=>2.0,3=>1.0)
7	Dict(2=>2.0,3=>1.0)
8	Dict(2=>1.0,3=>2.0)
9	Dict(2=>1.0,3=>2.0)
10	Dict(4=>1.0,3=>2.0)
11	Dict(4=>1.0,3=>1.0)
12	Dict(4=>1.0,3=>1.0)
13	Dict(4=>1.0)
14	Dict(4=>1.0)
15	Dict(5=>1.0)
julia> VERSION
v"1.1.0"

Auf diese Weise lassen sich die Häufigkeiten etwa 150-mal schneller berechnen als über die Auflistung aller Möglichkeiten (Test für ). Die Verteilungen für bis 200 erhält man innerhalb weniger Sekunden.

Links zu den Ergebnistabellen: Vorzeichen-Test, Rangsummen-Test.