Simulation Assuming equal covariance matrix .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_classification.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 359,
     "status": "ok",
     "timestamp": 1632587563990,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "JaKmXGFmdcRW",
    "outputId": "8d45e8a3-8ca3-4abe-b78d-ffa9e4a700db",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "nMC = 1000 #number of runs for MC simulation\n",
    "r20 = np.array([iter_i_error(.2) for i in range(nMC)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 40483,
     "status": "ok",
     "timestamp": 1632581918726,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "g7BKoRRddjyA",
    "outputId": "a25a41a2-2fa5-4e61-d5d9-8824b4121721"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 20% missing rate\n",
      "[0.15014167 0.1339625  0.13375417 0.70351667 0.13585   ]\n",
      "[0.02482887 0.023012   0.02271074 0.16174921 0.02287238]\n"
     ]
    }
   ],
   "source": [
    "    m20 = np.array([np.mean(r20[:,i]) for i in range(5)])\n",
    "    s20 = np.array([np.std(r20[:,i]) for i in range(5)])\n",
    "    print('mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 20% missing rate')\n",
    "    print(m20)\n",
    "    print(s20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "nMC = 1000 #number of runs for MC simulation\n",
    "r35 = np.array([iter_i_error(.35) for i in range(nMC)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 35% missing rate\n",
      "[0.15088333 0.13805417 0.13681667 0.71826667 0.14319167]\n",
      "[0.02437449 0.02278116 0.02301762 0.16132215 0.02400522]\n"
     ]
    }
   ],
   "source": [
    "    m35 = np.array([np.mean(r35[:,i]) for i in range(5)])\n",
    "    s35 = np.array([np.std(r35[:,i]) for i in range(5)])\n",
    "    print('mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 35% missing rate')\n",
    "    print(m35)\n",
    "    print(s35)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 427346,
     "status": "ok",
     "timestamp": 1632587999158,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "R4R6CLG-duSC",
    "outputId": "72523679-fc83-4b97-df0c-993aea5c7afb"
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "nMC = 1000 #number of runs for MC simulation\n",
    "r50 = np.array([iter_i_error(.5) for i in range(nMC)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 89722,
     "status": "ok",
     "timestamp": 1632588139241,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "T-s6pAC-dxxt",
    "outputId": "490dea9d-2bfa-4336-c632-6ae13a50b6d1"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 50% missing rate\n",
      "[0.152725   0.14746667 0.14435    0.72409583 0.1533375 ]\n",
      "[0.0253778  0.02490689 0.02421592 0.16116405 0.02732374]\n"
     ]
    }
   ],
   "source": [
    "    m50 = np.array([np.mean(r50[:,i]) for i in range(5)])\n",
    "    s50 = np.array([np.std(r50[:,i]) for i in range(5)])\n",
    "    print('mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 50% missing rate')\n",
    "    print(m50)\n",
    "    print(s50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "nMC = 1000 #number of runs for MC simulation\n",
    "r65 = np.array([iter_i_error(.65) for i in range(nMC)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "id": "HvmsGarG1Kra"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 65% missing rate\n",
      "[0.15224167 0.16002083 0.14684167 0.71433333 0.16245417]\n",
      "[0.02563044 0.03044946 0.02419662 0.16291232 0.03266594]\n"
     ]
    }
   ],
   "source": [
    "    m65 = np.array([np.mean(r65[:,i]) for i in range(5)])\n",
    "    s65 = np.array([np.std(r65[:,i]) for i in range(5)])\n",
    "    print('mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 65% missing rate')\n",
    "    print(m65)\n",
    "    print(s65)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "nMC = 1000\n",
    "#number of runs for MC simulation\n",
    "r80 = np.array([iter_i_error(.8) for i in range(nMC)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 80% missing rate\n",
      "[0.15300833 0.16718333 0.15113333 0.70077083 0.1691875 ]\n",
      "[0.02745473 0.03410681 0.02576874 0.16436033 0.03788253]\n"
     ]
    }
   ],
   "source": [
    "    m80 = np.array([np.mean(r80[:,i]) for i in range(5)])\n",
    "    s80 = np.array([np.std(r80[:,i]) for i in range(5)])\n",
    "    print('mean and sd of dper_err,mice_err, soft_err, resi_err, missf_err at 80% missing rate')\n",
    "    print(m80)\n",
    "    print(s80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.15 , 0.134, 0.134, 0.704, 0.136],\n",
       "       [0.151, 0.138, 0.137, 0.718, 0.143],\n",
       "       [0.153, 0.147, 0.144, 0.724, 0.153],\n",
       "       [0.152, 0.16 , 0.147, 0.714, 0.162],\n",
       "       [0.153, 0.167, 0.151, 0.701, 0.169]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assume_equal_cov_err = np.vstack((m20,m35,m50,m65,m80))\n",
    "assume_equal_cov_err.round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# column left to right: dper_err,mice_err, soft_err, resi_err, missf_err\n",
    "# row0 -> row5 denote the missing rate from 20% -> 80%\n",
    "assume_equal_cov_err = np.array([[0.15 , 0.134, 0.134, 0.704, 0.136],\n",
    "       [0.151, 0.138, 0.137, 0.718, 0.143],\n",
    "       [0.153, 0.147, 0.144, 0.724, 0.153],\n",
    "       [0.152, 0.16 , 0.147, 0.714, 0.162],\n",
    "       [0.153, 0.167, 0.151, 0.701, 0.169]])\n",
    "not_assume_equal_cov_err = np.array([[0.106, 0.107, 0.109, 0.697, 0.107],\n",
    "       [0.108, 0.112, 0.113, 0.716, 0.11 ],\n",
    "       [0.11 , 0.121, 0.117, 0.714, 0.113],\n",
    "       [0.113, 0.131, 0.121, 0.708, 0.118],\n",
    "       [0.128, 0.153, 0.128, 0.711, 0.126]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "    G=4\n",
    "    def iter_i_error(missing_rate):\n",
    "         # func to compute the misclassification rate for the i^{th} iteration\n",
    "         # based on the al_err_func compute error function\n",
    "         G, n_per_class = 4, 120\n",
    "         mu_1, mu_2 = np.array([2,-2]), np.array([2,2])\n",
    "         mu_3, mu_4 = np.array([-2,-2]), np.array([-2,2])\n",
    "         S= np.array([[2,.8],[.8,2]])\n",
    "         X_1 = np.random.multivariate_normal(mu_1,S,n_per_class)\n",
    "         X_2 = np.random.multivariate_normal(mu_2,S,n_per_class)\n",
    "         X_3 = np.random.multivariate_normal(mu_3,S,n_per_class)\n",
    "         X_4 = np.random.multivariate_normal(mu_4,S,n_per_class)        \n",
    "         X = np.vstack((X_1,X_2, X_3, X_4))\n",
    "         y = np.hstack((np.repeat(0,n_per_class), np.repeat(1,n_per_class), \n",
    "                        np.repeat(2,n_per_class),np.repeat(3,n_per_class)))   \n",
    "         Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.5)\n",
    "         dper_err = compute_err_MLE(Xtrain, ytrain, Xtest, ytest, G, missing_rate)\n",
    "         mice_err = compute_err_mice(Xtrain, ytrain, Xtest, ytest, G, missing_rate)\n",
    "         soft_err = compute_err_SOFT(Xtrain, ytrain, Xtest, ytest, G, missing_rate);\n",
    "         missf_err = compute_err_MissF(Xtrain, ytrain, Xtest, ytest, G, missing_rate);\n",
    "         resi_err = compute_err_resi(Xtrain, ytrain, Xtest, ytest, G, missing_rate)\n",
    "         return np.hstack((dper_err,mice_err, soft_err, resi_err, missf_err))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 25710,
     "status": "ok",
     "timestamp": 1632587563633,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "iEdO9dJPION8",
    "outputId": "c5425039-d3f5-4dcf-c147-b1d2c9fa94fc"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nIn case importing missingpy result in the error\\n\"No module named \\'sklearn.neighbors._base\\'\"\\nThen run the following\\nimport sklearn.neighbors._base\\nimport sys\\nsys.modules[\\'sklearn.neighbors.base\\'] = sklearn.neighbors._base\\n'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# !pip install impyute\n",
    "# !pip install fancyimpute\n",
    "from sklearn import datasets\n",
    "from sklearn.preprocessing import LabelEncoder, StandardScaler\n",
    "from sklearn.experimental import enable_iterative_imputer\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.impute import IterativeImputer\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as skLDA\n",
    "# import impyute as impy\n",
    "from fancyimpute import SoftImpute\n",
    "import pandas as pd\n",
    "import math\n",
    "import time\n",
    "# !pip install missingpy\n",
    "import sklearn.neighbors._base\n",
    "import sys\n",
    "sys.modules['sklearn.neighbors.base'] = sklearn.neighbors._base\n",
    "from missingpy import MissForest\n",
    "'''\n",
    "In case importing missingpy result in the error\n",
    "\"No module named 'sklearn.neighbors._base'\"\n",
    "Then run the following\n",
    "import sklearn.neighbors._base\n",
    "import sys\n",
    "sys.modules['sklearn.neighbors.base'] = sklearn.neighbors._base\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "executionInfo": {
     "elapsed": 15,
     "status": "ok",
     "timestamp": 1632587563640,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "m2pnqBY91dul"
   },
   "outputs": [],
   "source": [
    "def diag_term(i,X,y):\n",
    "  arr0 = X[:,i]\n",
    "  nar2 = 0\n",
    "  arr = arr0[~np.isnan(arr0)]\n",
    "  y_arr = y[~np.isnan(arr0)]\n",
    "\n",
    "  _, counts = np.unique(y_arr, return_counts=True)\n",
    "  ind = np.insert(np.cumsum(counts), 0, 0)\n",
    "  \n",
    "  return sum([(ind[g]-ind[g-1])*np.var(arr[ind[g-1]:ind[g]]) for \n",
    "                       g in range(1,G+1)])/len(y_arr)                       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "executionInfo": {
     "elapsed": 14,
     "status": "ok",
     "timestamp": 1632587563641,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "S7A2gflMad84"
   },
   "outputs": [],
   "source": [
    "def mle(X,y,G):\n",
    "    '''\n",
    "    X: input, should be a numpy array\n",
    "    y: label\n",
    "    G: number of classes\n",
    "    output:\n",
    "    - mus: each row is a class mean\n",
    "    - S: common covariance matrix of class 1,2,..., G \n",
    "    '''\n",
    "    epsilon = 1e-5 # define epsilon to put r down to 0 if r < epsilon\n",
    "    n,p = X.shape[0], X.shape[1]\n",
    "\n",
    "    # Estimating class means\n",
    "    mus = np.array([np.nanmean(X[y==g,:],axis=0) for g in range (G)]).T # so that each column is the mean of a class\n",
    " \n",
    "    S = np.diag([diag_term(i,X,y) for i in range(p)]) \n",
    "\n",
    "    for i in range(p):      \n",
    "      for j in range(i):\n",
    "        mat = X[:,[i,j]]\n",
    "\n",
    "        # drop rows with NA\n",
    "        idx = ~np.isnan(mat).any(axis=1)\n",
    "        mat, y_arr = mat[idx], y[idx]\n",
    "\n",
    "        _, counts = np.unique(y_arr, return_counts=True)\n",
    "        ind = np.insert(np.cumsum(counts), 0, 0)\n",
    "\n",
    "        m_g = counts\n",
    " \n",
    "        A = len(y_arr) \n",
    "        scaled_mat = [mat[ind[g-1]:ind[g],:]-mus[[i,j],g-1] for g in range(1,G+1)]   \n",
    "\n",
    "        q = lambda g: np.dot(scaled_mat[g][:,0],scaled_mat[g][:,0])\n",
    "        s11 = sum(map(q,range(G))) \n",
    "        q = lambda g: np.dot(scaled_mat[g][:,1],scaled_mat[g][:,1])\n",
    "        s22 = sum(map(q,range(G))) \n",
    "        d = lambda g: np.dot(scaled_mat[g][:,0],scaled_mat[g][:,1])\n",
    "        s12 = sum(map(d,range(G)))  \n",
    "\n",
    "        start_solve = time.time()\n",
    "        B = S[i,i]*S[j,j]*A - s22 * S[i,i] - s11 * S[j,j]\n",
    "        coefficient = [-A, s12, B, s12*S[i,i]*S[j,j]]\n",
    "        r = np.roots(coefficient)\n",
    "\n",
    "        r = r[abs(np.imag(r)) < epsilon]\n",
    "        r = np.real(r)\n",
    "        r[abs(r) < epsilon] = 0\n",
    "\n",
    "        if len(r)>1:\n",
    "          condi_var = S[j,j] - r**2/S[i,i]\n",
    "          eta = -A*np.log(condi_var)-(S[j,j]-2*r/S[i,i]*s12 +\n",
    "                                      r**2/S[i,i]**2*s11)/condi_var\n",
    "          # if condi_var <0 then eta = NA. in practice, it's impossible for cov to be negative\n",
    "          #  therefore, we drop NA elements of eta  \n",
    "          r = r[eta == max(eta[~np.isnan(eta)])]\n",
    "\n",
    "        if len(r) > 1:        \n",
    "            w = [m_g[g-1]*np.cov(mat[ind[g-1]:ind[g],], rowvar=False) for\n",
    "                 g in range(1,G+1)]\n",
    "            w = np.sum(w, axis = 0)    \n",
    "            r = r[np.abs(r-w[0,1]).argmin()] # select r that is closet to w[0,1] \n",
    "              \n",
    "        S[i,j] = S[j,i] = r\n",
    "    return [mus, S]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IwG2bsDOOxls"
   },
   "source": [
    "### compute_err function \n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "executionInfo": {
     "elapsed": 13,
     "status": "ok",
     "timestamp": 1632587563641,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "wbGqteoMenWb"
   },
   "outputs": [],
   "source": [
    "def lda_miss(mus, S, Xtest, ytrain, ytest, G):\n",
    "    # This function is from https://github.com/thunguyen177/EPEM\n",
    "    f = lambda g: np.log(np.mean(ytrain==g)) - np.matmul(\n",
    "                  np.matmul(mus[:,g].T, np.linalg.inv(S)), mus[:,g]/2)\n",
    "    last2 = [f(g) for g in np.arange(G)]\n",
    "    \n",
    "    h = lambda g,i: last2[g] + np.matmul(mus[:,g].T, np.matmul(\n",
    "                    np.linalg.inv(S), Xtest[i,:].T))\n",
    "    pred_label = [np.argmax([h(g,i) for g in np.arange(G)]) \n",
    "                  for i in np.arange(len(Xtest))]\n",
    "    pred_label = np.asarray(pred_label)\n",
    "    return np.mean(pred_label.flatten() != ytest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "executionInfo": {
     "elapsed": 14,
     "status": "ok",
     "timestamp": 1632587563642,
     "user": {
      "displayName": "Thu Nguyễn",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GjaRl8WgKOjLdrEa69IqFm__Rq7udZ763hhFDl9=s64",
      "userId": "04714273295929622355"
     },
     "user_tz": -120
    },
    "id": "fOIkN_O_eqdx"
   },
   "outputs": [],
   "source": [
    "def compute_err_MLE(Xtrain, ytrain, Xtest, ytest, G, missing_rate):    \n",
    "    # missing_rate: missing rate on training set\n",
    "    Xtr_nan,ytrain = generate_nan(Xtrain,ytrain, missing_rate)\n",
    "    scaler = StandardScaler()\n",
    "    scaler.fit(Xtr_nan)\n",
    "    Xtr_nan = scaler.transform(Xtr_nan)\n",
    "    Xtest = scaler.transform(Xtest)\n",
    "    \n",
    "    # MLEs approach\n",
    "    start = time.time()\n",
    "    mus, S = mle(Xtr_nan, ytrain, G)\n",
    "    mle_err = lda_miss(mus, S, Xtest, ytrain, ytest, G)\n",
    "    mle_time = time.time()-start\n",
    "    return mle_err#, mle_time\n",
    "    \n",
    "def compute_err_SOFT(Xtrain, ytr, Xtest, ytest, G, missing_rate):    \n",
    "    # missing_rate: missing rate on training set\n",
    "    Xtr_nan,ytr = generate_nan(Xtrain,ytr, missing_rate)\n",
    "    scaler = StandardScaler()\n",
    "    scaler.fit(Xtr_nan)\n",
    "    Xtr_nan = scaler.transform(Xtr_nan)\n",
    "    Xtest = scaler.transform(Xtest)\n",
    "\n",
    "    start = time.time()\n",
    "    Xtr_softimpute = SoftImpute(max_iters = 100).fit_transform(Xtr_nan)\n",
    "    clf_softimpute = skLDA().fit(Xtr_softimpute, ytr)\n",
    "    softimpute_err = np.mean(clf_softimpute.predict(Xtest).flatten() != ytest)\n",
    "    softimpute_time = time.time()-start\n",
    "\n",
    "    return softimpute_err#, softimpute_time\n",
    "def compute_err_MissF(Xtrain, ytr, Xtest, ytest, G, missing_rate):    \n",
    "    # missing_rate: missing rate on training set\n",
    "    Xtr_nan,ytr = generate_nan(Xtrain,ytr, missing_rate)\n",
    "    scaler = StandardScaler()\n",
    "    scaler.fit(Xtr_nan)\n",
    "    Xtr_nan = scaler.transform(Xtr_nan)\n",
    "    Xtest = scaler.transform(Xtest)\n",
    "\n",
    "    start = time.time()\n",
    "    Xtr_missf = MissForest(random_state=0).fit_transform(Xtr_nan)\n",
    "    clf_missf = skLDA().fit(Xtr_missf, ytr)\n",
    "    missf_err = np.mean(clf_missf.predict(Xtest).flatten() != ytest)\n",
    "    missf_time = time.time()-start\n",
    "\n",
    "    return missf_err#, missf_time\n",
    "\n",
    "from sklearn.experimental import enable_iterative_imputer\n",
    "from sklearn.impute import IterativeImputer\n",
    "def compute_err_mice(Xtrain, ytr, Xtest, ytest, G, missing_rate):    \n",
    "    # missing_rate: missing rate on training set\n",
    "    Xtr_nan,ytr = generate_nan(Xtrain,ytr, missing_rate)\n",
    "    scaler = StandardScaler()\n",
    "    scaler.fit(Xtr_nan)\n",
    "    Xtr_nan = scaler.transform(Xtr_nan)\n",
    "    Xtest = scaler.transform(Xtest)\n",
    "\n",
    "    start = time.time()\n",
    "    Xtr_mice= IterativeImputer(max_iter=100).fit(Xtr_nan).transform(Xtr_nan)\n",
    "    clf_mice= skLDA().fit(Xtr_mice, ytr)\n",
    "    mice_err = np.mean(clf_mice.predict(Xtest).flatten() != ytest)\n",
    "    mice_time = time.time()-start\n",
    "\n",
    "    return mice_err#, mice_time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_attribute_weights(A):\n",
    "  from sklearn.preprocessing import minmax_scale\n",
    "  y = minmax_scale(A, axis = 0)\n",
    "  p = y.copy()\n",
    "  ysum = y.sum(axis = 0)\n",
    "  for i in range(y.shape[1]):\n",
    "    if (ysum[i] == 0):\n",
    "      p[:, i] = 1\n",
    "    else:\n",
    "      p[:, i] = y[:, i] / ysum[i]\n",
    "  from scipy.special import xlogy\n",
    "  E = - xlogy(p,p).sum(axis=0) / math.log(p.shape[0])\n",
    "  w = (1 - E)/(len(E) - E.sum())\n",
    "  if np.isnan(w).any():\n",
    "    return None\n",
    "  else:\n",
    "    return w                                            \n",
    "def generate_tuple_partition(CT, ICT, m):\n",
    "  w = compute_attribute_weights(CT)\n",
    "  r = np.ones(ICT.shape[0])\n",
    "  for i in range(ICT.shape[0]):\n",
    "    for j in range(ICT.shape[1]):\n",
    "      if np.isnan(ICT[i,j]):\n",
    "        r[i] -= w[j] #If NoneType then insufficient CT set has been used\n",
    "  ICT = ICT[r.argsort()[::-1],:]\n",
    "  return np.array_split(ICT, m)\n",
    "def resi(Xtr_nan, m, n_neighbors):\n",
    "  CT = [Xtr_nan[~np.isnan(Xtr_nan).any(axis=1)]]\n",
    "  Tp = []\n",
    "  T = generate_tuple_partition(CT[0],Xtr_nan[np.isnan(Xtr_nan).any(axis=1)], m)\n",
    "  from sklearn.impute import KNNImputer\n",
    "  for i in range(m):\n",
    "    imputer = KNNImputer(n_neighbors=n_neighbors)\n",
    "    imputer.fit(CT[-1])\n",
    "    Tp.append(imputer.transform(T[i]))\n",
    "    CT.append(np.concatenate((CT[-1],Tp[-1])))\n",
    "  Tpp = []\n",
    "  for i in range(m):\n",
    "    train = CT[0]\n",
    "    for j in range(1,m):\n",
    "      if j != i:\n",
    "        train = np.concatenate((train,T[j]))\n",
    "    imputer = KNNImputer(n_neighbors=n_neighbors)\n",
    "    imputer.fit(train)\n",
    "    Tpp.append(imputer.transform(T[i]))\n",
    "  imputer = KNNImputer(n_neighbors=n_neighbors)\n",
    "  imputer.fit(CT[0])\n",
    "  Tpp.append(imputer.transform(T[m-1]))\n",
    "  CT = [CT[0]]\n",
    "  for i in range(m):\n",
    "    CT.append(np.concatenate((CT[-1],np.mean(np.array([Tp[i], Tpp[i]]), axis=0 ))))\n",
    "  return CT[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_nan(Xtrain,ytr, missing_rate):\n",
    "  minimum_complete_samples = 3\n",
    "  ct_id = np.random.choice(range(Xtrain.shape[0]), size = minimum_complete_samples, replace = False)\n",
    "  CT = Xtrain[ct_id]\n",
    "  ICT = Xtrain[[i for i in range(Xtrain.shape[0]) if i not in ct_id]]\n",
    "  ICTshape = ICT.shape\n",
    "  na_id = np.random.randint(0,ICT.size, round(missing_rate*ICT.size))\n",
    "  ICT = ICT.flatten()\n",
    "  ICT[na_id] = np.nan\n",
    "  xxx = np.concatenate((CT,ICT.reshape(ICTshape)))\n",
    "  ytrain = np.hstack((ytr[ct_id], np.array([ytr[i] for i in range(Xtrain.shape[0]) if i not in ct_id])))\n",
    "  return xxx, ytrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_err_resi(Xtrain, ytr, Xtest, ytest, G, missing_rate):    \n",
    "    # missing_rate: missing rate on training set\n",
    "    Xtr_nan,ytr = generate_nan(Xtrain,ytr, missing_rate)\n",
    "    scaler = StandardScaler()\n",
    "    scaler.fit(Xtr_nan)\n",
    "    Xtr_nan = scaler.transform(Xtr_nan)\n",
    "    Xtest = scaler.transform(Xtest)\n",
    "\n",
    "    start = time.time()\n",
    "    Xtr_resi = resi(Xtr_nan, 3, 3) #Parameters: (Dataset, m, n_neighbors) \n",
    "    clf_mice= skLDA().fit(Xtr_resi, ytr)\n",
    "    mice_err = np.mean(clf_mice.predict(Xtest).flatten() != ytest)\n",
    "    mice_time = time.time()-start\n",
    "\n",
    "    return mice_err#, mice_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "assuming equal covariance matrix Digits DPER - MICE-SOFT IMPUTE.ipynb",
   "provenance": [
    {
     "file_id": "1eazQvNk_WLTaG8heDoHQ3AgEg59y9f0v",
     "timestamp": 1608824112061
    },
    {
     "file_id": "126-leKywWg1oVxbfuH1_ORl8UMYHusn_",
     "timestamp": 1591364956450
    },
    {
     "file_id": "1Lb-WwE7STckIojvf1h3ADYX29-gfA6u-",
     "timestamp": 1590379038322
    },
    {
     "file_id": "1CnbOqwmJydQa8uCCpIR6HRVxw0YnveIP",
     "timestamp": 1590113304994
    }
   ],
   "toc_visible": true
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}