Module stikpetP.distributions.dist_kendall_tau
Expand source code
import pandas as pd
from math import exp
from math import factorial
from statistics import NormalDist
def di_kendall_tau(n, tau, method="as71", cc=False):
if method=="as71":
S = n * (n - 1) / 2 * abs(tau)
if cc:
S = abs(S) - 1
pvalue = 2 * he_AS71(S, n)
elif method=="kendall":
c = (tau + 1) * (n * (n - 1)) / 4
pvalue = he_kendall(n, c)
return pvalue
def he_AS71(S, n):
#Upper tail probabilities of tau
h = [0]*15
L = [[0]*16]*3
L = pd.DataFrame(L)
#Check validity of IS and N values
PRTAUS = 1
if (n < 1):
he_AS71 = PRTAUS
else:
m = n * (n - 1) / 2 - abs(S)
if ((m < 0) or (m > (m / 2) * 2)):
he_AS71 = PRTAUS
elif m == 0 and S <= 0:
he_AS71 = PRTAUS
else:
if (n > 8):
#GOTO 7
#7
#Calculation of Tchebycheff-Hermite polynomials
x = (S - 1) / ((6 + n * (5 - n * (3 + 2 * n))) / (-18))**0.5
h[0] = x
h[1] = x * x - 1
for i in range(2, 15):
#8
h[i] = x * h[i - 1] - (i) * h[i - 2]
#Probabilities calculated by modified Edgeworth series for n greater than 8
r = 1 / n
sc = r * (h[2] * (-0.09 + r * (0.045 + r * (-0.5325 + r * 0.506))) + r * (h[4] * (0.036735 + r * (-0.036735 + r * 0.3214)) + h[6] * (0.00405 + r * (-0.023336 + r * 0.07787)) + r * (h[8] * (-0.0033061 - r * 0.0065166) + h[10] * (-0.0001215 + r * 0.0025927) + r * (h[12] * 0.00014878 + h[14] * 0.0000027338))))
PRTAUS = 1 - NormalDist().cdf(x) + sc * 0.398942 * exp(-0.5 * x * x)
if (PRTAUS < 0):
PRTAUS = 0
elif (PRTAUS > 1):
PRTAUS = 1
he_AS71 = PRTAUS
else:
#Probabilities calculated by recurrence relation for N less than 9
if (S < 0):
m = m - 2
IM = m / 2 + 1
L.iloc[1, 1] = 1
L.iloc[2, 1] = 1
if (IM >= 2):
for i in range(2, int(IM)+1):
L.iloc[1, i] = 0
L.iloc[2, i] = 0
#2
IL = 1
i = 1
m = 1
j = 1
jj = 2
#3
GOTO3 = True
while (GOTO3):
if (i == n):
GOTO3 = False
else:
GOTO3 = False
IL = IL + i
i = i + 1
m = m * i
j = 3 - j
jj = 3 - jj
INN = 1
IO = 0
k = min(IM, IL)
#4
GOTO4 = True
while (GOTO4):
GOTO4 = False
INN = INN + 1
if (INN > k):
GOTO3 = True
else:
L.iloc[jj, INN] = L.iloc[jj, INN - 1] + L.iloc[j, INN]
if (INN <= i):
GOTO4 = True
else:
IO = IO + 1
L.iloc[jj, INN] = L.iloc[jj, INN] - L.iloc[j, IO]
GOTO4 = True
#5
k = 0
for i in range(1, int(IM)+1):
#6
k = k + L.iloc[jj, i]
PRTAUS = k / m
if (S < 0):
PRTAUS = 1 - PRTAUS
he_AS71 = PRTAUS
return he_AS71
def he_kendall(n, c):
c = int(min(c, int(n * (n - 1) / 2) - c))
newt = []
subt = []
if (n == 1):
Prob = 1
elif (n == 2):
Prob = 1
elif (c == 0):
if (n < 171):
Prob = 2 / factorial(n)
else:
Prob = 0
elif (c == 1):
if (n < 172):
Prob = 2 / factorial(n - 1)
else:
Prob = 0
elif (4 * c == n * (n - 1)):
Prob = 1
elif (n < 10):
newt = [0]*(c+1)
subt = [0]*(c+1)
newt[0] = 1
newt[1] = 1
for j in range(2, n):
cf = 0
for ff in range(0, c+1):
cf = cf + newt[ff]
newt[ff] = cf
if (j <= c):
for i in range(0, c - j + 1):
subt[i] = newt[i]
for i in range(j + 1, c + 1):
ff = 1
newt[i] = newt[i] - subt[ff]
ff = ff + 1
snew=0
for i in range(0, c + 1):
snew = snew + newt[i]
Prob = 2 * snew / factorial(n)
else:
newt = [0]*(c+1)
subt = [0]*(c+1)
newt[0] = 1
newt[1] = 1
for j in range(2, n):
cf = 0
for ff in range(0, c+1):
cf = cf + newt[ff] / j
newt[ff] = cf
if (j <= c):
for i in range(0, c - j + 1):
subt[i] = newt[i]
for i in range(j + 1, c + 1):
ff = 1
newt[i] = newt[i] - subt[ff]
ff = ff + 1
snew=0
for i in range(0, c + 1):
snew = snew + newt[i]
Prob = snew
return ProbFunctions
def di_kendall_tau(n, tau, method='as71', cc=False)-
Expand source code
def di_kendall_tau(n, tau, method="as71", cc=False): if method=="as71": S = n * (n - 1) / 2 * abs(tau) if cc: S = abs(S) - 1 pvalue = 2 * he_AS71(S, n) elif method=="kendall": c = (tau + 1) * (n * (n - 1)) / 4 pvalue = he_kendall(n, c) return pvalue def he_AS71(S, n)-
Expand source code
def he_AS71(S, n): #Upper tail probabilities of tau h = [0]*15 L = [[0]*16]*3 L = pd.DataFrame(L) #Check validity of IS and N values PRTAUS = 1 if (n < 1): he_AS71 = PRTAUS else: m = n * (n - 1) / 2 - abs(S) if ((m < 0) or (m > (m / 2) * 2)): he_AS71 = PRTAUS elif m == 0 and S <= 0: he_AS71 = PRTAUS else: if (n > 8): #GOTO 7 #7 #Calculation of Tchebycheff-Hermite polynomials x = (S - 1) / ((6 + n * (5 - n * (3 + 2 * n))) / (-18))**0.5 h[0] = x h[1] = x * x - 1 for i in range(2, 15): #8 h[i] = x * h[i - 1] - (i) * h[i - 2] #Probabilities calculated by modified Edgeworth series for n greater than 8 r = 1 / n sc = r * (h[2] * (-0.09 + r * (0.045 + r * (-0.5325 + r * 0.506))) + r * (h[4] * (0.036735 + r * (-0.036735 + r * 0.3214)) + h[6] * (0.00405 + r * (-0.023336 + r * 0.07787)) + r * (h[8] * (-0.0033061 - r * 0.0065166) + h[10] * (-0.0001215 + r * 0.0025927) + r * (h[12] * 0.00014878 + h[14] * 0.0000027338)))) PRTAUS = 1 - NormalDist().cdf(x) + sc * 0.398942 * exp(-0.5 * x * x) if (PRTAUS < 0): PRTAUS = 0 elif (PRTAUS > 1): PRTAUS = 1 he_AS71 = PRTAUS else: #Probabilities calculated by recurrence relation for N less than 9 if (S < 0): m = m - 2 IM = m / 2 + 1 L.iloc[1, 1] = 1 L.iloc[2, 1] = 1 if (IM >= 2): for i in range(2, int(IM)+1): L.iloc[1, i] = 0 L.iloc[2, i] = 0 #2 IL = 1 i = 1 m = 1 j = 1 jj = 2 #3 GOTO3 = True while (GOTO3): if (i == n): GOTO3 = False else: GOTO3 = False IL = IL + i i = i + 1 m = m * i j = 3 - j jj = 3 - jj INN = 1 IO = 0 k = min(IM, IL) #4 GOTO4 = True while (GOTO4): GOTO4 = False INN = INN + 1 if (INN > k): GOTO3 = True else: L.iloc[jj, INN] = L.iloc[jj, INN - 1] + L.iloc[j, INN] if (INN <= i): GOTO4 = True else: IO = IO + 1 L.iloc[jj, INN] = L.iloc[jj, INN] - L.iloc[j, IO] GOTO4 = True #5 k = 0 for i in range(1, int(IM)+1): #6 k = k + L.iloc[jj, i] PRTAUS = k / m if (S < 0): PRTAUS = 1 - PRTAUS he_AS71 = PRTAUS return he_AS71 def he_kendall(n, c)-
Expand source code
def he_kendall(n, c): c = int(min(c, int(n * (n - 1) / 2) - c)) newt = [] subt = [] if (n == 1): Prob = 1 elif (n == 2): Prob = 1 elif (c == 0): if (n < 171): Prob = 2 / factorial(n) else: Prob = 0 elif (c == 1): if (n < 172): Prob = 2 / factorial(n - 1) else: Prob = 0 elif (4 * c == n * (n - 1)): Prob = 1 elif (n < 10): newt = [0]*(c+1) subt = [0]*(c+1) newt[0] = 1 newt[1] = 1 for j in range(2, n): cf = 0 for ff in range(0, c+1): cf = cf + newt[ff] newt[ff] = cf if (j <= c): for i in range(0, c - j + 1): subt[i] = newt[i] for i in range(j + 1, c + 1): ff = 1 newt[i] = newt[i] - subt[ff] ff = ff + 1 snew=0 for i in range(0, c + 1): snew = snew + newt[i] Prob = 2 * snew / factorial(n) else: newt = [0]*(c+1) subt = [0]*(c+1) newt[0] = 1 newt[1] = 1 for j in range(2, n): cf = 0 for ff in range(0, c+1): cf = cf + newt[ff] / j newt[ff] = cf if (j <= c): for i in range(0, c - j + 1): subt[i] = newt[i] for i in range(j + 1, c + 1): ff = 1 newt[i] = newt[i] - subt[ff] ff = ff + 1 snew=0 for i in range(0, c + 1): snew = snew + newt[i] Prob = snew return Prob