{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "dfadef99-d54a-f7fd-9992-c14bcd04b9b6"
      },
      "source": [
        "I suspect that the pups will be difficult to count on the images because they are dark grey just like the rocks, they are usually very close to their mother, and they are pretty similar in shape, size, and color to the seals (see 5.jpg for example). However, all other sea lion types are light brown so they should be easier to locate and count. I explore in this notebook whether it is possible to predict the number of pups based on the number of other sea lions. That is, the pups are not directly counted but their number is inferred/estimated based on the number of males, females, juveniles, and subadult males using a regression model. This alternative approach could turn out to be more accurate than directly counting the pups and I hope some of you will find it useful to improve the overall accuracy of the counts."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a7703ccf-1493-a7d2-09a9-c6813981837f"
      },
      "outputs": [],
      "source": [
        "# load packages and read in the data\n",
        "\n",
        "from subprocess import check_output\n",
        "import numpy as np\n",
        "import xgboost as xgb\n",
        "from sklearn.model_selection import GridSearchCV\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.metrics import mean_squared_error\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib\n",
        "\n",
        "data = np.genfromtxt('../input/Train/train.csv', delimiter=',',skip_header=1,usecols=(1,2,3,4,5))\n",
        "X = data[:,:4]\n",
        "Y = data[:,4]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "2b744e4f-b395-e593-9e0f-d503c121bf79"
      },
      "source": [
        "**Step 1)** It is probably relatively easy to count the light brown blobs on the image (the sum of non-pup sea lions). So the first step is to check if the number of pups correlate with the number of non-pup sea lions. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "80768df3-cb7a-4356-1484-5c0033ea2a34"
      },
      "outputs": [],
      "source": [
        "blobs = np.sum(X,axis=1)\n",
        "\n",
        "plt.scatter(blobs,Y)\n",
        "plt.xlabel('#non-pup sea lions')\n",
        "plt.ylabel('#pups')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5cdcb197-dfd9-4dff-49bc-41f101b3e44e"
      },
      "source": [
        "The correlation looks pretty weak because often the sea lion colony has zero or very few pups. \n",
        "\n",
        "**Step 2)** Develop a regression model and use the number of adult_males, subadult_males, adult_females, and juveniles as features. The idea is that the family demography of sea lion colonies might help to improve the prediction.  \n",
        "\n",
        "The data is split into train (80%) and test (20%),  GridSearchCV is used on the training set to tune some xgboost parameters (learning rate and the number of trees only), and the test set is used to report RMSE and to save the feature importances. These steps (random splitting, grid search, RMSE, feature importances) are performed 10 times to assess the impact of random splitting on the RMSE. \n",
        "\n",
        "The prediction of the last model is visualized as well as the feature importances with their standard deviations.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a7547ad1-6b84-aac5-8477-418c4ecff88f"
      },
      "outputs": [],
      "source": [
        "# a function to do the training and prediction\n",
        "def train_pred(n_sims,X,Y,f_names,test_size):\n",
        "    RMSE = np.zeros(n_sims)\n",
        "    f_imp = np.zeros([n_sims,np.shape(X)[1]])\n",
        "\n",
        "    for i in range(n_sims):\n",
        "\n",
        "        # split the data\n",
        "        X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=test_size)\n",
        "        # initialize XGBRegressor\n",
        "        GB = xgb.XGBRegressor()\n",
        "\n",
        "        # the parameter grid below was too much on the kaggle kernel\n",
        "        #param_grid = {\"learning_rate\": [0.01,0.03,0.1],\n",
        "        #              \"objective\": ['reg:linear'],\n",
        "        #              \"n_estimators\": [300,1000,3000]}\n",
        "        # do GridSearch\n",
        "        #search_GB = GridSearchCV(GB,param_grid,cv=4,n_jobs=-1).fit(X_train,Y_train)\n",
        "        # the best parameters should not be on the edges of the parameter grid\n",
        "        #print('   ',search_GB.best_params_)\n",
        "        # train the best model\n",
        "        #xgb_pups = xgb.XGBRegressor(**search_GB.best_params_).fit(X_train, Y_train)\n",
        "\n",
        "        # preselected parameters\n",
        "        param_grid = {\"learning_rate\": 0.03,\n",
        "                      \"objective\": 'reg:linear',\n",
        "                      \"n_estimators\": 300}\n",
        "        xgb_pups = xgb.XGBRegressor(**param_grid).fit(X_train, Y_train)\n",
        "\n",
        "        # predict on the test set\n",
        "        preds = xgb_pups.predict(X_test)\n",
        "\n",
        "        # feature importance\n",
        "        b = xgb_pups.booster()\n",
        "        f_imp[i,:] = list(b.get_fscore().values())\n",
        "\n",
        "        # rmse of prediction\n",
        "        RMSE[i] = np.sqrt(mean_squared_error(Y_test, preds))\n",
        "    \n",
        "    # visualize the prediction of the last model\n",
        "    plt.scatter(Y_test,preds,label = 'regression model')\n",
        "    plt.plot(np.arange(np.max(Y_test)),np.arange(np.max(Y_test)),color='k',label='perfect prediction')\n",
        "    plt.title('predictions of the last model')\n",
        "    plt.legend(loc='best')\n",
        "    plt.xlabel('true #pups')\n",
        "    plt.ylabel('predicted #pups')\n",
        "    plt.show()\n",
        "    \n",
        "    return RMSE, f_imp\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "10119cdd-f4e4-64af-88cd-a90862476650"
      },
      "outputs": [],
      "source": [
        "f_names = ['adult males','subadult males','adult females','juveniles']\n",
        "RMSE, f_imp = train_pred(10,X,Y,f_names,test_size=0.2)\n",
        "\n",
        "print('RMSE = ',np.around(np.mean(RMSE),1),'+/-',np.around(np.std(RMSE),1))\n",
        "\n",
        "plt.bar(range(len(f_names)),np.mean(f_imp,axis=0),width=0.8,yerr = np.std(f_imp,axis=0))\n",
        "plt.ylabel('f score')\n",
        "plt.xticks(range(len(f_names)), f_names)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "db9008be-a9cb-35b1-860c-fa0f6f690d4a"
      },
      "source": [
        "**Step 3)** Let's do some feature engineering! The sea lion counts are normalized by the total count, and the total count is added as an additional feature. \n",
        "\n",
        "As before, the prediction of the last model is visualized as well as the feature importances with their standard deviations."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ba54545a-c465-33d2-5692-da8b07a279f6"
      },
      "outputs": [],
      "source": [
        "X_sum = np.sum(X,axis=1)\n",
        "X_new = np.hstack((X/X_sum[:,None],X_sum[:,None]))\n",
        "f_names_new = np.append(f_names,'sum')\n",
        "\n",
        "RMSE_new, f_imp_new = train_pred(10,X_new,Y,f_names_new,test_size=0.2)\n",
        "\n",
        "print('RMSE = ',np.around(np.mean(RMSE_new),1),'+/-',np.around(np.std(RMSE_new),1))\n",
        "\n",
        "plt.bar(range(len(f_names_new)),np.mean(f_imp_new,axis=0),width=0.8,yerr = np.std(f_imp_new,axis=0))\n",
        "plt.ylabel('f score')\n",
        "plt.xticks(range(len(f_names_new)), f_names_new,rotation='vertical')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "a867cb6b-2cb1-7848-2e77-8dbe767c3386"
      },
      "source": [
        "The RMSE did not improve significantly."
      ]
    }
  ],
  "metadata": {
    "_change_revision": 0,
    "_is_fork": false,
    "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.6.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}