From 48e10fc7f2464771e42cc39244cf3ddc4b8c8901 Mon Sep 17 00:00:00 2001 From: mp360288 Date: Tue, 27 Aug 2024 17:00:10 +0200 Subject: [PATCH] Add files via upload --- README.md | 3 ++ test.py | 135 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ voap.py | 112 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 250 insertions(+) create mode 100644 README.md create mode 100644 test.py create mode 100644 voap.py diff --git a/README.md b/README.md new file mode 100644 index 0000000..9a9f1da --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# VoAPython +Implementation of tools proposed in the paper "Validation of Association" by B. Ćmiel and T. Ledwina +https://arxiv.org/pdf/1904.06519 \ No newline at end of file diff --git a/test.py b/test.py new file mode 100644 index 0000000..147a61a --- /dev/null +++ b/test.py @@ -0,0 +1,135 @@ +import numpy as np +import voap + +if __name__ == "__main__": + true_copula = np.array( + [ + -0.1505847, + -0.2742042, + -0.3742406, + -0.4714045, + -0.5741693, + 0.69006556, + 0.5741693, + 0.4714045, + 0.3742406, + 0.2742042, + 0.1505847, + -0.2742042, + -0.4993070, + -0.6814663, + 1.1200012, + 0.8286734, + 1.25656172, + 1.0455226, + 0.8583951, + 0.6814663, + 0.4993070, + 0.2742042, + -0.3742406, + -0.6814663, + -0.9300817, + 0.4485395, + 0.1078143, + 0.20579830, + 1.4269545, + 1.1715584, + 0.9300817, + 0.6814663, + 0.3742406, + -0.4714045, + -0.8583951, + 0.4485395, + 1.4395893, + 0.9643376, + 0.80237742, + 1.7974341, + 1.4757296, + 1.1715584, + 0.8583951, + 0.4714045, + -0.5741693, + -1.0455226, + 0.1078143, + 0.9643376, + 0.4270426, + 0.05847053, + 0.8811133, + 0.4165482, + -0.1078143, + 1.0455226, + 0.5741693, + 0.6900656, + -0.5863955, + 0.2057983, + 0.8023774, + 1.3448223, + 0.63245553, + 1.3448223, + 0.8023774, + 0.2057983, + 1.2565617, + 0.6900656, + 0.5741693, + 1.0455226, + 1.4269545, + 1.7974341, + 2.1892691, + 1.34482230, + 1.7351985, + 0.9643376, + 0.1078143, + 0.8286734, + -0.5741693, + 0.4714045, + 0.8583951, + 1.1715584, + 1.4757296, + 1.7974341, + 0.80237742, + 0.9643376, + 1.4395893, + 0.4485395, + 1.1200012, + -0.4714045, + 0.3742406, + 0.6814663, + 0.9300817, + 1.1715584, + 1.4269545, + 1.71498585, + 1.6425832, + 2.0686374, + 0.8705564, + 1.5173982, + -0.3742406, + 0.2742042, + 0.4993070, + 0.6814663, + 0.8583951, + 1.0455226, + 1.25656172, + 0.8286734, + 1.1200012, + 1.5173982, + 2.1858551, + -0.2742042, + 0.1505847, + 0.2742042, + 0.3742406, + 0.4714045, + 0.5741693, + 0.69006556, + -0.5741693, + -0.4714045, + -0.3742406, + -0.2742042, + -0.1505847, + ] + ) + true_copula = true_copula.reshape((11, 11)) + + r = np.array(range(1, 11)) + s = np.array([3, 6, 2, 9, 4, 1, 7, 5, 8, 10]) + test_copula = voap.calculate_copula_grid(r, s) + print(np.sum(true_copula - test_copula)) diff --git a/voap.py b/voap.py new file mode 100644 index 0000000..affc173 --- /dev/null +++ b/voap.py @@ -0,0 +1,112 @@ +import numpy as np +import numba +# from tqdm.auto import tqdm + + +@numba.njit(parallel=True) +def Q_function(Q_grid, C_grid): + n = len(C_grid) - 1 + z = [] + + for i in numba.prange(len(Q_grid)): + row_idx = int(np.ceil(Q_grid[i, 0] * (n + 1)) - 1) + col_idx = int(np.ceil(Q_grid[i, 1] * (n + 1)) - 1) + if row_idx > n: + row_idx = n + if col_idx > n: + col_idx = n + z.append(C_grid[row_idx, col_idx]) + + return {"x": Q_grid[:, 0], "y": Q_grid[:, 1], "z": np.array(z)} + + +def create_Q_grid(n=10): + u = np.linspace(0.5 / (n + 1), (0.5 + n) / (n + 1), n + 1) + grid = np.stack(np.meshgrid(u, u)).T.reshape(-1, 2) + + return grid + +@numba.njit(parallel=True) +def calculate_copula_part(t, c): + n = len(t) + n_inv = 1 / n # n**-1 ## numba doesn't get n**-1 xD + c = np.zeros(c.shape) + for i in numba.prange(1, int(n / 2) + 2): + for j in numba.prange(int(n / 2) + 2): + # c[i, j] = c[i - 1, j] + (j >= t[i - 1]) * n_inv + c[i, j] = c[0, j] + n_inv * np.sum(j >= t[:i]) + return c + +@numba.njit +def wn(n, c, i, j): + u = (i + 0.5) / (n + 1) + v = (j + 0.5) / (n + 1) + return n**0.5 * (c - u * v) * (u * v * (1 - u) * (1 - v)) ** -0.5 + +@numba.njit +def fill_matrix(n, ks, c, x_prim, y_prim): + sign = (-1) ** (x_prim + y_prim) + + for i in range(int(np.floor((n + 1) / 2)) + 1): + for j in range(int(np.floor((n + 1) / 2)) + 1): + x = n - i if x_prim else i + y = n - j if y_prim else j + ks[x, y] = sign * wn(n, c[i, j], i, j) + return ks + +@numba.njit +def arma_copula(rx, ry): + n = len(rx) + + ctab = np.zeros((n + 1, n + 1)) + ctabs22 = np.zeros((n + 1, n + 1)) + ctabs12 = np.zeros((n + 1, n + 1)) + ctabs21 = np.zeros((n + 1, n + 1)) + + rsx = np.zeros(n) + rsy = np.zeros(n) + + ks = np.zeros((n + 1, n + 1)) + + rsx = len(rx) + 1 - rx + rsy = len(ry) + 1 - ry + + t = ry[np.argsort(rx)] + ts22 = rsy[np.argsort(rsx)] + ts12 = rsy[np.argsort(rx)] + ts21 = ry[np.argsort(rsx)] + + ctab = calculate_copula_part(t, ctab) + ctabs22 = calculate_copula_part(ts22, ctabs22) + ctabs12 = calculate_copula_part(ts12, ctabs12) + ctabs21 = calculate_copula_part(ts21, ctabs21) + + ks = fill_matrix(n, ks, ctab, False, False) + ks = fill_matrix(n, ks, ctabs22, True, True) + ks = fill_matrix(n, ks, ctabs12, False, True) + ks = fill_matrix(n, ks, ctabs21, True, False) + + return ks + +@numba.njit +def calculate_copula_grid(x, y): + return arma_copula( + np.argsort(np.argsort(x)) + 1, np.argsort(np.argsort(y)) + 1 + ) + +# @numba.njit(parallel=True) +def calculate_copula_mc_grid(x, y, mc=100, seed=0): + rng = np.random.RandomState(seed) + k = len(x) + g = np.zeros((k + 1, len(y) + 1)) + + for i in range(mc): + indx = rng.choice(k, size=k) + mat = calculate_copula_grid(x[indx], y[indx]) + g += mat / mc + + return g + + +# def create_Q_plot(x, y): +# return Q_function(create_Q_grid(), calculate_copula_mc_grid(x, y))