-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathtest.py
246 lines (223 loc) · 10.8 KB
/
test.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
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
import numpy as np
from preprocessing import compute_kernel_mats, centering_kernel_mats
from statistic import compute_streitberg_4way_statistics, compute_dHSIC_statistics, compute_streitberg_5way_statistics, \
compute_lancaster_4way_statistics, compute_3way_statistics
def permute_K(K, n_perms, seed, perm_var, composite=False):
"""
this function permutes the kernel matrices of the perm_var
:param composite: if True, the permutation is applied to the composite subhypotheses
:param K: list of kernel matrices
:param n_perms: number of permutations
:param seed: seed
:param perm_var: list of variable indices to be permuted,
the permutation is the same for all the indices in this list
:return: K_perms, list of permuted kernel matrices
"""
K_perms = []
for seed in np.random.SeedSequence(seed).spawn(n_perms):
K_perm = np.zeros_like(K)
n_var = K.shape[0]
if composite:
# this is used when composite independence is tested
rng = np.random.default_rng(seed)
# variable ordering after permutation
var_perm = rng.permutation(K.shape[1])
for i in range(n_var):
if i in perm_var:
K_perm[i, :, :] = K[i, var_perm, var_perm[:, None]]
else:
K_perm[i, :, :] = K[i, :, :]
K_perms.append(K_perm)
else:
# this is used when joint independence is tested
for d, sub_seed in enumerate(seed.spawn(n_var)):
rng = np.random.default_rng(sub_seed)
var_perm = rng.permutation(K.shape[1])
K_perm[d, :, :] = K[d, var_perm, var_perm[:, None]]
K_perms.append(K_perm)
return K_perms
def reject_H0(K0, K_perms, stat_fun, alpha=0.05):
"""
Approximates the null distribution by permutation. Using Monte Carlo approximation.
:param K0: kernel matrix of the original data
:param K_perms: list of kernel matrices of the permuted data
:param stat_fun: function to compute the test statistic
:param alpha: significance level
:return: 1 if H0 is rejected, 0 otherwise
"""
s0 = stat_fun(K0)
stats = list(map(stat_fun, K_perms))
return int((1 + sum(s >= s0 for s in stats)) / (1 + len(stats)) < alpha)
def pairwise_test(X, n_perms, seed=0, alpha=0.05):
"""
Tests pairwise independence of the data.
:param X: list of data matrices of shape
:param n_perms: number of permutations
:param seed: seed
:param alpha: significance level
:return: 1 if H0 is rejected, 0 otherwise
"""
K = compute_kernel_mats(X)
K_perms = permute_K(K, n_perms, seed, None)
reject = reject_H0(K, K_perms, compute_dHSIC_statistics, alpha=alpha)
return reject
def composite_4way_factorisability_test_dHSIC(X, n_perms, seed=0, alpha=0.05):
"""
Tests composite 4-way factorisability of the data.
:param X: list of data matrices of shape
:param n_perms: number of permutations
:param seed: seed
:param alpha: significance level
:return: 1 if H0 is rejected, 0 otherwise
"""
K = compute_kernel_mats(X)
reject7 = 0
# test1: p1234 = p1p234
K[1] = K[1] * K[2] * K[3]
K_perms = permute_K(K[:2], n_perms, seed, None, composite=False)
reject1 = reject_H0(K[:2], K_perms, compute_dHSIC_statistics, alpha=alpha)
# do next test only if the first test is rejected
if reject1:
# test2: p1234 = p2p134
K = compute_kernel_mats(X)
K[0] = K[0] * K[2] * K[3]
K_perms = permute_K(K[:2], n_perms, seed, None, composite=False)
reject2 = reject_H0(K[:2], K_perms, compute_dHSIC_statistics, alpha=alpha)
# do next test only if the second test is rejected
if reject2:
# test3: p1234 = p3p124
K = compute_kernel_mats(X)
K[3] = K[0] * K[1] * K[3]
K_perms = permute_K(K[2:], n_perms, seed, None, composite=False)
reject3 = reject_H0(K[2:], K_perms, compute_dHSIC_statistics, alpha=alpha)
# do next test only if the third test is rejected
if reject3:
# test4: p1234 = p4p123
K = compute_kernel_mats(X)
K[2] = K[0] * K[1] * K[2]
K_perms = permute_K(K[2:], n_perms, seed, None, composite=False)
reject4 = reject_H0(K[2:], K_perms, compute_dHSIC_statistics, alpha=alpha)
# do next test only if the fourth test is rejected
if reject4:
# test5: p1234 = p12p34
K = compute_kernel_mats(X)
K[0] = K[0] * K[1]
K[1] = K[2] * K[3]
K_perms = permute_K(K[:2], n_perms, seed, None, composite=False)
reject5 = reject_H0(K[:2], K_perms, compute_dHSIC_statistics, alpha=alpha)
# do next test only if the fifth test is rejected
if reject5:
# test6: p1234 = p13p24
K = compute_kernel_mats(X)
K[0] = K[0] * K[2]
K[1] = K[1] * K[3]
K_perms = permute_K(K[:2], n_perms, seed, None, composite=False)
reject6 = reject_H0(K[:2], K_perms, compute_dHSIC_statistics, alpha=alpha)
# do next test only if the sixth test is rejected
if reject6:
# test7: p1234 = p14p23
K = compute_kernel_mats(X)
K[0] = K[0] * K[3]
K[1] = K[1] * K[2]
K_perms = permute_K(K[:2], n_perms, seed, None, composite=False)
reject7 = reject_H0(K[:2], K_perms, compute_dHSIC_statistics, alpha=alpha)
return reject7
def composite_4way_factorisability_test_streitberg(X, n_perms, seed, alpha=0.05):
"""
Tests composite 4-way factorisability of the data.
:param X: list of data matrices of shape
:param n_perms: number of permutations
:param seed: seed
:param alpha: significance level
:return: 1 if H0 is rejected, 0 otherwise
"""
K = compute_kernel_mats(X)
K = centering_kernel_mats(K)
reject7 = 0
# test1: p1234 = p1p234
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[0], composite=True)
reject1 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# do next test only if the first test is rejected
if reject1:
# test2: p1234 = p2p134
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[1], composite=True)
reject2 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# do next test only if the second test is rejected
if reject2:
# test3: p1234 = p3p124
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[2], composite=True)
reject3 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# do next test only if the third test is rejected
if reject3:
# test4: p1234 = p4p123
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[3], composite=True)
reject4 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# do next test only if the fourth test is rejected
if reject4:
# test5: p1234 = p12p34
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[0, 1], composite=True)
reject5 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# do next test only if the fifth test is rejected
if reject5:
# test6: p1234 = p13p24
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[0, 2], composite=True)
reject6 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# do next test only if the sixth test is rejected
if reject6:
# test7: p1234 = p14p23
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[0, 3], composite=True)
reject7 = reject_H0(K, K_perms, compute_streitberg_4way_statistics, alpha=alpha)
# if reject7:
# print('Passed test')
return reject7
def composite_lancaster_4way(X, n_perms, seed, alpha=0.05):
"""
Composite test for Lancaster's 4-way interaction.
"""
# compute K
K = compute_kernel_mats(X)
# center K
K = centering_kernel_mats(K)
reject4 = 0
# test1: p1234 = p1p234
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[0], composite=True)
reject1 = reject_H0(K, K_perms, compute_lancaster_4way_statistics, alpha=alpha)
# do next test only if the first test is rejected
if reject1:
# test2: p1234 = p2p134
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[1], composite=True)
reject2 = reject_H0(K, K_perms, compute_lancaster_4way_statistics, alpha=alpha)
# do next test only if the second test is rejected
if reject2:
# test3: p1234 = p3p124
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[2], composite=True)
reject3 = reject_H0(K, K_perms, compute_lancaster_4way_statistics, alpha=alpha)
# do next test only if the third test is rejected
if reject3:
# test4: p1234 = p4p123
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[3], composite=True)
reject4 = reject_H0(K, K_perms, compute_lancaster_4way_statistics, alpha=alpha)
return reject4
def composite_3way(X, n_perms, seed, alpha=0.05):
"""
Composite test for 3-way interaction.
"""
# compute K
K = compute_kernel_mats(X)
# center K
K = centering_kernel_mats(K)
reject3 = 0
# test1: p123 = p12p3
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[0], composite=True)
reject1 = reject_H0(K, K_perms, compute_3way_statistics, alpha=alpha)
# do next test only if the first test is rejected
if reject1:
# test2: p123 = p13p2
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[1], composite=True)
reject2 = reject_H0(K, K_perms, compute_3way_statistics, alpha=alpha)
# do next test only if the second test is rejected
if reject2:
# test3: p123 = p23p1
K_perms = permute_K(K=K, n_perms=n_perms, seed=seed, perm_var=[2], composite=True)
reject3 = reject_H0(K, K_perms, compute_3way_statistics, alpha=alpha)
return reject3