Fig1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7d21a5af",
   "metadata": {},
   "source": [
    "# Defining distributional shift (DS)\n",
    "\n",
    "We define distributional shift (DS) as the concentration of frequencies towards the lowest discrete class, which we measure via the sum of cumulative frequencies.\n",
    "\n",
    "We begin by simply deriving DS as: \n",
    "\n",
    "$$DS = (\\sum{F}/n - 1)/(k-1)$$<br><br>\n",
    "\n",
    "We then refine DS to include exponentiated cumulative frequencies:\n",
    "\n",
    "$$DS = (\\sum{F^{z}}/n^{z} - 1)/(k-1)$$<br><br>\n",
    "\n",
    "Finally, we refine the exponent (*z*) to take fractional values:\n",
    "\n",
    "$$DS = (\\sum{F^{k+1/k}}/n^{k+1/k} - 1)/(k-1)$$<br><br>\n",
    "\n",
    "Below, we provide a function to allow users to explore the calculation of DS and to vary the number of observations (n_obs) and the number of bins (n_bins)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8b748901",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  from pandas import Int64Index as NumericIndex\n"
     ]
    }
   ],
   "source": [
    "\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "import scipy as sc\n",
    "from itertools import product\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import PIL\n",
    "from numpy import log10\n",
    "import random\n",
    "from math import factorial\n",
    "from scipy.stats import linregress, gaussian_kde, skew\n",
    "from scipy import stats\n",
    "from scipy.spatial import distance\n",
    "import warnings\n",
    "import pandas as pd\n",
    "import re\n",
    "import os\n",
    "import math\n",
    "from collections import Counter\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from statsmodels.stats.outliers_influence import summary_table\n",
    "import statsmodels.api as sm\n",
    "from scipy.optimize import linear_sum_assignment\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "%config InlineBackend.print_figure_kwargs={'facecolor' : \"w\"}\n",
    "\n",
    "pd.set_option('display.max_columns', None)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "abecec2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_cum_dists(n_obs, n_bins, cum=True):\n",
    "    \"\"\"Generate all possible discrete distributions for a given number of \n",
    "        observations distributed across a given number of bins\n",
    "    \"\"\"\n",
    "    \n",
    "    def partitions(n, k):\n",
    "        if k == 1:\n",
    "            yield (n,)\n",
    "        else:\n",
    "            for i in range(n + 1):\n",
    "                for result in partitions(n - i, k - 1):\n",
    "                    yield (i,) + result\n",
    "\n",
    "    for combo in partitions(n_obs, n_bins):\n",
    "        if cum:\n",
    "            cum_dist = [sum(combo[:i + 1]) for i in range(n_bins)]\n",
    "            yield cum_dist\n",
    "        else:\n",
    "            yield combo\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cd4b0bf5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11 (0, 1, 4, 0) (0, 1, 5, 5) 0.377916\n",
      "14 (0, 2, 2, 1) (0, 2, 4, 5) 0.358233\n",
      "16 (0, 3, 0, 2) (0, 3, 3, 5) 0.352044\n",
      "25 (1, 0, 3, 1) (1, 1, 4, 5) 0.341363\n",
      "28 (1, 1, 1, 2) (1, 2, 3, 5) 0.326641\n",
      "37 (2, 0, 0, 3) (2, 2, 2, 5) 0.318108\n",
      "\n",
      "\n",
      "56 members of the feasible set for n and k\n"
     ]
    }
   ],
   "source": [
    "n = 5\n",
    "k = 4\n",
    "\n",
    "# Generate cumulative distributions one at a time and add to a list\n",
    "_set = []\n",
    "for dist in generate_cum_dists(n, k, cum=False):\n",
    "    _set.append(dist)\n",
    "\n",
    "\n",
    "DS_ls = []\n",
    "for i, d in enumerate(_set):\n",
    "\n",
    "    # get cumulative distribution\n",
    "    cd = np.array([sum(d[:i+1]) for i in range(len(d))])\n",
    "    \n",
    "    # exponentiate cumulative frequencies\n",
    "    z = 1\n",
    "    z = (k + 1)/k\n",
    "    \n",
    "    DS = str(((sum(cd**z)/(n**z)) - 1) / (k - 1))\n",
    "    DS = DS[:8]\n",
    "    \n",
    "    if sum(cd)/(n) == 2.2:    \n",
    "        #print(i+1, d, tuple(cd), sum(cd**z)/(n**z))\n",
    "        print(i+1, d, tuple(cd), DS)\n",
    "\n",
    "    DS_ls.append(DS)\n",
    "    \n",
    "print('\\n')\n",
    "print(len(DS_ls), 'members of the feasible set for n and k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "385170d4",
   "metadata": {},
   "source": [
    "# Figure 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "996bffc4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n: 10 k: 5\n",
      "cardinality: 1001 | no. of unique ΣF/n: 60 | no. of unique ΣF^z/n^z: 355 \n",
      "\n",
      "n: 10 k: 5\n",
      "cardinality: 1001 | no. of unique ΣF/n: 60 | no. of unique ΣF^z/n^z: 956 \n",
      "\n",
      "n: 10 k: 5\n",
      "cardinality: 1001 | no. of unique ΣF/n: 60 | no. of unique ΣF^z/n^z: 1001 \n",
      "\n",
      "\n",
      "\n",
      "n: 10 k: 5\n",
      "cardinality: 1001 | no. of unique ΣF/n: 60 | no. of unique ΣF^z/n^z: 1001 \n",
      "\n",
      "n: 10 k: 10\n",
      "cardinality: 92378 | no. of unique ΣF/n: 208 | no. of unique ΣF^z/n^z: 92378 \n",
      "\n",
      "n: 20 k: 10\n",
      "cardinality: 10015005 | no. of unique ΣF/n: 491 | no. of unique ΣF^z/n^z: 10015005 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "n_obs = 10\n",
    "n_bins = 5\n",
    "fig = plt.figure(figsize=(10.5, 7))\n",
    "   \n",
    "################  TOP ROW  ###########\n",
    "\n",
    "_set = []\n",
    "for dist in generate_cum_dists(n_obs, n_bins, cum=False):\n",
    "    _set.append(dist)\n",
    "\n",
    "plot_num = 1\n",
    "for z in [2, 4, 7]:\n",
    "    \n",
    "    t1 = []\n",
    "    t2 = []\n",
    "    for d in _set:\n",
    "\n",
    "        cd = [sum(d[:i+1]) for i in range(len(d))]\n",
    "        ncd = np.array(cd)/max(cd)\n",
    "        t1.append(sum(ncd))\n",
    "\n",
    "        G = [sum(d[:i+1])**(z) for i in range(len(d))]\n",
    "        G_ = np.array(G)/(n_obs**z)\n",
    "\n",
    "        t2.append(sum(G_))\n",
    "\n",
    "    print('n:', n_obs, 'k:', n_bins)\n",
    "    print('cardinality:', len(t1), '| no. of unique ΣF/n:', \n",
    "          len(list(set(t1))), '| no. of unique ΣF^z/n^z:', len(list(set(t2))), '\\n')\n",
    "\n",
    "    ax = plt.subplot(2, 3, plot_num)\n",
    "    plt.scatter(t1, t2, s=1, c='k')\n",
    "    plt.xlabel(r'$\\sum{F/n}$', fontsize= 14)\n",
    "    plt.ylabel(r'$\\sum{F^{' + str(z) + '}/n^{' + str(z) + '}}$', fontsize= 14)\n",
    "    \n",
    "    s = '|A' + r'$_{n' + '=' + str(n_obs) + ', ' + 'k' + '=' + str(n_bins) + '}$' + '| = ' + str(len(_set)) + '\\n'\n",
    "    s += '\\nValues of ' + r'$\\sum{F^{' + str(z) + '}/n^{' + str(z) + '}}$' + '\\n = ' + str(len(list(set(t2))))\n",
    "    \n",
    "    plt.text(1.01, 3.65, s, fontsize=10)\n",
    "    plt.tick_params(axis='both', labelsize=10)\n",
    "    plot_num += 1\n",
    "    \n",
    "print('\\n')\n",
    "\n",
    "################  BOTTOM ROW  ###########  \n",
    "\n",
    "N_obs = [10, 10, 20]\n",
    "N_bins = [5, 10, 10]\n",
    "sets = []\n",
    "for i, n_obs in enumerate(N_obs):\n",
    "    n_bins = N_bins[i]\n",
    "    \n",
    "    _set = []\n",
    "    for dist in generate_cum_dists(n_obs, n_bins, cum=False):\n",
    "        _set.append(dist)\n",
    "    sets.append(_set)\n",
    "\n",
    "\n",
    "for i, n_obs in enumerate(N_obs):\n",
    "    n_bins = N_bins[i]\n",
    "    \n",
    "    _set = sets[i]\n",
    "    \n",
    "    t1 = []\n",
    "    t2 = []\n",
    "    for d in _set:\n",
    "\n",
    "        cd = [sum(d[:ii+1]) for ii in range(len(d))]\n",
    "        ncd = np.array(cd)/max(cd)\n",
    "        t1.append(sum(ncd))\n",
    "\n",
    "        z = (n_bins + 1)/n_bins\n",
    "        G = [sum(d[:ii+1])**(z) for ii in range(len(d))]\n",
    "        G_ = np.array(G)/(n_obs**z)\n",
    "        t2.append(sum(G_))\n",
    "\n",
    "    print('n:', n_obs, 'k:', n_bins)\n",
    "    print('cardinality:', len(t1), '| no. of unique ΣF/n:', \n",
    "          len(list(set(t1))), '| no. of unique ΣF^z/n^z:', len(list(set(t2))), '\\n')\n",
    "    \n",
    "    s = '|A' + r'$_{n' + '=' + str(n_obs) + ', ' + 'k' + '=' + str(n_bins) + '}$' + '| = ' + str(len(_set)) + '\\n'\n",
    "    s += '\\nValues of ' + r'$\\sum{F_{i}^{' + str(z) + '}/n^{' + str(z) + '}}$' + '\\n = ' + str(len(_set))\n",
    "    \n",
    "    if len(_set) > 10**5: \n",
    "        indices = np.random.choice(len(_set), 10**5, replace=False)\n",
    "        t1 = np.array(t1)\n",
    "        t2 = np.array(t2)\n",
    "        t1 = t1[indices]\n",
    "        t2 = t2[indices]\n",
    "        t1 = t1.tolist()\n",
    "        t2 = t2.tolist()\n",
    "    \n",
    "    ax = plt.subplot(2, 3, plot_num)\n",
    "    plt.scatter(t1, t2, s=1, c='k')\n",
    "    plt.xlabel(r'$\\sum{F/n}$', fontsize= 14)\n",
    "    plt.ylabel(r'$\\sum{F^{' + str(z) + '}/n^{' + str(z) + '}}$', fontsize= 14)\n",
    "    \n",
    "    if plot_num == 4:\n",
    "        plt.text(1. * min(t1), 0.73*max(t2), s, fontsize=10)\n",
    "    elif plot_num == 5:\n",
    "        plt.text(1. * min(t1), 0.7*max(t2), s, fontsize=10)\n",
    "    elif plot_num == 6:\n",
    "        plt.text(0.9 * min(t1), 0.71*max(t2), s, fontsize=10)\n",
    "        \n",
    "    plt.tick_params(axis='both', labelsize=10)\n",
    "    plot_num += 1\n",
    "    \n",
    "\n",
    "fig.patch.set_facecolor('white')\n",
    "plt.subplots_adjust(hspace=0.45, wspace=0.4)\n",
    "plt.savefig('Final_Figs/manuscript/Fig1.pdf', bbox_inches='tight', format='pdf', dpi=600)\n",
    "plt.savefig('Final_Figs/manuscript/Fig1.jpg', bbox_inches='tight', format='jpg', dpi=600)\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}