Skip to content

Statistical analysis

jmetal.lab.statistical_test implements frequentist and Bayesian tests for comparing algorithms across multiple problems and runs: the Friedman, Friedman aligned-rank, and Quade tests; post-hoc p-value adjustment procedures (Bonferroni-Dunn, Holm, Hochberg, Holland, Finner, Li, Shaffer, Nemenyi); Bayesian sign and signed-rank tests; and critical-distance plots. See Experiments for how to get from a set of algorithm runs to the tidy summary these functions expect as input.

Note

This module may be superseded by SAES, a dedicated statistical-analysis package under the same jMetal organization, once it becomes installable alongside jMetalPy -- currently blocked on a SAES release with a relaxed numpy pin.

API

functions

bonferroni_dunn(p_values, control)

Bonferroni-Dunn's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def bonferroni_dunn(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Bonferroni-Dunn's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        APVs[i, 0] = np.min([(k - 1) * p_values[0, argsorted_pvals[i]], 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Bonferroni"])

holland(p_values, control)

Holland's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def holland(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Holland's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    # --------------------------------------------------------------------------
    # ------------------------------- Procedure --------------------------------
    # --------------------------------------------------------------------------
    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        aux = k - 1 - np.arange(i + 1)
        v = np.max(1 - (1 - p_values[0, argsorted_pvals[: (i + 1)]]) ** aux)
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Holland"])

finner(p_values, control)

Finner's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def finner(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Finner's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        aux = float(k - 1) / (np.arange(i + 1) + 1)
        v = np.max(1 - (1 - p_values[0, argsorted_pvals[: (i + 1)]]) ** aux)
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Finner"])

hochberg(p_values, control)

Hochberg's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def hochberg(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Hochberg's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        aux = np.arange(k, i, -1).astype(np.uint8)
        v = np.max(p_values[0, argsorted_pvals[aux - 1]] * (k - aux))
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Hochberg"])

li(p_values, control)

Li's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm. If provided, control vs all comparisons are considered, else all vs all.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def li(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Li's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm. If provided, control vs all
            comparisons are considered, else all vs all.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        APVs[i, 0] = np.min(
            [
                p_values[0, argsorted_pvals[-2]],
                p_values[0, argsorted_pvals[i]]
                / (p_values[0, argsorted_pvals[i]] + 1 - p_values[0, argsorted_pvals[-2]]),
            ]
        )
    return pd.DataFrame(data=APVs, index=comparison, columns=["Li"])

holm(p_values, control=None)

Holm's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str | None

Index or name of the control algorithm. If provided, control vs all comparisons are considered, else all vs all.

 None

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def holm(p_values: pd.DataFrame | np.ndarray, control: int | str | None = None) -> pd.DataFrame:
    """Holm's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm. If provided, control vs all
            comparisons are considered, else all vs all.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])

    if type(control) == int:
        k = p_values.shape[1]

        # sort p-values p(0) <= p(1) <= ... <= p(k-1)
        argsorted_pvals = np.argsort(p_values[0, :])

        APVs = np.zeros((k - 1, 1))
        comparison = []
        for i in range(k - 1):
            aux = k - 1 - np.arange(i + 1)
            comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
            v = np.max(aux * p_values[0, argsorted_pvals[: (i + 1)]])
            APVs[i, 0] = np.min([v, 1])

    elif control is None:
        k = p_values.shape[1]
        m = int((k * (k - 1)) / 2.0)

        # sort p-values p(0) <= p(1) <= ... <= p(m-1)
        pairs_index = np.triu_indices(k, 1)
        pairs_pvals = p_values[pairs_index]
        pairs_sorted = np.argsort(pairs_pvals)

        APVs = np.zeros((m, 1))
        aux = pairs_pvals[pairs_sorted] * (m - np.arange(m))
        comparison = []
        for i in range(m):
            row = pairs_index[0][pairs_sorted[i]]
            col = pairs_index[1][pairs_sorted[i]]
            comparison.append(algorithms[row] + " vs " + algorithms[col])
            v = np.max(aux[: i + 1])
            APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Holm"])

shaffer(p_values)

Shaffer's procedure for adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
def shaffer(p_values: pd.DataFrame | np.ndarray) -> pd.DataFrame:
    """Shaffer's procedure for adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    def S(k):
        """Computes the set of possible numbers of true hypotheses.

        Args:
            k: Number of algorithms being compared.

        Returns:
            Set of true hypotheses.
        """

        from scipy.special import binom as binomial

        TrueHset = [0]
        if k > 1:
            for j in np.arange(k, 0, -1, dtype=int):
                TrueHset = list(set(TrueHset) | set([binomial(j, 2) + x for x in S(k - j)]))
        return TrueHset

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if p_values.ndim != 2 or p_values.shape[0] != p_values.shape[1]:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")

    # define parameters
    k = p_values.shape[0]
    m = int(k * (k - 1) / 2.0)
    s = np.array(S(k)[1:])

    # sort p-values p(0) <= p(1) <= ... <= p(m-1)
    pairs_index = np.triu_indices(k, 1)
    pairs_pvals = p_values[pairs_index]
    pairs_sorted = np.argsort(pairs_pvals)

    # compute ti: max number of hypotheses that can be true given that any
    # (i-1) hypotheses are false.
    t = np.sort(-np.repeat(s[:-1], (s[1:] - s[:-1]).astype(np.uint8)))
    t = np.insert(-t, 0, s[-1])

    # Adjust p-values
    APVs = np.zeros((m, 1))
    aux = pairs_pvals[pairs_sorted] * t
    comparison = []
    for i in range(m):
        row = pairs_index[0][pairs_sorted[i]]
        col = pairs_index[1][pairs_sorted[i]]
        comparison.append(algorithms[row] + " vs " + algorithms[col])
        v = np.max(aux[: i + 1])
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Shaffer"])

nemenyi(p_values)

Nemenyi's procedure for adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def nemenyi(p_values: pd.DataFrame | np.ndarray) -> pd.DataFrame:
    """Nemenyi's procedure for adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if p_values.ndim != 2 or p_values.shape[0] != p_values.shape[1]:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")

    # define parameters
    k = p_values.shape[0]
    m = int(k * (k - 1) / 2.0)

    # sort p-values p(0) <= p(1) <= ... <= p(m-1)
    pairs_index = np.triu_indices(k, 1)
    pairs_pvals = p_values[pairs_index]
    pairs_sorted = np.argsort(pairs_pvals)

    # Adjust p-values
    APVs = np.zeros((m, 1))
    comparison = []
    for i in range(m):
        row = pairs_index[0][pairs_sorted[i]]
        col = pairs_index[1][pairs_sorted[i]]
        comparison.append(algorithms[row] + " vs " + algorithms[col])
        APVs[i, 0] = np.min([pairs_pvals[pairs_sorted[i]] * m, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Nemenyi"])

ranks(data, descending=False)

Computes the rank of the elements in data.

Parameters:

Name Type Description Default
data array

2-D matrix.

 required
descending bool

If true, rank is sorted in descending order.

 False

Returns:

Type Description
array

A matrix of the same shape as data where entry (i, j) is the rank of the i-th row
array

with respect to the j-th column.
Source code in src/jmetal/lab/statistical_test/functions.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def ranks(data: np.array, descending: bool = False) -> np.array:
    """Computes the rank of the elements in data.

    Args:
        data: 2-D matrix.
        descending: If true, rank is sorted in descending order.

    Returns:
        A matrix of the same shape as ``data`` where entry `(i, j)` is the rank of the i-th row
        with respect to the j-th column.
    """
    s = 0 if (descending is False) else 1

    # Compute ranks. (ranks[i][j] == rank of the i-th treatment on the j-th sample.)
    if data.ndim == 2:
        ranks = np.ones(data.shape)
        for i in range(data.shape[0]):
            values, indices, rep = np.unique(
                (-1) ** s * np.sort((-1) ** s * data[i, :]),
                return_index=True,
                return_counts=True,
            )
            for j in range(data.shape[1]):
                ranks[i, j] += indices[values == data[i, j]] + 0.5 * (rep[values == data[i, j]] - 1)
        return ranks
    elif data.ndim == 1:
        ranks = np.ones((data.size,))
        values, indices, rep = np.unique(
            (-1) ** s * np.sort((-1) ** s * data),
            return_index=True,
            return_counts=True,
        )
        for i in range(data.size):
            ranks[i] += indices[values == data[i]] + 0.5 * (rep[values == data[i]] - 1)
        return ranks

sign_test(data)

Given the results drawn from two algorithms/methods X and Y, the sign test analyses if there is a difference between X and Y.

.. note:: Null Hypothesis: Pr(X<Y)= 0.5

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :return p_value: The associated p-value from the binomial distribution. :return bstat: Number of successes.

Source code in src/jmetal/lab/statistical_test/functions.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def sign_test(data):
    """Given the results drawn from two algorithms/methods X and Y, the sign test analyses if
    there is a difference between X and Y.

    .. note:: Null Hypothesis: Pr(X<Y)= 0.5

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :return p_value: The associated p-value from the binomial distribution.
    :return bstat: Number of successes.
    """

    if type(data) == pd.DataFrame:
        data = data.values

    if data.shape[1] == 2:
        X, Y = data[:, 0], data[:, 1]
        n_perf = data.shape[0]
    else:
        raise ValueError("Initialization ERROR. Incorrect number of dimensions for axis 1")

    # Compute the differences
    Z = X - Y
    # Compute the number of pairs Z<0
    Wminus = sum(Z < 0)
    # If H_0 is true ---> W follows Binomial(n,0.5)
    p_value_minus = 1 - binom.cdf(k=Wminus, p=0.5, n=n_perf)

    # Compute the number of pairs Z>0
    Wplus = sum(Z > 0)
    # If H_0 is true ---> W follows Binomial(n,0.5)
    p_value_plus = 1 - binom.cdf(k=Wplus, p=0.5, n=n_perf)

    p_value = 2 * min([p_value_minus, p_value_plus])

    return pd.DataFrame(
        data=np.array([Wminus, Wplus, p_value]),
        index=["Num X<Y", "Num X>Y", "p-value"],
        columns=["Results"],
    )

friedman_test(data)

Friedman ranking test.

..note:: Null Hypothesis: In a set of k (>=2) treaments (or tested algorithms), all the treatments are equivalent, so their average ranks should be equal.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :return p_value: The associated p-value. :return friedman_stat: Friedman's chi-square.

Source code in src/jmetal/lab/statistical_test/functions.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def friedman_test(data):
    """Friedman ranking test.

    ..note:: Null Hypothesis: In a set of k (>=2) treaments (or tested algorithms), all the treatments are equivalent, so their average ranks should be equal.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :return p_value: The associated p-value.
    :return friedman_stat: Friedman's chi-square.
    """

    # Initial Checking
    if type(data) == pd.DataFrame:
        data = data.values

    if data.ndim == 2:
        n_samples, k = data.shape
    else:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions")
    if k < 2:
        raise ValueError("Initialization Error. Incorrect number of dimensions for axis 1.")

    # Compute ranks.
    datarank = ranks(data)

    # Compute for each algorithm the ranking average.
    avranks = np.mean(datarank, axis=0)

    # Get Friedman statistics
    friedman_stat = (
        (12.0 * n_samples) / (k * (k + 1.0)) * (np.sum(avranks**2) - (k * (k + 1) ** 2) / 4.0)
    )

    # Compute p-value
    p_value = 1.0 - chi2.cdf(friedman_stat, df=(k - 1))

    return pd.DataFrame(
        data=np.array([friedman_stat, p_value]),
        index=["Friedman-statistic", "p-value"],
        columns=["Results"],
    )

friedman_aligned_rank_test(data)

Method of aligned ranks for the Friedman test.

..note:: Null Hypothesis: In a set of k (>=2) treaments (or tested algorithms), all the treatments are equivalent, so their average ranks should be equal.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :return p_value: The associated p-value. :return aligned_rank_stat: Friedman's aligned rank chi-square statistic.

Source code in src/jmetal/lab/statistical_test/functions.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
def friedman_aligned_rank_test(data):
    """Method of aligned ranks for the Friedman test.

    ..note:: Null Hypothesis: In a set of k (>=2) treaments (or tested algorithms), all the treatments are equivalent, so their average ranks should be equal.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :return p_value: The associated p-value.
    :return aligned_rank_stat: Friedman's aligned rank chi-square statistic.
    """

    # Initial Checking
    if type(data) == pd.DataFrame:
        data = data.values

    if data.ndim == 2:
        n_samples, k = data.shape
    else:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions")
    if k < 2:
        raise ValueError("Initialization Error. Incorrect number of dimensions for axis 1.")

    # Compute the average value achieved by all algorithms in each problem
    control = np.mean(data, axis=1)
    # Compute the difference between control an data
    diff = [data[:, j] - control for j in range(data.shape[1])]
    # rank diff
    alignedRanks = ranks(np.ravel(diff))
    alignedRanks = np.reshape(alignedRanks, newshape=(n_samples, k), order="F")

    # Compute statistic
    Rhat_i = np.sum(alignedRanks, axis=1)
    Rhat_j = np.sum(alignedRanks, axis=0)
    si, sj = np.sum(Rhat_i**2), np.sum(Rhat_j**2)

    A = sj - (k * n_samples**2 / 4.0) * (k * n_samples + 1) ** 2
    B1 = k * n_samples * (k * n_samples + 1) * (2 * k * n_samples + 1) / 6.0
    B2 = si / float(k)

    alignedRanks_stat = ((k - 1) * A) / (B1 - B2)

    p_value = 1 - chi2.cdf(alignedRanks_stat, df=k - 1)

    return pd.DataFrame(
        data=np.array([alignedRanks_stat, p_value]),
        index=["Aligned Rank stat", "p-value"],
        columns=["Results"],
    )

quade_test(data)

Quade test.

..note:: Null Hypothesis: In a set of k (>=2) treaments (or tested algorithms), all the treatments are equivalent, so their average ranks should be equal.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :return p_value: The associated p-value from the F-distribution. :return fq: Computed F-value.

Source code in src/jmetal/lab/statistical_test/functions.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
def quade_test(data):
    """Quade test.

    ..note:: Null Hypothesis: In a set of k (>=2) treaments (or tested algorithms), all the treatments are equivalent, so their average ranks should be equal.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :return p_value: The associated p-value from the F-distribution.
    :return fq: Computed F-value.
    """

    # Initial Checking
    if type(data) == pd.DataFrame:
        data = data.values

    if data.ndim == 2:
        n_samples, k = data.shape
    else:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions")
    if k < 2:
        raise ValueError("Initialization Error. Incorrect number of dimensions for axis 1.")

    # Compute ranks.
    datarank = ranks(data)
    # Compute the range of each problem
    problemRange = np.max(data, axis=1) - np.min(data, axis=1)
    # Compute problem rank
    problemRank = ranks(problemRange)

    # Compute S_stat: weight of each observation within the problem, adjusted to reflect
    # the significance of the problem when it appears.
    S_stat = np.zeros((n_samples, k))
    for i in range(n_samples):
        S_stat[i, :] = problemRank[i] * (datarank[i, :] - 0.5 * (k + 1))
    Salg = np.sum(S_stat, axis=0)

    # Compute Fq (Quade Test statistic) and associated p_value
    A = np.sum(S_stat**2)
    B = np.sum(Salg**2) / float(n_samples)

    if A == B:
        Fq = np.Inf
        p_value = (1 / (np.math.factorial(k))) ** (n_samples - 1)
    else:
        Fq = (n_samples - 1.0) * B / (A - B)
        p_value = 1 - f.cdf(Fq, k - 1, (k - 1) * (n_samples - 1))

    return pd.DataFrame(
        data=np.array([Fq, p_value]), index=["Quade Test statistic", "p-value"], columns=["Results"]
    )

friedman_ph_test(data, control=None, apv_procedure=None)

Friedman post-hoc test.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :param control: optional int or string. Default None. Index or Name of the control algorithm. If control = None all FriedmanPosHocTest considers all possible comparisons among algorithms. :param apv_procedure: optional string. Default None. Name of the procedure for computing adjusted p-values. If apv_procedure is None, adjusted p-values are not computed; otherwise the values are computed according to the specified procedure: for 1-vs-all comparisons, one of 'Bonferroni', 'Holm', 'Hochberg', 'Holland', 'Finner', 'Li'; for all-vs-all comparisons, one of 'Shaffer', 'Holm', 'Nemenyi'.

:return z_values: Test statistic. :return p_values: The p-value according to the Studentized range distribution.

Source code in src/jmetal/lab/statistical_test/functions.py 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
def friedman_ph_test(data, control=None, apv_procedure=None):
    """Friedman post-hoc test.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :param control: optional int or string. Default None. Index or Name of the control algorithm. If control = None all FriedmanPosHocTest considers all possible comparisons among algorithms.
    :param apv_procedure: optional string. Default None. Name of the procedure for computing
        adjusted p-values. If ``apv_procedure`` is None, adjusted p-values are not computed;
        otherwise the values are computed according to the specified procedure: for 1-vs-all
        comparisons, one of ``'Bonferroni'``, ``'Holm'``, ``'Hochberg'``, ``'Holland'``,
        ``'Finner'``, ``'Li'``; for all-vs-all comparisons, one of ``'Shaffer'``, ``'Holm'``,
        ``'Nemenyi'``.

    :return z_values: Test statistic.
    :return p_values: The p-value according to the Studentized range distribution.
    """

    # Initial Checking
    if type(data) == pd.DataFrame:
        algorithms = data.columns
        data = data.values
    elif type(data) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(data.shape[1])])

    if control is None:
        index = algorithms
    elif type(control) == int:
        index = [algorithms[control]]
    else:
        index = [control]

    if data.ndim == 2:
        n_samples, k = data.shape
    else:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")
    if k < 2:
        raise ValueError("Initialization Error. Incorrect number of dimensions for axis 1.")

    if control is not None:
        if type(control) == int and control >= data.shape[1]:
            raise ValueError("Initialization ERROR. control is out of bounds")
        if type(control) == str and control not in algorithms:
            raise ValueError("Initialization ERROR. %s is not a column name of data" % control)

    if apv_procedure is not None:
        if apv_procedure not in [
            "Bonferroni",
            "Holm",
            "Hochberg",
            "Hommel",
            "Holland",
            "Finner",
            "Li",
            "Shaffer",
            "Nemenyi",
        ]:
            raise ValueError("Initialization ERROR. Incorrect value for APVprocedure.")

    # Compute ranks.
    datarank = ranks(data)
    # Compute for each algorithm the ranking average.
    avranks = np.mean(datarank, axis=0)

    # Compute z-values
    aux = np.sqrt((k * (k + 1)) / (6.0 * n_samples))

    if control is None:
        z = np.zeros((k, k))
        for i in range(k):
            for j in range(i + 1, k):
                z[i, j] = abs(avranks[i] - avranks[j]) / aux
        z += z.T
    else:
        if type(control) == str:
            control = int(np.where(algorithms == control)[0])
        z = np.zeros((1, k))
        for j in range(k):
            z[0, j] = abs(avranks[control] - avranks[j]) / aux

    # Compute associated p-value
    p_value = 2 * (1.0 - norm.cdf(z))

    pvalues_df = pd.DataFrame(data=p_value, index=index, columns=algorithms)
    zvalues_df = pd.DataFrame(data=z, index=index, columns=algorithms)

    if apv_procedure is None:
        return zvalues_df, pvalues_df
    else:
        if apv_procedure == "Bonferroni":
            ap_vs_df = bonferroni_dunn(pvalues_df, control=control)
        elif apv_procedure == "Holm":
            ap_vs_df = holm(pvalues_df, control=control)
        elif apv_procedure == "Hochberg":
            ap_vs_df = hochberg(pvalues_df, control=control)
        elif apv_procedure == "Holland":
            ap_vs_df = holland(pvalues_df, control=control)
        elif apv_procedure == "Finner":
            ap_vs_df = finner(pvalues_df, control=control)
        elif apv_procedure == "Li":
            ap_vs_df = li(pvalues_df, control=control)
        elif apv_procedure == "Shaffer":
            ap_vs_df = shaffer(pvalues_df)
        elif apv_procedure == "Nemenyi":
            ap_vs_df = nemenyi(pvalues_df)

        return zvalues_df, pvalues_df, ap_vs_df

friedman_aligned_ph_test(data, control=None, apv_procedure=None)

Friedman Aligned Ranks post-hoc test.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :param control: optional int or string. Default None. Index or Name of the control algorithm. If control = None all FriedmanPosHocTest considers all possible comparisons among algorithms. :param apv_procedure: optional string. Default None. Name of the procedure for computing adjusted p-values. If apv_procedure is None, adjusted p-values are not computed; otherwise the values are computed according to the specified procedure: for 1-vs-all comparisons, one of 'Bonferroni', 'Holm', 'Hochberg', 'Holland', 'Finner', 'Li'; for all-vs-all comparisons, one of 'Shaffer', 'Holm', 'Nemenyi'.

:return z_values: Test statistic. :return p_values: The p-value according to the Studentized range distribution.

Source code in src/jmetal/lab/statistical_test/functions.py 
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
def friedman_aligned_ph_test(data, control=None, apv_procedure=None):
    """Friedman Aligned Ranks post-hoc test.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :param control: optional int or string. Default None. Index or Name of the control algorithm. If control = None all FriedmanPosHocTest considers all possible comparisons among algorithms.
    :param apv_procedure: optional string. Default None. Name of the procedure for computing
        adjusted p-values. If ``apv_procedure`` is None, adjusted p-values are not computed;
        otherwise the values are computed according to the specified procedure: for 1-vs-all
        comparisons, one of ``'Bonferroni'``, ``'Holm'``, ``'Hochberg'``, ``'Holland'``,
        ``'Finner'``, ``'Li'``; for all-vs-all comparisons, one of ``'Shaffer'``, ``'Holm'``,
        ``'Nemenyi'``.

    :return z_values: Test statistic.
    :return p_values: The p-value according to the Studentized range distribution.
    """

    # Initial Checking
    if type(data) == pd.DataFrame:
        algorithms = data.columns
        data = data.values
    elif type(data) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(data.shape[1])])

    if control is None:
        index = algorithms
    elif type(control) == int:
        index = [algorithms[control]]
    else:
        index = [control]

    if data.ndim == 2:
        n_samples, k = data.shape
    else:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")
    if k < 2:
        raise ValueError("Initialization Error. Incorrect number of dimensions for axis 1.")

    if control is not None:
        if type(control) == int and control >= data.shape[1]:
            raise ValueError("Initialization ERROR. control is out of bounds")
        if type(control) == str and control not in algorithms:
            raise ValueError("Initialization ERROR. %s is not a column name of data" % control)

    # Compute the average value achieved by all algorithms in each problem
    problemmean = np.mean(data, axis=1)
    # Compute the difference between control an data
    diff = np.zeros((n_samples, k))
    for j in range(k):
        diff[:, j] = data[:, j] - problemmean

    alignedRanks = ranks(np.ravel(diff))
    alignedRanks = np.reshape(alignedRanks, newshape=(n_samples, k))

    # Average ranks
    avranks = np.mean(alignedRanks, axis=0)

    # Compute test statistics
    aux = 1.0 / np.sqrt(k * (n_samples + 1) / 6.0)
    if control is None:
        z = np.zeros((k, k))
        for i in range(k):
            for j in range(i + 1, k):
                z[i, j] = abs(avranks[i] - avranks[j]) * aux
        z += z.T
    else:
        if type(control) == str:
            control = int(np.where(algorithms == control)[0])
        z = np.zeros((1, k))
        for j in range(k):
            z[0, j] = abs(avranks[control] - avranks[j]) * aux

    # Compute associated p-value
    p_value = 2 * (1.0 - norm.cdf(z))

    pvalues_df = pd.DataFrame(data=p_value, index=index, columns=algorithms)
    zvalues_df = pd.DataFrame(data=z, index=index, columns=algorithms)

    if apv_procedure is None:
        return zvalues_df, pvalues_df
    else:
        if apv_procedure == "Bonferroni":
            ap_vs_df = bonferroni_dunn(pvalues_df, control=control)
        elif apv_procedure == "Holm":
            ap_vs_df = holm(pvalues_df, control=control)
        elif apv_procedure == "Hochberg":
            ap_vs_df = hochberg(pvalues_df, control=control)
        elif apv_procedure == "Holland":
            ap_vs_df = holland(pvalues_df, control=control)
        elif apv_procedure == "Finner":
            ap_vs_df = finner(pvalues_df, control=control)
        elif apv_procedure == "Li":
            ap_vs_df = li(pvalues_df, control=control)
        elif apv_procedure == "Shaffer":
            ap_vs_df = shaffer(pvalues_df)
        elif apv_procedure == "Nemenyi":
            ap_vs_df = nemenyi(pvalues_df)

        return zvalues_df, pvalues_df, ap_vs_df

quade_ph_test(data, control=None, apv_procedure=None)

Quade post-hoc test.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :param control: optional int or string. Default None. Index or Name of the control algorithm. If control = None all FriedmanPosHocTest considers all possible comparisons among algorithms. :param apv_procedure: optional string. Default None. Name of the procedure for computing adjusted p-values. If apv_procedure is None, adjusted p-values are not computed; otherwise the values are computed according to the specified procedure: for 1-vs-all comparisons, one of 'Bonferroni', 'Holm', 'Hochberg', 'Holland', 'Finner', 'Li'; for all-vs-all comparisons, one of 'Shaffer', 'Holm', 'Nemenyi'.

:return z_values: Test statistic. :return p_values: The p-value according to the Studentized range distribution.

Source code in src/jmetal/lab/statistical_test/functions.py 
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
def quade_ph_test(data, control=None, apv_procedure=None):
    """Quade post-hoc test.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :param control: optional int or string. Default None. Index or Name of the control algorithm. If control = None all FriedmanPosHocTest considers all possible comparisons among algorithms.
    :param apv_procedure: optional string. Default None. Name of the procedure for computing
        adjusted p-values. If ``apv_procedure`` is None, adjusted p-values are not computed;
        otherwise the values are computed according to the specified procedure: for 1-vs-all
        comparisons, one of ``'Bonferroni'``, ``'Holm'``, ``'Hochberg'``, ``'Holland'``,
        ``'Finner'``, ``'Li'``; for all-vs-all comparisons, one of ``'Shaffer'``, ``'Holm'``,
        ``'Nemenyi'``.

    :return z_values: Test statistic.
    :return p_values: The p-value according to the Studentized range distribution.
    """

    # Initial Checking
    if type(data) == pd.DataFrame:
        algorithms = data.columns
        data = data.values
    elif type(data) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(data.shape[1])])

    if control is None:
        index = algorithms
    elif type(control) == int:
        index = [algorithms[control]]
    else:
        index = [control]

    if data.ndim == 2:
        n_samples, k = data.shape
    else:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")
    if k < 2:
        raise ValueError("Initialization Error. Incorrect number of dimensions for axis 1.")

    if control is not None:
        if type(control) == int and control >= data.shape[1]:
            raise ValueError("Initialization ERROR. control is out of bounds")
        if type(control) == str and control not in algorithms:
            raise ValueError("Initialization ERROR. %s is not a column name of data" % control)

    # Compute ranks.
    datarank = ranks(data)
    # Compute the range of each problem
    problemRange = np.max(data, axis=1) - np.min(data, axis=1)
    # Compute problem rank
    problemRank = ranks(problemRange)

    # Compute average rakings
    W = np.zeros((n_samples, k))
    for i in range(n_samples):
        W[i, :] = problemRank[i] * datarank[i, :]
    avranks = 2 * np.sum(W, axis=0) / (n_samples * (n_samples + 1))
    # Compute test statistics
    aux = 1.0 / np.sqrt(
        k * (k + 1) * (2 * n_samples + 1) * (k - 1) / (18.0 * n_samples * (n_samples + 1))
    )
    if control is None:
        z = np.zeros((k, k))
        for i in range(k):
            for j in range(i + 1, k):
                z[i, j] = abs(avranks[i] - avranks[j]) * aux
        z += z.T
    else:
        if type(control) == str:
            control = int(np.where(algorithms == control)[0])
        z = np.zeros((1, k))
        for j in range(k):
            z[0, j] = abs(avranks[control] - avranks[j]) * aux

    # Compute associated p-value
    p_value = 2 * (1.0 - norm.cdf(z))

    pvalues_df = pd.DataFrame(data=p_value, index=index, columns=algorithms)
    zvalues_df = pd.DataFrame(data=z, index=index, columns=algorithms)

    if apv_procedure is None:
        return zvalues_df, pvalues_df
    else:
        if apv_procedure == "Bonferroni":
            ap_vs_df = bonferroni_dunn(pvalues_df, control=control)
        elif apv_procedure == "Holm":
            ap_vs_df = holm(pvalues_df, control=control)
        elif apv_procedure == "Hochberg":
            ap_vs_df = hochberg(pvalues_df, control=control)
        elif apv_procedure == "Holland":
            ap_vs_df = holland(pvalues_df, control=control)
        elif apv_procedure == "Finner":
            ap_vs_df = finner(pvalues_df, control=control)
        elif apv_procedure == "Li":
            ap_vs_df = li(pvalues_df, control=control)
        elif apv_procedure == "Shaffer":
            ap_vs_df = shaffer(pvalues_df)
        elif apv_procedure == "Nemenyi":
            ap_vs_df = nemenyi(pvalues_df)

        return zvalues_df, pvalues_df, ap_vs_df

apv_procedures

bonferroni_dunn(p_values, control)

Bonferroni-Dunn's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def bonferroni_dunn(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Bonferroni-Dunn's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        APVs[i, 0] = np.min([(k - 1) * p_values[0, argsorted_pvals[i]], 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Bonferroni"])

holland(p_values, control)

Holland's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def holland(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Holland's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    # --------------------------------------------------------------------------
    # ------------------------------- Procedure --------------------------------
    # --------------------------------------------------------------------------
    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        aux = k - 1 - np.arange(i + 1)
        v = np.max(1 - (1 - p_values[0, argsorted_pvals[: (i + 1)]]) ** aux)
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Holland"])

finner(p_values, control)

Finner's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def finner(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Finner's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        aux = float(k - 1) / (np.arange(i + 1) + 1)
        v = np.max(1 - (1 - p_values[0, argsorted_pvals[: (i + 1)]]) ** aux)
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Finner"])

hochberg(p_values, control)

Hochberg's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def hochberg(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Hochberg's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        aux = np.arange(k, i, -1).astype(np.uint8)
        v = np.max(p_values[0, argsorted_pvals[aux - 1]] * (k - aux))
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Hochberg"])

li(p_values, control)

Li's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str

Index or name of the control algorithm. If provided, control vs all comparisons are considered, else all vs all.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def li(p_values: pd.DataFrame | np.ndarray, control: int | str) -> pd.DataFrame:
    """Li's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm. If provided, control vs all
            comparisons are considered, else all vs all.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])
    if control is None:
        raise ValueError("Initialization ERROR. Incorrect value for control.")

    k = p_values.shape[1]

    # sort p-values p(0) <= p(1) <= ... <= p(k-1)
    argsorted_pvals = np.argsort(p_values[0, :])

    APVs = np.zeros((k - 1, 1))
    comparison = []
    for i in range(k - 1):
        comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
        APVs[i, 0] = np.min(
            [
                p_values[0, argsorted_pvals[-2]],
                p_values[0, argsorted_pvals[i]]
                / (p_values[0, argsorted_pvals[i]] + 1 - p_values[0, argsorted_pvals[-2]]),
            ]
        )
    return pd.DataFrame(data=APVs, index=comparison, columns=["Li"])

holm(p_values, control=None)

Holm's procedure for the adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values obtained from a ranking test.

 required
control int | str | None

Index or name of the control algorithm. If provided, control vs all comparisons are considered, else all vs all.

 None

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def holm(p_values: pd.DataFrame | np.ndarray, control: int | str | None = None) -> pd.DataFrame:
    """Holm's procedure for the adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values obtained from a ranking test.
        control: Index or name of the control algorithm. If provided, control vs all
            comparisons are considered, else all vs all.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if type(control) == str:
        control = int(np.where(algorithms == control)[0])

    if type(control) == int:
        k = p_values.shape[1]

        # sort p-values p(0) <= p(1) <= ... <= p(k-1)
        argsorted_pvals = np.argsort(p_values[0, :])

        APVs = np.zeros((k - 1, 1))
        comparison = []
        for i in range(k - 1):
            aux = k - 1 - np.arange(i + 1)
            comparison.append(algorithms[control] + " vs " + algorithms[argsorted_pvals[i]])
            v = np.max(aux * p_values[0, argsorted_pvals[: (i + 1)]])
            APVs[i, 0] = np.min([v, 1])

    elif control is None:
        k = p_values.shape[1]
        m = int((k * (k - 1)) / 2.0)

        # sort p-values p(0) <= p(1) <= ... <= p(m-1)
        pairs_index = np.triu_indices(k, 1)
        pairs_pvals = p_values[pairs_index]
        pairs_sorted = np.argsort(pairs_pvals)

        APVs = np.zeros((m, 1))
        aux = pairs_pvals[pairs_sorted] * (m - np.arange(m))
        comparison = []
        for i in range(m):
            row = pairs_index[0][pairs_sorted[i]]
            col = pairs_index[1][pairs_sorted[i]]
            comparison.append(algorithms[row] + " vs " + algorithms[col])
            v = np.max(aux[: i + 1])
            APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Holm"])

shaffer(p_values)

Shaffer's procedure for adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
def shaffer(p_values: pd.DataFrame | np.ndarray) -> pd.DataFrame:
    """Shaffer's procedure for adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    def S(k):
        """Computes the set of possible numbers of true hypotheses.

        Args:
            k: Number of algorithms being compared.

        Returns:
            Set of true hypotheses.
        """

        from scipy.special import binom as binomial

        TrueHset = [0]
        if k > 1:
            for j in np.arange(k, 0, -1, dtype=int):
                TrueHset = list(set(TrueHset) | set([binomial(j, 2) + x for x in S(k - j)]))
        return TrueHset

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if p_values.ndim != 2 or p_values.shape[0] != p_values.shape[1]:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")

    # define parameters
    k = p_values.shape[0]
    m = int(k * (k - 1) / 2.0)
    s = np.array(S(k)[1:])

    # sort p-values p(0) <= p(1) <= ... <= p(m-1)
    pairs_index = np.triu_indices(k, 1)
    pairs_pvals = p_values[pairs_index]
    pairs_sorted = np.argsort(pairs_pvals)

    # compute ti: max number of hypotheses that can be true given that any
    # (i-1) hypotheses are false.
    t = np.sort(-np.repeat(s[:-1], (s[1:] - s[:-1]).astype(np.uint8)))
    t = np.insert(-t, 0, s[-1])

    # Adjust p-values
    APVs = np.zeros((m, 1))
    aux = pairs_pvals[pairs_sorted] * t
    comparison = []
    for i in range(m):
        row = pairs_index[0][pairs_sorted[i]]
        col = pairs_index[1][pairs_sorted[i]]
        comparison.append(algorithms[row] + " vs " + algorithms[col])
        v = np.max(aux[: i + 1])
        APVs[i, 0] = np.min([v, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Shaffer"])

nemenyi(p_values)

Nemenyi's procedure for adjusted p-value computation.

Parameters:

Name Type Description Default
p_values DataFrame | ndarray

2-D array or DataFrame containing the p-values.

 required

Returns:

Type Description
DataFrame

DataFrame containing the adjusted p-values.
Source code in src/jmetal/lab/statistical_test/apv_procedures.py 
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def nemenyi(p_values: pd.DataFrame | np.ndarray) -> pd.DataFrame:
    """Nemenyi's procedure for adjusted p-value computation.

    Args:
        p_values: 2-D array or DataFrame containing the p-values.

    Returns:
        DataFrame containing the adjusted p-values.
    """

    # Initial Checking
    if type(p_values) == pd.DataFrame:
        algorithms = p_values.columns
        p_values = p_values.values
    elif type(p_values) == np.ndarray:
        algorithms = np.array(["Alg%d" % alg for alg in range(p_values.shape[1])])

    if p_values.ndim != 2 or p_values.shape[0] != p_values.shape[1]:
        raise ValueError("Initialization ERROR. Incorrect number of array dimensions.")

    # define parameters
    k = p_values.shape[0]
    m = int(k * (k - 1) / 2.0)

    # sort p-values p(0) <= p(1) <= ... <= p(m-1)
    pairs_index = np.triu_indices(k, 1)
    pairs_pvals = p_values[pairs_index]
    pairs_sorted = np.argsort(pairs_pvals)

    # Adjust p-values
    APVs = np.zeros((m, 1))
    comparison = []
    for i in range(m):
        row = pairs_index[0][pairs_sorted[i]]
        col = pairs_index[1][pairs_sorted[i]]
        comparison.append(algorithms[row] + " vs " + algorithms[col])
        APVs[i, 0] = np.min([pairs_pvals[pairs_sorted[i]] * m, 1])
    return pd.DataFrame(data=APVs, index=comparison, columns=["Nemenyi"])

bayesian

bayesian_sign_test(data, rope_limits=[-0.01, 0.01], prior_strength=0.5, prior_place='rope', sample_size=50000, return_sample=False, rng=None)

Bayesian version of the sign test.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :param rope_limits: array_like. Default [-0.01, 0.01]. Limits of the practical equivalence. :param prior_strength: positive float. Default 0.5. Value of the prior strengt :param prior_place: string {left, rope, right}. Default 'left'. Place of the pseudo-observation z_0. :param sample_size: integer. Default 10000. Total number of random_search samples generated :param return_sample: boolean. Default False. If true, also return the samples drawn from the Dirichlet process. :param rng: Optional random generator for reproducible results. When None, a fresh np.random.default_rng() is used (i.e. results vary between calls, as before this parameter existed).

:return: List of posterior probabilities: [Pr(algorith_1 < algorithm_2), Pr(algorithm_1 equiv algorithm_2), Pr(algorithm_1 > algorithm_2)]

Source code in src/jmetal/lab/statistical_test/bayesian.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def bayesian_sign_test(
    data,
    rope_limits=[-0.01, 0.01],
    prior_strength=0.5,
    prior_place="rope",
    sample_size=50000,
    return_sample=False,
    rng: np.random.Generator | None = None,
):
    """Bayesian version of the sign test.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :param rope_limits: array_like. Default [-0.01, 0.01]. Limits of the practical equivalence.
    :param prior_strength: positive float. Default 0.5. Value of the prior strengt
    :param prior_place: string {left, rope, right}. Default 'left'. Place of the pseudo-observation z_0.
    :param sample_size: integer. Default 10000. Total number of random_search samples generated
    :param return_sample: boolean. Default False. If true, also return the samples drawn from the Dirichlet process.
    :param rng: Optional random generator for reproducible results. When None, a fresh
        np.random.default_rng() is used (i.e. results vary between calls, as before this
        parameter existed).

    :return: List of posterior probabilities:
        [Pr(algorith_1 < algorithm_2),
        Pr(algorithm_1 equiv algorithm_2),
        Pr(algorithm_1 > algorithm_2)]
    """
    rng = rng if rng is not None else np.random.default_rng()

    # Initial Checking
    if type(data) == pd.DataFrame:
        data = data.values

    if data.shape[1] == 2:
        sample1, sample2 = data[:, 0], data[:, 1]
        n = data.shape[0]
    else:
        raise ValueError("Initialization ERROR. Incorrect number of dimensions for axis 1")

    if prior_strength <= 0:
        raise ValueError("Initialization ERROR. prior_strength mustb be a positive float")

    if prior_place not in ["left", "rope", "right"]:
        raise ValueError("Initialization ERROR. Incorrect value fro prior_place")

    # Compute the differences
    Z = sample1 - sample2

    # Compute the number of pairs diff > right_limit
    Nright = sum(rope_limits[1] < Z)
    # Compute the number of pairs diff < right_lelft
    Nleft = sum(rope_limits[0] > Z)
    # Compute the number of pairs diff in rope_limits
    Nequiv = n - Nright - Nleft

    # compute the the probabilities that the mean difference of accuracy is in
    # the interval (−Inf, left), [left, right], or (ringth, Inf).

    # Parameters of the Dirichlet distribution
    alpha = np.array([Nleft, Nequiv, Nright], dtype=float) + 1e-6
    alpha[["left", "rope", "right"].index(prior_place)] += prior_strength
    # Simulate dirichlet process
    Dprocess = rng.dirichlet(alpha, sample_size)

    # Compute posterior probabilities
    winner_id = np.argmax(Dprocess, axis=1)
    win_left = sum(winner_id == 0)
    win_rifht = sum(winner_id == 2)
    win_rope = sample_size - win_left - win_rifht

    if return_sample is True:
        return np.array([win_left, win_rope, win_rifht]) / float(sample_size), Dprocess
    else:
        return np.array([win_left, win_rope, win_rifht]) / float(sample_size)

bayesian_signed_rank_test(data, rope_limits=[-0.01, 0.01], prior_strength=1.0, prior_place='rope', sample_size=10000, return_sample=False, rng=None)

Bayesian version of the signed rank test.

:param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem. :param rope_limits: array_like. Default [-0.01, 0.01]. Limits of the practical equivalence. :param prior_strength: positive float. Default 0.5. Value of the prior strengt :param prior_place: string {left, rope, right}. Default 'left'. Place of the pseudo-observation z_0. :param sample_size: integer. Default 10000. Total number of random_search samples generated :param return_sample: boolean. Default False. If true, also return the samples drawn from the Dirichlet process. :param rng: Optional random generator for reproducible results. When None, a fresh np.random.default_rng() is used (i.e. results vary between calls, as before this parameter existed).

:return: List of posterior probabilities: [Pr(algorith_1 < algorithm_2), Pr(algorithm_1 equiv algorithm_2), Pr(algorithm_1 > algorithm_2)]

Source code in src/jmetal/lab/statistical_test/bayesian.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def bayesian_signed_rank_test(
    data,
    rope_limits=[-0.01, 0.01],
    prior_strength=1.0,
    prior_place="rope",
    sample_size=10000,
    return_sample=False,
    rng: np.random.Generator | None = None,
):
    """Bayesian version of the signed rank test.

    :param data: An (n x 2) array or DataFrame contaning the results. In data, each column represents an algorithm and, and each row a problem.
    :param rope_limits: array_like. Default [-0.01, 0.01]. Limits of the practical equivalence.
    :param prior_strength: positive float. Default 0.5. Value of the prior strengt
    :param prior_place: string {left, rope, right}. Default 'left'. Place of the pseudo-observation z_0.
    :param sample_size: integer. Default 10000. Total number of random_search samples generated
    :param return_sample: boolean. Default False. If true, also return the samples drawn from the Dirichlet process.
    :param rng: Optional random generator for reproducible results. When None, a fresh
        np.random.default_rng() is used (i.e. results vary between calls, as before this
        parameter existed).

    :return: List of posterior probabilities:
        [Pr(algorith_1 < algorithm_2), Pr(algorithm_1 equiv algorithm_2), Pr(algorithm_1 > algorithm_2)]
    """
    rng = rng if rng is not None else np.random.default_rng()

    def weights(n, s):
        alpha = np.ones(n + 1)
        alpha[0] = s
        return rng.dirichlet(alpha, 1)[0]

    # Initial Checking
    if type(data) == pd.DataFrame:
        data = data.values

    if data.shape[1] == 2:
        sample1, sample2 = data[:, 0], data[:, 1]
        n = data.shape[0]
    else:
        raise ValueError("Initialization ERROR. Incorrect number of dimensions for axis 1")

    if prior_strength <= 0:
        raise ValueError("Initialization ERROR. prior_strength must be a positive float")

    if prior_place not in ["left", "rope", "right"]:
        raise ValueError("Initialization ERROR. Incorrect value for prior_place")

    # Compute the differences
    Z = sample1 - sample2
    Z0 = [-float("Inf"), 0.0, float("Inf")][["left", "rope", "right"].index(prior_place)]
    Z = np.concatenate(([Z0], Z), axis=None)

    # compute the the probabilities that the mean difference of accuracy is in
    # the interval (−Inf, left), [left, right], or (ringth, Inf).

    Dprocess = np.zeros((sample_size, 3))
    for mc in range(sample_size):
        W = weights(n, prior_strength)
        for i in range(n + 1):
            for j in range(i, n + 1):
                aux = Z[i] + Z[j]
                sumval = 2 * (W[i] * W[j]) if i != j else (W[i] * W[j])
                if aux < 2 * rope_limits[0]:
                    Dprocess[mc, 0] += sumval
                elif aux > 2 * rope_limits[1]:
                    Dprocess[mc, 2] += sumval
                else:
                    Dprocess[mc, 1] += sumval

    # Compute posterior probabilities
    winner_id = np.argmax(Dprocess, axis=1)
    win_left = sum(winner_id == 0)
    win_rifht = sum(winner_id == 2)
    win_rope = sample_size - win_left - win_rifht

    if return_sample is True:
        return np.array([win_left, win_rope, win_rifht]) / float(sample_size), Dprocess
    else:
        return np.array([win_left, win_rope, win_rifht]) / float(sample_size)

critical_distance

NemenyiCD(alpha, num_alg, num_dataset)

Computes Nemenyi's critical difference: * CD = q_alpha * sqrt(num_alg(num_alg + 1)/(6num_prob)) where q_alpha is the critical value, of the Studentized range statistic divided by sqrt(2). :param alpha: {0.1, 0.999}. Significance level. :param num_alg: number of tested algorithms. :param num_dataset: Number of problems/datasets where the algorithms have been tested.

Source code in src/jmetal/lab/statistical_test/critical_distance.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def NemenyiCD(alpha: float, num_alg, num_dataset):
    """Computes Nemenyi's critical difference:
    * CD = q_alpha * sqrt(num_alg*(num_alg + 1)/(6*num_prob))
    where q_alpha is the critical value, of the Studentized range statistic divided by sqrt(2).
    :param alpha: {0.1, 0.999}. Significance level.
    :param num_alg: number of tested algorithms.
    :param num_dataset: Number of problems/datasets where the algorithms have been tested.
    """

    # get critical value
    q_alpha = qsturng(p=1 - alpha, r=num_alg, v=num_alg * (num_dataset - 1)) / np.sqrt(2)

    # compute the critical difference
    cd = q_alpha * np.sqrt(num_alg * (num_alg + 1) / (6.0 * num_dataset))

    return cd

CDplot(results, alpha=0.05, higher_is_better=False, alg_names=None, output_filename='cdplot.eps')

CDgraph plots the critical difference graph show in Janez Demsar's 2006 work: * Statistical Comparisons of Classifiers over Multiple Data Sets. :param results: A 2-D array containing results from each algorithm. Each row of 'results' represents an algorithm, and each column a dataset. :param alpha: {0.1, 0.999}. Significance level for the critical difference. :param alg_names: Names of the tested algorithms.

Source code in src/jmetal/lab/statistical_test/critical_distance.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def CDplot(
    results,
    alpha: float = 0.05,
    higher_is_better: bool = False,
    alg_names: list = None,
    output_filename: str = "cdplot.eps",
):
    """CDgraph plots the critical difference graph show in Janez Demsar's 2006 work:
    * Statistical Comparisons of Classifiers over Multiple Data Sets.
    :param results: A 2-D array containing results from each algorithm. Each row of 'results' represents an algorithm, and each column a dataset.
    :param alpha: {0.1, 0.999}. Significance level for the critical difference.
    :param alg_names: Names of the tested algorithms.
    """

    def _join_alg(avranks, num_alg, cd):
        """
        join_alg returns the set of non significant methods
        """

        # get all pairs
        sets = (-1) * np.ones((num_alg, 2))
        for i in range(num_alg):
            elements = np.where(
                np.logical_and(avranks - avranks[i] > 0, avranks - avranks[i] < cd)
            )[0]
            if elements.size > 0:
                sets[i, :] = [avranks[i], avranks[elements[-1]]]
        sets = np.delete(sets, np.where(sets[:, 0] < 0)[0], axis=0)

        # group pairs
        group = sets[0, :]
        for i in range(1, sets.shape[0]):
            if sets[i - 1, 1] < sets[i, 1]:
                group = np.vstack((group, sets[i, :]))

        return group

    # Initial Checking
    if type(results) == pd.DataFrame:
        alg_names = results.index
        results = results.values
    elif type(results) == np.ndarray and alg_names is None:
        alg_names = np.array(["Alg%d" % alg for alg in range(results.shape[1])])

    if results.ndim == 2:
        num_alg, num_dataset = results.shape
    else:
        raise ValueError("Initialization ERROR: In CDplot(...) results must be 2-D array")

    # Get the critical difference
    cd = NemenyiCD(alpha, num_alg, num_dataset)

    # Compute ranks. (ranks[i][j] rank of the i-th algorithm on the j-th problem.)

    rranks = ranks(results.T, descending=higher_is_better)

    # Compute for each algorithm the ranking averages.
    avranks = np.transpose(np.mean(rranks, axis=0))
    indices = np.argsort(avranks).astype(np.uint8)
    avranks = avranks[indices]

    # Split algorithms.
    spoint = np.round(num_alg / 2.0).astype(np.uint8)
    leftalg = avranks[:spoint]
    rightalg = avranks[spoint:]
    rows = np.ceil(num_alg / 2.0).astype(np.uint8)

    # Figure settings.
    highest = np.ceil(np.max(avranks)).astype(np.uint8)  # highest shown rank
    lowest = np.floor(np.min(avranks)).astype(np.uint8)  # lowest shown rank
    width = 6  # default figure width (in inches)
    height = 0.575 * (rows + 1)  # figure height

    """
                        FIGURE
      (1,0)
        +-----+---------------------------+-------+
        |     |                           |       |
        |     |                           |       |
        |     |                           |       |
        +-----+---------------------------+-------+ stop
        |     |                           |       |
        |     |                           |       |
        |     |                           |       |
        |     |                           |       |
        |     |                           |       |
        |     |                           |       |
        +-----+---------------------------+-------+ sbottom
        |     |                           |       |
        +-----+---------------------------+-------+
            sleft                       sright     (0,1)
    """

    stop, sbottom, sleft, sright = 0.65, 0.1, 0.15, 0.85

    # main horizontal axis length
    lline = sright - sleft

    # Initialize figure
    fig = plt.figure(figsize=(width, height), facecolor="white")
    ax = fig.add_axes([0, 0, 1, 1])
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.set_axis_off()

    # Main horizontal axis
    ax.hlines(stop, sleft, sright, color="black", linewidth=0.7)
    for xi in range(highest - lowest + 1):
        # Plot mayor ticks
        ax.vlines(
            x=sleft + (lline * xi) / (highest - lowest),
            ymin=stop,
            ymax=stop + 0.05,
            color="black",
            linewidth=0.7,
        )
        # Mayor ticks labels
        ax.text(
            x=sleft + (lline * xi) / (highest - lowest),
            y=stop + 0.06,
            s=str(lowest + xi),
            ha="center",
            va="bottom",
        )
        # Minor ticks
        if xi < highest - lowest:
            ax.vlines(
                x=sleft + (lline * (xi + 0.5)) / (highest - lowest),
                ymin=stop,
                ymax=stop + 0.025,
                color="black",
                linewidth=0.7,
            )

    # Plot lines/names for left models
    vspace = 0.5 * (stop - sbottom) / (spoint + 1)
    for i in range(spoint):
        ax.vlines(
            x=sleft + (lline * (leftalg[i] - lowest)) / (highest - lowest),
            ymin=sbottom + (spoint - 1 - i) * vspace,
            ymax=stop,
            color="black",
            linewidth=0.7,
        )
        ax.hlines(
            y=sbottom + (spoint - 1 - i) * vspace,
            xmin=sleft,
            xmax=sleft + (lline * (leftalg[i] - lowest)) / (highest - lowest),
            color="black",
            linewidth=0.7,
        )
        ax.text(
            x=sleft - 0.01,
            y=sbottom + (spoint - 1 - i) * vspace,
            s=alg_names[indices][i],
            ha="right",
            va="center",
        )

    # Plot lines/names for right models
    vspace = 0.5 * (stop - sbottom) / (num_alg - spoint + 1)
    for i in range(num_alg - spoint):
        ax.vlines(
            x=sleft + (lline * (rightalg[i] - lowest)) / (highest - lowest),
            ymin=sbottom + i * vspace,
            ymax=stop,
            color="black",
            linewidth=0.7,
        )
        ax.hlines(
            y=sbottom + i * vspace,
            xmin=sleft + (lline * (rightalg[i] - lowest)) / (highest - lowest),
            xmax=sright,
            color="black",
            linewidth=0.7,
        )
        ax.text(
            x=sright + 0.01,
            y=sbottom + i * vspace,
            s=alg_names[indices][spoint + i],
            ha="left",
            va="center",
        )

    # Plot critical difference rule
    if sleft + (cd * lline) / (highest - lowest) <= sright:
        ax.hlines(
            y=stop + 0.2, xmin=sleft, xmax=sleft + (cd * lline) / (highest - lowest), linewidth=1.5
        )
        ax.text(
            x=sleft + 0.5 * (cd * lline) / (highest - lowest),
            y=stop + 0.21,
            s="CD=%.3f" % cd,
            ha="center",
            va="bottom",
        )
    else:
        ax.text(x=(sleft + sright) / 2, y=stop + 0.2, s="CD=%.3f" % cd, ha="center", va="bottom")

    # Get pair of non-significant methods
    nonsig = _join_alg(avranks, num_alg, cd)
    if nonsig.ndim == 2:
        if nonsig.shape[0] == 2:
            left_lines = np.reshape(nonsig[0, :], (1, 2))
            right_lines = np.reshape(nonsig[1, :], (1, 2))
        else:
            left_lines = nonsig[: np.round(nonsig.shape[0] / 2.0).astype(np.uint8), :]
            right_lines = nonsig[np.round(nonsig.shape[0] / 2.0).astype(np.uint8) :, :]
    else:
        left_lines = np.reshape(nonsig, (1, nonsig.shape[0]))

    # plot from the left
    vspace = 0.5 * (stop - sbottom) / (left_lines.shape[0] + 1)
    for i in range(left_lines.shape[0]):
        ax.hlines(
            y=stop - (i + 1) * vspace,
            xmin=sleft + lline * (left_lines[i, 0] - lowest - 0.025) / (highest - lowest),
            xmax=sleft + lline * (left_lines[i, 1] - lowest + 0.025) / (highest - lowest),
            linewidth=2,
        )

    # plot from the rigth
    if nonsig.ndim == 2:
        vspace = 0.5 * (stop - sbottom) / (left_lines.shape[0])
        for i in range(right_lines.shape[0]):
            ax.hlines(
                y=stop - (i + 1) * vspace,
                xmin=sleft + lline * (right_lines[i, 0] - lowest - 0.025) / (highest - lowest),
                xmax=sleft + lline * (right_lines[i, 1] - lowest + 0.025) / (highest - lowest),
                linewidth=2,
            )

    plt.savefig(output_filename, bbox_inches="tight")
    plt.show()