incubator/multiclass-classification-neural-network.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import theano.tensor as tt\n",
    "import theano\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from pymc3.backends import HDF5, text\n",
    "\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theano.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "I'd like to try doing multiclass classification with a deep neural network. The network architecture is a feed-forward network with one hidden layer in between the input and output; `n_hidden` units is a hyperparameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Specification\n",
    "\n",
    "Because neural networks are the highlight here, I will first do the model specification up-front."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_nn(ann_input, ann_output, n_hidden):\n",
    "    \"\"\"\n",
    "    Makes a feed forward neural network with n_hidden layers for doing multi-\n",
    "    class classification.\n",
    "    \n",
    "    Feed-forward networks are easy to define, so I have not relied on any \n",
    "    other Deep Learning frameworks to define the neural network here.\n",
    "    \"\"\"\n",
    "    init_1 = np.random.randn(ann_input.shape[1], n_hidden)\n",
    "    init_2 = np.random.randn(n_hidden, n_hidden)\n",
    "    init_out = np.random.randn(n_hidden, ann_output.shape[1])\n",
    "\n",
    "    with pm.Model() as nn_model:\n",
    "        # Define weights\n",
    "        weights_1 = pm.Normal(\n",
    "            \"w_1\", mu=0, sd=1, shape=(ann_input.shape[1], n_hidden), testval=init_1\n",
    "        )\n",
    "        weights_2 = pm.Normal(\n",
    "            \"w_2\", mu=0, sd=1, shape=(n_hidden, n_hidden), testval=init_2\n",
    "        )\n",
    "        weights_out = pm.Normal(\n",
    "            \"w_out\", mu=0, sd=1, shape=(n_hidden, ann_output.shape[1]), testval=init_out\n",
    "        )\n",
    "\n",
    "        # Define activations\n",
    "        acts_1 = pm.Deterministic(\n",
    "            \"activations_1\", tt.tanh(tt.dot(ann_input, weights_1))\n",
    "        )\n",
    "        acts_2 = pm.Deterministic(\"activations_2\", tt.tanh(tt.dot(acts_1, weights_2)))\n",
    "        acts_out = pm.Deterministic(\n",
    "            \"activations_out\", tt.nnet.softmax(tt.dot(acts_2, weights_out))\n",
    "        )  # noqa\n",
    "\n",
    "        # Define likelihood\n",
    "        out = pm.Multinomial(\"likelihood\", n=1, p=acts_out, observed=ann_output)\n",
    "\n",
    "    return nn_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocess Data\n",
    "\n",
    "## Basic Cleaning\n",
    "\n",
    "Now, let's read in the dataset. There's a bunch of preprocessing that has to happen. I happened to have this code written for the IACS 2017 data science bootcamp, and copied it over from there. It's commented out because it takes some time to execute, and if executed once, it needn't be executed again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df = pd.read_csv('datasets/covtype.data', header=None)\n",
    "\n",
    "# columns = [\n",
    "#     'Elevation',\n",
    "#     'Aspect',\n",
    "#     'Slope',\n",
    "#     'Horizontal_Distance_To_Hydrology',\n",
    "#     'Vertical_Distance_To_Hydrology',\n",
    "#     'Horizontal_Distance_To_Roadways',\n",
    "#     'Hillshade_9am',\n",
    "#     'Hillshade_Noon',\n",
    "#     'Hillshade_3pm',\n",
    "#     'Horizontal_Distance_To_Fire_Points'\n",
    "# ]\n",
    "\n",
    "# # Add in wilderness area data (binary)\n",
    "# for i in range(1, 5):\n",
    "#     columns.append('Wilderness_Area_{0}'.format(i))\n",
    "\n",
    "# # Add in soil type data (binary)\n",
    "# for i in range(1, 41):\n",
    "#     columns.append('Soil_Type_{0}'.format(i))\n",
    "\n",
    "# # Add in soil cover type\n",
    "# columns.append('Cover_Type')\n",
    "\n",
    "# df.columns = columns\n",
    "\n",
    "# # Add in soil codes. These were downloaded from the UCI repository.\n",
    "# soil_codes = pd.read_csv('datasets/climatic_geologic_zone.csv')\n",
    "# soil_dict = soil_codes.set_index('soil_type').to_dict()\n",
    "\n",
    "# # Add geologic and climatic zone code to soil type\n",
    "# for i in range(1, 41):\n",
    "#     df.loc[df['Soil_Type_{i}'.format(i=i)] == 1, 'Climatic_Zone'] = soil_dict['climatic_zone'][i]\n",
    "#     df.loc[df['Soil_Type_{i}'.format(i=i)] == 1, 'Geologic_Zone'] = soil_dict['geologic_zone'][i]\n",
    "\n",
    "# # Encode one-of-K for the geologic zone, climatic zone, and cover_type encodings.\n",
    "# # This is important because the geologic and climatic zones aren't ordinal - they are strictly categorical.\n",
    "# enc = OneHotEncoder()\n",
    "# clm_zone_enc = enc.fit_transform(df['Climatic_Zone'].values.reshape(-1, 1)).toarray()\n",
    "# geo_zone_enc = enc.fit_transform(df['Geologic_Zone'].values.reshape(-1, 1)).toarray()\n",
    "# cov_type_enc = enc.fit_transform(df['Cover_Type'].values.reshape(-1, 1)).toarray()\n",
    "\n",
    "# for i in range(clm_zone_enc.shape[1]):\n",
    "#     df['Climatic_Zone_{i}'.format(i=i)] = clm_zone_enc[:, i]\n",
    "\n",
    "# for i in range(geo_zone_enc.shape[1]):\n",
    "#     df['Geologic_Zone_{i}'.format(i=i)] = geo_zone_enc[:, i]\n",
    "\n",
    "# del df['Climatic_Zone']\n",
    "# del df['Geologic_Zone']\n",
    "\n",
    "# df.to_csv('datasets/covtype_preprocess.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"../datasets/covtype_preprocess.csv\", index_col=0)\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now make the X and Y matrices. We use patsy to give us a quick and fast way to turn categorical variables into one-of-K columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"Cover_Type\"] = df[\"Cover_Type\"].apply(lambda x: str(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_col = \"Cover_Type\"\n",
    "input_cols = [c for c in df.columns if c != output_col]\n",
    "input_formula = \"\".join(c + \" + \" for c in input_cols)\n",
    "input_formula = input_formula + \"-1\"\n",
    "\n",
    "import patsy\n",
    "from sklearn.preprocessing import scale, normalize\n",
    "\n",
    "X = patsy.dmatrix(formula_like=input_formula, data=df, return_type=\"dataframe\")\n",
    "# X = normalize(X)\n",
    "\n",
    "Y = patsy.dmatrix(formula_like=\"Cover_Type -1\", data=df, return_type=\"dataframe\")\n",
    "print(X.shape, Y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Balance Classes\n",
    "\n",
    "We will balance out the classes to make them evenly distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "downsampled_targets = []\n",
    "\n",
    "for i in range(1, 7 + 1):\n",
    "    # print(f'target[{i}]')\n",
    "    target = Y[Y[\"Cover_Type[{i}]\".format(i=i)] == 1]\n",
    "    # print(len(target))\n",
    "    downsampled_targets.append(target.sample(2747))\n",
    "\n",
    "mms = MinMaxScaler()\n",
    "X_tfm = pm.floatX(mms.fit_transform(X[input_cols]))\n",
    "\n",
    "Y_downsampled = pd.concat(downsampled_targets)\n",
    "Y_downsamp = pm.floatX(Y_downsampled)\n",
    "\n",
    "X_downsampled = X_tfm[Y_downsampled.index]\n",
    "X_downsamp = pm.floatX(X_downsampled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_downsamp.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from models.feedforward import ForestCoverModel\n",
    "\n",
    "fcm = ForestCoverModel(n_hidden=20,)\n",
    "fcm.fit(X_downsamp, Y_downsamp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_downsampled.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_downsamp.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_hidden = 20  # define the number of hidden units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "variables": {
     "n_hidden": "20"
    }
   },
   "source": [
    "# Model Execution\n",
    "\n",
    "We now make the model with {{n_hidden}} hidden units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = make_nn(X_downsamp, Y_downsamp, n_hidden=n_hidden)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    # s = theano.shared(pm.floatX(1.1))\n",
    "    # inference = pm.ADVI(cost_part_grad_scale=s)\n",
    "    approx = pm.fit(\n",
    "        500000, callbacks=[pm.callbacks.CheckParametersConvergence(tolerance=1e-1)]\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(approx.hist)\n",
    "plt.yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    trace = approx.sample(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pm.traceplot(trace, varnames=['w_1'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pm.traceplot(trace, varnames=['w_2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pm.traceplot(trace, varnames=['w_out'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    samp_ppc = pm.sample_ppc(trace, samples=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds_proba = samp_ppc[\"likelihood\"].mean(axis=0)\n",
    "preds = (preds_proba == np.max(preds_proba, axis=1, keepdims=True)) * 1\n",
    "plt.pcolor(preds)\n",
    "plt.savefig(\"figures/class_predictions.png\", dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolor(preds_proba)\n",
    "plt.savefig(\"figures/class_probabilities.png\", dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolor(samp_ppc[\"likelihood\"].std(axis=0))\n",
    "plt.savefig(\"figures/class_uncertainties.png\", dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import classification_report\n",
    "\n",
    "print(classification_report(Y_downsamp, preds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to the logistic regression notebook, we have higher performance!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit `sklearn`-like estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bayesian",
   "language": "python",
   "name": "bayesian"
  },
  "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.6.6"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "138px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_position": {
    "height": "593px",
    "left": "0px",
    "right": "1122px",
    "top": "106px",
    "width": "158px"
   },
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}