{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "xicorr.ipynb",
"provenance": [],
"collapsed_sections": [],
"mount_file_id": "1ij5R4zzazK9Ji3J2SAa2RUSZ2b1hkyI4",
"authorship_tag": "ABX9TyPzAIixJDAhWWQYfGcpa942",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/github/morteza/notebooks/blob/master/xicorr.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"source": [
"import numpy as np\n",
"\n",
"import seaborn as sns; sns.set()"
],
"metadata": {
"id": "-4vX69aFhGlb"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "-1bLlqB6hDjQ"
},
"outputs": [],
"source": [
"\n",
"def xicorr(X,Y):\n",
" n = X.size\n",
" xi = np.argsort(X,kind='quicksort')\n",
" Y = Y[xi]\n",
" _,b,c = np.unique(Y,return_counts=True,return_inverse=True)\n",
" r = np.cumsum(c)[b]\n",
" _,b,c = np.unique(-Y,return_counts=True,return_inverse=True)\n",
" l = np.cumsum(c)[b]\n",
" return 1 - n*np.abs(np.diff(r)).sum() / (2*(l*(n-l)).sum())\n",
"\n",
"\n",
"def xicorr2(x: 'ArrayLike', y: 'ArrayLike') -> float:\n",
" \"\"\"xi correlation coefficient (alternative implementation)\n",
"\n",
" Written by github/atarashansky, modified by github/escherba\n",
" https://github.com/czbiohub/xicor/issues/17#issue-965635013\n",
" \"\"\"\n",
" x = np.asarray(x)\n",
" y = np.asarray(y)\n",
" n = x.size\n",
" assert y.size == n, \"arrays must be of the same size\"\n",
" y = y[np.argsort(x, kind='quicksort')]\n",
" _, inverse, counts = np.unique(y, return_inverse=True, return_counts=True)\n",
" right = np.cumsum(counts)[inverse]\n",
" left = np.cumsum(np.flip(counts))[(counts.size - 1) - inverse]\n",
" return 1. - 0.5 * float(np.abs(np.diff(right)).sum() / np.mean(left * (n - left)))\n",
"\n",
"\n",
"def xicorr_orig(x: 'ArrayLike', y: 'ArrayLike') -> float:\n",
" \"\"\"Original R implementation of Xi translated into Python\n",
"\n",
" R implementation:\n",
" https://github.com/cran/XICOR/blob/master/R/calculateXI.R\n",
" \"\"\"\n",
"\n",
" from scipy.stats import rankdata\n",
"\n",
" x = np.asarray(x)\n",
" y = np.asarray(y)\n",
" n = x.size\n",
" assert y.size == n, \"arrays must be of the same size\"\n",
" PI = rankdata(x, method=\"average\")\n",
" fr = rankdata(y, method=\"average\")\n",
" CU = np.mean((float(n + 1) - fr) * (fr - 1.))\n",
" A1 = np.abs(np.diff(fr[np.argsort(PI, kind=\"quicksort\")])).sum()\n",
" return 1. - 0.5 * float(A1 / CU)"
]
},
{
"cell_type": "code",
"source": [
"X = np.load('/content/drive/MyDrive/dosenbach2007_timeseries.npz')['arr_0']"
],
"metadata": {
"id": "AxqcY-XSqNqL"
},
"execution_count": 88,
"outputs": []
},
{
"cell_type": "code",
"source": [
"X[0][0]"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "_rqDenM3qclG",
"outputId": "9e05be78-6093-4c03-9937-300e896483ff"
},
"execution_count": 99,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"32"
]
},
"metadata": {},
"execution_count": 99
}
]
},
{
"cell_type": "code",
"source": [
"x = np.random.rand(125) + 2\n",
"y = x**2 + x + np.sqrt(x) + 12 + np.random.randn(125)\n",
"\n",
"\n",
"for i in range(len(X)):\n",
" x = X[i][1]\n",
" y = X[i][2]\n",
"\n",
" from scipy.stats import pearsonr\n",
"\n",
" print(f'sub_{i}:', xicorr(x, y), xicorr2(x, y), xicorr_orig(x, y), pearsonr(x, y)[0])\n"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "_lUdB0FxhSGw",
"outputId": "4ea2457d-2f7a-4570-a28d-47c6c6cd9304"
},
"execution_count": 103,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"sub_0: 0.05971582181259605 0.05971582181259605 0.05971582181259605 0.007057805036315024\n",
"sub_1: 0.02188940092165903 0.02188940092165903 0.02188940092165903 0.05709108235559736\n",
"sub_2: 0.11847158218125964 0.11847158218125964 0.11847158218125964 0.35156058295897297\n",
"sub_3: -0.004992319508448473 -0.004992319508448473 -0.004992319508448473 0.007753525326876064\n",
"sub_4: 0.3160522273425499 0.3160522273425499 0.3160522273425499 0.688780863414386\n",
"sub_5: -0.08870967741935476 -0.08870967741935476 -0.08870967741935476 0.04074056728328548\n",
"sub_6: 0.030721966205837226 0.030721966205837226 0.030721966205837226 0.31886534096959906\n",
"sub_7: 0.18682795698924726 0.18682795698924726 0.18682795698924726 0.5295077968729054\n",
"sub_8: 0.081989247311828 0.081989247311828 0.081989247311828 -0.4073567278311721\n",
"sub_9: 0.0643241167434716 0.0643241167434716 0.0643241167434716 0.16035001468148263\n",
"sub_10: 0.004224270353302639 0.004224270353302639 0.004224270353302639 -0.0017708462342737451\n",
"sub_11: 0.035522273425499185 0.035522273425499185 0.035522273425499185 0.0884169650086374\n",
"sub_12: 0.046658986175115214 0.046658986175115214 0.046658986175115214 -0.00015848439682089626\n",
"sub_13: 0.03590629800307221 0.03590629800307221 0.03590629800307221 -0.07561658397244023\n",
"sub_14: -0.005952380952380931 -0.005952380952380931 -0.005952380952380931 0.027196366379864764\n",
"sub_15: 0.037634408602150504 0.037634408602150504 0.037634408602150504 -0.023365355634358476\n",
"sub_16: -0.004992319508448473 -0.004992319508448473 -0.004992319508448473 0.04638615068057758\n",
"sub_17: -0.043010752688172005 -0.043010752688172005 -0.043010752688172005 -0.03483487350881016\n",
"sub_18: 0.049347158218125964 0.049347158218125964 0.049347158218125964 0.24801276531894492\n",
"sub_19: 0.19854070660522272 0.19854070660522272 0.19854070660522272 0.49034516733068445\n",
"sub_20: 0.0038402457757296116 0.0038402457757296116 0.0038402457757296116 0.18363832525018858\n",
"sub_21: 0.028033794162826475 0.028033794162826475 0.028033794162826475 0.32805666138666145\n",
"sub_22: 0.022849462365591378 0.022849462365591378 0.022849462365591378 0.3428984291846277\n",
"sub_23: 0.0005760368663594306 0.0005760368663594306 0.0005760368663594306 0.37827926330590994\n",
"sub_24: 0.019009216589861766 0.019009216589861766 0.019009216589861766 -0.014662758255547994\n",
"sub_25: 0.0071044546850997925 0.0071044546850997925 0.0071044546850997925 0.31079563307701513\n",
"sub_26: 0.092741935483871 0.092741935483871 0.092741935483871 0.34952238435895855\n",
"sub_27: 0.20852534562211977 0.20852534562211977 0.20852534562211977 0.6566733607025322\n",
"sub_28: -0.033986175115207296 -0.033986175115207296 -0.033986175115207296 -0.13475542316468594\n",
"sub_29: 0.06682027649769584 0.06682027649769584 0.06682027649769584 0.44675728779363655\n",
"sub_30: 0.05030721966205842 0.05030721966205842 0.05030721966205842 0.4287767220184844\n",
"sub_31: 0.006144393241167445 0.006144393241167445 0.006144393241167445 0.16001212766044476\n"
]
}
]
}
]
}
