{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analysis of the clusters\n", "\n", "Here is a streamlined code that shows just analysis without blind alleys I encountered before. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "#general\n", "import numpy as np\n", "import scipy\n", "from matplotlib import pyplot as plt\n", "%pylab inline\n", "import pandas as pd\n", "import MySQLdb\n", "import os\n", "import sys\n", "sys.setrecursionlimit(3000)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "data loaded\n" ] } ], "source": [ "con=MySQLdb.connect(user=user, passwd=passwd, db=dbname, host=host)\n", "df = pd.read_sql(\"SELECT * FROM business WHERE state ='AZ' \", con)\n", "print \"data loaded\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The low activity of some businesses in their particular cluster caused their Location Quotient to be marked as inf because of the computing roundup.
\n", "The trick is with the LQ equation. Here you can see the equation again.
\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$LQ_{ij}=\\frac{\\frac{E_{ij}}{E_i}}{\\frac{\\sum_i E_{ij}}{\\sum_i E_i}}$$\n", "Where \n", "$E_{ij}$ is economic activity in subarea i, department j \n", "$E_i$ is total economic activity in subarea i \n", "$\\sum_i E_{ij}$ is economic activity of department j in the whole area \n", "$\\sum_i E_i$ is total economic activity in the whole area" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So when a business has low activity, the real value of nominator is of an order of $10^{-6}$. That is small enough number that my Python code on computer decides that the value is so close to the zero and thus can be zero; causing division to produce infinity. This creates a problem to me since, in reality, those businesses are not the most popular, they are the least popular in the cluster.
\n", "Therefore, I decided to replace all infinity values with the minimum value of Location Quotient. Yes, this will assign a higher value of LQ to those unpopular businesses, but in clusters that have more than four business each, such low visited businesses do not even come into future calculations. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/Lexa/anaconda/lib/python2.7/site-packages/pandas/core/indexing.py:118: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", " self._setitem_with_indexer(indexer, value)\n" ] } ], "source": [ "i=0\n", "for t in df.LQ:\n", " if t == np.inf:\n", " df.LQ.loc[i]=df.LQ.min()\n", " i=i+1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us see individual cluster." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAJ1CAYAAABdKx0tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcXmV99/HPb5bse0gISQghkrDKvgSREDZZVBDRorXi\nUvdH9LFqrX2eCt2sWrRKXUotFGy11CogKogiRBHZ9yWsIWQlQPY9mZnf88e58zBOZiYzyZ2558x8\n3q/XvGbuc67rnN89SSbfuc65rhOZiSRJksqnrtYFSJIkaecY5CRJkkrKICdJklRSBjlJkqSSMshJ\nkiSVlEFOkiSppAxykiRJJWWQk9SnRcT8iHi+1nUARMTsiGiJiItrXYukvsEgJ2m3qYSWlm72mRER\n34qIJyNibUSsq3z9rYiYsRNlZOWjN+mxeiLiqsqfw5SeOqeknmOQk7S7dTm0RMQngCeAjwCLge8A\n36p8/RHg8Yi4qJvnPwU4tZt9+preFmQlVUlDrQuQJICIuBD4OrAcOC8zf9dm/+uB64FvRMTKzPzP\nrhw3M3vFZdUai8qHpD7GETlJNRcRwylCXAJ/3DbEAVS2vavy8usRMayLx97uHrmIeG/lcuN7IuLk\niJgTEWsiYnVE/CwiDtiJ9/CGiPhpRLwUEZsiYkFEXB8ROxwN7Ow+voi4pFLrrDbbT6ycb1HlfEsj\n4s6I+EKrNi3AhZWXz2+71N3O92NMRPxDRMyNiA0RsSoibomI09upp/X37szK9251dy+hS6oOR+Qk\n9QZvA0YBd2fmrzpqlJk3R8S9wDHA+cDVXTx+R5cW3wScC9xIcRn3YOBs4JiIOCgzl3fl4BHx18Bf\nAWspRg0XApOA11GEz1/vQo3tne9M4OfAKuAGikvPY4CDgI8Cf1Np+tfAW4DDKILyqsr2Va2OtQ8w\nB9gH+C3F92IYxffmFxHx4cz8t3bKeBtwZqX9tyv9JfUwg5yk3uD1lc+3dKHtryiC3Al0Pch15Fzg\njMy8bduGiPgi8BfA+4F/3NEBIuINFCFuHnBiZi5ts3/SLtbYng9SXCqdnZmPtjnfmG1fZ+ZfR8S+\nVIJcZi5o51hXA3sD78jMH7Y6zkiKgHdZRNyQmS+16XcWcHZm/rIab0jSzvHSqqTeYK/K54VdaLuo\nTZ9dcU3rEFfxr5XPx3TxGNsmX3y6bYgDyMzFO1tcF2xq53wruto5Ig4DZgE/bh3iKsdZDVwCDKIY\n/WzrJ4Y4qfYckWtHRPwtcA7FpY7lwHszc2GbNoOA3wADgQEUP9Q+31n/iHgX8JlWhzkUOCIzH4mI\nXwATgEbgLuAjmbm1kxoPAP4dOAL4P5n51V1/51Kp1FfhGPe1s21bUBzdxWPMBFqAX1Shnq76T+A8\n4O6I+G+KkbM7MnNRp722d3zl86iIuKSd/eMqnw9sZ9893TyXpN2g3we5iJgNvCcz39dq81cy868q\n+y8CLgY+0LpfZm6KiJMzc0NENAC/i4jXV27Ibrd/Zn4f+H5l+yHAdZn5SOWQb8vMdZV9PwIuoPhh\n3ZHlFCMBb9mFty/1Fi9WPndlrbO9K5+rMdK1qu2GzGyKCOh6UBwFrMzMzVWop0sy87qIeBPwaYpL\nwB8GiIj7gc9nZlcuUQOMrXw+vfLR7umAoe1sf7GdbZJ6mJdW27nBODPXtno5DHil3Y6ZGypfDqD4\nob+iG/3/GLim1bG2hbjGyvFeqbweFxE/ioh7Kh+vq7R/OTPvAzoctZNK5PbK59O60HZbmzt2Uy3d\ntQoYXRml31ktdPyL9aj2NmbmjZl5amX/qcA/UUzW+FlEtDeC1p7Vlc+fyMy6Dj7qM/NP2yuhi+eQ\ntBsZ5DpYWyki/j4iFgDvAb7UQZu6iHgIWAbclplPdKP/HwH/1eZ4N1eOtTEzt12m+QbwT5l5LMUs\nsfZmj0ll9z/ASuDYiOgwzFWWwziG4pemH/VQbTtyJ8XP0jN34RgrgT0ro/ttHd1Zx8zcmJm3Zean\ngS9S/CJ4VqsmzZXP7Y0w3ln5PKudfZJKoN8GuYi4KyIeBL4LnBMRD1Y+TgfIzP+TmVOAqyh+091O\nZrZk5uHAZGBW5TItO+ofEccBG1oHv0qfMyhu4B4YEe+pbD4N+Gal1p8AwyNiyK69e6l3qYxIf6ry\n8gfbRp5bq2z7AcXo1Qe2jWL3Av9c+fzViJjYdmd729pxN8X9sa1v8SAi3kuxhEm22T4rItoLZhMq\nn9e32rZtCZXtlgfJzPspRkPfGhHva7u/cq7XRsS49vZJqr1+e49cZs4EiIiTKCYjtPtDjOI/jht3\ncKzVEfFzit+c53Sh/zsq29s71uaI+DFwHMWyAAEcl5lbOqtB6sUiIq7qYF8CH6uMKn0vIkYBXwVu\nj4g5wAOVNkcBJ1Os0/ahzLx+95fdNZn5q4j4O+D/AnMj4nqKCRN7UiyrcidtAlo7/rnS5juVBYQX\nAYdTTKT4GcWabq1dBkyMiDuAF4AtvPo9mk+r2zYolnT5DPDdiLiW4nu4MjO/Vdn/x8CtwBVRPCLt\nHorLxZMpJmQdXKnj5S5+SyT1oH4b5FrZ7tJqREzPzGcqL88FHmynzR5AU2auiojBFDcK//WO+kdE\nHfB2Xl03i4gYCozIzKWVSytvArZN6/8l8Ang0krbwzPzoc7ql3qZ5NWnC7TeFpXPnwQ2AmTmZZVb\nDD5Bcd/XTGBwpc8a4LVtZ5B38fztbavaPV6Z+YWIuJOi7jdRTA5YRjErdodr3WXm3Mol5S8Cb6a4\n9/V2ivd/PvDGNl3+nmLW6tEUo/YtFIHu7ynWi9t27xuZ+cuI+DTF2nOfpLj0Op/iGbZk5uKIOIpi\n8tT5FMGuHlhK8dzbbwCPtS4X74+Teo3I7N//Hisjcu/JzPe32vYjYH+Ke0ueAz6amS9VLpF8NzPf\nGBGHUlw2rat8/Edm/mNn/Sv7ZgNfzMzXtTrfeIrfugdS/Od2M/DnmZkRMZbiB+6BFMH7N5n5sYiY\nANwLjKD4Ib4WOKgXXW6SqqJyCfFGil+WrsjMD9a4JEnqNfp9kJPU+0XECIpZqgdTrJv4DzUuSZJ6\nBS+tSur1MnNNRLyR4j6y+ogYlZnbrQEnSf2NI3KSJEkl1W9H5CLCBCtJkkojM7eb4Nhvgxy0/w1p\nKyIuycxLeqAc9QL+efc//pn3P/6Z9z994c+8owGofrsgsCRJUtkZ5CRJkkrKILdjc2pdgHrUnFoX\noB43p9YFqMfNqXUB6nFzal3A7tJvZ61GRHblHjlJkqRa6yi3OCInSZJUUgY5SZKkkjLISZIklZRB\nTpIkqaQMcpIkSSVlkJMkSSopg5wkSVJJGeQkSZJKyiAnSZJUUgY5SZKkkjLISZIklZRBTpIkqaQM\ncpIkSSVlkJMkSSopg5wkSVJJGeQkSZJKyiAnSZJUUgY5SZKkkjLISZIklZRBTpIkqaQMcpIkSSVl\nkJMkSSopg5wkSVJJGeQkSZJKyiAnSZJUUgY5SZKkkjLISZIklZRBTpIkqaQMcpIkSSVlkJMkSSop\ng5wkSVJJGeQkSZJKyiAnSZJUUgY5SZKkkjLISZIklZRBTpIkqaQMcpIkSSVV0yAXEWdGxJMR8UxE\nfK6DNpdV9j8cEUfsqG9EHBsR90TEgxFxb0Qc0xPvRZIkqafVLMhFRD3wTeBM4CDgnRFxYJs2ZwP7\nZeZ04EPAd7rQ9yvAX2XmEcAXKq8lSZL6nFqOyB0LPJuZ8zNzK3ANcG6bNucAVwNk5t3AqIiYsIO+\nS4GRla9HAYt379uQJEmqjYYannsSsLDV60XAcV1oMwmY2EnfvwB+FxGXUgTV46tYsyRJUq9RyyCX\nXWwX3TzuFcAnMvO6iHg7cCVwersHjrik1cs5mTmnm+eSJEmquoiYDczeUbtaBrnFwN6tXu9NMbLW\nWZvJlTaNnfQ9NjNPq3z9I+DfOiogMy/pdtWSJEm7WWVwac621xFxcXvtanmP3H3A9IiYGhEDgAuA\nG9q0uQG4ECAiZgKrMnPZDvo+GxEnVb4+BXh6N78PSZKkmqjZiFxmNkXEx4GbgXrgisycGxEfruy/\nPDNvjIizI+JZYD3wvs76Vg79IeBbETEQ2Fh5LUmS1OdEZldvVetbIiIzs7v330mSJPW4jnKLT3aQ\nJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOc\nJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOcJElSSRnk\nJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOcJElSSRnkJEmSSsog\nJ0mSVFIGOUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIG\nOUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIy\nyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWU\nQU6SJKmkDHKSJEklZZCTJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWUQU6SJKmk\nahrkIuLMiHgyIp6JiM910Oayyv6HI+KIHfWNiGsi4sHKx/MR8WBPvBdJkqSe1lCrE0dEPfBN4DRg\nMXBvRNyQmXNbtTkb2C8zp0fEccB3gJmd9c3Md7TqfymwqufelSRJUs+p5YjcscCzmTk/M7cC1wDn\ntmlzDnA1QGbeDYyKiAld6RsRAfwR8F+7921IkiTVRi2D3CRgYavXiyrbutJmYhf6nggsy8znqlKt\nJElSL1OzS6tAdrFd7OTx3wn8oNMDR1zS6uWczJyzk+eSJEmqmoiYDczeUbtaBrnFwN6tXu9NMbLW\nWZvJlTaNnfWNiAbgPODIzgrIzEu6W7QkSdLuVhlcmrPtdURc3F67Wl5avQ+YHhFTI2IAcAFwQ5s2\nNwAXAkTETGBVZi7rQt/TgLmZuWR3vwlJkqRaqdmIXGY2RcTHgZuBeuCKzJwbER+u7L88M2+MiLMj\n4llgPfC+zvq2OvwFOMlBkiT1cZHZ1VvV+paIyMzc2fvvJEmSekxHucUnO0iSJJWUQU6SJKmkDHKS\nJEklZZCTJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCT\nJEkqKYOcJElSSRnkJEmSSsogJ0mSVFIGOUmSpJIyyEmSJJWUQU6SpF0UhdERMajWtah/aah1AZIk\nlVlENMI+H4H9j4SVmyIGXJa5ZW6t61L/YJCTJGnXHAZvPhaGj4WJr8DmDwCf3rYzIgIYB6zLzA01\nq1J9kpdWJUnaNQ0wFFg1DLbWQ/3AP9y9x/lw0C2wz39ExJiaVKg+yyAnSdKueQSunwdbnoGbNsOz\n10TEkIiYHhHjYPhbYfJ+sOfrgCm1LlZ9i5dWJUnaBZm5ISK+CE9NAdYAW2Gff4exR8LGV+D5qyD3\nhqYngGdrW636GoOcJEm7KDM3A88ARAx/J0yZBfsPhQV7wLoTYf6pQFNmttS2UvU1BjlJkqoqm4CE\nzQnZArk1M7fUuir1TQY5SZKqav3N8PyNsGYmbHgJFl1a64rUd0Vm1rqGmoiIzMyodR2SpL6nWFuO\n8cDqzFxX63pUfh3lFoOcJElSL9dRbnH5EUmSpJIyyEmSJJWUQU6SJKmkDHKSJEklZZCTJEkqKYOc\nJElSSRnkJEmSSsogJ0mSVFI+okuSREQ0AHtQ/L+wJjPX1LgkSV1gkJOkfiwihsOQ82D6n8C4KVBX\nBxvWR+x5C7x0VWY+XusaJXXMR3RJUj8VEeNh6nfhdUfBmQPg8BYYlLCwHn7bDDevgkcvzlz3w1rX\nKvV3Pmu1DYOcpP4sIgbAlP+EC06EjwXstQUGthR7m4G1A+C39fCVVXDX+zOb7qhpwVI/11Fu8dKq\nJPVPx8FBx8KHAqZsgqX1cMsw2BSw3xaYvRlOboTnR8LCj0fE77O//uYv9WLOWpWkfmnShXDqIJjQ\nDMvr4MoRMG0gzBwAzwyHGwfDsK0wC9j7aGBarSuWtD2DnCT1SyOnw/R6GLIVHm6Ew+tgygAYOxDO\nTnhgMAQwOWHiQGBKrSuWtD2DnCT1ewG0ACRk5evWt+K01KIoSV1gkJOkfmnVXHiqGTY0wuFb4dEW\nmLcVXtoCPwOO2VCEuhcClm4Cnq9xwZLaYZCTpH5pyffg15tgcT2MboEPrIalG+GRzXD4WjhjE6xt\nhN8GLLwLeKHWFUvansuPSFI/FBGNsPe/w1tPgYvqYNJWGNRc7G0OWDMAfl0Hl66A+y7MbLq3thVL\n/ZvryLVhkJPU30XEGNjncjhmJpwxEA5PGJiwsA5+0wS3roAnPp+5/qe1rrUMIqIOxnwSBh0BSy7O\nTC9Hq2oMcm0Y5CQJImIIDHwjTHw3jH8N1NfB+rWw6CZY/p+Z+UytayyLiJgIh/wG9hoJD16W+fLf\n1bom9R0GuTYMcpL0qmI0idFAPbAuMzfUuKTSKZ6WsdfXYdBhsOATmU3317om9R0GuTYMcpLU9xXP\nk2VzZq7uofMFUJeZzT1xPvUfHeUWZ61KkvqkiFEnwkmXw1HfjIh9euKcWTDEqcf4rFVJUh+113Fw\nzmtg+UZ4ZD9cQkV9kEFOktRHPXsz/McRsHkVbH201tVIu4P3yEmS+qxiAgLNXu5U2XWUWxyRkyT1\nWZm5pdY1SLuTkx0kSZJKyiAnSZJUUgY5SZKkkjLISZIklZRBTpIkqaQMcpIkSSVlkJMkSSqpmga5\niDgzIp6MiGci4nMdtLmssv/hiDiiK30j4qKImBsRj0XEl3f3+5AkSaqFmi0IHBH1wDeB04DFwL0R\ncUNmzm3V5mxgv8ycHhHHAd8BZnbWNyJOBs4BDs3MrRExroffmiRJUo+o5YjcscCzmTk/M7cC1wDn\ntmlzDnA1QGbeDYyKiAk76PtR4B8q28nMl3f/W5EkSep5tQxyk4CFrV4vqmzrSpuJnfSdDsyKiLsi\nYk5EHF3VqiVJknqJWj5rNbvYrrsPtm8ARmfmzIg4BvghMK3dA0dc0urlnMyc081zSZLU70XEMBg4\nEyaeADEAVj0FK27JzCW1rq2sImI2MHtH7WoZ5BYDe7d6vTfFyFpnbSZX2jR20ncRcC1AZt4bES0R\nMTYzl7ctIDMv2ZU3IElSfxYRAUNOgSM/B4dNhbGNxfjLxoTHPxYx8VpYellmbqp1rWVTGVyas+11\nRFzcXrtaXlq9D5geEVMjYgBwAXBDmzY3ABcCRMRMYFVmLttB3+uBUyp9ZgAD2gtxkiRpVw05BU7+\nRzjrAJg0Ao4bCrOGwdQhcMwUeOMHYNLnI6KWA0d9Ws2+sZnZFBEfB24G6oErKrNOP1zZf3lm3hgR\nZ0fEs8B64H2d9a0c+krgyoh4FNhCJQhKkqTqiYghcMSfw+F7wPSB8KaB0NJY3Dn1uhZ4aiNcNwxm\nng8/vg54qNY190WR2dVb1fqWiMjM7O79d5IkCYgYcBK867vFSNynBsHGQfDKANgKjG6C0VvgFxvg\nwfVw7bWZ8z5d65rLrKPc4lCnJEnaCRNPKO6Je21dMRL3SiM8EzAA2NgIQ5rhkDp4uhGGHxYRdZnZ\n0vYoETESpr0HRh0MW1fCo1e3XlNWnTPISZKknRD1xecBldfNUcSKYcBmIAMGVj7X19POKhTFAv/T\nL4LPTIWpG6FuKHzxMxFxcWa2nQCpdvisVUmStBNWzIX1LfBkC9Q1w6gm2KOyb2IzDGiGZ1qgeSus\nX5CZze0cZDTMmAbDW+Cm2fCrGfD2Ohh0QA++kVJzRE6SJO2ENXPg0aUwcio8sgkOjeJyagKNLbB2\nI9zVAi9ugRf+o4ODNMOWgFFbYG0T7LsB1gds3dJz76PcHJGTJEndlpmvwLP/BcvWwy+a4CcbYMk6\nWL4Ofr8e/q0JmjbC7++DzXd2cJhVMPcOuG04nH8H7PUKfH8VND/So2+mxJy1KkmSdkpENMKET8Fx\n74LXjIBhjcWeLU3FSNyd98BTn8nMFzs5RgOMPh32OgjWLYcFPytColrrKLcY5CRJ0k6LiDrgAJjy\nDhhxDEQDbFoA8/8Ttt6dmRtqXWNfYJBrwyAnSVJ1FY/sguyv4WI3ch05SZK0Wxngep6THSRJkkrK\nICdJklRSBjlJkqSSMshJkiSVlEFOkiSppAxykiRJJWWQkyRJKimDnCRJUkm5ILAkSVVWPLZq7Nkw\n5VTYvBqeuDozn691Xep7HJGTJKnqBh4L510A3wYuHQdH/VlEDKl1Vep7DHKSJFXd+Glwxga4ZRq8\nNAhmDAbG1Loq9T0GOUmSqm75fPjVEDhtHkzYCM9sAlbUuir1PdFfn28bEZmZUes6JEl9T3GP3Lg3\nwaTTYOtKePzqzJxX67pUXh3lFoOcJElSL9dRbvHSqiRJUkkZ5CRJkkrKICdJklRSBjlJkqSSMshJ\nkiSVlEFOkiSppAxykiRJJWWQkyRJKimDnCRJUkkZ5CRJkkrKICdJklRSBjlJkqSSMshJkiSVlEFO\nkiSppAxykiRJJWWQkyRJKimDnCRJUkkZ5CRJkkqqodYFSFJfFhF1wAzY67Uw5kCI/aFlKDRvgA2P\nwtLfQtNDmbm+1rVKKp/IzFrXUBMRkZkZta5DUt8UEQENh8NB74Qjx8GRo2HIZNi7AcYGrA94YgA8\nvgnuegoW/Rhe/Glmbqx17ZJ6n45yi0FOkqosIgbDlAth1uvgz16B3ANWHQ5HroNRW19tmcDLg+DO\nQXDvPPjpPHjk25k5r2bFS+qVOsot3iMnSVVUhLj9PgUXHQdXzIdJCS+/Fl635g9DHEAA4zfBSevh\n+L3hGwPhpL+MiP1rUbuk8jHISVKVFJdT93k/fGQ/+LMFMCDh2SkwFRjU/GrLdQ2wsf7V16O2woQG\nGDME/mk1HP+piBjX0/VLKh+DnCRVzcCj4eTj4JMLX/3xunJvmLTh1TYLh8Jv94fbZsArA1/dPmkL\nvDIJjlgLFwHTLqxMlJCkDvlDQpKqICIa4MB3wWeWQUOrm49bGmBAy6uvXxkKYwbD8MGwYtCr2xtb\noHlA8fUFL8Ks1wIH9Ejxkkprp4JcRAyMiEkRMXDHrSWpXzgQjh8JB7dZRqRxE6xvtdTTPqtg7Upo\nWgGTWrXd0ACNlZG7OuBtG2D66bu9akml1q0gFxFHRcRtwDpgAXBCZfueEXFrRJy2G2qUpBLY92g4\nffP228fNgwWDoRn40Rj429fAdUPgwWZY26rdCw0w6YVXX5/+Mow7LCIGtT2iJG3T5SAXEYcDvwWm\nAd+jmG4FQGYuAwYD76l2gZJUDiP2h4PWbL99xmJYtBWungBPT4P3joTPDoOJE+Fr02BzwMIhsGYt\nTHvl1X4DEvZLYEKPvQVJpdOdEbm/AZYChwCfa2f/r4Fjq1GUJJVJZfHfPeE17SzmO3wLzLgb5uwJ\n5wyA6XUwCTi1HsaMhOvHwyPNcPTdUN9mYc99Acb0wFuQVFLdCXInAt/NzLUd7F9A8dNJkvqbuuLH\naUMHK6yPWw0Na2F4S7EIcAIRMKQe5i+B1/8GxrfziK5GgPrtt0tSoTtBbhCwqpP9I3axFkkqqxZo\n3gKrOnh+9YgmGLakeBzXpi2wqQlWboanX4GzHoSRm9rvtxKgg32SBB380GnXPOCoTvafDDyxa+VI\nUvlkZkYcMg8e3wtOaOcX3jrg3b+Hy0fBM+NhRMBDG2DqnXDouo6P/GQAL+6uuiWVX3eC3PeBL0TE\n/wAPbNtY3BvCnwFnAZ+sbnmSVBZLHoMHZ7Qf5ACOWQOTr4Wb94S1A+GdL8Oxqzs+3vxBsHgdsHy3\nlCt1Q0TjwXDAn8Lqx2Dh1ZnZvONe6gndubT6VeBO4GaK2asAXwMWA/8I/BL4dlWrk6TSWHkv/LQO\ntmz3UOtX7bUF3rsQLnoWZq7u/EfwT8bB/F9mZksnjaQecuBb4YRDYd9zgIk7ah2x5xsjpn8vYvDh\nPVBcv9blIJeZm4E3AJ+muGdjE7A/8DLwWeBNJnRJ/VVmvgTP3Af/tdeuH23hQPhxC6z6/a4fS6qG\nJXfB/cth6aMU/+/vwIjTYfpBMP6Y3V5aPxeZHUyy6uMiIjOzk9+cJal7ImIMHP9FuHoNTN+w4x7t\naQE+tS9cdXXm6lurWqC0CyJiBLAhM5u60HZfGHc0vDwnM7sQ/LQjHeUWg5wkVVHE4CPhjE/CZUth\nSjdnnLYAX9oH/vVBeOFbXuWQtE23g1xEvIdisaNuyczvdb+8nmeQk7S7RAw9HmZ9EP58HZy8omu9\nlg6Af5wE1z0A8/+lcjtLvxMRg4GxQBOwLPvraIPUxs4EuZ26wTYzu/X81loxyEmqhoioo1iBJNts\nnwYHfxDeOBHOXwFHdzC5YekAuGE8XAM89ENYdWtXLl31NRExDCaeA1Nmw7Q6WF8H85bB09fC5vsM\ndOrvdibIzW6zqRH4MsXjYv4FmFvZfhDwEeAV4HOZeUuVat6tDHKSdkVENMLEC2Cv2ZAtMP+nsOLn\nrWeZRsQAGHAUzDgbJk+CAxKm1hUPa1iRMDfh2S2w8BZYdnsxYaL/iYihcMDn4COT4V1LYY+txZ4H\nh8PXxsHP/ytzxU21rVKqrV2+Ry4i/gZ4GzAzM9e02TcCuBv4YWZe3I2izgS+TvFT7d8y88vttLmM\nYo26DcB7M/PBzvpGxCXAB3h1Vs3nM/MX7RzXICdpp0VMOB/e8xZ43Wioa4Ib1sD/XJG56jfbt40A\nRlMs2zCKYmhuI8Xzq5dl5tYeLb6XiRj3Fvjzt8Bn58PmelgwBoZsgUmr4ZVGeMde8OvP9ZWb5it/\nHyYC4yn+r1rsiKN2pKPc0p0Fgd8L/HPbEAeQmWsi4krg40CXglxE1APfBE6jWIvu3oi4ITPntmpz\nNrBfZk6PiOOA7wAzd9A3ga9l5te68d4kqZv2OhaOHQYPzIaWFjjlRvj9kcB2Qa7yn/SKyodaiYgG\nOOo0eMcpe4j4AAAgAElEQVRS2FoHvz4RYhps2QLLfwOHLoS3JDx0HPCzWte7q4qR3L0/CvufB9OH\nwrwN8MQNEfHNzNxS6/pUPt25n20cnT+8uR7YsxvHOxZ4NjPnV34bvQY4t02bc4CrATLzbmBUREzo\nQl9H2iTtZk3roGUDNGyGoRuL56du6uRxW+rAUBg7GPbeDC8Nh+YpMGMcTB4Pi/Yvmuy/EcZOqW2Z\n1TL8dDjrvfCV0fD5RvjyKHjTu2HkGbWuTOXUnSD3JPCBYp2kPxQRY4EP8up9c10xCVjY6vWiyrau\ntJm4g74XRcTDEXFFRIzqRk2S1EWP/Qgu3wx73Qgjfgn/uhbmeR9X922B9cCmOhi2CZrWw7LNsGoT\nDKk8nmxVI2zuIyF58jlweiOMHQrr94NRw+CMepj0llpXpnLqzqXVS4DrgCcj4t8pgh3AgcD7KCZB\nvK0bx+vq/QDdHV37DvA3la//luLRYn/azWNIUqcyc25x7/AjR0JzE6y4JzN9wH03ZebGiP0ehF8c\nDG9ZBkfcBvP2g0Hr4cinirX1fj4AXrin1rVWR9RBQwtkHdTVFZ8boXsDK9L/1+Ugl5k/iYjzgW9Q\nPJKrtUXAH2Xmdd0492Jg71av964cp7M2kyttGjvq23rWV0T8G/DTjgqoTIzYZk5mzuly9ZL6vcyc\nD8yvcRm7rLhPbeKHgBZY8t1qLkRcHJuDYNgkWPci8Nj2kzue+xl888jinrGDV8DUSmhrAf5lMtzx\nNPB0tWqqrSW/gDnHwYw1MGoNrAZuqYNl203KU/9WWT1k9g7bdXeiTGWiwVHAtMqm54D7u/tg58o/\n7qeAU4ElwD3AO9uZ7PDxzDw7ImYCX8/MmZ31jYi9MnNppf+ngGMy84/bOb+zViUJiIg94Mh/h+YW\nePg9mbmqSsdtgKkfh2NnwpBRsHk13PMAPPf1tgseRzQeAod/FE4dAodvhfX18Is6eOAxmHd5ZvaJ\nS6vFgsdTPwtHngEzhsGz6+H+W+D5L2XmTj7WTf1Br3xEV0ScxatLiFyRmf8QER8GyMzLK22+CZxJ\ncRPF+zLzgY76VrZ/Dzic4tLt88CHM3NZO+c2yEkS25bDaDgcaMnc+nAVj3skvPNz8NypMGg4bFwP\nB9wK//21zM2/b6f9IKg/FMZPgY2bYdWjwAt9bWmOyiLS0ykmEb4CPN3dwRD1P70yyNWSQU6Sdq+I\niW+Hkz8EL86C4+vggYQhv4cHv5f53JW1rq+vqqxTN4liJYmXgYV9LQz3R7u8jlzlkV3J9pMPtv3l\nCIrlkjpbokSS1G8sXwKNq2HdZnh0EKzcDONXw8sLal1ZX1Xc/jT53XDkSXBUCzxcB/fcGRFX9sdH\nv/UH3Zm1+r0O+k8DZgIPAw9VoyhJUs+JiJEw7nxYeUfm1qeqd+QtD8Ddc2FmIywfBwcsh/uehLV3\nV+8cauMweOMp8L82wkv7wnnz4d9eD996BLir1sWp+qpyaTUiXgfcALwpM0vxF8VLq5JUiBj8Jjj6\na7Dwt5nzP1DdY8dQGH4CTJoBy56Hlbe394QgVUfEfh+Arx4NdcfAkc3wYAPU3wWfejLzyctqXZ92\nXjUe0dWhzPx9RFwFfBk4qRrHlCT1lE33wfwfwqqbq33kzFwP/LLyod1u/VpYWQ/j18PiERDrYFUD\nbFxd68q0e1RzAcJngKOreDxJUg/IzBczF/7fzLW317oW7aoXfwf/UQcDH4Ohd8GwR+DqBliw3TOA\n1TdUbdZqRPwAeENm7lGVA+5mXlqVJPVFEQMOhP3fCUOmwIbF8OR/ZW59rNZ1adfs8vIjEfEe2n+s\n1hjgdOAsivXcPrgrhfYUg5ykWouIKbDvWTD6MGjZCAvmwIrb+srit6qdyhIk9UCzS4/0DdUIcp0t\nVtgEXAV8qnI/RK9nkJNUSxFxAJz8GThvIGwcDHXNsHkzfP8FmPulsvwsldQzqjHZ4ZR2tiWwAnje\n3yAlqWuKlf1f+1744wb43Rlw0AhYm7B0HrwD+Nos4KYalympBLoc5HygvCRVzUTYfzw8egjMGg1H\nJ6wJ+MU0GPooTJ2NQU5SF3R51mpEPB8R53Sy/80RMa86ZUlSn9YAQxKaB8DAhJFZvB4cQB3UDah1\ngZLKoTvLj+wDDOtk/1Bg6i5VI0n9w1J4YjNMexLu2gK3B9wHPL0S6jfDSw/WukBJ5VCVBYErxgMb\nqng8SeqTMnNzxOjr4LF3w2t/DvdMh8ZNcOKTcNUWWPyrWtcoqRw6DXIRcRLFkxq2zZJ4a0Ts107T\nsRR36PqsVUnqklW3wPUt8PB5MGMdrAu46QWY+73MXFrr6qSdVVn6pBHY6tInu1+ny49ExCXAF7p4\nrGeBd2XmvVWoa7dz+RFJvUFENAJ7AluAl/2PT2VUzMRmP5h2Kow6GgbUwZZmWHEPzL8VeM6/27tm\np9aRi4iRwOjKy3nAp4CftGmWwLrMXF6lWnuEQU6S+peIqAf2oLg/fFVmbqxxSX1CRAyGqR+C1x8B\nb90EZ7wEQ1pgXT38YjxcOwDuuBcWXJmZm2pdb1lVY0Hg2cATmflSlWurCYOcJJVHRAwAplEsQP98\nZjZ3r+/wWfCas2DqKGhMWNhSPH90yc2Z+cou1nYwxWTA+7pTV19QjChP+yR87CD41IIiI7cAaxpg\nRFPxuingK1PgXx+EF76dmU01LruUdnlBYNeRk6TtVa5czACagbmO8lRfcW/2EZ+A44bCxoB7lkbE\nN7oysBARA2HaRXDBofDuZXDgwmLPK41ww8nw7eMi4suZuXgna5sER30dRo6EWy8C7t6Z45TXgGPh\nrYfCp+ZBS8CVU+GOg2HzMBi8Fk58DC58Af7iBVh6FHzncIop2qqSDoNcRFxMcdn07zOzudXrTmXm\n31SxPknqtSJiPBz+l3DGBNic8ItnK6HAJ91USTHic/hF8PU62DAJGjfDmU1w8fuBL+34CBPOhfce\nAv9nXjE6tLUOmupgj63w/kUwZSz8709ExF+2HU0r/nwZCCzq5P6uTbBhVeXrdZV+dTDhbTD6ZFh9\nJyz5QV8chSomNbz2jfCOV4rv7T8dAC+cALNHwaQoRj3v2xPW3g4XPQMXrIJbz46I+71frno6G5G7\nuPL5SxS/aV7cSdvWDHKS+ol9z4X3T4YZr4HmFhjTAJfOAm6sdWV9yD5w1DDYuBfUnwAbWmDq9TBp\nRkSMyMw1HXUs7t067lT48OIiaNy/LzxxMrQMgHGPwhvuhNOWwwn7wOMHAo+16vtaOOp7MHQIPPQF\n4L/bO0dmLo+IjwCDWo3qHQQzPgqD9oO9DoMlc4FSTATspkkwYwIcsQAWD4QHDoW3jYLXRzHZ4VBg\n2Bi4/nB453x43UrYdyo8MQ7oE7dp9QadBblpAJm5pfVrSdI2Q0bCqEaYUA8b62FsPYwcVeuq+pit\nxdIsgzbAyuZiJmTjFthcR3G/XGdeA8c3wPitsHwIPH4mjBlXDLK9cgI8sBSOnQdnbIZfHk2rIAeD\nj4JDxsOYgIWz6CDIQRHm2mxqgGiAgXXQ3EB112ztTYbCXi1FSH5sOIwfAlProKEeBgMtdbBvM4we\nCo8Ph5NWwJ4tFPcTGuSqpMO/XJk5v7PXkqQlD8Jvj4LhS6C5GeZshKWPVvMMlTW5GoCmfno5agE8\n9AI8tzccei1EE9w0ABb+NjN3tAj9QBhZ+XLdQIghMAEYF7C2EdZUQveIrTBoyB923Xgz3Hk+DBgF\nC67qZs2PwbwfwJBTYNOdwAPd7F8WTbBtEuqIJljXDOuBaCl+sYnm4mrzhiYYubVot6XSr/qK+1Xr\nD4GBw2DDMuDxzNy6O87Vm3T5t4SIuA34u8z8dQf7Twb+KjNPqVZxktS7rbwNrh8G978BsgmevRaa\nHttxvx0rAtzIWXDYuTB4NKxbGNH4w8ytVTl+WWRmRsRl8MULi8uUTXWw6DZY1OEIWSurYX5llt8e\n66DuRZg/FRbUwfq1cMiiYt+CwbByWZvzLo2Ic4G67t7fVmn/nYi4PDNbutO3ZJbBEwmrGuCY1XDl\nUnhgNAwaWNyD+HLA/Ztg8GI4dC281AhPN1Hl0bji38ros2DW+XBWHYxLeLQOfrM6Iv45M5+t5vl6\nm+4M954EfLeT/XsCs3epGkkqkcp/0j9h+/U1q2DM6XDen8Bb1sL8epg6Bi7/bER8KTPnVv98vVdm\nrgS+ERFDgObM3NzFrvPgoRXw2DA4ZB2c/hO4+2hoHgyHPQozXiqWxri+HpZtN9u08ue700Gsj4c4\nMnNdxN63w89OhD9ZDP/rLvj6QJg/EcYPhGUbYfUS+PQ9xeXXG/aEBb+s/lpyw0+Ac94BX1sAY1qF\n7vtGwCc/GxH/NzNfru45e49qXrcfCXT1H5ckqQPFTM2jz4Vz1sGlfwKThsC1y+GT18Jz5wL9Ksht\n04VLqW3bt0QM/RH8w8fgG5th/Dp485xXW7QAX50Cj9wH7NTyI1r0a7hiFhw3BA5fC5fdBDftCS8N\ngRPWw5kvwbDmIkxf3Qwv/qaaZy8WeT7sLfDZZUWIu3U0fHU2fOx38MaX4d0j4ZlZwI+red7eZEfP\nWj0MOIxXn7V6YkS012cs8DHgieqWJ0n90nAYNwjuGgsHDIGT6uGa0bC5HgZPqXVx5bLhbvj5SFhx\nAZzfUtxwP7AFHh4B1w2BOQ/BC1f29vsPK0uh7Am82JtGlzJzccTgb8LHPgF/vgFOXg4XLHm1RVPA\njePg0kFw11czc1nHR9spe8C+o+HghUUwv/x4WD0V/uclWDQEbt8XGj4UEQ9k5vNVPnevsKMRufP4\nw2etfrjy0Z61wCeqUZQk9XNr4eVNMHMJ/NMGWDsE1qyEgc2wcUGtiyuTSkC7OSIegYdmwYQjIRph\n3RPwzK3AU735aQzF/V8T/gjecCYc1gIP1UVM+Bksu7a3hM/MjQ9FxN/DovNh2oFwWsCIFlhVB7ck\nPP8YPPXjnpk0OW4zTF4Adx8AQ46EWfUwowV+e1PEqG/D6n/uLd+3atnRs1anAlMrL28Fvgjc0qZZ\nUkxLebxMz1DzEV1S71d5EHefv9eoPRFj3wBvfTectwaeHQHTVsPlw+BnX+5v98h1pLgETSOwsa/9\n57xN8fivd30OLhwFK14DY5+GK9fANV/MzKdqXV9bETEReA0MGgybNgJP74ZRuNbnq4fDvgLfb4GD\n18PCgfCRs+DUKXBUPWwZCFs2w9Bl8A3g+vdm5u93Vz270049oquSnudXDvB+4Dd9dWhSUq/UJ/9z\n7poVv4IfboF7z4HBTbBuJTx5uSEOImI0TH4bHHYWNAyE9XMjGv8jc+vDta6t+iYeAic1wAuHwdNj\nYcZQmH0H/OogoNcFucxcAizZYcPqna85Yvj1cOkH4asLoC6hbiIckzBoELymDpoaYN5YeNMqeOBP\ngFIGuY5051mrV+3GOtRNETEGmAgsz8ylta5H2h366ihLV1Te+5yI+A39ex25PxARw2D/v4UDzoGW\n0dBYB5tnwoRjIoZ9NnPdPbWusbrWroDlzTBlNYweDsPWwHPNsH7Vjvv2F+t+Bz8dAc+/FY4YXix/\nsqkeRgZM2AgbgWUNML4ZRkyudbXV1u1ZqxFxDHAsMJpiPvEf8Fmru1dxv8S4N8PJ58KRwJN1EVPu\ngIVX9cVn+Un9XSW89flFTbtu6PEw/Q0waRzMqitupl9VDzcdClO33dTeh34Wrr0XfvRm+NAjxZIp\nixN+vAI23V/rynqLyr+Rn0fE7+G3R8IRR8LWgbB2BLxYVwS5NZthYT0s73WjmLuqOwsCDwauA96w\ng6YGud1rOsw+H766DF4eB3uugEtPgm8/Cfyu1sVJ0u61z0kQe8BxdbCoAfZJGBIwehBsmA6PTwAW\n1brKasnMVRHxRfiHN8Ko/WDl0/DCzzt7xmx/VVlv8NcRe/w73PKncM4r8Mjw4kkTdavhx82w9Kpa\n11lt3RmR+wJwOvB3wK+B24D3UqzQ/BfAEODCKten7ex5KJy7BZ49GIbuAw+vhPMegF8ej0FOUp+X\nzRT3TmYxULk+oKHydUszu7CAb2+VmS8CV9S6jvJY/o9wzRh45I1w7CpYDdy1Dh69ODMfr3V11dad\nIPc24EeZ+YWI2KOybVFm3hoRtwD3UQS7v6hyjfoDWzbC6jrYd33xAOkB64tnBm5dX+vKJGn3m/9r\n2PftMGdfOL0J6gKeb4b162HZE8CLta5QtZWZW4BPF492+/URwAbg9szcWOPSdovt7nHrxN7AnMrX\n29bcGQD//7lyPwAuqFpl6sDKe+GaZmheCofcAuMfhyuGwjO31boySdr9Nt4DT14HyxfB9zfCf26G\nO9bAK3fDs9/pj0vVqH2Z+UJmXp+Zv+yrIQ66NyK3tlX7tRTD1xNb7V8D7FWlutSBzHwpIi6Fj74b\nRk2GdSvgmW+7JIGk/iAzNxWLz664CyadC/XDYN09MO/6zJxX6/qkntbpgsB/0DDibuCezLyo8voR\nYHFmnlVZtPMmYFpmTt9t1VZR2RcELmavMhDY4m+gkiT1bR3llu5cWv0V8LZiFWUA/gU4IyKeA56h\nmAjhzZg9JAubDHGSJPVf3RmRGwZMBp7LzK2VbX8GvBtoAn4EfKUsC1aWfUROkiT1Hx3lli4HuS6c\n4CPAJzLzoKoccDczyEmSpLKoxqXVHdkDOKCKx5MkSVInqhnkJEmS1IMMcpIkSSVlkJMkSSopg5wk\nSVJJdfpkh4j4NMXDibvidd1oK0mSpF3U6fIjEdHtxWYzsxSjfC4/IkmSyqKj3LKjZ62e0s3zOCIn\nSZLUQ6q2IHDZOCInSZLKoicWBJYkSVIPMshJkiSVlEFOkiSppHY02UFSLxERdcBooDkzV9W6HklS\n7TnZQSqBiMGHw4x3wOTxsAVY9Aw8+f3MXFDr2iRJu5+THaSSimg8BE7/33DlAPj5ArhpIXx1Msz8\nfESMr3V9kqTa8dKq1ItFRMDBb4fPrYKj1sA1E2H8Jjj7ZXhxMrxwCnBNreuUJNWGI3JS7zYYxu4N\nx6+Cm8bBVW+Fv3k7bAqYtRwmHFHrAiVJteOInNS7NcHmFthQD9PWQ/NaGLkOGhJWNcLW1bUuUJJU\nOwY5qRfLzC0RU++C/zkO3rcYfnV1sacF+NFYeO4nNS1QklRTBjmp13vhWvjaDHhxHzh5FWysh5+O\nhOsfgY131ro6SaqWiKiH+qPggFNhwChY8Qy8cEtmzq91bb2Vy49IJRARw2D48bDP0dCyFZ66HZof\nyMytta5NkqqhCHGTPwhHnANjxkNLI8RamLcQ7r00c9P9ta6xljrKLQY5SZJUcxFxOJz1t3DkFHj/\nepi8Fe4aClcn3P8QPHxRZm6qdZ214jpykiSpF5txIowfD+/cCNO2wICEWevgxAaYshewf60r7I1q\nGuQi4syIeDIinomIz3XQ5rLK/ocj4oiu9o2IT0dES0SM2Z3vQZIkVcOAYUAD7NH0h9vHJgxqBAbV\noqrermZBrrgWzjeBM4GDgHdGxIFt2pwN7JeZ04EPAd/pSt+I2Bs4HXihB96KJEnaZUsehy1r4Pah\nr27bFHBXHSxdAyyqWWm9WC1nrR4LPLttJkpEXAOcC8xt1eYc4GqAzLw7IkZFxARg3x30/Rrw54BL\nM0iSVAor7oCn3gw3HAovjISJwP0BC1fAwt9k5uJaV9gb1TLITQIWtnq9CDiuC20mUfzptts3Is4F\nFmXmIxHOZZAkqQwyc3lEfBE2fgSWTIchjfDKBlh8Kyz4Xq3r661qGeS6Ol22y2ksIgYDf0lxWbXb\n/SVJUu1k5vMR8XmYOxUYAryYma/UuKxerZZBbjGwd6vXe7P99e+2bSZX2jR20Pc1wFTg4cpo3GTg\n/og4NjNfaltARFzS6uWczJyzE+9DkiRVSWa2APNqXUetRcRsYPYO29VqHbmIaACeAk4FlgD3AO/M\nzLmt2pwNfDwzz46ImcDXM3NmV/pW+j8PHJWZK9o5v+vISZKkUugot9RsRC4zmyLi48DNQD3w/9q7\n92i/yvrO4+8PCeFuwmUgAkFQLhKsiK1I1WIqVmNYwqjjOIytAnbpGku7HC9Dkc6Utna66mV1qi6t\nTtUytWDVQUUFJVNNrY4FUWQoJIYIkXCtIiACkpB854/fjuvnj985OZCcs8+T836ttdc5+9nPs/ez\n82SHD/v6kapaneQN3fIPVdVlSVYkWQc8AJw1Wdtxm5mRnZEkSeqBX3aQJEma5fyygyRJ0k7GICdJ\nktQog5wkSVKjDHKSJEmNMshJkiQ1yiAnSZLUKIOcJElSowxykiRJjTLISZIkNcogJ0mS1CiDnCRJ\nUqMMcpIkSY0yyEmSJDXKICdJktQog5wkSVKjDHKSJEmNMshJkiQ1yiAnSZLUKIOcJElSowxykiRJ\njTLISZIkNWp+3x2QJEltSbILcDQcfDzstgf8+E6471tVdXfffZtrUlV996EXSaqq0nc/JElqSZLF\n8NTfgWccDgfvAynY8ghcez+sXgl3fKqqHum7nzubiXKLZ+QkSdKUJNkXTjgXXnQ43HkkLNwH9poH\n6x+Cp94GR58Gl84DPt53X+cKg5wkSZqig06BXz8MfrIUXrkQnrYF9tgCtz8BVu0FNxUcf0qSr1TV\n7X33di7wYQdJkrRNSRbAYafAI4vgmQvh5I2wZU/40f5w+BZ4/jzYfCj8yh5w0HP67u9cYZCTJElT\n8QQ4aAHcsz88pWDLfHjwCbDX7nDbfrB4MyzZbXCxb/8n9d3ZucIgJ0mSpuIReDiw+8Nw9y6w62Z4\nZDPcX7D7Rngo8NPNML9g86a+OztXGOQkSdJU3Ad33AYH3gY3PAw37wpH3gmH3AX73gfXBm79Cdyz\nCW6+uu/OzhU+7CBJkrapqirZ/Ytw15NhwU3w6SPgmN1hYcH6+fCde+AZ18Hf3Q8br+m7v3OF75GT\nJElTMngR8KFnw6kvgEPnwe2HwsO7wRPvhP3uhS88DF9/d9XG1X33dWczUW4xyEmSpClLMg/2Ohme\nciocvS/sU/CDXeCW78K6S6tqfd993BkZ5EYY5CRJevwGgY7FwALg3qq6p+cu7dQMciMMcpIkqRUT\n5RafWpUkSWqUQU6SJKlRBjlJkqRGGeQkSZIaZZCTJElqlEFOkiSpUQY5SZKkRhnkJEmSGmWQkyRJ\napRBTpIkqVEGOUmSpEYZ5CRJkhplkJMkSWqUQU6SJKlRBjlJkqRGGeQkSWpUBvZNMr/vvqgfBjlJ\nkpq1/ynwjL+BJ725756oHwY5SZKatfeTYOlC2O3YJOm7N5p5qaq++9CLJFVV/qWXJDUrySJY9Dy4\n94aquqnv/mj6TJRbDHKSJEmz3ES5xUurkiRJjTLISZIkNcogJ0mS1CiDnCRJUqMMcpIkSY0yyEmS\nJDXKICdJktQog5wkSVKjDHKSJEmN6jXIJVmeZE2SG5OcO0Gd93bLr01ywrbaJvmTru53k/xDkiUz\nsS+SJEkzrbdPdCWZB3wPeCFwG/At4IyqWj1UZwVwTlWtSPJs4C+r6qTJ2ibZp6ru79r/LnB8Vf32\nmO37iS5JkjqD77byVGA+g/+2rq+5+h3PWWii3DK/j850TgTWVdV6gCSfAE4HVg/VOQ24EKCqrkyy\nKMli4IiJ2m4NcZ29gR9N835IktSswcmRxa+A578YnrYPzJsHGx6A69cm+WBV3d13HzWxPoPcIcCG\noflbgWdPoc4hwMGTtU3yp8BvAQ8CJ+24LkuStLM58KXwkpfDnoth/4WwT8HDm+GYPeDStyS5oKo2\n9t1LjdfnPXJTPV37mC9/VtX5VXUY8DfAXzzW9pIkzQVJ9oInnwr7HAjPWwSP7A93LYYzgAefAs89\nHOY9ve9+amJ9npG7DRh+EGEJgzNrk9U5tKuz6xTaAlwEXDZRB5JcMDS7qqpWbavTkiTtRI6ApXvA\ngn1hw95wwGJYGPj8LvC8O+DKveDJzwKu7rujc02SZcCybdXrM8hdDRyV5HDgduBVDP4XYNilwDnA\nJ5KcBNxbVXcluXuitkmOqqobu/anA9dM1IGqumAH7YskSS3aBXbZBeYX7LYF7uuugi3YAgu2VvFV\nZT3oTi6t2jqf5A/H1estyFXVI0nOAb4MzAM+0j11+oZu+Yeq6rIkK5KsAx4AzpqsbbfqP0tyDLAZ\n+D7wn2Z2zyRJasYtsPYh2O2ncMoW+NIdsGEB/Pvb4Yp58OBDcMt1fXdSE+vt9SN98/UjkiRBsuS1\n8JLT4aZnwQl7wwGBqzbDLqvhpm/Dd95eVQ/03c+5bqLc4ulSSZLmtFv/Hj53J/xSwX4L4O494AVb\n4JaF8J2/MsTNbgY5SZLmsKr6Gey3Fo68GW6eB9kF7tgMx1wH+NqRWc4gJ0nSnHfnOvjXe+Dwu2DL\n/XDELXDzz4Af9t0zTa7Pp1YlSdKscO8/wid/FV60CRbvDp+9D274ZFX9pO+eaXI+7CBJkkiyO/BL\nsOtC2HQjcIvfWp09JsotBjlJkqRZzqdWJUmSdjIGOUmSpEYZ5CRJkhplkJMkSWqUQU6SJKlRBjlJ\nkqRGGeQkSZIaZZCTJElqlEFOkiSpUQY5SZKkRhnkJEmSGmWQkyRJapRBTpIkqVEGOUmSpEYZ5CRJ\nkhplkJMkSWqUQU6SJKlRBjlJkqRGGeQkSZIaZZCTJElqlEFOkiSpUQY5SZKkRhnkJEmSGmWQkyRJ\napRBTpIkqVEGOUmSpEYZ5CRJkhplkJMkSWqUQU6SJKlRBjlJkqRGGeQkSZIaZZCTJElqlEFOkiSp\nUQY5SZKkRhnkJEmSGmWQkyRJapRBTpIkqVEGOUmSpEYZ5CRJkhplkJMkSWqUQU6SJKlRBjlJkqRG\nGeQkSZIaZZCTJElqlEFOkiSpUQY5SZKkRhnkJEmSGmWQkyRJapRBTpIkqVEGOUmSpEYZ5CRJkhpl\nkJMkSWqUQU6SJKlRBjlJkqRGGeQkSZIaZZCTJElqlEFOkiSpUb0HuSTLk6xJcmOScyeo895u+bVJ\nTgSv/04AAA6HSURBVNhW2yTvSrK6q39JkoUzsS+SJEkzqdcgl2Qe8H5gObAUOCPJsSN1VgBHVtVR\nwOuBD06h7RXAcVV1PLAWOG8GdkeSJGlG9X1G7kRgXVWtr6pNwCeA00fqnAZcCFBVVwKLkiyerG1V\nrayqLV37K4FDp39XJEmSZlbfQe4QYMPQ/K1d2VTqHDyFtgBnA5dtd08lSZJmmb6DXE2xXh7PypOc\nD2ysqoseT3tJkqTZbH7P278NWDI0v4TBmbXJ6hza1dl1srZJzgRWAKdMtPEkFwzNrqqqVVPuuSRJ\n0jRJsgxYts16VVM9KbbjJZkPfI9B2LoduAo4o6pWD9VZAZxTVSuSnAT8j6o6abK2SZYD7wGeX1U/\nmmDbVVWP60yfJEnSTJoot/R6Rq6qHklyDvBlYB7wkS6IvaFb/qGquizJiiTrgAeAsyZr2636fcAC\nYGUSgG9W1RtndOckSZKmWa9n5PrkGTlJktSKiXJL3w87SJIk6XEyyEmSJDXKICdJktQog5wkSVKj\nDHKSJEmNMshJkiQ1yiAnSZLUKIOcJElSowxykiRJjTLISZIkNcogJ0mS1CiDnCRJUqMMcpIkSY0y\nyEmSJDXKICdJktQog5wkSVKjDHKSJEmNMshJkiQ1yiAnSZLUKIOcJElSowxykiRJjTLISZIkNcog\nJ0mS1CiDnCRJUqMMcpIkSY0yyEmSJDXKICdJktQog5wkSVKjDHKSJEmNMshJkiQ1yiAnSZLUKIOc\nJElSowxykiRJjTLISZIkNcogJ0mS1CiDnCRJUqMMcpIkSY0yyEmSJDXKICdJktQog5wkSVKjDHKS\nJEmNMshJkiQ1yiAnSZLUKIOcJElSowxykiRJjTLISZIkNcogJ0mS1CiDnCRJUqMMcpIkSY0yyEmS\nJDXKICdJktQog5wkSVKjDHKSJEmNMshJkiQ1yiAnSZLUKIOcJElSowxykiRJjTLISZIkNcogJ0mS\n1CiDnCRJUqMMcpIkSY0yyEmSJDXKICdJktQog5wkSVKjeg9ySZYnWZPkxiTnTlDnvd3ya5OcsK22\nSV6Z5Pokm5M8cyb2Q5Ikaab1GuSSzAPeDywHlgJnJDl2pM4K4MiqOgp4PfDBKbS9DngZ8LWZ2A9J\nkqQ+9H1G7kRgXVWtr6pNwCeA00fqnAZcCFBVVwKLkiyerG1VramqtTO1E5IkSX3oO8gdAmwYmr+1\nK5tKnYOn0FaSJGmnNb/n7dcU62U6Np7kgqHZVVW1ajq2I0mS9FgkWQYs21a9voPcbcCSofklDM6s\nTVbn0K7OrlNoO6mquuCx1JckSZoJ3cmlVVvnk/zhuHp9X1q9GjgqyeFJFgCvAi4dqXMp8BqAJCcB\n91bVXVNsC9t5Nq9LxJojHO+5xzGfexzzuWdnHvNeg1xVPQKcA3wZuAH4+6paneQNSd7Q1bkMuCnJ\nOuBDwBsnawuQ5GVJNgAnAV9Mcvl2dHPZdrRVe5b13QHNuGV9d0AzblnfHdCMW9Z3B6ZL35dWqarL\ngctHyj40Mn/OVNt25Z8BPrMDuylJkjTr9H1pVZIkSY9Tqqb64OjOJcnc3HFJktSkqnrUff9zNshJ\nkiS1zkurkiRJjTLISZIkNWpOBrkk+yVZmWRtkiuSLJqg3vIka5LcmOTcofJ3JVmd5NoklyRZ2JUf\nnuShJNd00wdmap80ueka827ZeV39NUleNBP7o23bAWP+yiTXJ9mc5JlD5R7ns9B0jXe3zGN8FtoB\nYz62fWvH+JwMcsDvAyur6mjgH7r5X5BkHvB+YDmwFDgjybHd4iuA46rqeGAtcN5Q03VVdUI3vXE6\nd0KPybSMeZKlDF5GvbRr94Ekc/W4mm22d8yvA14GfG3Muj3OZ59pGW+P8Vlte8d8svbNHONz9S/j\nacCF3e8XAv92TJ0TGQzk+qraBHwCOB2gqlZW1Zau3pUMPhum2W26xvx04OKq2lRV64F13XrUv+0d\n8zVVtXZGeqodYbrG22N89tquMZ9i+1lvrga5g7rPfAHcBRw0ps4hwIah+Vu7slFnA5cNzR/RnYpd\nleR5O6S32hGma8wP5he/8TtRG828HTnmozzOZ5/pGm+P8dlre8d8svbNHOO9f9lhuiRZCSwes+j8\n4ZmqqgneKbfN97IkOR/YWFUXdUW3A0uq6p7uHovPJjmuqu5/jN3X49DTmI/jO31myEyM+Rge5z3p\nabzH8RifIdMw5hlTNtq+qWN8pw1yVfUbEy1LcleSxVV1Z5InAv86ptptwJKh+SUM/V9ZkjOBFcAp\nQ9vcCGzsfv9Oku8DRwHf2Y5d0RT1MeZj2hzalWkGTPeYT7BNj/Oe9DHeY9p4jM+gaRjz4fEb2761\nY3yuXlq9FHht9/trgc+OqXM1cFT39MoCBje7XgqDJ2CAtwGnV9XPtjZIckB3YyVJnsxg4G+atr3Q\nYzEtY94t/w9JFiQ5gsGYXzVN+6DHZrvGfMTP36bucT5rTct44zE+m23vmI9t39wxXlVzbgL2A/4P\ng6cPrwAWdeUHA18cqvcS4HsMbm49b6j8RuAHwDXd9IGu/BXAv3Rl3wZO7XtfnaZ3zLtlb+/qrwFe\n3Pe+Ou2wMX8Zg3trHgLuBC7vyj3OZ+E0XePdLfMYn4XTDhjzidq/vKVj3E90SZIkNWquXlqVJElq\nnkFOkiSpUQY5SZKkRhnkJEmSGmWQkyRJc0KSVya5Psnm7mW/E9X7aPeeuutGyt+VZHWSa5NckmRh\nV75fkq8muT/J+0baLEjy4STf69q+fBt9fHW3/v+X5BtJnj5ZfYOcJEmaK65j8KqZr22j3seA5WPK\nrwCOq6rjGby25Lyu/GfAHwBvHdPmfODOqjqmqo4F/nEb274JOLmqng78CfDhySob5CRpB0myPslX\ne9juliQfm+ntSq2pqjVVtXYK9f4JuGdM+cqq2tLNXsngSxFU1YNV9Q3g4TGrOwv4s6F13A2Q5N8k\n+XSSq7rpOd3yb1bVfaPbmIhBTtKslmRZF1Te0ndfpqAY+Y5jkjclee0E9Xf0tiXNnLOBy0bKRo//\nRd2v70jy7SSfTHJgV/aXwF9U1YnAvwP+esw2XjdmG79gp/3WqqSdTgtB5Wge3c83MbhUcuHMd0ea\ne5KsBBaPWfT2qvr8DtrG+cDGqrpoG1XnMzij9o2qekuS/wy8G3gN8ELg2OTnX4TbJ8meVfVgt41f\nZxAWn7utDUiSdoCq2tR3H6S5rqp+YzrXn+RMYAVwyhSq3w08WFWXdPOfZnCWDQbf9H12VW0cs42n\nA/8TWF5Vj7rEO8xLq5J2CklOTrIyyb1JHuwuY5w9pt6qJDcneWKSi5P8OMkDSb6U5Kgx9Q9P8r+T\n/CTJfUk+25U96n640bIkW4DDgK2Xh7dOh21dPu7etiRndstOHik/ruvnT5PcneTjQ5dpxv2ZvCrJ\n17u+P5Dkn5O8Yip/ntIckG1XGWmQLAfeBpxeVT/b1jpr8B3Uz3dn12AQ/q7vfr8C+L2hdT+j+3kY\ncAnwm1W1blt9MshJal6SlwJfAY5hcNniPGAT8NdJ3jFSvYC9GDy1tqmr+35gGfC5JD//dzHJ/sA/\nAacCHwX+C/AA8FVgTx59GXX0HrnfAn4ErAZ+c2j60UibqezjEV1fngu8D/ivwAHAlyao/w7gYuA+\nBk/TnQs8CHwqyRunsk1pZ5PkZUk2ACcBX0xyeVd+cJIvDtW7GPi/wNFJNiQ5q1v0PmBvYGWSa5J8\nYKjNeuA9wJlJbkny1G7RucAFSa4FXg1svd/394Bf6V41cj3w+q78vwH7Ah/stnHVpDtVVU5OTk6z\ndmIQsLYAb55g+TzgB8CPgcVD5bsCXwceAY4cKl/Vre+tI+t5a1f+oqGyd3ZlZ4zU/fOu/Csj5eun\nUja0bAvw0THlZ3bLTh4qu6gre/5I3UtG1wM8syt7x5h1f4ZBuNu777F1cnLa/skzcpJa98vAEgZB\n5s6thTW4X+2dDK48nD7SZjPw3pGyrZdEjxwqeylwe1VdPFL33dvb6ceiO0v4UuBbVTX6Dqp3jmny\nagZn+v5XkgOGJ+DzwD7Ar05rpyXNCB92kNS6I7qf149ZdsNIna1ur0ffYHx393P/kXX/8+hKq+qH\nSe4bLZ9GBzK4HLxmzLLVY8qOZXCvzrj6MAh5E95bJ6kdBjlJc9HmSZY95hugp8H2/tscBmFtORPv\n6w0TlEtqiEFOUuu+3/182phlS7ufNz3Oda8HjkqSqvr5Qwndk6ILp7iOyR5m+DGw35jyJ4/M/xD4\nKfDUMXWXjilbC7wY2FBVE52Vk7QT8B45Sa27BrgFOCvJQVsLk+zK4DUBW4DPPc51Xwo8EThjpHzc\n9xQn8lN+8XLtsLXAc5LssbUgyb4MPunz8wBYVZuBLwDPSrJsqG4YPEk76m+7n/99+CncoXYHjZZJ\napNn5CS14oVJ9hxT/kPgHAZPY34ryYcZhKdXAc8G/rSqvj/SZqqXT/8c+I/Ax5KcCHwP+DXgOQxe\nITKVV4d8E3hdkj9mcM/aFuDSGry9/f3Ax4GvJPk4sAj4bQZnAkfD1h8ALwG+kOR9wG0MHoA4YHSD\nVXV1kguAC4DvJvkUcAeDUPrL3Xp2m9KfgKRZzSAnabbbGpZezOCer1FrqmppklMYhJ23AQsY3AP2\nuqoafeHuo76HOuGGq+5O8jwG74Y6u2u3CngBcBXw0AR9HXY+g8unv8MgqMHgIYpbquqiJAczCKLv\nYXCZ+I+69Zw40pebkvxaV+93GXyc+zIG76W7a0zf/zjJ1QzeVfUmBg9L3AX8S9de0k4gQ7d9SJKm\noHtR8A+Bv6oqX64rqTfeIydJkxi+f23I73c/V85kXyRplGfkJGkS3bdT1zN4qGIXBt9KPBX4BoMv\nL/iPqKTeGOQkaRJJ3gy8Bjgc2APYwOCzWH9UVQ/02DVJMshJkiS1ynvkJEmSGmWQkyRJapRBTpIk\nqVEGOUmSpEYZ5CRJkhplkJMkSWrU/wcnX3KBFekCKQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "lqCluster=df[df.cluster ==1]\n", "plt.figure(figsize=(10, 10), dpi=100)\n", "df_scatter = plt.scatter(lqCluster['longitude'], lqCluster['latitude'], c='b', alpha=.5, s=lqCluster['LQ']*10)\n", "plt.title( 'LQ in cluster', fontsize=20)\n", "plt.xlabel('Longitude', fontsize=18)\n", "plt.ylabel('Latitude', fontsize =18)\n", "plt.xlim(lqCluster.longitude.min()-0.002,lqCluster.longitude.max()+0.002)\n", "plt.ylim(lqCluster.latitude.min()-0.002,lqCluster.latitude.max()+0.002)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Determining the enviroment\n", "\n", "\n", "Now is time to determine which businesses are carrying economic activity in each cluster.
\n", "I also wish to see how those categories relate to the most common category of the cluster. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I used multiple comparisons, so that needs to be corrected. I used a method called: \n", "\n", "##Sidak correction\n", "\n", "It is a simple method to control the familywise error rate that is probabilistically exact when the individual tests are independent of each other but is conservative otherwise. \n", "The test is used if the test statistics are independent of each other then testing each of m hypotheses at level $$\\alpha_{SID} = 1-(1-\\alpha)^\\frac{1}{m}$$ is Sidak's multiple testing procedure. \n", "This test is more powerful than Bonferroni, but the gain is small: for $\\alpha_{SID} = 0.05$ and $m= 10$ and $10^{12}$, Bonferroni vs Sidak give 0.005 and 5 $10^{-14}$ vs 0.005116 and $5.129 10^{-14}$, respectively. The main merit of the correction is that it is exact probabilistically when the tests are independent of each other. Bonferroni is an easier approximate way to calculate the Sidak correction. \n", "\n", "The Šidák correction is derived by assuming that the individual tests are independent. Let the significance threshold for each test be $\\alpha_1$; then the probability that at least one of the tests is significant under this threshold is (1 - the probability that none of them is significant). Since it is assumed that they are independent, the probability that all of them are not significant is the product of the probabilities that each of them are not significant, or $1 - (1 - \\alpha_1)^n$. Our intention is for this probability to equal \\alpha, the significance level for the entire series of tests. By solving for $\\alpha_1$, we obtain $\\alpha_1 = 1 - (1 - \\alpha)^{1/n}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Helper functions" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def getTopCategories(newcat,listCat,pomDF):\n", " busi= pomDF.name.ix[pomDF.categories == newcat].values\n", " if len(busi)>0:\n", " for b in busi:\n", " rcat = pom.categories.ix[pom.name == b].values\n", " if len(rcat)>0:\n", " for t in rcat:\n", " if t in listCat:\n", " pass\n", " else:\n", " listCat.append(t)\n", " return listCat" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def SortingBusinessCategories(categ,pomDF):\n", " busR = {}\n", " for rc in categ:\n", " rcat=rc.flatten().tolist()[0]\n", " busR.update({rcat: pomDF.name[pomDF.categories == rcat].count()})\n", " sortedBus = sorted(busR.items(), key=operator.itemgetter(1), reverse=True)\n", " return sortedBus" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def getCategoryStats(pomDF, categ):\n", " RcatList= arrayToList(categ)\n", " pom2 = pomDF.loc[pomDF['categories'].isin(RcatList)]\n", " Large = pom2.LQ.max()\n", " Bus = pom2.name.ix[pom2.LQ == Large].values[0]\n", " BusCat = pom2.categories.ix[(pom2.name == Bus) & (pom2.LQ == Large)].values[0]\n", " CatNum = pom2.name[pom2.categories == BusCat].count()\n", " answer=[]\n", " answer.append(Large)\n", " answer.append(BusCat)\n", " answer.append(CatNum)\n", " return answer" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def arrayToList(array):\n", " i=0\n", " newlist=[]\n", " for r in array:\n", " if type(r)==np.ndarray:\n", " st=r[0]\n", " newlist.append(st)\n", " return newlist" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Performing analysis for some other business in the cluster\n", "\n", "However, I cannot answer the question where to put stand-alone business, not with just Yelp data. To answer that question, I would need an additional data source.
\n", "The question I can answer is can some business play well together.
\n", "So, now, I'll concentrate on business that cluster around top business. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I'll pick experimental category and run with it.
\n", "You can see other examples I tried. I have to give up on those, mostly because there was no more than five different clusters where that business appeared. In that case, I would recommend not to use this analysis as determination, but different factors, different data. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "yeah\n" ] } ], "source": [ "#print df.categories.unique()\n", "targetCategory =\"Coffee & Tea\"#\"Bookstores\" # \"Taxis\" #\"Rugs\"#\n", "if targetCategory in df.categories.unique():\n", " print 'yeah'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, so there is something connected with Coffee.
\n", "Next step let us see how many clusters have this category." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "139\n" ] } ], "source": [ "pom=df[df['categories']== targetCategory]\n", "pom2=pom.cluster.unique().tolist()\n", "#print pom2\n", "clustersDF = df.loc[df['cluster'].isin(pom2)]\n", "targetCategoryClusters=clustersDF.cluster.unique()\n", "print len(targetCategoryClusters)\n", "#print clustersDF.cluster.unique()" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "Neat, over 100 clusters. \n", "Next step is to determine the top four businesses in each category and get LQ of targeted category for each cluster. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "df3=pd.DataFrame(df.cluster.unique())\n", "df3['BusNum']=0\n", "df3['topBusCat']=''\n", "df3['topCatNum']=0\n", "df3['LQmax']=0\n", "df3['cat1']=''\n", "df3['LQ2']=0\n", "df3['cat1num']=0\n", "df3['cat2']=''\n", "df3['LQ3']=0\n", "df3['cat2num']=0\n", "df3['cat3']=''\n", "df3['LQ4']=0\n", "df3['cat3num']=0\n", "df3['targetLQ']=0\n", "df3 = df3.rename(columns={0: 'cluster'})" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:22: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:29: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:36: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n" ] }, { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
clusterBusNumtopBusCattopCatNumLQmaxcat1LQ2cat1numcat2LQ3cat2numcat3LQ4cat3numtargetLQ
038637Beer Bar11.770652Religious Schools1.7706521Horse Racing1.7706521Hearing Aid Providers1.77065210.076894
10853Fireplace Services139.935685Fishing39.9356851Fireworks39.9356851Chimney Sweeps32.67465210.141500
2828000.0000000.00000000.00000000.00000000.000000
3919000.0000000.00000000.00000000.00000000.000000
418995Ski Resorts121.027324Bistros21.0273241Rugs21.0273241Cooking Classes21.02732410.147098
\n", "
" ], "text/plain": [ " cluster BusNum topBusCat topCatNum LQmax \\\n", "0 3 8637 Beer Bar 1 1.770652 \n", "1 0 853 Fireplace Services 1 39.935685 \n", "2 828 0 0 0.000000 \n", "3 919 0 0 0.000000 \n", "4 18 995 Ski Resorts 1 21.027324 \n", "\n", " cat1 LQ2 cat1num cat2 LQ3 cat2num \\\n", "0 Religious Schools 1.770652 1 Horse Racing 1.770652 1 \n", "1 Fishing 39.935685 1 Fireworks 39.935685 1 \n", "2 0.000000 0 0.000000 0 \n", "3 0.000000 0 0.000000 0 \n", "4 Bistros 21.027324 1 Rugs 21.027324 1 \n", "\n", " cat3 LQ4 cat3num targetLQ \n", "0 Hearing Aid Providers 1.770652 1 0.076894 \n", "1 Chimney Sweeps 32.674652 1 0.141500 \n", "2 0.000000 0 0.000000 \n", "3 0.000000 0 0.000000 \n", "4 Cooking Classes 21.027324 1 0.147098 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import operator\n", "for c in clustersDF.cluster.unique():\n", " pom = clustersDF[clustersDF.cluster == c]\n", " large = pom.LQ.max()\n", " topBus = pom.name.ix[pom.LQ == large].values[0]\n", " topBusCat = pom.categories.ix[(pom.name == topBus) & (pom.LQ == large)].values[0]\n", " df3.topBusCat.loc[df3.cluster == c] = topBusCat\n", " df3.LQmax.loc[df3.cluster == c]=large\n", " bNum=pom.name.unique()\n", " df3.BusNum.loc[df3.cluster == c]=len(bNum)\n", " topCatStart=[]\n", " topCat = getTopCategories(topBusCat, topCatStart, pom)\n", " topCatNum =pom.name[pom.categories == topBusCat].count()\n", " df3.topCatNum.loc[df3.cluster ==c]=topCatNum\n", " df3.targetLQ.loc[df3.cluster ==c]=pom.LQ[(pom.categories == targetCategory) & (pom.name != topBus)].max()\n", " cat = pom.categories.unique()\n", " Rcat = cat[np.argwhere(np.in1d(cat,np.intersect1d(cat,topCat))==False)]\n", " if len(Rcat)>0:\n", " ans=getCategoryStats(pom,Rcat)\n", " df3.cat1.loc[df3.cluster == c] = ans[1]\n", " df3.cat1num.loc[df3.cluster == c] = ans[2]\n", " df3.LQ2[df3.cluster == c] = ans[0]\n", " topCate = getTopCategories(ans[1], topCat, pom)\n", " Rcat1 = cat[np.argwhere(np.in1d(cat,np.intersect1d(cat,topCate))==False)]\n", " if len(Rcat1)>0:\n", " ans2=getCategoryStats(pom,Rcat1)\n", " df3.cat2.loc[df3.cluster == c] = ans2[1]\n", " df3.cat2num.loc[df3.cluster == c] = ans2[2]\n", " df3.LQ3[df3.cluster == c] = ans2[0]\n", " topCateg = getTopCategories(ans2[1], topCate, pom)\n", " Rcat2 = cat[np.argwhere(np.in1d(cat,np.intersect1d(cat,topCateg))==False)]\n", " if len(Rcat2)>0:\n", " ans3=getCategoryStats(pom,Rcat2)\n", " df3.cat3.loc[df3.cluster == c] = ans3[1]\n", " df3.cat3num.loc[df3.cluster == c] = ans3[2]\n", " df3.LQ4[df3.cluster == c] = ans3[0]\n", "df3.head(5) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us determine is our category in one of the 4 four here. If it is, then we can proceed with the analysis. " ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "yeah\n" ] } ], "source": [ "allCatIndf3=df3.topBusCat.unique().tolist()+df3.cat1.unique().tolist()+df3.cat2.unique().tolist()+df3.cat3.unique().tolist()\n", "uniqCatDF3= set(allCatIndf3)\n", "if targetCategory in df.categories.unique():\n", " print 'yeah'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The goal is to find business that preform well. When business is in one off the four most popular in the cluster, it si more probable that I will find a positive influence in that cluster. Taking all clusters with targeted business in consideration will drown the signal with the noise. " ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['Coffee & Tea', 'Gas & Service Stations', 'Gastropubs']\n" ] } ], "source": [ "import operator\n", "from scipy.stats import kendalltau\n", "par=[]\n", "for c in df3.topBusCat.unique():\n", " if c == targetCategory:\n", " par.append(df3.topBusCat[df3['topBusCat']==c].tolist()[0])\n", "for c in df3.cat1.unique():\n", " if c == targetCategory:\n", " par.append(df3.topBusCat[df3['cat1']==c].tolist()[0])\n", "for c in df3.cat2.unique():\n", " if c == targetCategory:\n", " par.append(df3.topBusCat[df3['cat2']==c].tolist()[0])\n", "for c in df3.cat3.unique():\n", " if c == targetCategory:\n", " par.append(df3.topBusCat[df3['cat3']==c].tolist()[0])\n", "print par" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now with these categories, we'll find clusters where there are target category and one of the listed categories and tests with Spearman correlation how they influence coffee shops. \n", "\n", "#Spearman's rank correlation coefficient\n", "\n", "It is often denoted by the Greek letter $\\rho$ (rho) or as $r_s$, is a nonparametric measure of statistical dependence between two variables. It assesses how well the relationship between two variables can be described using a monotonic function. If there are no repeated data values, a perfect Spearman correlation of +1 or −1 occurs when each of the variables is a perfect monotone function of the other.\n", "\n", "Spearman's coefficient, like any correlation calculation, is appropriate for both continuous and discrete variables, including ordinal variables. Spearman's $\\rho$ and Kendall's $\\tau$ can be formulated as special cases of a more a general correlation coefficient. \n", "\n", "The Spearman correlation coefficient is defined as the Pearson correlation coefficient between the ranked variables.For a sample of size n, the n raw scores $X_i$, $Y_i$ are converted to ranks $x_i$, $y_i$, and $\\rho$ is computed from:\n", "\n", " $$\\rho = {1- \\frac {6 \\sum d_i^2}{n(n^2 - 1)}}$$\n", "where $d_i = x_i - y_i$, is the difference between ranks.\n", "\n", "And for those who do not know:\n", "\n", "##Pearson product-moment correlation coefficient\n", "\n", "The Pearson product-moment correlation coefficient (sometimes referred to as the PPMCC or PCC or Pearson's r) is a measure of the linear correlation (dependence) between two variables X and Y, giving a value between +1 and −1 inclusive, where 1 is total positive correlation, 0 is no correlation, and −1 is total negative correlation. It is widely used in the sciences as a measure of the degree of linear dependence between two variables. \n", "Pearson's correlation coefficient is the covariance of the two variables divided by the product of their standard deviations. The form of the definition involves a \"product moment\", that is, the mean (the first moment about the origin) of the product of the mean-adjusted random variables; hence the modifier product-moment in the name. \n", "Pearson's correlation coefficient when applied to a population is commonly represented by the Greek letter $\\rho$ (rho) and may be referred to as the population correlation coefficient or the population Pearson correlation coefficient. The formula is:\n", "\n", " $$\\rho_{X,Y}= \\frac{\\operatorname{cov}(X,Y)}{\\sigma_X \\sigma_Y} $$\n", "where: \n", " $\\operatorname{cov}$ is the covariance \n", " $\\sigma_X$ is the standard deviation of X " ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "resul=pd.DataFrame(range(100))\n", "resul['category']=''\n", "resul['Spearman']=0.\n", "resul['P']=0.\n", "resul['sidak']=0." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:17: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:19: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Coffee & Tea -> 1.0 p= 0.0 Sidak correction: 0.0\n", "Gas & Service Stations -> 0.541412272587 p= 0.00010216046362 Sidak correction: 0.000306450081646\n", "Gastropubs -> 0.928571428571 p= 0.00251947240379 Sidak correction: 0.00753938998076\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/Lexa/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:20: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame\n", "\n", "See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n" ] } ], "source": [ "import operator\n", "from scipy.stats import spearmanr\n", "i=0\n", "for ca in par:\n", " tarClustDF=clustersDF[clustersDF['categories']==ca]\n", " pomList=tarClustDF.cluster.unique().tolist()\n", " clustersDF2 = clustersDF.loc[clustersDF['cluster'].isin(pomList)]\n", " if len(tarClustDF)>5:\n", " para=[]\n", " resu=[]\n", " for c in pomList:\n", " para.append(clustersDF2.LQ[(clustersDF2['categories']==ca) & (clustersDF2['cluster']==c) ].max())\n", " resu.append(clustersDF2.LQ[(clustersDF2['categories']==targetCategory) & (clustersDF2['cluster']==c)].max())\n", " cou=len(par)\n", " spear=spearmanr(para,resu)\n", " sida=1.-(1.-spear[1])**cou\n", " resul.category[i]=ca\n", " resul.Spearman[i]=spear[0]\n", " resul.P[i]=spear[1]\n", " resul.sidak[i]=sida\n", " i+=1\n", " print ca, '->', spear[0],'p=',spear[1],'Sidak correction:', sida " ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "result=resul[resul['sidak'] < 0.1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Validation\n", "\n", "Using the same random pick of the categories" ] }, { "cell_type": "code", "execution_count": 106, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0857142857143 0.761334126191\n" ] } ], "source": [ "import random\n", "category = random.sample(df.LQ,15)\n", "LQrandom = random.sample(df.LQ, 15)\n", "spear = spearmanr(category, LQrandom)\n", "print spear[0], spear[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Shows there is no correlation between randomly chosen categories. Correlation is close to 0, and the P value is way too high. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##let's make a map\n", "\n", "First we need to sort out which Spearman coefficient is significant, larger than 0.5 and then sort which one is negative and which one is positive." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ "negCat=result.category[result['Spearman']<-0.5].tolist()\n", "posCat=result.category[result['Spearman']>0.5].tolist()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we have to get clusters that have a positive influence business and do not have negative influence business.
\n", "In case that one of the categories is empty, it is skipped. The rest of the clusters are treated as neutral. Meaning, there, a personal performance of a business, together with unaccounted for factors will have the greatest influence. " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "allClust=df.cluster.unique()\n", "goodList=[]\n", "neutralList=[]\n", "if len(posCat)>0:\n", " goodf= df.loc[df['categories'].isin(posCat)]\n", " good=goodf.cluster.unique()#.tolist()\n", " withGood =allClust[np.argwhere(np.in1d(allClust,np.intersect1d(allClust,good))==True)]\n", " goodList=arrayToList(withGood)\n", "\n", "if len(negCat)>0:\n", " badf= df.loc[df['categories'].isin(negCat)]\n", " bad=badf.cluster.unique()#.tolist()\n", " withoutBad = allClust[np.argwhere(np.in1d(allClust,np.intersect1d(allClust,bad))==False)]\n", " neutralList=arrayToList(withoutBad)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we get coordinates for the good and neutral clusters." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if (len(goodList)==0) & (len(neutralList)==0):\n", " print 'No good places for business. Try put it as stand alone business.'\n", "if len(goodList)>0:\n", " coordGood = pd.DataFrame(goodList)\n", " coordGood = coordGood.rename(columns={0: 'cluster'})\n", " coordGood['lati']=0\n", " coordGood['longi']=0\n", " coordGood['ratioNbus']=0\n", " coordGood['LQmean']=0\n", " coordGood['scale']=0\n", "if len(neutralList)>1:\n", " coordNeutral = pd.DataFrame(neutralList)\n", " coordNeutral = coordNeutral.rename(columns={0: 'cluster'})\n", " coordNeutral['lati']=0\n", " coordNeutral['longi']=0\n", " coordNeutral['ratioNbus']=0\n", " coordNeutral['LQmean']=0\n", " coordNeutral['scale']=0" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import division\n", "if len(goodList)>0:\n", " for c in goodList:\n", " pom4= df[df.cluster == c]\n", " NumBus=len(pom4.name.unique())\n", " NumGoodBus = pom4.name[pom4.categories == cat[0]].count()\n", " ratio=NumGoodBus/NumBus\n", " meanLQ = pom4.LQ[pom4.categories==cat[0]].mean()\n", " lati =pom4.latitude.mean()\n", " longi =pom4.longitude.mean()\n", " scal = ratio*meanLQ\n", " coordGood.lati.loc[coordGood.cluster == c] = lati\n", " coordGood.longi.loc[coordGood.cluster == c] = longi\n", " coordGood.ratioNbus.loc[coordGood.cluster == c] = ratio\n", " coordGood.LQmean.loc[coordGood.cluster == c] = meanLQ\n", " coordGood.scale.loc[coordGood.cluster == c] = scal\n", "if len(neutralList)>1:\n", " for c in neutralList:\n", " pom4= df[df.cluster == c]\n", " NumBus=len(pom4.name.unique())\n", " NumGoodBus = pom4.name[pom4.categories == cat[0]].count()\n", " ratio=NumGoodBus/NumBus\n", " meanLQ = pom4.LQ[pom4.categories==cat[0]].mean()\n", " lati =pom4.latitude.mean()\n", " longi =pom4.longitude.mean()\n", " scal = ratio*meanLQ\n", " coordNeutral.lati.loc[coordNeutral.cluster == c] = lati\n", " coordNeutral.longi.loc[coordNeutral.cluster == c] = longi\n", " coordNeutral.ratioNbus.loc[coordNeutral.cluster == c] = ratio\n", " coordNeutral.LQmean.loc[coordNeutral.cluster == c] = meanLQ\n", " coordNeutral.scale.loc[coordNeutral.cluster == c] = scal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot the result. Let us first check the sanity of our coordinates. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3, 0, 18, 25, 9, 236, 248, 21, 48, 4, 49, 80, 541, 39, 52, 167, 55, 100, 110, 10, 47, 5, 254, 75, 35, 16, 19, 507, 189, 78, 269, 6, 28, 30, 23, 390, 206, 322, 130, 99, 298, 125, 59, 94, 149, 232, 200, 74, 304, 77, 491, 427, 136, 27, 355, 51, 283, 7, 98, 161, 82, 58, 152, 92, 267, 11, 44, 24, 119, 228, 29, 205, 181, 135, 243, 33, 286, 45, 112, 61, 565, 40, 102, 230, 2, 291, 151, 109, 188, 302, 290, 20, 62, 178, 116, 281, 498, 14, 159, 142, 321, 199, 327, 72, 106, 81, 22, 218, 68, 50, 1, 115, 12, 13, 144, 154, 477, 192, 104, 103, 174, 15, 65, 414, 191, 141, 96, 584, 655, 225, 122, 209, 26, 426, 86, 128, 198, 255, 702, 163, 168, 471, 306, 31, 37, 8, 264, 299, 220, 244, 107, 706, 242, 145, 34, 180, 182, 138, 56, 36, 233, 190, 312, 308, 195, 217, 87, 261, 359, 193, 921, 554, 644, 902, 621, 666, 97, 449] []\n" ] } ], "source": [ "print goodList, neutralList" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def intersect(a,b):\n", " return list(set(a) & set(b))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAADSCAYAAADaDs5CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHiJJREFUeJzt3Xn4HFWd7/H3B5KwRUBENiEGBYdlEPRh0Cszkos6AhdZ\nBkV90GFTmXsddQSVRdSoKKBXRVScqwZkURwUBmFQBJUojrghkQBGZCQQlkTZ12FJvvePc5pUOr13\nVf/696vP63l+T7q7qs73VNU59a06Vd1RRGBmZlYna0x0BczMzEbNyc/MzGrHyc/MzGrHyc/MzGrH\nyc/MzGrHyc/MzGpnwpOfpC9LOrGksmZJeliS8vv5ko4so+xc3vckvbWs8vqIe5Kkv0i6a9SxO5G0\nQtILRhDnBkmvLLG8xZJeVVZ548p9q6e4Y9G3xq1Nlr0/ym4vZZhWZeGSFgObAE8Dy4GbgHOAr0T+\ngmFE/O8+yjoiIn7cbp6IuB14VvGj/DdI3ecCL4yIZxpAROwzSFnDkDQLOBrYKiLuHXX8QUmaD5wb\nEfP6XO7rwJKI+FDjs4j468L0uTTtlwEM3C7GhfvW8Masb7Xdnq36ROWVKewPSYcBR0bE3w1TJGPW\n56q+8gtg34hYH5gFnAIcC/R1QCyUpXYTJVWayCfQLODeMeic/Rq7xj7FuG8Nb7L0LfelKkREZX/A\nrcCeTZ/9DelMdYf8/uvAx/PrjYH/AO4H7gV+SuqU5+ZlHgMeBt4HzAZWAEcAtwHzgefnz9bI5V0F\nfBL4JfAgcDHw7DxtDulsqli3xcCrgL2AJ4Anc7zr8vT5pDMgcr1OzMssA84G1s/TGnX7x1y3vwAn\ndNhOG5DO2v+cy/tgLv/VeZ2X53qc2cM2XwP4TI75J+Cfm7bJFsAlefv+EXhbYdm1gNOAO/Pf54AZ\nhenvB+4C7sjbfQXwgjb1uIp0NdFq2reBu4EHgJ8U2sI78jZ/Iq/vd3vcL4uBVxXKn0u66my8f2ve\nD/cAJ1Bol3k7Hwfckqf/W6GNrA2clz+/H/gVsEmVfcZ9a7z6FvABVrb5t1Fo8+3K7rYO3dpkizqc\n1diPLaa9ndSP7wW+C2xemLYCOAq4Oe/3L/ZxnJgPHAlsB/w3aYThYeC+5v2V3x8GXF14/xpgEamP\nf6HF/EeQRivuAy4HZhWmfS5vsweB64Edq+hDI7/nFxG/JjWkxiV08azmGGAJqaNuAhwfyVuB20ln\nus+KiP9bKPKVpB30WlY/exWpkxwObE7agad3ql6qYlxO6tjfyvFe0qKuhwOHkjr6C4CZwBebytsd\neBGp039Y0nZt4n6BNKS0NbBHo84R8UNgb+CuXI8jOtS94R2kA8zOwEuBA1j1rPFbpG25OfB64JOS\n/mee9kFgt7zszvn1iQCS9iLtn1fndXp1D3Vp5zJgG+C5wG+BbwBExFfy61Pz+u6f5+9nv1B8LWkH\n4AzgEFLifw6wZWHedwP7kdrR5qSDxJfytEOB9fP8G5EOJI8Psd6Vct9qaeC+ldv8e3OMbXN9upbd\nbR16aJM9kbQnaVu+gbQPbiP176L/BewKvBg4WNJr8+fdjhON/bWI1O6vydtpo+L0NvXaGLiQlNSf\nA/wXaX9Fnr4/cDxwIKk9Xg2cn6e9ltR+t42IDfK6VXJlPlEPvNxFOpg0e5K0E2dHxPKI+M8eypob\nEY9HxBMtpgVwTkTcFBGPAR8iNYC2QzwFosNQEKnhfiYiFkfEo6Sd+SZJxW360Yh4IiKuB35Hamir\nBpHWBN5IOhg9GhG3kc7IGvdDeqlr0cHAaRFxV0Q8AJzcKEPSVsArgGMj4smI+B3wNVKnbazTxyLi\nnoi4B/hooR4Hk86OG9vyI33W6xkR8fW8rk/lGDtLKt5P6rTO3fZL8/KvBy6NiJ9FxJOkNrCiMP0o\n4MS8vRr1eX3eL0+SOu+2OVFcFxEP97SSE8d9qxFk+L7VaPO/j4jHKbT5Hsputw5r0r1N9uoQYF5E\nLMjlHA/8j3wvs+GUiHgoIpaQrtYb26ntcaKFfo9B+wA3RMRFua2dBiwtTP8n4OSI+ENErMixd8n1\nfpJ0QrG9pDXyPEtXi1CCiUp+W5IudxsaG/fTpOGnKyT9l6RjeyhrSR/Tbwemk842htU40yqWPQ3Y\ntPBZcac9BqzXopyNc52ay3reEPUqrvMdhddbkIYtHm2KtUVh2eZ6FKc1b8u+SVpD0imSbpH0IGm4\nB8rZJ61sQWEb5AN18UxyNvDvku6XdD9pKOZp0tXRucAPgG9JulPSqZPg/pf71krD9q1Ofalb2Z3W\nYXM6t8lerRIj9+t7WXX9mrfTzMKy7dZtWKv0uawY6/nA5wt9rrHuW0TEVaQr5C8ByyT9v6YT49KM\nPPlJ+hvSxvlZ87SIeCQi3hcRLyQNRR1dGJJrd8O3243gWU2vnyKNsz8KrFuo15qkYbhey72LdOAs\nlv00aay6H/fkOjWXNWhjvBvYqvC++PouYCNJMwufzSLd32tMb65HY9rdrL4tB3EIad++Kg9rbJ0/\nbxyku233VtMfZdWD32aF+e6isA0krUu6mmu4HdgrIp5d+Fs3Iu6OiKcj4mMRsSPpinlfVl4ljx33\nrdUM27c69aVuZbdbh6XN5bZok6202marxJC0Xi7nzhbzNuu0br3EbtXnivUqrp+ayr8deEdTn1sv\nIn4BEBFfiIhdgR1IQ9vv72F9+jaK5NcYcltf0r6ksd1zI+LG4vQ8z76Stskb6yHSzejGcMAy4IUD\nxH6LpO1zA/sY8O2ICNJN4LUl7SNpOune1lqFZZcCszsM45wPvFfS7JxMGvcxOg1frFZWRCwHLgA+\nIWmmpOeT7jOc19+qPuMC4D2StpC0IekJwMaj70uAnwMnS1pL0otJN54bsc4HTpS0cR63/3Bh2gXA\nYYVt2cuw53RJaxf+ppPOPJ8A7sud9ZNNyywj3SNpp9V+WUAaUpomaVfgoMK0C4F9Je0uaQapDRTb\n/b+S7nvOApD0XEn75ddzJO2UD94Pkw52y3tY71Fx31q1PqsooW9dABwuabu8jsWv33Qru9M6dGuT\nrdZtWlNfmpFjHC5pZ0lr5Ri/iPS1lHblNLZT2+NEC8uALfO+bFgA/IOkdSRtQ3o4puF7wI6SDswj\nJe9m1eT4r8AJSvc+kbSBpDfk17tKelmO9RjpYZtK+twokt+lkh4iZfvjSePihxemF2+cbgNcSTrQ\n/Bz4UkT8JE87mXRgvl/S0YVlmzXftD2H9NTb3cAM0o4gIh4E/g/pntcdwCOsemn+7fzvvZJ+0yLO\nmaRhsZ+SnpZ6DHhXm3p0+oy83KO5nKtJD32c1cNyrXwVuIL0lNS1pIdLlhcOHG8mnS3eBVwEfDhW\nfr/rJOA3ednr8+uTAPKDCqcBPyYd3H7UQ72+TNoujb95pP1xG+ns9AbgmqZy5gE75P18UYsyW+2X\nD5EO3veTnvT8RmPmnAjeCXwzr/N9rLqfP096+vWK3E6vIT3oA6nDfpv01NlNpCfWzu2yzqPkvtX5\nMxiib+U2fzrpXtnNpLYB6eStW9lt16GHNtlq3Y5j1b70w4j4EantX5jL2Rp4U4d1K7aHbseJoh8B\nNwJLJf05f/Y50v25ZXmdz2uUnZ8XeAPp6zf3kNreM6MREXExcCrpdsKDwELSQ1WQHjD7St4mi/Py\nn+6wbQbWeCy39cT0gMQ5pPsfQfoC7elKX1J9G+kxWUg3fS+vooI2HEl7A1+OiNkTXRdrz31t/Ena\nnnSgntHlKnTSqeNxolvy2wzYLCIW5Ev3a0mPxB4MPBwRnx1NNa1XktYG9iSd1W1KOiv8eUQc3XFB\nm1Dua+NJ0oGkYbx1Sd/Vezoi/mFiazU8Hye6DHtGxNKIWJBfPwL8npVPEvX7+KuNhkhDf/eRvkN3\nI+nenY0x97Wx9Q7S0N4tpHu+Pf1k3CRQ++NExyu/VWaUZpN+jWNH0hdmDyfdC/kNcEz+roiZDcl9\nzax6PSW/PAwzHzgpIi6WtAkr70F8nPSTOkc2LdPPQxpmk0ZEVHYl5r5mtlKVfa2X3xCcTvqi77+0\nmT4bWNji8+hWdhl/pF+hmPJxgIBo89f/tm7EKbvccdtuFcQZept0KNt9zXEcZ2WcqLL8jvf88vdw\n5gE3RfqJmsbnmxdmO5D0BJSZDch9zWy0uv1U0+7AW4DrJV2XPzsBeLOkXUiPZN9K+n1EMxuc+5rZ\nCHVMfhHxM1o/Efr9aqozkPmO4zgjjFMJ9zXHcZzR6vlpz74LliKqvFlZM+mhhnb7Sgy6rasqd6oa\nx3Y9jnUyG1bV7Xqi/lcHMzOzCTPu/z2LDcmPwZuZrc7JrxY65T+PlplZ/XjY08zMasfJz8zMasfJ\nz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zM\nasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMasfJz8zMamfaRFegTiRFp+kRoVHVxczGT7dj\nBPg4URYnv5Fr17bdns0M2h8jwMeJ8njY08zMaqdj8pO0laSrJN0o6QZJ786fbyTpSkk3S7pC0oaj\nqa7Z1OS+ZjZaimh/iS1pM2CziFggaSZwLXAAcDhwT0R8StKxwLMj4rimZcNj06tK4/nthz07ba9B\nl+28XFp20DrVUVXt2n3NoLf+Wpd9XXW77njlFxFLI2JBfv0I8HvgecB+wNl5trNJndTMBuS+ZjZa\nPd/zkzQbeAnwS2DTiFiWJy0DNi29ZmY15b5mVr2envbMwzAXAu+JiIellVeiERHtHs+VNLfwdn5E\nzB+8qmajJ2kOMGeE8dzXRqyXrxeAv2JQtZH3tU73/AAkTQf+A/h+RJyWP1sEzImIpZI2B66KiO2a\nlvN9iCa+5zf5Vdmu3dcmRvc+AqPqC77nt9KE3vNTOu2cB9zU6IzZJcCh+fWhwMXVVM+sHtzXzEar\n29Oefwv8FLielacjxwO/Ai4AZgGLgYMj4oGmZX022sRXfpNfhU97uq9NEF/5jaeq23XXYc+BC3aH\nXI2T3+Q3ju16HOs0mTj5jacJHfY0MzObipz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8\nzMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdqZNdAVs\npfR/eZnZVOT+PV6c/MZKt/901swmt+7/aa6Nhoc9zcysdpz8zMysdpz8zMysdpz8zMysdpz8zMys\ndpz8zMysdpz8zMysdromP0lnSlomaWHhs7mS7pB0Xf7bq9pqmk1t7mdmo9XLld9ZQHOnC+CzEfGS\n/Hd5+VUzqxX3M7MR6pr8IuJq4P4Wk/xTBGYlcT8zG61h7vm9S9LvJM2TtGFpNTKzIvczswoM+tue\nXwY+ll9/HPgMcGTzTJLmFt7Oj4j5A8abNKbaj9d2Wp+ImPJXJZLmAHMmKHxP/Qzq2dd6MdX641Q2\n6r6miO5tQ9Js4NKI2KnXaZKiDgfHZqmztdumYrBp3Zdtt60712eYOrWPOZVV2a4H6WdV12my697+\noXvfS/MMu43LqUt9+l3V7XqgYU9JmxfeHggsbDevmQ3G/cysOl2HPSWdD+wBbCxpCfARYI6kXUin\nKLcCR1VaS7Mpzv3MbLR6GvYcqOCaDsV42HNqG8d2PY51Ghce9py8xnLY08zMbDJz8jMzs9oZ9KsO\nZmYTbpRfZegWa1yGI3vdJuNS34ni5Gdmk1y3e9qTLU4Zerm3WG8e9jQzs9px8jMzs9px8jMzs9px\n8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMz\ns9px8jMzs9px8jMzs9px8jMzs9px8jMzs9px8jMzs9qZNtEVsHJIiomug1mdjaoPuq+Xw8lvymjX\nHzTSWpjVV7ecVFZfHFWcqc3DnmZmVjtdk5+kMyUtk7Sw8NlGkq6UdLOkKyRtWG01zaY29zOz0erl\nyu8sYK+mz44DroyIFwE/yu/NbHDuZ2Yj1DX5RcTVwP1NH+8HnJ1fnw0cUHK9zGrF/cxstAa957dp\nRCzLr5cBm5ZUHzNbyf3MrCJDP+0ZEdHu0VtJcwtv50fE/GHjmY2SpDnAnAmuRsd+Bu5rNvmNuq8p\novtXRiTNBi6NiJ3y+0XAnIhYKmlz4KqI2K5pmYiI2j1zmw5Qnb52MMi0YZatrty67t+q1nuQflZ1\nncZd5/4G3dt/WfNMvjjj3maqbteDDnteAhyaXx8KXFxOdcyswP3MrCJdr/wknQ/sAWxMuu/wYeC7\nwAXALGAxcHBEPNC0XC3PRn3lN7VV1a4H7WdV1mky8JXf4POMe5upul33NOw5UME17ZBOflPbOLbr\ncazTqDj5DT7PuLeZcR32NDMzm7Sc/MzMrHac/MzMrHac/MzMrHac/MzMrHac/MzMrHac/MzMrHac\n/MzMrHac/MzMrHac/MzMrHac/MzMrHac/MzMrHac/MzMrHac/MzMrHac/MzMrHamTXQFJqP0f4iZ\nWZXcz6rVbfuO+//3Nywnv4F1+g9gzawcvfzHrTaYbv9p7tTmYU8zM6sdJz8zM6sdJz8zM6sdJz8z\nM6sdJz8zM6sdJz8zM6sdJz8zM6udob7nJ2kx8BCwHHgqInYro1JmtpL7mVn5hv2SewBzIuK+Mipj\nZi25n5mVrIxhz6n/UwBmE8/9zKxEwya/AH4o6TeS3l5GhcxsNe5nZiUbdthz94i4W9JzgSslLYqI\nqxsTJc0tzDs/IuYPGc8mianyo7mS5gBzJrgaHfsZTL6+5h+tHn+97KMy+/Go+5oiymmDkj4CPBIR\nn8nvY7Ic4PqVGkWnH7Yue9p4lttp/3bbRpO1bUx0u27uZ+NQp0F0bh/PzMXw85RRhuO0m6fKdld1\nux542FPSupKelV+vB/w9sLCsipmZ+5lZVYYZ9twU+HdJjXK+ERFXlFIrM2twPzOrQGnDnqsVPAmH\nYnrlYc80zcOe42Ec69SNhz2nRpxaDnuamZlNVv6f3M2mAEkHAH/dw6wnR8TyqutjNu487DkAD3um\naR72HA9pW290Cez+Otipw5wnC2KtiHhyNHXqZnyG7xxn0DK6G7S/VN3XfOVnNjUI3ig4pMMsp65I\nPw86Kt0Orjb5lZMgJ4Lv+ZmZWe04+ZmZWe04+ZmZWe04+ZmZWe04+ZmZWe34aU8bmH+Z38wmKyc/\nG4IfZTezycnDnmZmVjtOfmZmVjtOfmZmVjsTfs9P0kbAjA6zPBYRD42qPmZmNvVNePKDmZfCil1h\nWosfHXxiGvBV4J2jrpVVq9OTou1+zLaXp0vH7UenpyI/5Wv96NZeJqrPjsGw51oz4Lsz4MF1Vv87\ndTpMX3Oia2hViDZ/gy7n4/FoddoP3hdWNJ7tZAySn5mZ2Wg5+ZmZWe04+ZmZWe04+ZmZWe04+ZmZ\nWe2MwVcdunnkKElHVVHyMI/Um01FbvtWF5Mg+UH7R2LVYVq36d2+WuIfbba66pb/3P5t8vOwp5mZ\n1c7AyU/SXpIWSfqjpGPLrFR/5juO44wwzui5rzlOPeNUa6DkJ2lN4IvAXsAOwJslbV9mxXo333Ec\nZ4RxRst9zXHqG6dag1757QbcEhGLI+Ip4FvA/uVVy8wy9zWzCgz6wMvzgCWF93cALxusqOXL4ZhH\n4TlPrz7t9rWAtQcr12xK6LGvPb0cTnoc5j3ZvqjlG5RcN7NJSxH9P9ks6SBgr4h4e37/FuBlEfGu\nwjx+ZNqmpFH+Cr37mtVZlX1t0Cu/O4GtCu+3Ip2RPsP/tYxZKdzXzCow6D2/3wDbSpotaQbwRuCS\n8qplZpn7mlkFBrryi4inJf0z8ANgTWBeRPy+1JqZmfuaWUUGuudnZmY2qUVEX3+kpy9/CSwAbgJO\nzp9/Gvg98DvgImCDNssfD9wILAS+CazVZ5yP5xgLgB8BW7VZfi9gEfBH4NgB1qdrHNL9l6vy+twA\nvLuKOIUy1gSuAy6tKg6wIfCdvC9vAl5eUZyh2kFh+jHACmCjKtpBL3H6aQd99LM35PKWAy8tfL5R\njvUw8IXC5+sAl+X9dkNz/ZvW83zg+rye36giTp73xcA1eb7bcrzS4+T5ZwGPAGdXtN1eQxqCvj7/\nO7fC7XZ8bq+LSP2p5zh52ieA24GHu7T3gdtBr3GGbQf9xGlqB8d0nXfAjrlu/nca8Avgb3PjWCN/\nfgpwSovlZgN/Ih/ogH8DDu0zzrMK098FfK3FcmsCt+R400kHtO0riLMZsEt+PRP4QxVxCtOPzo30\nkgH2T09xSAePIwrLtzyJGXK7Dd0O8vutgMuBW2mdlIZuBz3G6asd9NjHtgNeRDogFA8S6wK7A0ex\n+sF1j/x6OvBT0lOizeUeBpxfWOYOYI8K4kwjnQTtlN/vBvxV2XEK838nt6NTK9puuwCb5dc7Aksr\nirNDbqfTc7u9rZ/tVtjWm9E5+Q3VDvqIM1Q76DVOi3bQNfkN9MBLRDyWX84gHWDui4grI2JF/vyX\nwJYtFn0IeApYV9K0vMJ39hnn4cIsM4F7Wiza1xeDB40TEUsjYkF+/QjprG6LCtYHSVsC+wBfo8sv\nCw8aR9IGwN9FxJm5nKcj4sEK1mfodpDffxb4QLvlKKEd9BKn33bQi4hYFBE3t6pjRPwn8ETT549H\nxE/y66eA35K+I9jsbmC9/Msx6wGPkg5OZcf5e+D6iFiY5/1VRPyhgjhIOoB0MnUT8OcqtltELIiI\npfntTaTkdGsF67M/KSk9FRGLSW3p2b3GydN+VahrO0O1gz7iDNUO+ojT3A66GvTnzdaQtABYBlwV\nEc3BjgC+17xcRNwHfIZ0CXsX8EBE/LDfOJI+Iel24FDSVWazVl8MbtlxhoxTLGM28BJS4q8izueA\n95OG3joaIs7WwF8knSXpt5K+KmndsuOU0Q4k7Q/cERHXd9gUQ7eDHuMUy5hNl3ZQkuhQhw2B15GG\nnVddKOIHpJOPu4HFwKcj4oGy4wDbAiHpcknXSnp/hxgDx5E0k3RiMrdL+UPFaXIQcG1OYmXH2YJV\nv8rSsc12itNJme2gi9LaQScDtIOBr/xWRMQupKu7V0qaU6jEB4EnI+KbLSr4QuBfSJfzWwAzJR3S\nb5yI+GBEzAK+TkoKqy1axvr0EKexXjNJl9vvyWf+pcaRtC/pjPY6evj/ZIZYn2nAS4EzIuKlpLPB\n4ypYn2HbwT6k+yIfKRbbatF2ZZYcJ03osR0U5r9S0sIWf6/rp95NZU4j3cv5fL5yaI6zhHTgXUYa\nYnqfpK3LjkNqN28i3YP5JHCgpD0riHMrKTn8kjScNpAetlvj72bgdNK2Kz0O8Gbg5KZ2MGjiaY5d\nejvoFoeS2kEP5gKfy6M3PX3vdaj/zy8iHpR0GbArMF/SYaShuVe1WWRX4OcRcS+ApIuAV5DuY/Uc\npzDpm7S4wqSHLwaXFAdJ04ELgfMi4uJuMQaM8wpgv3wwXhtYX9I5EfGPJce5g3SV8+v8/jt0SH5D\nxBm2HbyUdJX6O0mQktW1knaLiD8XFhm2HfQaZ9B28Jpe5uvTV4A/RMTpreJIOoO07c/L7/ck7Y+y\n47wR2DsiDsvvtydtz7Lj/JS0X2cCewOvlfR42XFyrC1JV2sHRcRqQ55lxJF0XP7slPz+nXS4JdCP\nitpBtzhltYNudgMOkvQp0kN7KyQ9HhFntFug7ys/SRvny3YkrUN60OU6SXuRhuX2j4j/brP4IuDl\nktZROpq8mjbjsx3ibFOYbX/S04/Nev5i8DBx8jrMA26KiNParPPQcSLihIjYKiK2Jp1F/bhd4hsy\nzlJgiaQX5Y9eTXoCrNQ4DN8OromITSNi67xN7iDdPP9z0+LDtoOe4vTTDgbU6kx2tc8knQSsD7y3\nQ1mLgD3z/OsBLyfdVyo7zg+AnfI+nkZ6mKLRlkqLExGvLOyf04BPFA54pcXJbeMy0hPD13Qrc9A4\npPb5Jkkz8pXYtsCv+onTo6HbQY+Gbge96NIO2i7U1x+wE+lm7QLSY7Lvz5//kfRk0nX574z8+RbA\nZYXlP8DKR9zPBqb3Gec7edkFpDPtTdrE2Zv01N0twPEDrE/XOKSnG1fkeRrr3fKJtGHXp1DOHnR4\n2rOE7bYz8Gu6f2Vl2DhDtYOmef5Efgqz7HbQS5x+2kEf/exA0v3Kx0lPFn6/MG0xcC/psfAlpCdD\nt8x1uLFQh8ZTu68DPppfrwWcl7f7jcBZVcTJ7w8hPd6+MLeFSuIUyvkIcGZF2+1E0iP0jXn+RLoi\nq2K7nUBqr4uAj/azPvnzT+X3T+d/P1x2O+g1zrDtoJ84Te3g6G59zF9yNzOz2hn4f3I3MzObrJz8\nzMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdpz8zMysdv4/bRYG3bB3CGAAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "figsize(15, 3)\n", "if len(goodList)>0:\n", " tornTupleG = coordGood.cluster.tolist()\n", " latsG = coordGood.lati.tolist()\n", " lonsG = coordGood.longi.tolist()\n", " #scalesG = coordGood.scale.tolist()\n", "\n", " subplot(141)\n", " title(\"Distribution of good Latitudes\");\n", " hist(latsG, bins=20);\n", "\n", " subplot(142)\n", " title(\"Distribution of good Longitudes\");\n", " hist(lonsG, bins=20);\n", "\n", " \n", "if len(neutralList)>0:\n", " tornTupleN = coordNeutral.cluster.tolist()\n", " latsN = coordNeutral.lati.tolist()\n", " lonsN = coordNeutral.longi.tolist()\n", " #scalesN = coordNeutral.scale.tolist()\n", "\n", " subplot(141)\n", " title(\"Distribution of neutral Latitudes\");\n", " hist(latsN, bins=20);\n", "\n", " subplot(142)\n", " title(\"Distribution of neutral Longitudes\");\n", " hist(lonsN, bins=20);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Let us plot the result using Folium. " ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from IPython.display import HTML\n", "import folium\n", "def inline_map(map):\n", " \"\"\"\n", " Embeds the HTML source of the map directly into the IPython notebook.\n", " \n", " This method will not work if the map depends on any files (json data). Also this uses\n", " the HTML5 srcdoc attribute, which may not be supported in all browsers.\n", " \"\"\"\n", " map._build_map()\n", " return HTML(''.format(srcdoc=map.HTML.replace('\"', '"')))\n" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "meanlat=df.latitude.mean()\n", "meanlong=df.longitude.mean()\n", "\n", "map = folium.Map(width=600,height=600,location=[meanlat,meanlong], zoom_start=10)\n", "\n", "if len(goodList)>0:\n", " for i in range(len(tornTupleG)):\n", " map.simple_marker([latsG[i], lonsG[i]], popup=str(tornTupleG[i])+' Reccomended',marker_color='green',marker_icon='ok-sign')\n", "\n", "if len(neutralList)>0:\n", " for i in range(len(tornTupleG)):\n", " map.simple_marker([latsN[i], lonsN[i]], popup=str(tornTupleN[i])+' Neutral',marker_color='blue',marker_icon='ok-sign') \n", " \n", "inline_map(map)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "92 Skate Parks\n", "['Local Services' 'Dry Cleaning & Laundry' 'Optometrists'\n", " 'Health & Medical' 'Beauty & Spas' 'Nail Salons' 'Home Services'\n", " 'Real Estate' 'Apartments' 'Food' 'Grocery' 'Shipping Centers'\n", " 'Printing Services' 'Notaries' 'Chiropractors' 'Fast Food' 'Sandwiches'\n", " 'Restaurants' 'Sushi Bars' 'Japanese' 'Donuts' 'Web Design' 'Marketing'\n", " 'Graphic Design' 'Professional Services' 'Nightlife' 'Irish' 'Bars'\n", " 'Sports Bars' 'Libraries' 'Public Services & Government' 'Landscaping'\n", " 'Veterinarians' 'Pets' 'Pizza' 'Automotive' 'Gas & Service Stations'\n", " 'Plumbing' 'Active Life' 'Skate Parks' 'Swimming Pools' 'Parks'\n", " 'Fitness & Instruction' 'Fashion' 'Shopping' 'Hair Salons'\n", " \"Women's Clothing\" 'Day Spas' 'Laser Hair Removal' 'Hair Removal'\n", " 'Medical Spas' 'Real Estate Agents' 'Juice Bars & Smoothies' 'Gyms'\n", " 'Trainers' 'Lounges']\n" ] } ], "source": [ "i=92\n", "hpom=df[df.cluster == i]\n", "large = hpom.LQ.max()\n", "topBus = hpom.name.ix[hpom.LQ == large].values[0]\n", "topBusCat = hpom.categories.ix[(hpom.name == topBus) & (hpom.LQ == large)].values[0]\n", "bus=hpom.categories.unique()\n", "print i, topBusCat\n", "print bus" ] }, { "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.9" } }, "nbformat": 4, "nbformat_minor": 0 }