notebooks/log_likelihood.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pickle\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# %load '../predictor/util.py'\n",
    "'''\n",
    "Utility Functions.\n",
    "TODO(nautilik): These functions have currently only been tested with the Boston\n",
    "data set!\n",
    "\n",
    "Authors:\n",
    "    Alex Wang (alexwang@college.harvard.edu)\n",
    "    Luis Perez (luis.perez.live@gmail.com)\n",
    "Copyright 2015, Harvard University\n",
    "'''\n",
    "\n",
    "# The mapping of data matrix columns to indexes.\n",
    "columns = {'t': 0, 'x': 1, 'y': 2, 'count': 3}\n",
    "\n",
    "\n",
    "def split(X, tr_size):\n",
    "    '''\n",
    "    Splits input matrix X. If tr_size = 0, splits by final year.\n",
    "    Note that the ratio for the final year is currently hard-coded!\n",
    "    '''\n",
    "    n_col = np.shape(X)[1]\n",
    "    if tr_size != 0:\n",
    "        Y = np.copy(X)\n",
    "        np.random.shuffle(Y)\n",
    "        break_pt = tr_size * np.shape(X)[0]\n",
    "        train, test = Y[:break_pt, :], Y[break_pt:, :]\n",
    "    else:\n",
    "        break_pt = (3500. / 4400.) * np.shape(X)[0]\n",
    "        train, test = X[:break_pt, :], X[break_pt:, :]\n",
    "\n",
    "    tr_t, te_t = train[:, n_col - 1], test[:, n_col - 1]\n",
    "    tr, te = train[:, range(n_col - 1)], test[:, range(n_col - 1)]\n",
    "    return tr, tr_t, te, te_t\n",
    "\n",
    "\n",
    "def normalize_features(X_train):\n",
    "    '''\n",
    "    Implementation notes: set NaN to mean.\n",
    "    Generally normalizes X_train across all columns.\n",
    "    '''\n",
    "    mean_X_train = np.nanmean(X_train, 0)\n",
    "    for i in xrange(np.shape(X_train)[1]):\n",
    "        col = X_train[:, i]\n",
    "        col[np.isnan(col)] = mean_X_train[i]\n",
    "    std_X_train = np.std(X_train, 0)\n",
    "    std_X_train[std_X_train == 0] = 1\n",
    "    X_train_normalized = (X_train - mean_X_train) / std_X_train\n",
    "    return X_train_normalized\n",
    "\n",
    "\n",
    "def bucket(X, cols, num_buckets):\n",
    "    '''\n",
    "    Note: bucket edits in place\n",
    "    '''\n",
    "    Y = np.copy(X)\n",
    "    for col in cols:\n",
    "        buckets = np.linspace(np.min(X[:, col]), np.max(\n",
    "            X[:, col]), num=num_buckets + 1)\n",
    "        for i in xrange(num_buckets):\n",
    "            X_col = Y[:, col]\n",
    "            X_col[(buckets[i] <= X_col) & (X_col <= buckets[i + 1])] = i\n",
    "            Y[:, col] = X_col\n",
    "    return Y\n",
    "\n",
    "\n",
    "def rmse(predict, true):\n",
    "    '''\n",
    "    Root mean square error between the predictions and the true value.\n",
    "    '''\n",
    "    return np.sqrt(1.0 / np.shape(predict)[0] * np.sum(np.square(predict - true)))\n",
    "\n",
    "\n",
    "def createBuckets(good_data, n_buckets=15, logSpace=True):\n",
    "    '''\n",
    "    Count data for each cell. If logSpace is true, returns log values.\n",
    "    '''\n",
    "\n",
    "    data_b = bucket(good_data, [1, 2], n_buckets)\n",
    "\n",
    "    n_time = int(data_b[np.argmax(data_b[:, 0])][0])\n",
    "\n",
    "    # buckets = np.zeros((n_time, n_buckets, n_buckets))\n",
    "    buckets2 = np.zeros((n_buckets * n_buckets * n_time, 4))\n",
    "\n",
    "    # divide the data up by year and month\n",
    "    for i in xrange(n_time):\n",
    "        for j in xrange(n_buckets):\n",
    "            for k in xrange(n_buckets):\n",
    "                count = data_b[(data_b[:, 0] == i + 1) &\n",
    "                               (data_b[:, 1] == j) &\n",
    "                               (data_b[:, 2] == k)]\n",
    "                # buckets[i][j][k] = np.size(count,0)\n",
    "                buckets2[i * (n_buckets * n_buckets) +\n",
    "                         j * (n_buckets) + k][0] = i\n",
    "                buckets2[i * (n_buckets * n_buckets) +\n",
    "                         j * (n_buckets) + k][1] = j\n",
    "                buckets2[i * (n_buckets * n_buckets) +\n",
    "                         j * (n_buckets) + k][2] = k\n",
    "                buckets2[i * (n_buckets * n_buckets) + j *\n",
    "                         (n_buckets) + k][3] = np.size(count, 0)\n",
    "    print np.shape(buckets2)\n",
    "\n",
    "    if logSpace:\n",
    "        buckets2[:, 3] = np.log(np.add(sys.float_info.epsilon, buckets2[:, 3]))\n",
    "\n",
    "    return buckets2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# %load '../predictor/plot.py'\n",
    "'''\n",
    "Plotting utilities.\n",
    "\n",
    "Authors:\n",
    "    Alex Wang (alexwang@college.harvard.edu)\n",
    "    Luis Perez (luis.perez.live@gmail.com)\n",
    "Copyright 2015, Harvard University\n",
    "'''\n",
    "\n",
    "def createHeatMap(X):\n",
    "    '''\n",
    "    Given a data set, creates a heatmap of it based on x,y coordinates.\n",
    "    Ignore the temporal feature. You should subset the data before passing\n",
    "    it into this function if you'd like a heatmap for a specific time period.\n",
    "    '''\n",
    "    n = X[:, columns['x']].astype(int).max()\n",
    "    m = X[:, columns['y']].astype(int).max()\n",
    "    heatmap = np.zeros((n, m))\n",
    "    for i in xrange(n):\n",
    "        for j in xrange(m):\n",
    "            total = X[(X[:, columns['x']] == i) &\n",
    "                      (X[:, columns['y']] == j), columns['count']].sum()\n",
    "            if total > 0:\n",
    "                heatmap[i, j] = total\n",
    "    heatmap = heatmap / float(heatmap.sum())\n",
    "    return heatmap\n",
    "\n",
    "\n",
    "def plotDistribution(predict, true, city, n, process='GP'):\n",
    "    '''\n",
    "    Make some plots for n = 15 for GP process\n",
    "    '''\n",
    "    minValue = min(len(predict), 100)\n",
    "    yPred = predict[-minValue:]\n",
    "    yTrue = true[-minValue:]\n",
    "\n",
    "    # Plot Crime for Final Time Period\n",
    "    plt.clf()\n",
    "    plt.plot(yPred, label=\"Predictions\")\n",
    "    plt.plot(yTrue, label=\"Actual Data\")\n",
    "    plt.title('Crimes using {} in Last Time Period'.format(process))\n",
    "    plt.xlabel('Final {}} Regions'.format(minValue))\n",
    "    plt.ylabel('Crime Count')\n",
    "    plt.legend()\n",
    "    savefile = os.path.abspath('figures/{}_results/{}_crime_n={}_periods=last.png'.format(\n",
    "        city, process, n))\n",
    "    plt.savefig(savefile)\n",
    "    plt.close()\n",
    "\n",
    "    print \"Crimes for final period saved to {}\".format(savefile)\n",
    "\n",
    "    # Plot crime distribution for final period\n",
    "    yPred = yPred / float(np.sum(yPred))\n",
    "    yTrue = yTrue / float(np.sum(yTrue))\n",
    "    plt.clf()\n",
    "    plt.plot(yPred, label=\"Predictions\")\n",
    "    plt.plot(yTrue, label=\"Actual Data\")\n",
    "    plt.title('Predictive Distribution using {} in Last Time Period'.format(process))\n",
    "    plt.xlabel('Final {}} Regions'.format(minValue))\n",
    "    plt.ylabel('Probability')\n",
    "    plt.legend()\n",
    "    savefile = os.path.abspath('figures/{}_results/{}_dist_n={}_periods=last.png'.format(\n",
    "        city, process, n))\n",
    "    plt.savefig(savefile)\n",
    "    plt.close()\n",
    "\n",
    "    print \"Distribution saved to {}\".format(savefile)\n",
    "\n",
    "    yPred = predict[:minValue]\n",
    "    yTrue = true[:minValue]\n",
    "\n",
    "    # Plot Crime for First Time Period\n",
    "    plt.clf()\n",
    "    plt.plot(yPred, label=\"Predictions\")\n",
    "    plt.plot(yTrue, label=\"Actual Data\")\n",
    "    plt.title('Crimes using {} in First Time Period'.format(process))\n",
    "    plt.xlabel('Final {}} Regions'.format(minValue))\n",
    "    plt.ylabel('Crime Count')\n",
    "    plt.legend()\n",
    "    savefile = os.path.abspath('figures/{}_results/{}_crime_n={}_periods=first.png'.format(\n",
    "        city, process, n))\n",
    "    plt.savefig(savefile)\n",
    "    plt.close()\n",
    "\n",
    "    print \"Crimes for first period saved to {}\".format(savefile)\n",
    "\n",
    "    yPred = yPred / float(np.sum(yPred))\n",
    "    yTrue = yTrue / float(np.sum(yTrue))\n",
    "    plt.clf()\n",
    "    plt.plot(yPred, label=\"Predictions\")\n",
    "    plt.plot(yTrue, label=\"Actual Data\")\n",
    "    plt.title(\n",
    "        'Predictive Distribution using {} in First Time Period'.format(process))\n",
    "    plt.xlabel('Final 100 Regions')\n",
    "    plt.ylabel('Probability')\n",
    "    plt.legend()\n",
    "    savefile = os.path.abspath('figures/{}_results/{}_dist_n={}_period=first.png'.format(\n",
    "        city, process, n))\n",
    "    plt.savefig(savefile)\n",
    "    plt.close()\n",
    "\n",
    "    print \"Distribution saved to {}\".format(savefile)\n",
    "\n",
    "\n",
    "def plotHeatMaps(X_test, predict, city, n, process='GP'):\n",
    "    '''\n",
    "    Plots the heatmap based on the X_test data matrix and the predictions.\n",
    "    Note that X_test must have a final column with the true values of the input\n",
    "    vectors.\n",
    "    '''\n",
    "    # Attach the predictions to the data\n",
    "    trueValues = np.copy(X_test)\n",
    "    predictedValues = np.copy(X_test)\n",
    "    predictedValues[:, columns['count']] = predict.reshape((predict.shape[0]))\n",
    "\n",
    "    # Now we want to plot the heatmaps for the predictions/actual data\n",
    "    # by time period\n",
    "    months = np.unique(X_test[:, columns['t']])\n",
    "    for month in months:\n",
    "        # Create the heatmaps\n",
    "        selected = (X_test[:, columns['t']] == month)\n",
    "        if selected.sum() > 0:\n",
    "            plt.clf()\n",
    "            m1 = createHeatMap(trueValues[selected, :])\n",
    "            m2 = createHeatMap(predictedValues[selected, :])\n",
    "\n",
    "            # Make a plot only if both have data available :)\n",
    "            if m1.sum() > 0 and m2.sum > 0:\n",
    "                sns.heatmap(m1)\n",
    "                plt.title('True Density Distribution in Month {}'.format(month))\n",
    "                savefile = os.path.abspath('figures/{}_results/{}_heatmap_true_n={}_t={}.png'.format(\n",
    "                    city, process, n, month))\n",
    "                plt.savefig(savefile)\n",
    "                plt.close()\n",
    "                print \"True heatmap saved to {}\".format(savefile)\n",
    "\n",
    "                sns.heatmap(m2)\n",
    "                plt.title(\n",
    "                    'Predicted Density Distribution in Month {}'.format(month))\n",
    "                savefile = os.path.abspath('figures/{}_results/{}_heatmap_pred_n={}_t={}.png'.format(\n",
    "                    city, process, n, month))\n",
    "                plt.savefig(savefile)\n",
    "                plt.close()\n",
    "                print \"Predictions heatmap saved to {}\".format(savefile)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def KL(p,q):\n",
    "    logs = np.copy(p + sys.float_info.epsilon)\n",
    "    logs = p*np.log(logs/q)\n",
    "    return np.sum(logs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Warnings raised:', [<warnings.WarningMessage object at 0x7f791b239e10>])\n",
      "('Warning message:', 'Columns (10) have mixed types. Specify dtype option on import or set low_memory=False.')\n",
      "(\"Applying <type 'str'> dtype to columns:\", [10])\n"
     ]
    }
   ],
   "source": [
    "# Load in the base data for boston\n",
    "'''\n",
    "Read in data\n",
    "'''\n",
    "\n",
    "import pandas as pd\n",
    "import re\n",
    "import warnings\n",
    "\n",
    "bos_file = '../data/boston.csv'\n",
    "target_type = str  # The desired output type\n",
    "\n",
    "with warnings.catch_warnings(record=True) as ws:\n",
    "    warnings.simplefilter(\"always\")\n",
    "\n",
    "    bos_data = pd.read_csv(bos_file, sep=\",\", header=0)\n",
    "    print(\"Warnings raised:\", ws)\n",
    "    # We have an error on specific columns, try and load them as string\n",
    "    for w in ws:\n",
    "        s = str(w.message)\n",
    "        print(\"Warning message:\", s)\n",
    "        match = re.search(r\"Columns \\(([0-9,]+)\\) have mixed types\\.\", s)\n",
    "        if match:\n",
    "            columns = match.group(1).split(',') # Get columns as a list\n",
    "            columns = [int(c) for c in columns]\n",
    "            print(\"Applying %s dtype to columns:\" % target_type, columns)\n",
    "            bos_data.iloc[:,columns] = bos_data.iloc[:,columns].astype(target_type)\n",
    "\n",
    "'''\n",
    "Featurize data\n",
    "'''\n",
    "# temporal features\n",
    "# day of week\n",
    "day = np.array(bos_data.DAY_WEEK)\n",
    "day[ day == \"Sunday\"] = 0\n",
    "day[ day == \"Monday\"] = 1\n",
    "day[ day == \"Tuesday\"] = 2\n",
    "day[ day == \"Wednesday\"] = 3\n",
    "day[ day == \"Thursday\"] = 4\n",
    "day[ day == \"Friday\"] = 5\n",
    "day[ day == \"Saturday\"] = 6\n",
    "\n",
    "# Split mm/dd/yyyy xx:yy:zz AM/PM into components\n",
    "date_time = np.array([x.split() for x in bos_data.FROMDATE])\n",
    "date = date_time[:,0]\n",
    "time = date_time[:,1]\n",
    "tod = date_time[:,2]\n",
    "\n",
    "# month, day, year\n",
    "date = np.array([x.split('/') for x in date])\n",
    "month = [int(x) for x in date[:,0]]\n",
    "dom = [int(x) for x in date[:,1]]\n",
    "year = [int(x) for x in date[:,2]]\n",
    "# months since Jan 2012\n",
    "time_feat = np.subtract(year, 2012)*12 + month\n",
    "\n",
    "# time of day\n",
    "time_c = [x.split(':') for x in time]\n",
    "time = [int(x[1]) if (y == 'AM' and int(x[0]) == 12) else 60*int(x[0])+int(x[1]) \n",
    "        if (y =='AM' and int(x[0]) != 12) or (int(x[0]) == 12 and y == 'PM') else 12*60+60*int(x[0])+int(x[1]) \n",
    "        for x,y in zip(time_c, tod)]\n",
    "\n",
    "# grab the features we want\n",
    "data_unnorm = np.transpose(np.vstack((time_feat, bos_data.X, bos_data.Y))).astype(float)\n",
    "# remove NaNs\n",
    "good_data = data_unnorm[~(np.isnan(data_unnorm[:,1]))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "good_data = good_data[(good_data[:,0] != good_data[:,0].max())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1636"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(train_t > 0).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4300, 4)\n",
      "CPU times: user 10.9 s, sys: 12 ms, total: 10.9 s\n",
      "Wall time: 10.9 s\n",
      "CPU times: user 4 ms, sys: 0 ns, total: 4 ms\n",
      "Wall time: 1.75 ms\n"
     ]
    }
   ],
   "source": [
    "log = False\n",
    "%time data = createBuckets(good_data, n_buckets=10, logSpace=log)\n",
    "%time train, train_t, test, test_t = split(good_data, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import GPy as gp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Collect likelihoods for different n values\n",
    "testN = range(2,7) + range(8,11,2)\n",
    "likelihoods = []\n",
    "optimize = True\n",
    "likelihoods = {}\n",
    "KLs = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(172, 4)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 66 µs\n",
      "CPU times: user 52 ms, sys: 0 ns, total: 52 ms\n",
      "Wall time: 52.7 ms\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      "    00s15  007   1.723476e+04   6.333227e+11 \n",
      "    00s29  015   1.018111e+04   2.487903e+09 \n",
      "    00s60  037   3.050891e+03   4.307706e+08 \n",
      "    01s22  080   1.243379e+03   5.520190e+09 \n",
      "    01s82  126   1.079361e+03   1.672625e+07 \n",
      "    02s18  152   1.026266e+03   1.578293e+03 \n",
      "Runtime:     02s18\n",
      "Optimization status: Maximum number of f evaluations reached\n",
      "\n",
      "KL 334.144929202\n",
      "(387, 4)\n",
      "CPU times: user 4 ms, sys: 0 ns, total: 4 ms\n",
      "Wall time: 60.1 µs\n",
      "CPU times: user 308 ms, sys: 0 ns, total: 308 ms\n",
      "Wall time: 307 ms\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      "    00s61  005   1.322473e+05   2.729948e+07 \n",
      "    01s53  017   2.571787e+03   9.403852e+05 \n",
      "    02s46  026   2.540849e+03   2.420951e+05 \n",
      "    04s73  051   2.481302e+03   4.044285e+05 \n",
      "    06s02  068   2.444315e+03   7.940125e+03 \n",
      "    12s43  152   2.443585e+03   6.918906e+02 \n",
      "Runtime:     12s43\n",
      "Optimization status: Maximum number of f evaluations reached\n",
      "\n",
      "KL 459.765828487\n",
      "(688, 4)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 99.9 µs\n",
      "CPU times: user 1.07 s, sys: 0 ns, total: 1.07 s\n",
      "Wall time: 1.07 s\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      "    03s48  009   1.304406e+05   6.188847e+10 \n",
      "    04s71  013   8.529696e+04   1.062896e+11 \n",
      "    05s65  016   2.227780e+04   6.218604e+09 \n",
      "    07s57  022   1.534771e+04   5.949242e+09 \n",
      "    09s71  029   1.069103e+04   1.728876e+07 \n",
      "    12s51  039   5.734223e+03   7.094386e+10 \n",
      "    19s63  063   4.791008e+03   2.433776e+06 \n",
      "    32s28  104   4.648944e+03   3.990199e+08 \n",
      "    34s34  111   4.552994e+03   3.463919e+07 \n",
      "    39s61  128   4.181903e+03   1.642305e+07 \n",
      "    41s13  133   4.179345e+03   8.754646e+06 \n",
      "    46s70  152   4.078065e+03   7.381075e+07 \n",
      "Runtime:     46s70\n",
      "Optimization status: Maximum number of f evaluations reached\n",
      "\n",
      "KL 104.389940851\n",
      "(1075, 4)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 91.1 µs\n",
      "CPU times: user 2.86 s, sys: 0 ns, total: 2.86 s\n",
      "Wall time: 2.86 s\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      "    03s86  002   1.794714e+05   1.187873e+05 \n",
      "    13s55  012   7.023699e+04   2.743041e+12 \n",
      "    18s32  017   8.399863e+03   3.226173e+24 \n",
      "    28s39  027   8.495260e+03   7.821173e+05 \n",
      "    45s98  045   8.366149e+03   5.266016e+24 \n",
      "Runtime:     45s98\n",
      "Optimization status: Converged\n",
      "\n",
      "KL 450.514446422\n",
      "(1548, 4)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 78.2 µs\n",
      "CPU times: user 6.65 s, sys: 0 ns, total: 6.65 s\n",
      "Wall time: 6.65 s\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      "    32s04  011   2.671914e+04   8.519395e+10 \n",
      "    56s56  021   1.437191e+04   5.630286e+06 \n",
      " 01m01s73  023   1.437157e+04   5.627959e+06 \n",
      " 01m18s40  031   1.437053e+04   5.621521e+06 \n",
      " 01m27s05  035   1.436983e+04   5.686845e+06 \n",
      " 01m29s43  036   1.436983e+04   5.617659e+06 \n",
      " 01m36s81  039   1.436913e+04   5.614126e+06 \n",
      " 02m01s83  051   1.436771e+04   5.608005e+06 \n",
      " 02m07s33  053   1.436735e+04   5.606663e+06 \n",
      " 02m24s05  061   1.436628e+04   5.603069e+06 \n",
      " 02m31s33  064   1.436556e+04   5.601020e+06 \n",
      " 02m44s35  071   1.436484e+04   5.599237e+06 \n",
      " 03m02s12  078   1.436374e+04   5.597042e+06 \n",
      " 03m27s66  091   1.436190e+04   5.594589e+06 \n",
      " 03m35s56  095   1.436116e+04   5.665674e+06 \n",
      " 03m47s77  101   1.436041e+04   5.593639e+06 \n",
      " 03m55s67  105   1.435966e+04   5.665739e+06 \n",
      " 04m08s17  111   1.435890e+04   5.593526e+06 \n",
      " 04m26s38  119   1.435738e+04   5.594194e+06 \n",
      " 04m29s24  121   1.435738e+04   5.594194e+06 \n",
      " 04m47s72  129   1.435584e+04   5.595592e+06 \n",
      " 05m10s99  141   1.435428e+04   5.597673e+06 \n",
      " 05m26s10  148   1.435310e+04   5.599658e+06 \n",
      " 05m33s91  152   1.435270e+04   5.600393e+06 \n",
      "Runtime:  05m33s91\n",
      "Optimization status: Maximum number of f evaluations reached\n",
      "\n",
      "KL 292.48975078\n",
      "(2752, 4)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 93 µs\n",
      "CPU times: user 27.1 s, sys: 268 ms, total: 27.4 s\n",
      "Wall time: 27.5 s\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      "    24s60  001   1.887743e+05   4.188282e+03 \n",
      " 01m31s26  007   1.816467e+05   7.608504e+08 \n",
      " 02m41s02  013   4.938900e+04   2.289498e+11 \n",
      " 03m51s74  019   3.441288e+04   1.100598e+12 \n",
      " 04m47s16  024   1.965860e+04   2.388761e+08 \n",
      " 06m30s32  033   1.904593e+04   7.729795e+07 \n",
      " 15m53s15  047   1.387074e+04   4.277675e+07 \n",
      " 18m08s52  059   1.367532e+04   7.741339e+04 \n",
      " 23m54s35  090   1.336880e+04   3.903092e+07 \n",
      " 25m20s92  098   1.315253e+04   3.337416e+07 \n",
      " 26m04s48  102   1.306470e+04   1.379396e+07 \n",
      " 26m37s02  105   1.321097e+04   1.475578e+09 \n",
      " 28m30s89  115   1.298314e+04   3.951632e+06 \n",
      " 30m05s30  123   1.297400e+04   4.191305e-01 \n",
      "Runtime:  30m05s30\n",
      "Optimization status: Converged\n",
      "\n",
      "KL 5860.27785014\n",
      "(4300, 4)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 96.1 µs\n",
      "CPU times: user 1min 12s, sys: 14.7 s, total: 1min 26s\n",
      "Wall time: 1min 35s\n",
      "WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n",
      "Running L-BFGS-B (Scipy implementation) Code:\n",
      "  runtime   i     f              |g|        \n",
      " 05m17s65  007   1.622823e+05   4.723162e+13 \n",
      " 11m44s76  017   6.707304e+04   7.615566e+11 \n",
      " 13m43s07  020   2.673694e+05   9.452543e-43 \n",
      " 18m24s06  027   2.085235e+04   9.295752e+09 \n",
      " 22m35s54  033   1.831159e+04   9.251310e+07 \n",
      " 33m39s29  050   1.795366e+04   5.426638e+08 \n",
      " 36m56s76  055   1.785071e+04   7.158137e+06 \n",
      " 42m54s10  064   1.962258e+04   1.018331e+11 \n",
      " 48m23s38  072   1.779754e+04   8.586540e+05 \n",
      " 01h02m41  094   1.727396e+04   3.246583e+05 \n",
      " 01h12m17  107   1.724344e+04   3.775361e+03 \n",
      " 01h19m09  118   1.724324e+04   1.697708e+04 \n",
      " 03h06m33  142   1.724321e+04   3.043073e+02 \n",
      " 03h13m11  152   1.724320e+04   1.470197e+02 \n",
      "Runtime:  03h13m11\n",
      "Optimization status: Maximum number of f evaluations reached\n",
      "\n",
      "KL 8802.04737447\n"
     ]
    }
   ],
   "source": [
    "for n in testN:\n",
    "    # Bucketize the data as specified! By default, does Boston data.\n",
    "    data = createBuckets(good_data, n, logSpace=False)\n",
    "\n",
    "    # Split for latest year.\n",
    "    %time train, train_t, test, test_t = split(data, 0)\n",
    "\n",
    "    # Train the model and optimize it too\n",
    "    k_lin = gp.kern.Linear(input_dim = 3)\n",
    "    k_per = gp.kern.StdPeriodic(input_dim = 3, ARD1= True, ARD2=True)\n",
    "    \n",
    "    # Linear times periodic\n",
    "    k = k_lin * k_per\n",
    "    \n",
    "    # Optmize the values\n",
    "    train_t = train_t.reshape((train_t.shape[0], 1))\n",
    "    %time m = gp.models.GPRegression(train, train_t, k)\n",
    "    \n",
    "    # Fix the variance\n",
    "    if optimize:\n",
    "        m.Gaussian_noise.variance.constrain_fixed(train_t.std())\n",
    "\n",
    "        # Optimize\n",
    "        m.optimize(messages=True, max_iters=150)\n",
    "\n",
    "    # Get likelihood\n",
    "    likelihoods[n] = m.log_likelihood()\n",
    "    \n",
    "    \n",
    "    predictions_optimal = m.predict(test)\n",
    "    preds2 = predictions_optimal[0]\n",
    "    preds_p = preds2.clip(min=sys.float_info.epsilon)\n",
    "    n_test_t = test_t / np.sum(test_t)\n",
    "    n_preds = preds_p / np.sum(preds_p)\n",
    "    \n",
    "    KLs[n] = KL(n_test_t, n_preds)\n",
    "    print \"KL {}\".format(KLs[n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{2: 334.14492920210171,\n",
       " 3: 459.76582848657102,\n",
       " 4: 104.38994085124376,\n",
       " 5: 450.51444642156696,\n",
       " 6: 292.48975078031987,\n",
       " 8: 5860.2778501384018,\n",
       " 10: 8802.0473744740884}"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "KLs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.clf()\n",
    "x = KLs.keys()\n",
    "kls = KLs.values()\n",
    "likelihood = likelihoods.values()\n",
    "plt.plot(x, likelihood)\n",
    "plt.title(\"Linear*Periodic Kernel on Boston City Data\")\n",
    "plt.xlabel(\"Dimension of Bucketization Grid\")\n",
    "plt.ylabel(\"Log Likelihood\")\n",
    "plt.savefig(\"kl_likelihood_periodic_linear_periodic.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.clf()\n",
    "x = KLs.keys()\n",
    "kls = KLs.values()\n",
    "likelihood = likelihoods.values()\n",
    "plt.plot(x, kls)\n",
    "plt.title(\"Linear*Periodic Kernel on Boston City Data\")\n",
    "plt.xlabel(\"Dimension of Bucketization Grid\")\n",
    "plt.ylabel(\"Measured KL Divergence\")\n",
    "plt.savefig(\"kl_periodic_linear_periodic.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}