-
Notifications
You must be signed in to change notification settings - Fork 1
/
statistic.py
89 lines (78 loc) · 2.8 KB
/
statistic.py
1
2
3
4
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
78
79
80
81
82
83
84
85
86
87
88
89
import numpy as np
def compute_dHSIC_statistics(k_list):
"""
Computes the dHSIC statistic
:param k_list: list of kernel matrices for each variable
:return: single value for dHSIC test statistic
"""
n_variables = len(k_list)
n_samples = k_list[0].shape[0]
term1 = k_list[0]
term2 = np.sum(k_list[0]) / (n_samples ** 2)
term3 = 2 * np.sum(k_list[0], axis=0) / (n_samples ** 2)
for j in range(1, n_variables):
term1 = term1 * k_list[j]
term2 = term2 * np.sum(k_list[j]) / (n_samples ** 2)
term3 = term3 * np.sum(k_list[j], axis=0) / n_samples
term1_sum = np.sum(term1)
term3_sum = np.sum(term3)
dHSIC = term1_sum / (n_samples ** 2) + term2 - term3_sum
return dHSIC
def compute_3way_statistics(k_list):
"""
Computes the 3-way Lancaster/Streitberg test statistic
:param k_list: list of kernel matrices for each variable
:return: single value for 3-way interaction test statistic
"""
K = k_list[0]
L = k_list[1]
M = k_list[2]
m = np.shape(K)[0]
H = np.eye(m) - 1 / m * np.ones(m)
Kc = H @ K @ H
Lc = H @ L @ H
Mc = H @ M @ H
statMatrix = Kc * Lc * Mc
threeway = 1 / (m ** 2) * np.sum(statMatrix)
return threeway
def compute_lancaster_4way_statistics(K):
"""
Computes the 4-way Lancaster test statistic
:param K: list of kernel matrices for each variable
:return: single value for 4-way Lancaster test statistic
"""
n = K[0].shape[0]
statMatrix = np.einsum('ij,ij,ij,ij->', K[0], K[1], K[2], K[3])
stat = 1 / (n ** 2) * np.sum(statMatrix)
return stat
def compute_streitberg_4way_statistics(K):
"""
Computes the 4-way Streitberg test statistic
:param K: list of kernel matrices for each variable
:return: single value for 4-way Streitberg test statistic
"""
n = K[0].shape[0]
# p1234,p1234
ip1 = np.einsum('ij,ij,ij,ij->', K[0], K[1], K[2], K[3])
# p1234, p12p34
ip2 = np.einsum('ij,ij,ik,ik->', K[0], K[1], K[2], K[3])
# p1234, p13p24
ip3 = np.einsum('ij,ij,ik,ik->', K[0], K[2], K[1], K[3])
# p1234, p14p23
ip4 = np.einsum('ij,ij,ik,ik->', K[0], K[3], K[1], K[2])
# p12p34, p12p34
ip5 = np.einsum('ij,ij,kl,kl->', K[0], K[1], K[2], K[3])
# p12p34, p13p24
ip6 = np.einsum('ij,il,kj,kl->', K[0], K[1], K[2], K[3])
# p12p34, p14p23
ip7 = np.einsum('ij,il,kl,kj->', K[0], K[1], K[2], K[3])
# p13p24, p13p24
ip8 = np.einsum('ij,ij,kl,kl->', K[0], K[2], K[1], K[3])
# p13p24, p14p23
ip9 = np.einsum('ij,il,kl,kj->', K[0], K[2], K[1], K[3])
# p14p23, p14p23
ip10 = np.einsum('ij,ij,kl,kl->', K[0], K[3], K[1], K[2])
# sum up
stat = 1 / (n ** 4) * (n ** 2 * ip1 - 2 * n * (ip2 + ip3 + ip4)
+ ip5 + ip8 + ip10 + 2 * (ip6 + ip7 + ip9))
return stat