import numpy as np
import pandas as pd
[docs]def bayesian_sign_test(data, rope_limits=[-0.01, 0.01], prior_strength=0.5, prior_place='rope', sample_size=50000,
return_sample=False):
""" 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.
:return: List of posterior probabilities:
[Pr(algorith_1 < algorithm_2),
Pr(algorithm_1 equiv algorithm_2),
Pr(algorithm_1 > algorithm_2)]
"""
# 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(Z > rope_limits[1])
# Compute the number of pairs diff < right_lelft
Nleft = sum(Z < rope_limits[0])
# 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 = np.random.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)
[docs]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):
""" 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.
:return: List of posterior probabilities:
[Pr(algorith_1 < algorithm_2), Pr(algorithm_1 equiv algorithm_2), Pr(algorithm_1 > algorithm_2)]
"""
def weights(n, s):
alpha = np.ones(n + 1)
alpha[0] = s
return np.random.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)
