5SSG2060_Practical_5_PPA_II.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5SSG2060_Week5_Point Pattern Analysis (PPA)  II\n",
    "<a href=\"#This Week's Overview\">This Week's Overview</a>\n",
    "\n",
    "<a href=\"#Learning Outcomes\">Learning Outcomes</a>\n",
    "\n",
    "<a href=\"#Get prepared\">Get prepared</a>\n",
    "\n",
    "- <a href=\"#Density Estimates\">Density Estimates</a>\n",
    "\n",
    "<a href=\"#Quadrat Based Statistics\">Quadrat Based Statistics</a>\n",
    "- <a href=\"#Create Windows\">Create Windows</a>\n",
    "- <a href=\"#Window attributes & methods\">Window attributes & methods</a>\n",
    "- <a href=\"#Point Intensity Estimates\">Point Intensity Estimates</a>\n",
    "  - <a href=\"#Minimum bounding box\">Minimum bounding box</a>\n",
    "  - <a href=\"#Convex hull\">Convex hull</a>\n",
    "- <a href=\"#Two window testing shapes\">Two window testing shapes</a>\n",
    "\n",
    "<a href=\"#Heatmap of London Pubs\">Heatmap of London Pubs</a>\n",
    "\n",
    "<a href=\"#Nearest Neighbors Classification\">Nearest Neighbors Classification</a>\n",
    "- <a href=\"#Finding Nearest Neighbors\">Finding Nearest Neighbors</a> \n",
    "- <a href=\"#Classification Accuracy\">Classification Accuracy</a>\n",
    "\n",
    "\n",
    "<a href=\"#Kernel Density Estimation (KDE)\">Kernel Density Estimation (KDE)</a>\n",
    "- <a href=\"#Univariate Distribution in 1D\">Univariate Distribution in 1D</a>\n",
    "- <a href=\"#Bivariate Distribution in 2D\">Bivariate Distribution in 2D</a> \n",
    "- <a href=\"#Others channels to understand KDE\">Others channels to understand KDE</a> \n",
    "- <a href=\"#Multivariate KDE\">Multivariate KDE<a/> \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "- <a href=\"#Task 1\">Task 1</a>\n",
    "- <a href=\"#Task 2\">Task 2</a>\n",
    "- <a href=\"#Task 3\">Task 3</a>\n",
    "- <a href=\"#Task 4\">Task 4</a>\n",
    "- <a href=\"#Task 5\">Task 5</a>\n",
    "- <a href=\"#Task 6\">Task 6</a>\n",
    "- <a href=\"#Task 7\">Task 7</a> (Optional)\n",
    "- <a href=\"#Task 8\">Task 8</a>\n",
    "- <a href=\"#Task 9\">Task 9</a>\n",
    "- <a href=\"#Task 10\">Task 10</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"This Week's Overview\">This Week's Overview</a>\n",
    "Last week, we had practiced point pattern analysis with london pubs data and got the preliminary idea that, various methods for point pattern analysis could be categorized into **density-based methods** and **distance-based methods**.\n",
    "\n",
    "- *Distance-based methods* are commonly used for interpreting the second order effects – where and what kind of clusters are formed; thus suitable for identifying local patterns.\n",
    "\n",
    "We had tried Nearest Neighbours Analysis from theoratical view a bit last week, and this week we will continue realize its classification through `Bokeh` plot; in addtion, we will move to Quadrat density estimates, and Kernel Density Estimates (KDE), which falls in \n",
    "- *Density-based methods*, they are often used for understanding the first order effects – how many and how often points can be found, and where their spatial mean is; thus good at capturing the global tendency.\n",
    "\n",
    "These include looking at the density-based methods on measuring point spatial patterns, the corresponding visualizations and how to interpret their utilization. The London pubs data on Week 4 will be used for interpreting **Minimum Bounding Box** and **K-Nearest Neighbours** methods; and the London Airbnb listings data will be the basis for our understanding of **Nearest Neighbours Classification** and **Kernel Density Estimates**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Learning Outcomes\">Learning Outcomes</a> \n",
    "Upon the completion of this practical, you will get comparative understanding on density-based methods for point pattern analysis, especially the following:\n",
    "1. Quadrat based statistics and point density estimates\n",
    "2. Density based visualization by interactive heatmap \n",
    "3. Nearest Neighbours Classification\n",
    "4. Kernel Density Estimates\n",
    "\n",
    "New libraries such as `Bokeh` and `statsmodels` will be introduced in this practical as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Get prepared\">Get prepared</a>\n",
    "Set up a \"**data**\" folder in your directory for this week's practical; \n",
    "\n",
    "copy the **London pubs data** we've used for Week 4, and the **Airbnb listings data** for Week 3 into the \"**data**\" folder for this practical. \n",
    "\n",
    "We will import the general libraries at this stage, but import other specific libraries by section to help you better understand their functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import pysal as ps\n",
    "import libpysal as lps\n",
    "import matplotlib.pyplot as plt\n",
    "import pointpats \n",
    "from pointpats import PointPattern\n",
    "import seaborn as sns\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter('ignore') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Density Estimates\">Density Estimates</a>\n",
    "Through your completion of Geocomputation module, you are already familiar with one density estimation technique:**histogram**, by visualizing the definition of data bins, and the tallied number of data points within each bin; however, its choice of binning can have a disproportionate effect on visualization, hence smoother density estimates using various kernels will be powerful in modeling the distribution of points. \n",
    "\n",
    "The most commonly used kernel density estimate is Gaussian kernel density estimate, in which each point contributes a Gaussian curve to the total. These algorithms could be found both in [scikit-learn](http://scikit-learn.org/stable/) and in [PySAL](http://pysal.readthedocs.io/en/latest/). In scikit-learn it is implemented in the `classification` section at `sklearn.neighbors.KernelDensity estimator`, which uses the Ball Tree or KD Tree for efficient queries. I had included a map here to sklearn's algorithms and how to navigate them:\n",
    "\n",
    "<a href=\"http://scikit-learn.org/stable/tutorial/machine_learning_map/\"><img alt=\"SciKit-Learn Algorithm Map\" src=\"http://scikit-learn.org/stable/_static/ml_map.png\"></a>\n",
    "\n",
    "While in PySAL it will be called from the Python Spatial Analysis Library. PySAL is similarly complex and _also_ has a map to help you navigate its complexities:\n",
    "\n",
    "![PySAL Map](http://darribas.org/gds_scipy16/content/figs/pysal.png)\n",
    "\n",
    "But I can't tell you which is the 'right' approach (as I said above) so far, as it all depends on what you're trying to accomplish and how you're _reasoning_ about your problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Quadrat Based Statistics\">Quadrat Based Statistics</a>\n",
    "\n",
    "As introduced in lecture, `CRS` process has two prominent characteristics: `Uniform` and `Independent`, which will help us to tell a `non-CSR` process by using Quadrat-based statistics. Against the null hypothesis of a `CRS` process, the expected point intensity $\\lambda$ is uniform across the whole study area $\\Re$, so the number of points within a sub-region $|A|$ should be $\\lambda |A|$. \n",
    "\n",
    "If we further define the window $|A|$ as a $m \\times n$ rectangular, then for a `CRS` process, the observed number of points within $|A|$ could be compared against the expected point counts using $\\chi^2$ test statistic:\n",
    "\n",
    "$$\\chi^2 = \\sum^m_{i=1} \\sum^n_{j=1} \\frac{[x_{i,j}-E(x_{i,j})]^2}{\\lambda |A_{i,j}|}$$\n",
    "\n",
    "* For a $\\chi^2$ distribution with $m \\times n -1$ degree of freedom, $p$-value could be derived from the $\\chi^2$ distribution table; and the the null hypothesis will be rejected at the $95\\%$ confidence level if $p$-value is smaller than $0.05$.\n",
    "\n",
    "* For sampling distribution simulations from a large number of $\\chi^2$ test statistics, under the null hypothesis of `CRS` process, if the $\\chi^2$ test statistic for the observed point pattern is among the largest $5%$ test statistics, then it is very unlikely to be result of a `CSR` process at $95\\%$ confidence level. Hence the null hypothesis is rejected followed by a pseudo $p$-value calculation. \n",
    "$$p(\\chi^2) = \\frac{1+\\sum^{nsim}_{i=1}\\phi_i}{nsim+1}$$\n",
    "where \n",
    "$$ \n",
    "\\phi_i =\n",
    " \\begin{cases}\n",
    "    1       & \\quad \\text{if } \\psi_i^2 \\geq \\chi^2 \\\\\n",
    "    0       & \\quad \\text{otherwise } \\\\\n",
    "  \\end{cases}\n",
    "$$\n",
    "\n",
    "$nsim$ is the number of simulations, $\\psi_i^2$ is the $\\chi^2$ test statistic for each simulated point pattern, $\\chi^2$ is the $\\chi^2$ test statistic for the observed point pattern, $\\phi_i$ is an indicator variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary libraries for this section\n",
    "from libpysal.cg import shapely_ext\n",
    "from pointpats import window\n",
    "from pointpats.window import poly_from_bbox, as_window, Window, to_ccf\n",
    "import pointpats.quadrat_statistics as qs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Create Windows\">Create Windows</a>\n",
    "For Quadrat based statistics, **windows** play important roles in  of point patterns analysis. A window can define the domain for the point pattern, can support corrections for **edge effects**, and can also be used for point density analysis. In the following sections, we will use London pubs data as the basis for point pattern analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get point data ready for use\n",
    "f = ps.open('data/London_Pubs.shp')\n",
    "pp_pubs = PointPattern(np.asarray([pnt for pnt in f]))\n",
    "f.close()\n",
    "# general information on london pubs dataset\n",
    "pp_pubs.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may find the summary reported the **Bounding Rectangle** and the **Area of window** for the London Pubs point pattern. As we passed it to the `PointPatterns` constructor the 3337 points' data for London pubs, PySAL finds the [minimum bounding box](https://en.wikipedia.org/wiki/Minimum_bounding_rectangle) and uses it as the window. So the area of window is simply the area of the bounding rectangle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.window.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.window.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.window.bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may get a bounding box been given from left to bottom, then from right to top."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `parts` for `window` is a list of polygons. Now let's check the point containment for the `window`. For example, I want to check whether The Edgar Wallace (with lat and lon at -0.112875 and 51.512704) is contained in this `window`, the test will return a boolean output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.window.contains_point((-0.112875, 51.512704)) # The Edgar Wallace contained in the window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.window.parts # We will introduce the concept of parts in later section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total rows add up to 3337, indicating 3337 pubs as our point targets. What else information you can derive from the summary of london pubs point data? Can you write down your answers to following questions:\n",
    "1. What are the coordinates for the vertices of bounding rectangle?\n",
    "\n",
    "2. What is the area of bounding rectangle window?\n",
    "\n",
    "3. What is the estimated intensity of pubs for this rectangle window?\n",
    "\n",
    "**Hint**: these will be used to compare with point intensity estimates results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the attributes of point patterns\n",
    "pp_pubs.points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Point Intensity Estimates\">Point Intensity Estimates</a>\n",
    "The Point intensity is the mean number of event points per unit of area at given point $p_j$, and can be defined as:\n",
    "\n",
    "$$\\lambda(p_j) = \\lim \\limits_{|\\mathbf{A}p_j| \\to 0} \\left \\{ \\frac{E(Y(\\mathbf{A}p_j)}{|\\mathbf{A}p_j|} \\right \\}   $$\n",
    "\n",
    "where $\\mathbf{A}p_j$ is a region surrounding location $p_j$ with area $|\\mathbf{A}p_j|$, and $E(Y(\\mathbf{A}p_j)$ is the expected number of event points in $\\mathbf{A}p_j$. \n",
    "\n",
    "We may recall one of the implications of CSR process last week, that its point intensity is constant in study area $\\Re$, which could be expressed as:\n",
    "\n",
    "$\\lambda(p_j) = \\lambda(p_{j+1}) = \\ldots = \\lambda(p_n) = \\lambda \\ \\forall p_j \\in \\Re$. Thus, if the area of $\\Re$ = $|\\Re|$, the expected number of event points in the study region is: $E(Y(\\Re)) = \\lambda |\\Re|.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last week, you've already been introduced the concept of `Convex Hull` (the smallest convex set that contains a point pattern) and `Minimum Bounding Rectangle (Box)` (the minimum bounding Rectangle, bigger than convex hull) for Shape Analysis. They could fit into the general concept of `window`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In PySAL, it is estimated by using a geometric object to encode the study region, which is the so-called \"window\", $W$, as introduced in last practical. The reason for distinguishing between $\\Re$ and $W$ is that the latter permits alternative definitions of the bounding object as:\n",
    "\n",
    "$$\\hat{\\lambda} = \\frac{n}{|W|}$$\n",
    "\n",
    "where $n$ is the number of points in the *window* $W$, and $|W|$ is the area of $W$. \n",
    "\n",
    "You may recap on what we've practiced last week covered a bit on the function of minimum bounding box:\n",
    "\n",
    "$$\\hat{\\lambda}_{mbb} = \\frac{n}{|W_{mbb}|}$$\n",
    "\n",
    "where $W_{mbb}$ is the minimum bounding box for the point pattern. Now let's try to understand the definition by using London pubs point dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <a id=\"Minimum bounding box\">Minimum bounding box</a> \n",
    "Now we will call the minimum bounding box function as `lambda_mbb`, to get the estimated density of pub points within the defined window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of pubs if within minimum bounding box window\n",
    "pp_pubs.lambda_mbb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What you've got? Is it exactly same with the above summary result at 10834.35...? YES, you will find the minimum bounding box for london pubs is exactly the rectangle window of greater London! It means based on current pubs density in London, if we take the whole area of the rectangle bounding window, then the estimated pubs in Greater London will reach at 10834 (or 1 more!), which is more than 3337 pubs so far, and a much densier calculated distribution of pubs will be expected.\n",
    "\n",
    "Let's further check the coordinates for bounding rectangle from summary by using Google Map. What have you returned with this time? If you find them as $Ripley$ and $Stondon Massey$, then you had successfully locate the left-bottom and right-top for the bounding rectangle, which is also the boundary for Greater London. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <a id=\"Convex hull\">Convex hull</a> \n",
    "Besides of the rectangle window, we want a shapely closer window. So we can recall the concept of `convex hull` introduced last week, and express the point intensity as:\n",
    "\n",
    "$$\\hat{\\lambda}_{hull} = \\frac{n}{|W_{hull}|}$$\n",
    "\n",
    "where $W_{hull}$ is the convex hull for the point pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of pubs if within convex hull window\n",
    "pp_pubs.lambda_hull"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the area of the window decreased by its getting closer to the London boundary, the estimated pubs number reduced to 14802, but still approximately 3 times more than the real value. \n",
    "\n",
    "We can conclude that, point density estimation is very sensitive to the **area of study region** (or we term it as \"window\" here) as introduced in lecture, hence geographical **scale** is cucial when we conduct point pattern analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Window attributes & methods\">Window attributes & methods</a>\n",
    "Besides of the simple bounding box, windows could more complicated in realistic scenarios. For example, we need a window with 2 parts and 1 hole in it. How you are going to build it? What challenges for this type of window will bring to related point pattern analysis? Let's start with drawing a simple rectangle to understand the concepts of **window, parts and holes**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the parts and holes within the window by their vertices\n",
    "parts = [[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)], # define the vertices for part 1\n",
    "         [(1.5, 1.5), (1.5, 2.0), (2.0, 2.0), (2.0, 1.5)]] # define the vertices for part 2\n",
    "holes = [[(0.5, 0.5), (0.8, 0.5), (0.8, 0.8), (0.5, 0.8)]] # define the vertices for hole\n",
    "\n",
    "# get part 1 plotted as a trial\n",
    "p0 = np.asarray(to_ccf(parts[0]))\n",
    "plt.plot(p0[:,0], p0[:,1])\n",
    "plt.xlim(-1,2)\n",
    "plt.ylim(-1,2) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot part 1 and part 2 as exterior rings in blue\n",
    "for part in parts:\n",
    "    part = np.asarray(to_ccf(part))\n",
    "    plt.plot(part[:,0], part[:,1], 'b') \n",
    "# plot hole associated with the left ring in red\n",
    "for hole in holes:   \n",
    "    hole = np.asarray(to_ccf(hole))\n",
    "    plt.plot(hole[:,0], hole[:,1], 'r') \n",
    "plt.xlim(-1,3)\n",
    "plt.ylim(-1,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have 6 points with varied locations: 2 outside of the exterior rings, 1 point in the hole ring, and 3 points are contained in either exterior ring, then how to evaluate these containments?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up your working framework\n",
    "from pointpats import Window\n",
    "window = Window(parts, holes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pnts = [(1.6,1.6), (0.4,0.4), (0.7, 0.7), (0.2,0.2), (2.5,0.1), (0.5,2.0)]\n",
    "for pnt in pnts:\n",
    "    plt.plot(pnt[0], pnt[1], 'g.') # plot the six points in green\n",
    "\n",
    "for part in parts:\n",
    "    part = np.asarray(to_ccf(part))\n",
    "    plt.plot(part[:,0], part[:,1], 'b') # plot \"parts\" in blue\n",
    "for hole in holes:\n",
    "    hole = np.asarray(to_ccf(hole))\n",
    "    plt.plot(hole[:,0], hole[:,1], 'r') # plot \"hole\" in red \n",
    "\n",
    "from pointpats.window import poly_from_bbox\n",
    "poly = np.asarray(poly_from_bbox(window.bbox).vertices)\n",
    "plt.plot(poly[:,0], poly[:,1], 'm-.') # plot the minimum bounding box in magenta\n",
    "\n",
    "plt.xlim(-1,3)\n",
    "plt.ylim(-1,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to check which points are contained in the window, we can check whether it will be filtered out by calling `window.filter_contained` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check which points are contained in the window\n",
    "pin = window.filter_contained(pnts)\n",
    "pin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you are comfortable with the definitions of **window**, **part** and **hole**, let's fit london pubs data into LSOA boundary and play with them. Firstly, we still need to call `ps.open` function to open the London pubs shapefile. However, you may think that we've already opened it at the beginning of this practical and define it as $\"f\"$, why we can't simply use $\"f\"$ to get the coordinates of all points? Try replace the newly defined $\"pubs\"$ with $\"f\"$, what will you get then?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pubs = ps.open ('data/London_Pubs.shp') # open London pubs polygon shapefile\n",
    "points = [shp for shp in pubs] # get x,y coordinates for all the points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do you still remember how we had loaded the `LDN-LSOAs.shp` data as $\"losas\"$ last week? We also defined them as `polygon` type and used *cascaded union* to get a \"boundary\" window. Recall your experience and write your code below:\n",
    "### <a id=\"Task 1\">Task 1</a>\n",
    "Load the LDN-LSOAs data and get the boundary window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "# open London lsoa polygon shapefile\n",
    "???\n",
    "# define the polygon shapes from London lsoa data\n",
    "???\n",
    "# Create the exterior polygons for Greater London from the union of the polygon shapes\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(boundary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(???) # the count of lsoas in London\n",
    "print(???) # the count of pubs in London"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check if there is any hole in the boundary window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = as_window(boundary)\n",
    "w.holes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the output returns empty, there is no hole in lsoas boundary window. Then we want to know the counts of parts and whether the centroid point in London contained in the boundary window (although we know it should be true as we choose whole London as the window)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(w.parts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.contains_point(w.centroid) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to create a sequence of windows to carry out a containment test, in order to find which lsoa contains the window centroid. \n",
    "\n",
    "We will firstly create a window for each of the individual lsoas, then check each one for the containment of window's centroid. You may find the window (count) with index 4662 is the only one that contains the centroid point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "windows = [as_window(lsoa) for lsoa in polys] # create a window for each of the individual lsoas \n",
    "cent_poly = [ (i, lsoa) for i,lsoa in enumerate(windows) if lsoa.contains_point(w.centroid)]\n",
    "cent_poly # check each lsoa for containment of the window's centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i, cent_poly = cent_poly[0]\n",
    "cent_poly.bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Two window testing shapes\">Two window testing shapes</a> (Optional)\n",
    "You may observed that the shapes of window very by definition; besides of the real \"window\" from geographical shapefile, there are mainly two shapes to conduct statistics for window testing: Rectangle and Hexagon. We are familiar with the former shape as it is set by default, so let's start with the **rectangle** firstly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.plot(window= True, title= \"London Pubs Point pattern\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_r_pubs = qs.QStatistic(pp_pubs,shape= \"rectangle\",nx = 4, ny = 4)\n",
    "q_r_pubs.plot() # plot out the quadrat count figure with 4*4 windows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you figure out the meanings of number in each window? It shows directly the counts of points in each window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Chi-squared test statistic for the observed point pattern is: '+ str(q_r_pubs.chi2)) \n",
    "print('Degree of freedom is: '+str(q_r_pubs.df)) \n",
    "print('P-valus for Chi-squared test statistic is: '+str.format('{0:.6f}', q_r_pubs.chi2_pvalue)) # 6 decimals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the $p$-value based on the analytical $\\chi^2$ distribution (degree of freedom = 15) is 0.000000, much smaller than 0.05. We might reject the null hypothesis of being a `CSR` process.\n",
    "\n",
    "Seconday, let's use hexagon quadrats as windows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_h_pubs = qs.QStatistic(pp_pubs,shape= \"hexagon\",lh =0.08) # lh is the length of hexagon edge, adjustable\n",
    "q_h_pubs.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 2\">Task 2</a>\n",
    "Get the statistics testing result for Hexagon window printed below, and conclude whether we should reject the `CSR` hypothesis? Is the distribution of pubs in London random or not?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the p-value is 0.00000 much smaller than 0.05, we reject the null hypothesis of CSR, and think the point pattern for London pubs is not random. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These varied shapes help to build up the sub windows for statistical purpose, and count the numbers of pubs falling in each sub window. This is actually the elementary basis for generating **heatmap**. So has it reminded you any idea of a high quality interactive heatmap for london pubs?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Heatmap of London Pubs\">Heatmap of London Pubs</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 3\">Task 3</a>\n",
    "Last week, we've tried to get a simple heatmap by calling `scipy`, but it looked not good enough. As we are familiar with using `folium` in the first three weeks' practicals, let's try to work it out in a much fancy way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use geopandas to read london pubs shapefile data\n",
    "pubs_gdf=???('data/london_pubs.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We call folium to realize interactive heatmap of London pubs\n",
    "import folium\n",
    "from folium import plugins\n",
    "from folium.plugins import HeatMap, MarkerCluster, FastMarkerCluster\n",
    "# Ensure you're handing it floats\n",
    "pubs_gdf['lat'] = ???.astype(float)\n",
    "pubs_gdf['lon'] = ???.astype(float)\n",
    "\n",
    "# Filter the DF for rows, then columns, then remove NaNs\n",
    "heat_df = pubs_gdf[['lat', 'lon']]\n",
    "heat_df = heat_df.dropna(axis=0, subset=['lat','lon'])\n",
    "\n",
    "# List comprehension to make out list of lists\n",
    "heat_data = [[row['lat'],row['lon']] for index, row in heat_df.iterrows()]\n",
    "\n",
    "# define the base map\n",
    "heatmap_map = folium.???([51.50632, -0.1271448], zoom_start=12)\n",
    "\n",
    "# Plot data on the map\n",
    "hm=plugins.HeatMap(heat_data)\n",
    "heatmap_map.add_child(hm)\n",
    "# get the map shown below \n",
    "# if it is blank for browser reason, please save it as html file\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to zoom in/out with the heatmap you've got, and find the densiest part of pubs in London? In addition, what's the rationale of `HeatMap` function in `folium`, especially its density basis? Try google it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is straightforward and fancinating, however if we want to know more about the pubs' density in each losa, how we are going to deal with it? \n",
    "\n",
    "Think about the \"Window\" concept, each LSOA could be taken as one sub window defined by administrative boundary, so the measure for density of pubs in each losa is actually to count the number of points (pubs) in each polygon (lsoa) with varied shapes. The starting point should be joining the pubs and lsoas datasets together, by recalling the joining function introduced in week 2 and week 3. Hint: We need to get the columns names listed for both datasets and find the shared column (or rename relevant column for indexing convenience) to realize the data joining."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here\n",
    "pubs_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use geopandas to read london lsoas shapefile data\n",
    "lsoas_gdf=gpd.read_file('data/LDN-LSOAs.shp')\n",
    "lsoas_gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have you decided which column should be used for joining? Yes, the \"lsoa\" and \"lsoa11nm\". But it is better to rename the columns for convenience. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames = ['code','lsoa','nmw','id','area','length','geometry']\n",
    "lsoas_gdf.columns = colnames\n",
    "lsoas_gdf.head(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 4\">Task 4</a>\n",
    "We want to realize the following purposes and you need to replace the question marks with reasonable codes.\n",
    "1. count the total number of pubs within each lsoa, ensure the data type is integer, and reset the index.\n",
    "2. join the lsoas and newly generated count datasets together, and save it as shapefile, get it plotted with raw count number.\n",
    "3. calculate the density of pubs in each lsoa and get it plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating total number of pubs per lsoa\n",
    "lsoas_pubs_count = pd.DataFrame(pubs_gdf['lsoa'].value_counts().astype(int)).reset_index()\n",
    "# only leave the 2 columns for plot generating\n",
    "lsoas_pubs_count.columns = ['lsoa','Numbers']\n",
    "lsoas_pubs_count.head(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add a new column of pubs number in each lsoa \n",
    "# by join lsoa and pubs through attribute join\n",
    "join_gdf = lsoas_gdf.???(lsoas_pubs_count, on='???')\n",
    "# save this newly joined .csv file into .shp file 'lsoa_numbers.shp'\n",
    "???\n",
    "# Make a Choropleth map on pubs per lsoa.\n",
    "join_gdf.plot(column='Numbers', cmap='coolwarm', scheme='quantiles', alpha=0.7, legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output looks quite different from the point pattern presented in heatmap, what is the problem?\n",
    "The output shows us the raw count data in each lsoa, without taking into account of the area of each window, which means if we want the density map, we need to divide the **numbers** of pubs by **area** of the lsoa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the data format, and think of why we need multiple 100000\n",
    "join_gdf['density'] = join_gdf.apply(lambda row: 100000*row.??? / row.???, axis=1)\n",
    "join_gdf.plot(column='density', cmap='coolwarm', scheme='quantiles', alpha=0.7, legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is it much similar with the heatmap point pattern now? \n",
    "\n",
    "Since we've spent a lot of time with London pubs data, now let's change to Airbnb listings data for following realizations on: Nearest Neighbors Classification and Kernel Density Estimation (KDE)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  <a id=\"Nearest Neighbors Classification\">Nearest Neighbors Classification</a>\n",
    "\n",
    "If we recall the theoretical introduction of nearest neighbors classification last week, it realizes classification by training defined sample data using specified weights. So the configuration of **weights** are crucial in neareast neighbours classification, and we normally use default *uniform weights*, where the value setting of weights keywords is:**weights = 'uniform'** (assigns uniform weights to each query neighbor from a simple majority \"vote\" of the nearest neighbors). You may also found another definition of weights by distance as **weights = 'distance'** (assigns weights proportional to the inversed distance from the query point). \n",
    "\n",
    "Upon the your clear setting of weights, we can call two types of nearest neighbors classifiers from `Scikit-learn`:\n",
    "1. **K Neighbors Classifier**: based on the $k$ (integer value, specified by you) nearest neighbors of each query point in the training sample, it is the most commonly used technique. In general, larger $k$ can help to suppress the effects from noise, whilst smaller $k$ can help to make the classification boundaries more distinct.\n",
    "2. **Radius Neighbors Classifier**: based on the number of neighbors within a radius $r$ (fixed, floating-point data, specified by you) of each point in the training sample. Rarely used and specialized for not uniformly sampled data, so you should get the \"curse of dimensionality\" problem in mind if you need to use it.\n",
    "\n",
    "For a better understanding and visualization of the classification results, we will hereby use London Airbnb data as the basis for analysis, call `K neighbors classifier` and try to use `Bokeh` for visulization. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "terrain = sns.color_palette(palette='terrain',n_colors=10)\n",
    "plasma = sns.color_palette(palette='plasma',n_colors=10)\n",
    "rainbow = sns.color_palette(palette='rainbow',n_colors=6)\n",
    "\n",
    "from bokeh.io import output_notebook\n",
    "from bokeh.layouts import gridplot,row,column\n",
    "from bokeh.plotting import figure,show\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bokeh is a powerful data visualization library, you may want to explore its wide application at [Bokeh website](https://bokeh.pydata.org/en/latest/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the dataset for airbnb listings as training dataframe\n",
    "traindf = pd.read_csv('data/airbnb_listings.csv')\n",
    "traindf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you figure out the room types for Airbnb in London from this dataset?\n",
    "As you may find, there are 3 types of room types, Entire home/apt, Private room and Shared room, so we will try to plot them by different colors below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up the plot frame, \n",
    "# define x axis and y axis scales by the upper and lower limits of corresponding values\n",
    "p = figure(title=\"London airbnb room types\",y_range=(51.28,51.7),x_range=(-0.6,0.4))\n",
    "p.xaxis.axis_label = 'longitude'\n",
    "p.yaxis.axis_label = 'latitude'\n",
    "# set up latitude and longitude variables for type of Private room\n",
    "privateLat=traindf['latitude'][traindf['room_type']=='Private room']\n",
    "privateLong=traindf['longitude'][traindf['room_type']=='Private room']\n",
    "# set up latitude and longitude variables for type of Shared room\n",
    "sharedLat=traindf['latitude'][traindf['room_type']=='Shared room']\n",
    "sharedLong=traindf['longitude'][traindf['room_type']=='Shared room']\n",
    "# set up latitude and longitude variables for type of Entire home/apt\n",
    "entireLat=traindf['latitude'][traindf['room_type']=='Entire home/apt']\n",
    "entireLong=traindf['longitude'][traindf['room_type']=='Entire home/apt']\n",
    "# define different colors for each type\n",
    "p.circle(privateLong,privateLat,size=3,color=terrain.as_hex()[1],fill_alpha=0.1,line_alpha=0.4,legend='Private room')\n",
    "p.circle(sharedLong,sharedLat,size=3,color=plasma.as_hex()[9],fill_alpha=0.1,line_alpha=0.4,legend='Shared room')\n",
    "p.circle(entireLong,entireLat,size=3,color=plasma.as_hex()[5],fill_alpha=0.1,line_alpha=0.4,legend='Entire home/apt')\n",
    "show(p, notebook_handle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may already noticed that some buttons generated together with the plot on the right column. Try click on the buttons, and get the right way to zoom in certain areas in London. Of course, you also need to find your way back with another button. Once you zoomed in, you may find that areas e.g. near East Village, Chelsea, Hell's Kitchen, and Upper East Side, look brighter than the neighbors, indicating higher interest or more requests from Airbnb \"hunters\". In order to figure out which specific type of Airbnb really contribute to the distribution, we can plot them out by type separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot private room\n",
    "p1 = figure(width=500, height=300, title='Private airbnb rooms distribution in London',y_range=(51.28,51.7),x_range=(-0.6,0.4))\n",
    "p1.circle(privateLong,privateLat,size=3,color=terrain.as_hex()[1],fill_alpha=0.2,line_alpha=0.1,legend='Private room')\n",
    "# plot shared room\n",
    "p2 = figure(width=500, height=300, title='Shared airbnb rooms distribution in London',y_range=(51.28,51.7),x_range=(-0.6,0.4))\n",
    "p2.circle(sharedLong,sharedLat,size=3,color=plasma.as_hex()[9],fill_alpha=0.2,line_alpha=0.1,legend='Shared room')\n",
    "# plot entire home/apt\n",
    "p3 = figure(width=500, height=300, title='Entire airbnb homes/apartments distribution in London',y_range=(51.28,51.7),x_range=(-0.6,0.4))\n",
    "p3.circle(entireLong,entireLat,size=3,color=plasma.as_hex()[5],fill_alpha=0.2,line_alpha=0.1,legend='Entire home/apt')\n",
    "# get plots show\n",
    "show(column(p1,p2,p3), notebook_handle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you may get some rough ideas on where are the hot airbnbs accummulating, are they normally around some popular venues, siteseeing places, etc.? Does an airbnb receiving high popularity in booking indicate the well acceptance from the customers? If so, how are other neighboring airbnb accommodations? To get the answers, we can use KNN as one of the optimal choices, and starting with creating dataframes from point coordinates data the dependent variable, which is **room_type**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Finding Nearest Neighbors\">Finding Nearest Neighbors</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loading the libraries\n",
    "from matplotlib.colors import ListedColormap\n",
    "from sklearn import neighbors, datasets\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "import scipy as sp\n",
    "\n",
    "# define independant variable X with the points' coordinates\n",
    "# define dependant variable y with room type\n",
    "X = pd.concat([traindf['latitude'],traindf['longitude']],axis=1)\n",
    "y = traindf['room_type']\n",
    "h = .02  # step size in the mesh\n",
    "\n",
    "# we'll split the training data into testing and training sets\n",
    "# test_size is the proportion of dataset to include in the test split\n",
    "# random_state is the count of seeds used by random number generator\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=45) \n",
    "# execute the classification, define k as 9\n",
    "neigh = KNeighborsClassifier(n_neighbors=9)\n",
    "neigh.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Upon classification, we can further predict the \"unseen\" room type, and check the capability of predication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the predicted y value for tested X\n",
    "predVal=neigh.predict(X_test)\n",
    "# include both the predicted y value and real tested y value\n",
    "mat=[predVal,y_test]\n",
    "# read both into a single dataframe\n",
    "testdf=pd.DataFrame(mat).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will get a dataframe with two columns \"0\" and \"1\", let's made them more understandable by renaming the column names into 'prey' and 'y'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testdf.columns=('prey','y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Classification Accuracy\">Classification Accuracy</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the accuracy in predicting the right room type\n",
    "# compare predicted y value and the real y value\n",
    "testdf['diff']=np.where(testdf.prey==testdf.y,1,0)\n",
    "# calculate the proportion of matched pairs\n",
    "print('% correct =',sum(testdf['diff'])/len(testdf['diff'])*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we had the predicting accuracy at 60.69%, we will further build out the [log loss function](http://wiki.fast.ai/index.php/Log_Loss) to see how much error in the predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PredProb=neigh.predict_proba(X_test)\n",
    "pred=np.asmatrix(PredProb)\n",
    "pred.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred.columns=('Entire home/apt','Private room','Shared room')\n",
    "s=np.asmatrix(pd.get_dummies(y_test))\n",
    "def f(x):\n",
    "    return sp.log(sp.maximum(sp.minimum(x,1-10**-5),10**-5))\n",
    "f=np.vectorize(f)\n",
    "predf=f(pred)\n",
    "mult=np.multiply(predf,s)\n",
    "print('log loss =',np.sum(mult)/-len(y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This log loss is quite high so let's see if we can improve this by increasing our k value. Since, it would be annoying to change the value and run it, I figure it'll be faster to run a for loop through values of $k$ from odd numbers between 3 to 39 (represented by $j$). I also wanted to have at least 5 samples in each $k$ to give us a good average (represented by $i$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accbig=[]\n",
    "loglossbig=[]\n",
    "\n",
    "def f(x):\n",
    "    return sp.log(sp.maximum(sp.minimum(x,1-10**-5),10**-5))\n",
    "f=np.vectorize(f)\n",
    "\n",
    "for j in range(3,40,2):\n",
    "    logloss=[]\n",
    "    acc=[]\n",
    "    for i in range(5):\n",
    "        #split data\n",
    "        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=i)\n",
    "        neigh = KNeighborsClassifier(n_neighbors=j)\n",
    "        #train classifier\n",
    "        neigh.fit(X_train, y_train)\n",
    "        \n",
    "        #find % predicted correctly for this k\n",
    "        predVal=neigh.predict(X_test)\n",
    "        mat=[predVal,y_test]\n",
    "        testdf=pd.DataFrame(mat).transpose()\n",
    "        testdf.columns=('prey','y')\n",
    "        testdf['diff']=np.where(testdf.prey==testdf.y,1,0)\n",
    "        acc.append(sum(testdf['diff'])/len(testdf['diff']))\n",
    "        \n",
    "        #find the logloss for this k\n",
    "        PredProb=neigh.predict_proba(X_test)\n",
    "        pred=np.asmatrix(PredProb)\n",
    "        pred.columns=('Entire home/apt','Private room','Shared room')\n",
    "        s=np.asmatrix(pd.get_dummies(y_test))\n",
    "        predf=f(pred)\n",
    "        mult=np.multiply(predf,s)\n",
    "        logloss.append(np.sum(mult)/-len(y_test))\n",
    "    loglossbig.append(np.mean(logloss))\n",
    "    accbig.append(np.mean(acc))\n",
    "print(accbig)\n",
    "print(loglossbig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot this against every K to see the decrease\n",
    "plt.plot(range(3,40,2),loglossbig)\n",
    "plt.ylabel('logloss')\n",
    "plt.xlabel('k value')\n",
    "plt.title('KNN logloss on longitude and latitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(range(3,40,2),accbig)\n",
    "plt.ylabel('% predicted correctly')\n",
    "plt.xlabel('k value')\n",
    "plt.title('KNN prediction on longitude and latitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though the performance for this isn't amazing, we can see that there are some predictive value that we can use from longtitude and latitude data. We could further check the accuracy of predication at specific $k$ value. For example, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "# instantiate learning model (k = 25) \n",
    "knn = KNeighborsClassifier(n_neighbors=25)\n",
    "# fitting the model\n",
    "knn.fit(X_train, y_train)\n",
    "# predict the response\n",
    "pred = knn.predict(X_test)\n",
    "# evaluate accuracy\n",
    "print (accuracy_score(y_test, pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Task 6'>Task 6</a> \n",
    "Try to change the value of $k$, have you received different predicted accuracy scores? Please list out the results below:\n",
    "1. k=3, accuracy score=???\n",
    "2. k=5, accuracy score=???\n",
    "3. k=10, accuracy score=???\n",
    "4. k=15, accuracy score=???\n",
    "5. k=20, accuracy score=???\n",
    "6. k=30, accuracy score=???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Task 7'>Task 7</a> (Optional)\n",
    "\n",
    "Price distribution for Airbnb in London (purely statistics), 3 missions to complete.\n",
    "1. Airbnb price distribution in London from west to east. \n",
    "\n",
    "   **Hint**: Scatter plotted the price against longitude calling `seaborn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to interpret the distribution in a geographcial language then? Your answers below:\n",
    "\n",
    "**---------------------------------------------------------**\n",
    "\n",
    "2. Airbnb price distribution in London from south to north.\n",
    "\n",
    "   **Hint:** Scatter plotted the price against latitude calling `seaborn`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to interpret the distribution in a geographcial language then? Your answers below:\n",
    "\n",
    "**---------------------------------------------------------**\n",
    "\n",
    "3. The average price and price range for each room type.\n",
    "   **Hint:** Price statistics by room type, using `groupby`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to interpret the Output? Your answers below:\n",
    "\n",
    "**---------------------------------------------------------**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, some unbelievable prices pop up, and we want to get rid of these outliers for further analysis. So let's firstly get a general idea of the spatial distribution among defined price ranges, say, below £100, £100-£200, £200-£400, £400-£600 and above £600, and call `Bokeh` to plot again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define your points' coordinates for each price range\n",
    "Lat100=traindf['latitude'][traindf['price']<100]\n",
    "Long100=traindf['longitude'][traindf['price']<100]\n",
    "Lat200=traindf['latitude'][(traindf['price']<200)&(traindf['price']>=100)]\n",
    "Long200=traindf['longitude'][(traindf['price']<200)&(traindf['price']>=100)]\n",
    "Lat400=traindf['latitude'][(traindf['price']<400)&(traindf['price']>=200)]\n",
    "Long400=traindf['longitude'][(traindf['price']<400)&(traindf['price']>=200)]\n",
    "Lat600=traindf['latitude'][(traindf['price']<600)&(traindf['price']>=400)]\n",
    "Long600=traindf['longitude'][(traindf['price']<600)&(traindf['price']>=400)]\n",
    "Latup=traindf['latitude'][(traindf['price']>=600)]\n",
    "Longup=traindf['longitude'][(traindf['price']>=600)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = figure(title=\"Airbnb prices\",y_range=(51.28,51.7),x_range=(-0.6,0.4))\n",
    "p.xaxis.axis_label = 'latitude'\n",
    "p.yaxis.axis_label = 'longitude'\n",
    "\n",
    "# plot points in each price range with different colors\n",
    "p.circle(Long100,Lat100,size=3,color=rainbow.as_hex()[0],fill_alpha=0.6,line_alpha=0.6,legend='<£100')\n",
    "p.circle(Long200,Lat200,size=3,color=rainbow.as_hex()[1],fill_alpha=0.6,line_alpha=0.6,legend='£200')\n",
    "p.circle(Long400,Lat400,size=3,color=rainbow.as_hex()[4],fill_alpha=0.6,line_alpha=0.6,legend='£400')\n",
    "p.circle(Long600,Lat600,size=3,color=rainbow.as_hex()[5],fill_alpha=0.6,line_alpha=0.6,legend='£600')\n",
    "p.circle(Latup,Longup,size=3,color=rainbow.as_hex()[3],fill_alpha=0.9,line_alpha=1,legend='up')\n",
    "p.legend.location = 'top_right'\n",
    "show(p, notebook_handle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I intentially highlighted the points with price above £600 by higher definition of its point transparency, however, we still can't spotted too much in that category. It indicates that majority of the renting prices are below £600 roughly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Kernel Density Estimation (KDE)\">Kernel Density Estimation (KDE)</a> \n",
    "\n",
    "A density plot is a smoothed, continuous description of the data distribution, and [kernel density estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation) could be the most typical density plot receiving wide popularity. It draws a continuous curve (the kernel) at each data point and adds all of these curves together to make a smooth density estimation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most often used kernel is a Gaussian, within which the kernel is mathmetically a positive function for `k` and `h`:\n",
    "<img  height=\"100\" src=\"kernel function.png\">\n",
    "\n",
    "where `h` as the bandwidth parameter, controlling the tradeoff between bias and variance in the result. Larger bandwidth leads to smoother (i.e. high-bias) density distribution, while smaller bandwidth leads to an unsmooth (i.e. high-variance) density distribution, and finally produces a Gaussian bell curve at each data point. More detailed explaination could be explored from [Kernel Density Estimation](https://scikit-learn.org/stable/modules/density.html#kernel-density-estimation). \n",
    "\n",
    "Besides of the common Gaussian Kernel, there are also other 5 kernel forms used in sklearn:\n",
    "<img  height=\"100\" src=\"kernel formats .png\">\n",
    "\n",
    "The most common use of KDE is in graphically representing distributions of points. There are many ways to visualize the results of Kernal Density Estimates, e.g. the `seaborn` plot we used last semester and in this practical. In the `seaborn` visualization library, KDE is built in and automatically used to help visualize points in one and two dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Univariate Distribution in 1D\">Univariate Distribution in 1D</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# histogram of price below £600\n",
    "# your code here\n",
    "sns.distplot(traindf['price'][traindf['price']<600])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(color_codes=True)\n",
    "kdeprice=traindf.price[traindf['price']<600]\n",
    "# call KDE function in seaborn\n",
    "sns.kdeplot(kdeprice, shade=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you check the parameters for [seaborn.kdeplot](https://seaborn.pydata.org/generated/seaborn.kdeplot.html), you will find the \"bw\" (bandwidth) parameter controlling how tightly the estimation is fit to the data, similar to the bin size in a histogram. As corresponding to the width of the kernels we plotted above, we may try larger or smaller values instead, to compare with the default behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(kdeprice)\n",
    "sns.kdeplot(kdeprice, bw=.2, label='bw: 0.2')\n",
    "sns.kdeplot(kdeprice, bw=2, label='bw: 2')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far for this Airbnb price data, it seems that bandwidth at 2 is much better than that bw at 0.2. The nature of the Gaussian KDE process means that estimation extends past the extreme (largest and smallest) values in the dataset; so we can control where we want to \"start\" or \"stop\" the past by defining the \"cut\" parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(kdeprice, shade=True, cut=0)\n",
    "sns.rugplot(kdeprice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot cut at \"0\" as defined, however, you may also spotted that only the drawn curve is affected, but the fit line doesn't change at all."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 8\">Task 8</a>\n",
    "Pandas dataframe plot we've introduced last semester can also plot KDE curve for univariate distribution. Could you please try to get it done? Further to compare it with the seaborn's kdeplot output.\n",
    "\n",
    "**Hint**: The sub function under 'plot' is $kde$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Bivariate Distribution in 2D\">Bivariate Distribution in 2D</a> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.jointplot(x='longitude', y='latitude', data=traindf, kind=\"kde\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also draw a two-dimensional kernel density plot with the `kdeplot` function. This allows you to draw this kind of plot onto a specific (and possibly already existing) matplotlib axes, whereas the `jointplot` function manages its own figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(6, 6))\n",
    "sns.kdeplot(traindf.longitude, traindf.latitude, ax=ax)\n",
    "sns.rugplot(traindf.longitude, color=\"g\", ax=ax)\n",
    "sns.rugplot(traindf.latitude, vertical=True, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish to show the bivariate density more continuously, you can simply increase the number of contour levels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(6, 6))\n",
    "cmap = sns.cubehelix_palette(as_cmap=True, dark=0, light=1, reverse=True)\n",
    "sns.kdeplot(traindf.longitude, traindf.latitude, cmap=cmap, n_levels=60, shade=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `jointplot` function uses a JointGrid to manage the figure. For more flexibility, you may want to draw your figure by using JointGrid directly. `jointplot` returns the JointGrid object after plotting, which you can use to add more layers or to tweak other aspects of the visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.jointplot(x='longitude', y='latitude', data=traindf, kind=\"kde\", color=\"m\")\n",
    "g.plot_joint(plt.scatter, c=\"w\", s=30, linewidth=1, marker=\"+\")\n",
    "g.ax_joint.collections[0].set_alpha(0)\n",
    "g.set_axis_labels(\"$longitude$\", \"$latitude$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 9\">Task 9</a> \n",
    "Plot bivariate distribution between Airbnb price and corresponding distance to City Center.\n",
    "1. Calculate the distances between airbnb and city center, which is Charing Cross at (51.5081, -0.1248).\n",
    "2. Apply what we've covered before to plot your own bivariate distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "# define the coordinates of each point\n",
    "x1 = traindf['latitude']\n",
    "y1 = traindf['longitude']\n",
    "\n",
    "# an empty list for distance\n",
    "distance=[]\n",
    "# loop to calculate the distance between each point and city center\n",
    "for i in range(len(x1)):\n",
    "    xa, ya=x1[i], y1[i]\n",
    "    distance.append((((xa-51.5081)**2)+((ya+0.1284)**2))**(1/2))\n",
    "# list the distanc results\n",
    "distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.jointplot(x=distance, y=traindf['price'], data=traindf, kind=\"kde\", color=\"m\")\n",
    "g.plot_joint(plt.scatter, c=\"r\", s=50, linewidth=1, marker=\"+\")\n",
    "g.ax_joint.collections[0].set_alpha(0)\n",
    "g.set_axis_labels(\"$distance$\", \"$price$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to interpret the result above?\n",
    "\n",
    "**----------------------------------------------------**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Others channels to understand KDE\">Others channels to understand KDE</a> \n",
    "\n",
    "Just as stated in previous practicals, there are always multiple ways to reach at your goals and get problems solved; so does KDE. Upon your understanding of the statistical rationale of KDE, you may find other methods to fulfill your requests, for example: \n",
    "1. `gaussian_kde` in `SciPy`: contains only the common Gaussian Kernel; \n",
    "2. `KernelDensity` in `Scikit-learn`: contains six kernels, each of which can be used with one of about a dozen distance metrics, resulting in a very flexible range of effective kernel shapes; \n",
    "3. `KDEUnivariate` & `KDEMultivariate` in Statsmodels: contains seven kernels.\n",
    "\n",
    "We will make use of the aforementioned methods in order for our airbnb data, to get you more understanding of the parameters of KDE: **kernel**, which specifies the shape of the distribution placed at each point; and **bandwidth**, which controls the size of the kernel at each point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Univariate KDE in 1D\n",
    "We will keep using the London Airbnb price data, which are less than £600, as defined in previous section as $kdeprice$, to demonstrate the principles of Kernel Density Estimation in one dimension. Firstly, let's get a rough idea of the data by creating histogram simply with the `plt.hist()` function, and with the density parameter specification to get the height of the bins reflecting probability density. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=np.array(kdeprice)\n",
    "mean=np.mean(X)\n",
    "var=np.var(X)\n",
    "std=np.sqrt(var)\n",
    "X_plot = np.linspace(min(X),max(X), 1000) # set up a sample size at 1000\n",
    "print(mean,var,std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can read the mean value for price from the output and use it as benchmark for this section:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute histogram of Samples\n",
    "n, bins,patches=plt.hist(X,bins=100,density=True,facecolor='g',edgecolor='red', alpha=0.8)\n",
    "# try to change the value of bins, figure out the variations\n",
    "plt.axvline(label='Mean=£92',x=92,linestyle='dashed')\n",
    "plt.axvline(label='Mean-2sigma=-£53',x=-53,linestyle='dashed')\n",
    "plt.axvline(label='Mean+2sigma=£237',x=237, linestyle='dashed')\n",
    "plt.xlabel('Prices(£)')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Histogram of Airbnb prices distribution')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "X = np.array(kdeprice) # load the data\n",
    "count = sum(norm(Xi).pdf(X_plot) for Xi in X)  # pdf stands for probability distribution function\n",
    "\n",
    "plt.fill_between(X_plot, count, alpha=0.75)\n",
    "plt.plot(X, np.full_like(X, -0.1), '|k', markeredgewidth=1)\n",
    "plt.axis([-20, 620, 0, 800]) # set the axis limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that scipy weights its bandwidth by the covariance of the input data.\n",
    "# To make the results comparable to the other methods, we use the bandwidth at 0.2\n",
    "from scipy.stats import gaussian_kde\n",
    "kde_scipy = gaussian_kde(X, bw_method=0.2)\n",
    "log_dens_scipy = kde_scipy.evaluate(X_plot)\n",
    "\n",
    "plt.fill_between(X_plot, log_dens_scipy, alpha=0.5)\n",
    "plt.plot(X, np.full_like(X, -0.1), '|k', markeredgewidth=1)\n",
    "plt.ylim(-0.001, 0.014)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import kurtosis, skew\n",
    "print(\"Excess kurtosis of normal distribution ( should be 0):{}\".format(kurtosis(X)))\n",
    "print(\"Skewness of normal distribution ( should be 0):{}\".format(skew(X)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's turn from `Scipy` to `Scikit-learn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors.kde import KernelDensity\n",
    "# instantiate and fit the Gaussian KDE model\n",
    "kde_sklearn = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X[:,None])\n",
    "# score_samples returns the log of the probability density\n",
    "log_dens_sklearn = kde_sklearn.score_samples(X_plot[:, None])\n",
    "\n",
    "plt.fill_between(X_plot, np.exp(log_dens_sklearn), alpha=0.75)\n",
    "plt.plot(X, np.full_like(X, -0.1), '|k', markeredgewidth=1)\n",
    "plt.ylim(-0.001, 0.014)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about the 3rd method in `statsmodels`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.nonparametric.kde import KDEUnivariate\n",
    "X_new=X.astype(np.double)  # Try to explain do we need this code? what will happen if you still use variable X?\n",
    "kde_statsmu = KDEUnivariate(X_new)\n",
    "kde_statsmu.fit() # Estimate the densities\n",
    "\n",
    "fig = plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.hist(X_new, bins=100, normed=True, color='red')\n",
    "ax.plot(kde_statsmu.support, kde_statsmu.density, lw=-2, color='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "# Plot the histrogram\n",
    "ax.hist(X_new, bins=20, normed=True, label='Histogram from samples', \n",
    "        zorder=5, edgecolor='k', alpha=0.5)\n",
    "# Plot the KDE as fitted using the default arguments\n",
    "ax.plot(kde_statsmu.support, kde_statsmu.density, lw=3, label='KDE from samples', zorder=10)\n",
    "ax.legend(loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "# Plot the true distribution\n",
    "true_values = (norm.pdf(loc=-1, scale=0.5, x=kde_statsmu.support)*0.25\n",
    "              + norm.pdf(loc=1, scale=0.5, x=kde_statsmu.support)*0.75)\n",
    "\n",
    "# different defination of location, scale and weights for 2 normal distribution dataset\n",
    "ax.plot(kde_statsmu.support, true_values, lw=3, label='True distribution', zorder=15)\n",
    "\n",
    "# Plot the samples\n",
    "ax.scatter(X_new, np.abs(np.random.randn(X_new.size))/40, \n",
    "           marker='o', color='k', zorder=20, label='Samples', alpha=0.5)\n",
    "ax.legend(loc='best')\n",
    "ax.grid(True, zorder=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 10\">Task 10</a>\n",
    "1. Try change the value of bandwidth and get your own plots;\n",
    "2. Try change the kernel and get your own plots;\n",
    "3. Could you summarize what are the differences among these methods?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional\n",
    "Be careful! It takes time!\n",
    "\n",
    "Since the KDE is a distribution, we can access attributes and methods such as:\n",
    "`entropy`, `cdf`, `icdf`, `sf` and `cumhazard`. Try to explore the meaning of them only if you want to include in your own project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 5))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(kde_statsmu.support, kde_statsmu.cdf, lw=3, label='CDF')\n",
    "ax.plot(np.linspace(0, 1, num = kde_statsmu.icdf.size), kde_statsmu.icdf, lw=3, label='Inverse CDF')\n",
    "ax.plot(kde_statsmu.support, kde_statsmu.sf, lw=3, label='Survival function')\n",
    "ax.legend(loc = 'best')\n",
    "ax.grid(True, zorder=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 5))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(kde_statsmu.support, kde_statsmu.cumhazard, lw=3, label='Cumulative Hazard Function')\n",
    "ax.legend(loc = 'best')\n",
    "ax.grid(True, zorder=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Multivariate KDE\">Multivariate KDE<a/> \n",
    "Now we will look at a slightly more sophisticated use of KDE's visualization for the distribution of different  airbnb room types in London distributions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the AirBnB data into a GeoDataFrame called `sdf` (spatial sample  dataframe)\n",
    "import os\n",
    "import shapely\n",
    "from shapely.geometry import Point\n",
    "geometry = [Point(xy) for xy in zip(traindf.longitude, traindf.latitude)]\n",
    "crs      = {'init': 'epsg:4326'} # What projection is lat/long?\n",
    "sdf      = gpd.GeoDataFrame(traindf, crs=crs, geometry=geometry)\n",
    "sdf      = sdf.to_crs({'init': 'epsg:27700'}) # Reproject into OSGB\n",
    "\n",
    "# Check the output\n",
    "sdf.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## airbnb listing data read csv file 3 types of airbnb\n",
    "sdf_group = sdf.groupby('room_type').size().reset_index(name='count')\n",
    "sdf_sort = sdf_group.sort_values(['count'], ascending=False).reset_index(drop=True)\n",
    "sdf_sort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_type=sdf_sort['room_type']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdf['longitude'].notnull().count() # check whether they all have longitude data\n",
    "# Your code here to replace the ???\n",
    "sdf['latitude'].notnull().count() # check whether they all have latitude data\n",
    "# If you get the same results on rows, that'll be great!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the mean value of coordinates for each room type\n",
    "sdf_group = sdf.groupby('room_type')\n",
    "loc_sdf_group = sdf_group[['latitude', 'longitude']].mean()\n",
    "x = loc_sdf_group.longitude.values\n",
    "y = loc_sdf_group.latitude.values\n",
    "loc_sdf_group.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "# Set up the plotting base\n",
    "new_gdb = gpd.GeoSeries(sdf[['longitude', 'latitude']].apply(Point, axis=1), crs=\"+init=epsg:4326\")\n",
    "bbox = new_gdb.total_bounds\n",
    "titles=[\"Kernel Density: \"+in_type[i] for i in range(3)]\n",
    "fig.tight_layout(pad = 0.4, w_pad = 4.0, h_pad = 4.0)\n",
    "\n",
    "fig, ax = plt.subplots(1,3,figsize=(15,6))\n",
    "\n",
    "# sets the x and y limits for each function \n",
    "# otherwise the density maps will be very small\n",
    "ax0 = plt.subplot2grid((8,8), (0,0),rowspan=3, colspan=3) \n",
    "ax0.set_title(titles[0], fontsize =16)\n",
    "ax0.set_xlim(bbox[0], bbox[2])\n",
    "ax0.set_ylim(bbox[1], bbox[3]) \n",
    "\n",
    "ax1 = plt.subplot2grid((8,8), (4,0),rowspan=3, colspan=3)\n",
    "ax1.set_title(titles[1], fontsize =16)\n",
    "ax1.set_xlim(bbox[0], bbox[2])\n",
    "ax1.set_ylim(bbox[1], bbox[3]) \n",
    "\n",
    "ax2 = plt.subplot2grid((8,8), (0,4),rowspan=3, colspan=3)\n",
    "ax2.set_title(titles[2], fontsize =16)\n",
    "ax2.set_xlim(bbox[0], bbox[2])\n",
    "ax2.set_ylim(bbox[1], bbox[3]) \n",
    "\n",
    "gdfnew0 = sdf[sdf['room_type']==in_type[0]]\n",
    "sns.kdeplot(gdfnew0.longitude, gdfnew0.latitude, shade = True, cmap = 'Reds', ax=ax0) \n",
    "gdfnew1 = sdf[sdf['room_type']==in_type[1]]\n",
    "sns.kdeplot(gdfnew1.longitude, gdfnew1.latitude, shade = True, cmap = 'Blues', ax=ax1)\n",
    "gdfnew2 = sdf[sdf['room_type']==in_type[2]]\n",
    "sns.kdeplot(gdfnew2.longitude, gdfnew2.latitude, shade = True, cmap = 'Greens', ax=ax2)\n",
    "\n",
    "sns.set(style = 'whitegrid') # aesthetics\n",
    "sns.despine(left=True) # aesthetics\n",
    "sns.set_context('paper') # aesthetics\n",
    "plt.axis('equal')\n",
    "plt.show()   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Credits!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take home work -- GIS realization (Optional)\n",
    "\n",
    "**Task**: Nearest Neighbour Distances analysis\n",
    "\n",
    "For ArcGIS user, hints:\n",
    "- in ArcMap, activate the Spatial Analyst extension (i.e. check that the checkbox next to the Spatial Analysis extension is checked in the “Extension” menu);\n",
    "- in ArcToolBox, open Spatial Statistics Tools > Analyzing Patterns > Average Nearest Neighbor to carry out NN analysis.\n",
    "- a report will be produced to explain the distribution of the nearest neighbor distances.\n",
    "\n",
    "For QGIS user, hints:\n",
    "- Load the London Pub data (CRS: WGS 84, Authority ID: EPSG:4326). \n",
    "- Select from the main menu, Vector -> Analysis Tools -> Distance Matrix….\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Geocomp",
   "language": "python",
   "name": "gsa2018"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}