ksample/ksamp_power_sampsize.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "import multiprocessing as mp\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from power_ksamp import power_ksamp_sample\n",
    "from hyppo.independence import CCA, Dcorr, HHG, Hsic, RV, MGC, KMERF\n",
    "from hyppo.tools import *\n",
    "\n",
    "sys.path.append(os.path.realpath('..'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "sns.set(color_codes=True, style='white', context='talk', font_scale=1.5)\n",
    "PALETTE = sns.color_palette(\"Set1\")\n",
    "sns.set_palette(PALETTE[1:5] + PALETTE[6:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are some constants that are used in this notebook. If running these notebook, please only manipulate these constants if you are not running more tests. They define the sample sizes tested upon and the number of replications. The simulations tested over and the independence tests tested over are defined also."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rpy2.robjects import Formula, numpy2ri\n",
    "from rpy2.robjects.packages import importr\n",
    "\n",
    "\n",
    "class Manova:\n",
    "    r\"\"\"\n",
    "    Wrapper of R MANOVA\n",
    "    \"\"\"\n",
    "    def __init__(self):\n",
    "        self.stats = importr('stats')\n",
    "        self.r_base = importr('base')\n",
    "        \n",
    "        numpy2ri.activate()\n",
    "\n",
    "        self.formula = Formula('X ~ Y')\n",
    "        self.env = self.formula.environment\n",
    "\n",
    "    def statistic(self, x, y):\n",
    "        r\"\"\"\n",
    "        Helper function to calculate the test statistic\n",
    "        \"\"\"\n",
    "        self.env['Y'] = y\n",
    "        self.env['X'] = x\n",
    "\n",
    "        stat = self.r_base.summary(self.stats.manova(self.formula), test=\"Pillai\")[3][0, 1]\n",
    "\n",
    "        return stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "MAX_SAMPLE_SIZE = 100\n",
    "STEP_SIZE = 5\n",
    "SAMP_SIZES = range(5, MAX_SAMPLE_SIZE + STEP_SIZE, STEP_SIZE)\n",
    "POWER_REPS = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIMULATIONS = {\n",
    "    \"linear\": \"Linear\",\n",
    "    \"exponential\": \"Exponential\",\n",
    "    \"cubic\": \"Cubic\",\n",
    "    \"joint_normal\": \"Joint Normal\",\n",
    "    \"step\": \"Step\",\n",
    "    \"quadratic\": \"Quadratic\",\n",
    "    \"w_shaped\": \"W-Shaped\",\n",
    "    \"spiral\": \"Spiral\",\n",
    "    \"uncorrelated_bernoulli\": \"Bernoulli\",\n",
    "    \"logarithmic\": \"Logarithmic\",\n",
    "    \"fourth_root\": \"Fourth Root\",\n",
    "    \"sin_four_pi\": \"Sine 4\\u03C0\",\n",
    "    \"sin_sixteen_pi\": \"Sine 16\\u03C0\",\n",
    "    \"square\": \"Square\",\n",
    "    \"two_parabolas\": \"Two Parabolas\",\n",
    "    \"circle\": \"Circle\",\n",
    "    \"ellipse\": \"Ellipse\",\n",
    "    \"diamond\": \"Diamond\",\n",
    "    \"multiplicative_noise\": \"Multiplicative\",\n",
    "    \"multimodal_independence\": \"Independence\"\n",
    "}\n",
    "\n",
    "TESTS = [\n",
    "    KMERF,\n",
    "    MGC,\n",
    "    Dcorr,\n",
    "    Hsic,\n",
    "    Manova,\n",
    "    HHG,\n",
    "    CCA,\n",
    "    RV,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function calculates the estimated power ``POWER_REPS`` number off times and averages them. It does this iterating over the number of sample sizes.\n",
    "\n",
    "**Note: We only recommend running this code if running the next 2 cells ONCE to generate the csv files used to visualize the plots. This code takes a very long time to run and if running, we recommend using a machine with many cores.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "RENAME_TESTS = {\n",
    "    \"Dcorr\": \"Energy\",\n",
    "    \"Hsic\": \"MMD\"\n",
    "}\n",
    "\n",
    "def estimate_power(sim, test):\n",
    "    est_power = np.array([np.mean([power_ksamp_sample(test, rot_ksamp, sim, n=i, noise=True) for _ in range(POWER_REPS)])\n",
    "                          for i in SAMP_SIZES])\n",
    "    np.savetxt('../ksample/ksamp_vs_samplesize/{}_{}.csv'.format(sim, test.__name__),\n",
    "               est_power, delimiter=',')\n",
    "    \n",
    "    return est_power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "outputs = Parallel(n_jobs=-1, verbose=100)(\n",
    "    [delayed(estimate_power)(sim, test) for sim in SIMULATIONS.keys() for test in TESTS]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code loops over each saved independence test file and generates absolute power curves for each test and for each simulation modality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_power():\n",
    "    fig, ax = plt.subplots(nrows=4, ncols=5, figsize=(25,20))\n",
    "    \n",
    "    plt.suptitle(\"Multivariate Three-Sample Testing Increasing Sample Size\", y=0.93, va='baseline')\n",
    "    \n",
    "    for i, row in enumerate(ax):\n",
    "        for j, col in enumerate(row):\n",
    "            count = 5*i + j\n",
    "            sim = list(SIMULATIONS.keys())[count]\n",
    "            \n",
    "            for test in TESTS:\n",
    "                test_name = test.__name__\n",
    "                power = np.genfromtxt('../ksample/ksamp_vs_samplesize/{}_{}.csv'.format(sim, test_name),\n",
    "                                      delimiter=',')\n",
    "#                 manova_power = np.genfromtxt('../ksample/ksamp_vs_samplesize/{}_Manova.csv'.format(sim), delimiter=',')\n",
    "\n",
    "                if test_name == \"MGC\":\n",
    "                    col.plot(SAMP_SIZES, power, color=\"#e41a1c\", label=test_name, lw=4)\n",
    "                elif test_name == \"KMERF\":\n",
    "                    col.plot(SAMP_SIZES, power, color=\"#e41a1c\", label=test_name, lw=4, linestyle='dashed')\n",
    "                elif test_name == \"Manova\":\n",
    "                    col.plot(SAMP_SIZES, power, color=\"#000000\", label=test_name, lw=4)\n",
    "                else:\n",
    "                    col.plot(SAMP_SIZES, power, label=test_name, lw=2)\n",
    "                col.set_xticks([])\n",
    "                if i == 3:\n",
    "                    col.set_xticks([SAMP_SIZES[0], SAMP_SIZES[-1]])\n",
    "                col.set_ylim(-0.05, 1.05)\n",
    "                col.set_yticks([])\n",
    "                if j == 0:\n",
    "                    col.set_yticks([0, 1])\n",
    "                col.set_title(SIMULATIONS[sim])\n",
    "    \n",
    "    fig.text(0.5, 0.08, 'Sample Size', ha='center')\n",
    "    fig.text(0.08, 0.5, 'Statistical Power Relative to Manova', va='center', rotation='vertical')\n",
    "    leg = plt.legend(bbox_to_anchor=(0.5, 0.07), bbox_transform=plt.gcf().transFigure,\n",
    "                     ncol=len(TESTS), loc='upper center')\n",
    "    leg.get_frame().set_linewidth(0.0)\n",
    "    for legobj in leg.legendHandles:\n",
    "        legobj.set_linewidth(5.0)\n",
    "    plt.subplots_adjust(hspace=.50)\n",
    "    plt.savefig('../ksample/figs/ksamp_power_sampsize.pdf', transparent=True, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_power()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}