Unverified Commit bb2ae0b6 authored by Simon Bowly's avatar Simon Bowly
Browse files

ADS1002 logistic regression.

parent 627e6093
%% Cell type:markdown id: tags:
# Logistic Regression
%% Cell type:markdown id: tags:
* Topic: Supervised machine learning
* Unit: ADS1002
* Level: Beginner
* Authors: Simon Clarke, Simon Bowly
* Version: 3.1
In this lesson we will discuss the logistic regression classifier model. This is closely related to linear regression, but predicts probabilities for classification, instead of continuous values. We will use the [Iris dataset](https://en.wikipedia.org/wiki/Iris_flower_data_set), which is built in to `seaborn` and `sklearn`.
We will firstly explain how logistic regression can be used to classify binary data, by considering just one species of iris. We will then explain how the accuracy of logistic regression can be evaluated using a range of methods. The exercise will then consider multiclass regression.
First, import various libraries which we require.
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns;
import pandas as pd
```
%% Cell type:markdown id: tags:
## Exploratory Data Analysis
%% Cell type:markdown id: tags:
We now import the iris data set and undertake some exploratory data analysis. The data frame has 150 rows and 5 columns. Viewing the header shows there are four features: sepal_length, sepal_width, petal_length and petal_width. These are all in cm. There is one target variable, which is the species of iris.
%% Cell type:code id: tags:
``` python
iris = sns.load_dataset('iris') # load the dataset from seaborn
print('Shape of the iris dataset is',iris.shape) # display the shape of the data
iris.head() # display the first few lines
```
%%%% Output: execute_result
sepal_length sepal_width petal_length petal_width species
0 5.1 3.5 1.4 0.2 setosa
1 4.9 3.0 1.4 0.2 setosa
2 4.7 3.2 1.3 0.2 setosa
3 4.6 3.1 1.5 0.2 setosa
4 5.0 3.6 1.4 0.2 setosa
%% Cell type:markdown id: tags:
By viewing the statistics we see that there are no missing entries and the standard deviation of each column is fairly similar. However, we will normalise the data before modelling.
%% Cell type:code id: tags:
``` python
iris.describe() # show the statistics of the numerical columns
```
%%%% Output: execute_result
sepal_length sepal_width petal_length petal_width
count 150.000000 150.000000 150.000000 150.000000
mean 5.843333 3.057333 3.758000 1.199333
std 0.828066 0.435866 1.765298 0.762238
min 4.300000 2.000000 1.000000 0.100000
25% 5.100000 2.800000 1.600000 0.300000
50% 5.800000 3.000000 4.350000 1.300000
75% 6.400000 3.300000 5.100000 1.800000
max 7.900000 4.400000 6.900000 2.500000
%% Cell type:markdown id: tags:
The target field is a string, and has three unique values.
%% Cell type:code id: tags:
``` python
iris['species'].unique() # show the unique values of the species column
```
%%%% Output: execute_result
array(['setosa', 'versicolor', 'virginica'], dtype=object)
%% Cell type:markdown id: tags:
Plotting the classification against the first two features, it is apparent that setosa is fairly well separated in feature space, whereas versicolor and virginica have a blurred boundary. Note that the boundary may differ using other features. We would expect that classifiers should be able to easily classify setosa, though the other two species may be a bit more difficult. You can experiment boundaries by changing the features or using `sns.pairplot`.
%% Cell type:code id: tags:
``` python
# plot the first two fields of the iris dataset and classify the points based on the species
# fit_reg=False indicates don't fit a regression model
# hue='species' indicates to classify the points based on the species field
sns.lmplot(x="sepal_length",y="sepal_width",hue='species',data=iris,fit_reg=False);
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAa8AAAFgCAYAAAAM1fZgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1lklEQVR4nO3df5xcdX3v8ddnZyfZXZKQKIniJhIoIhD5IaSC/EhTEAWl8fbKLfShwdUiEVHQXrTKtUjTNgWlV4lSDKBdiLaiKLdRfkYwJkhJDRgDIZHLDT+SFc0iCSTsbDKZ+dw/zlmYHXZnZufsmTNn9/3ksY/ZPT/mfOYsySfnzHfeX3N3RERE0qQl6QJERERGSs1LRERSR81LRERSR81LRERSR81LRERSR81LRERSpyHNy8wyZvYrM/vJEOu6zKzXzNaHXxc0oiYREUmv1gYd51JgEzBlmPW3uvsna32yM8880+++++5RKUxEZBRZ0gWMF7FfeZnZTOB9wE2j9ZzPP//8aD2ViIikUCNuG34N+BxQrLDNB8xsg5ndZmazhtrAzC40s3Vmtq63tzeOOkVEJCVibV5mdjaw3d0frrDZj4HZ7n40sBK4eaiN3P0Gd5/r7nOnT58eQ7UiIpIWcV95nQwsMLOnge8Bp5nZd0o3cPc/uPue8MebgONjrklERFIu1ubl7l9w95nuPhs4D7jf3T9Uuo2ZHVjy4wKCgR0iIiLDatRow0HMbDGwzt1XAJeY2QJgH/AC0JVETSIikh6WxilR5s6d6+vWrUu6DBGRchoq3yBK2BARkdRR8xIRkdRR8xIRkdRJZMCGSKlVm7ezbPUWtu7oY9a0DhbNO4T5h89IuiwRaWK68pJErdq8nStWbGT7rn6mtmfZvqufK1ZsZNXm7UmXJiJNTM1LErVs9RayGaNjQitmwWM2YyxbvSXp0kSkial5SaK27uijPZsZtKw9m2Hbjr6EKhKRNFDzkkTNmtZBLl8YtCyXLzBzWkdCFYlIGqh5SaIWzTuEfMHp27sP9+AxX3AWzTsk6dJEpImpeUmi5h8+g8UL5jBjchsv5vLMmNzG4gVzNNpQRCrSUHlJ3PzDZ6hZiciI6MpLRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR81LRERSR/N5yahYtXk7y1ZvYeuOPmZN62DRvEM0R5eIxEZXXhLZqs3buWLFRrbv6mdqe5btu/q5YsVGVm3ennRpIjJGqXlJZMtWbyGbMTomtGIWPGYzxrLVW5IuTUTGKDUviWzrjj7as5lBy9qzGbbt6EuoIhEZ69S8JLJZ0zrI5QuDluXyBWZO60ioIhEZ69S8JLJF8w4hX3D69u7DPXjMF5xF8w5JujQRGaPUvCSy+YfPYPGCOcyY3MaLuTwzJrexeMEcjTYUkdhoqLyMivmHz1CzEpGG0ZWXiIikjpqXiIikjpqXiIikjpqXiIikjgZsjDPKIBSRsUBXXuOIMghFZKxQ8xpHlEEoImOFmtc4ogxCERkr1LzGEWUQishYoeY1jiiDUETGCjWvcUQZhCIyVmio/DijDEIRGQt05SUiIqmj5iUiIqmj5iUiIqnTkPe8zCwDrAN63P3ssnUTgVuA44E/AOe6+9ONqEuan+KsRGQojbryuhTYNMy6vwJ2uPuhwFeBqxtUkzQ5xVmJyHBib15mNhN4H3DTMJu8H7g5/P424HQzs7jrkuanOCsRGU4jrry+BnwOKA6zvhPYCuDu+4AXgdeXb2RmF5rZOjNb19vbG1Op0kwUZyUiw4m1eZnZ2cB2d3846nO5+w3uPtfd506fPn0UqpNmpzgrERlO3FdeJwMLzOxp4HvAaWb2nbJteoBZAGbWCuxPMHBDxjnFWYnIcGJtXu7+BXef6e6zgfOA+939Q2WbrQA+HH5/TriNx1mXpIPirERkOInEQ5nZYmCdu68AvgUsN7MngRcImpwIoDgrERmapfEiZ+7cub5u3bqkyxARKaeR0g2ihA0REUkdNS8REUkdNS8REUkdzeclNVv60ye46YGneHlvgf0mZLjglIO55F2HJV2WiIxDal5Sk6U/fYJr73+SFoPWluDDwtfe/ySAGpiINJxuG0pNbnrgqbBxtdBiLeFjsFxEpNHUvKQmL+8t0FI2CLjFguUiIo2m5iU12W9ChmLZRwKLHiwXEWk0NS+pyQWnHEzRYV+xSNGL4WOwXESk0TRgQ2oyMChDow1FpBkoHkpEZPQoHqpBdNtQRERSR81LRERSR81LRERSRwM2xpjPfO8RVmz4HYWik2kxFhz9Rr563nFJl1XRqs3bWbZ6C1t39DFrWgeL5h2iObwkkjXb1tC9sZue3T10Tuqka04Xp848NfZ9pXF05TWGfOZ7j3D7+ucohB/IKhSd29c/x2e+90jClQ1v1ebtXLFiI9t39TO1Pcv2Xf1csWIjqzZvT7o0Sak129awZO0SenO9TJkwhd5cL0vWLmHNtjWx7iuNpeY1hqzY8DsAzF79Kl3ejJat3kI2Y3RMaMUseMxmjGWrtyRdmqRU98Zuspks7a3tmBntre1kM1m6N3bHuq80lprXGFIoj8CosrwZbN3RR3t2cEpHezbDth19CVUkadezu4e2TNugZW2ZNnp298S6rzSWmtcYkikPH6yyvBnMmtZBLj84HzGXLzBzWkdCFUnadU7qpL/QP2hZf6Gfzkmdse4rjaXmNYYsOPqNALi/+lW6vBktmncI+YLTt3cf7sFjvuAsmndI0qVJSnXN6SJfyJPbl8Pdye3LkS/k6ZrTFeu+0lhK2Bhj0jzacNuOPmZqtKGMggRHGzbvbY4xRs1LRGT0qHk1iG4biohI6qh5iYhI6qh5iYhI6qh5iYhI6ijbcIyJkhNYbV9lEIpIs9CV1xgSJSew2r7KIBSRZqLmNYZEyQmstq8yCEWkmah5jSFRcgKr7asMQhFpJmpeY0iUnMBq+yqDUESaiZrXGBIlJ7DavsogFJFmouY1hsw/fAaLF8xhxuQ2XszlmTG5jcUL5tQ0IrDavlGeW0RktCnbUERk9CjbsEF05SUiIqmj5iUiIqmj5iUiIqmjeKgI4opLUkyTpEnEyRtF6qIrrzrFFZekmCZJkzXb1rBk7RJ6c71MmTCF3lwvS9YuYc22NUmXJmOcmled4opLUkyTpEn3xm6ymSztre2YGe2t7WQzWbo3diddmoxxal51iisuSTFNkiY9u3toy7QNWtaWaaNnd09CFcl4oeZVp7jikhTTJGnSOamT/kL/oGX9hX46J3UmVJGMF2pedYorLkkxTZImXXO6yBfy5PblcHdy+3LkC3m65nQlXZqMcUrYiGBg1N+2HX3MjGG04XDPG9dxReqh0YaDKGGjQdS8RERGj5pXg+i2oYiIpI6al4iIpI6al4iIpI6al4iIpE6s2YZm1gasBiaGx7rN3b9Utk0X8BVg4FON33D3m+KsK+2W/vQJbnrgKV7eW2C/CRkuOOVgLnnXYVXXQXJ5jCIioynuYN49wGnuvtvMssADZnaXuz9Utt2t7v7JmGsZE5b+9Amuvf9JWgxaW4IPKF97/5OvrB9u3SXvOuyVXMRsxgblIi6GSI0mrucVERlOrLcNPbA7/DEbfqVvbH4TuemBp8Lm1EKLtYSPwfJK6yC5PEYRkdEW+3teZpYxs/XAdmClu68dYrMPmNkGM7vNzGYN8zwXmtk6M1vX29sbZ8lN7eW9BVrKPknSYsHySusguTxGEZHRFnvzcveCux8LzATeYWZvK9vkx8Bsdz8aWAncPMzz3ODuc9197vTp02OtuZntNyFDsezatejB8krrILk8RhGR0daw0YbuvhP4GXBm2fI/uPue8MebgOMbVVMaXXDKwRQd9hWLFL0YPgbLK62D5PIYRURGW9yjDacDeXffaWbtwBnA1WXbHOjuz4U/LgA2xVlT2g2MHKw0onC4dfMPn8FiGPVcxLieV0RkOLFmG5rZ0QS3ATMEV3nfd/fFZrYYWOfuK8zsnwia1j7gBeAid99c6XmVbSgiTUrZhg2iYF4RkdGj5tUgStgQEZHUUfMSEZHUGdGADTM7CZhdup+73zLKNaVGlEikajFO9e5braakak7EEyvhwWth5zMw9SA46VI47Iyadq02waImYBRJVs3veZnZcuCPgPXAwId63N0viae04TXDe16lkUjt2Qy5fIF8wVm8YE7VZlAa8dRiwWexig6XnnZo1WZQad+jZ06tWFNSNSfiiZVw12XQMgGy7ZDPQXEvnHVN1Qa2ZtsalqxdQjaTpS3TRn+hn3whz+UnXM6pM0+tul7GNb3n1SAjuW04FzjZ3T/h7p8KvxreuJpFlEikajFO9e5braakak7Eg9cGjWtCB5gFjy0TguVVdG/sJpvJ0t7ajpnR3tpONpOle2N3TetFJH4jaV6PAW+Mq5C0iRKJVC3Gqd59q9WUVM2J2PlMcMVVKtsOO5+tumvP7h7aMm2DlrVl2ujZ3VPTehGJX9XmZWY/NrMVwAHA42Z2j5mtGPiKv8TmFCUSqVqMU737VqspqZoTMfWg4FZhqXwOpr656q6dkzrpL/QPWtZf6KdzUmdN60UkfrVceV0D/DNwJfDfgCXhzwNf41KUSKRqMU717lutpqRqTsRJlwbvce3tA/fgsbg3WF5F15wu8oU8uX053J3cvhz5Qp6uOV01rReR+I1kwMbV7v431ZY1QjMM2IBXR+7VE4kU92jD4WpKquZEvDLa8NngikujDSV+GrDRICNpXo+4+3FlyzaEafAN1SzNS0SkjJpXg1T9nJeZXQR8AjjEzDaUrJoM/CKuwkRERIZTy4eU/w24C/gn4PMly3e5+wuxVCUiIlJBLc0rA7wEXFy+wsxepwYmIiKNVkvzehhwgnu5bwZ2hN9PBZ4FmnS4mYiIjFVVm5e7HwxgZjcCt7v7neHPZxEMnZcGq5RPGCW7UGq3Zu1SujfdQk9xD50tE+k64nxOPSH+wJnr11/P8k3L6cv30ZHtYOERC7no2ItiP65IsxnJaMNH3f2oassaYTyPNqyUTwjUnV0otVuzdilLNt5IFmizFvq9SB64fM7HYm1g16+/nmUblmFmZMhQoIC7s+joRWpgzUOjDRtkJPFQvzWzL5rZ7PDrfwG/jaswGVqlfMIo2YVSu+5Nt5AF2i2DYbRbhmy4PE7LNy3HzGi11kGPyzctj/W4Is1oJM3rL4HpwO3h14xwmTRQpXzCKNmFUrue4h7abPAfnTZroae4J9bj9uX7yDD495shQ19ev18ZPWZ2p5lNTbqOamqezyscVVg9W0diNWtaB9t39dMx4dVfXWk+YaV1Mjo6WybSW9xDu73aSPq9SGfLxFiP25HtILcvR2vJH9sCBTqy+v3K6HH39yZdQy1qCeb9Wvj449JA3vEezJuUSvmEUbILpXZdR5xPHsh5AcfJeYF8uDxOC49YiLuzz/cNelx4xMJYjyvNx8z2M7M7zOzXZvaYmZ1rZk+b2ZfN7FEz+y8zOzTcdrqZ/dDMfhl+nRwun2Rm/xpuv8HMPhAuf9rMDgi//1D4XOvNbJmZZcKv7vC4j5rZZ5I4B7VceQ3cUL8mzkKkNvMPn8FiGDafsNI6GR2nnnAJl0PDRxsODMrQaEMBzgR+6+7vAzCz/YGrgRfd/SgzOx/4GnA2cC3wVXd/wMzeDNwDHAH87cD24XNMKz2AmR0BnEswj2PezP4F+CCwEeh097eF202N+8UOZSSjDU8HHnT3XNWNYzaeRxuKSFNryGhDMzsMuBe4FfiJu68xs6eB09x9i5llgd+5++vNbDuDB9dNB94K/Bw4z93/b9lzP00w+fB5wOXA9nBVO/DvBM1wHXAncAdwr7sXY3mhFdT8nhdwPnC9mb0ArAFWAw+4+45YKhMRkSG5+xNmdhzwXuAfzOy+gVWlm4WPLcCJ7j5oEjqzqn3WgJvd/QuvWWF2DPAe4OPAXwAfHfGLiKjm0Ybu/mF3Pwz478BW4DqgN67CRERkaGb2JqDP3b8DfAUYmPHj3JLH/wy/vxf4VMm+x4bfrqQk9q/8tiFwH3COmc0I17/OzA4K3w9rcfcfAl8sOXZD1XzlZWYfAk4FjgKeB75BcAUmIiKNdRTwFTMrAnngIuA2YFo4+8ceXv0o0yXAdeHyVoK7Zh8H/iFc/hhQAP4O+NHAAdz9cTP7InCvmbWEx7kYyAH/Gi4DeM2VWSOM5D2v54H/B3wT+Jm7Px1jXRXpPS8RaVKJJWwMvFfl7s8nVUMjjeRzXgeY2RxgHvCPZvYW4DfunupxulGyACvtG+esw8ovLPHKbMnPwNSDRjRbclyq5R7GNQtzpeeteswI51GzSksSRnLlNQU4GfgTgtuHBwAPufuH4ytvaKN15VUpJ7BaM6i074ZtO7n2/idpMWgxKHrwdelph0ZuYFFqHnOeWAl3XQYtEyDbDvkcFPfCWdck1sCq5R6u2baGJWuXkM1kacu00V/oJ1/Ic/kJl0f6C7/S8wKVjxnhPMb1elJM2YYNMpJ4qAeAPwM2AOe6+1uTaFyjKUoWYKV9b3rgKVoMWltaaLGW8BFueuCpRGsecx68NvgLd0IHmAWPLROC5QmplnvYvbGbbCZLe2s7ZkZ7azvZTJbujd3RjlvheaseM8J5jOv1iFQzktuGR1dab2Zfd/dPVdqm2Wzd0cfU9uygZbVmAVba9+W9BVrL/lnQYvDy3kKiNY85O5+BtrIBUtl22PlsMvUQ5B5OqZB72LO7hykTpgxen2mjZ3dPtONWed6Kx4xwHuN6PSLVjOTKq5qTR/G5GmLWtA5y+cENpdYswEr77jchQ7HsbmzRYb8Jg0NVG13zmDP1oOAWV6l8Dqa+OZl6CHIP+8s+r1mae9g5qZP+Qv/g9YV+Oid1RjtuheeteswI5zGu1yNSzWg2r9SJkgVYad8LTjmYosO+YpGiF8NHuOCU6JNOK7+wxEmXBu/N7O0D9+CxuDdYnpBquYddc7rIF/Lk9uVwd3L7cuQLebrmdEU7boXnrXrMCOcxrtcjUk3NAzaqPpHZI+7ekA+rjeZQ+YGRe/VkAVbatxGjDZVfSMkouWeDKwWNNow42nDk51GjDQdJ5YANM+siiHlKzRyNo9m8fuXubx+VJ6tCn/MSkSZVV/Oa/fk7zgQ+CxwMPAV85emr3nf3aBZWiZmtAi5z99T8xTqatw2TG+IlIpJSYeO6DjgQeCF8vC5cXrdhpk053sx+bmYPm9k9ZnagmZ1DEMT73XDqk3YzO93MfhVOefJtM5sYPudVZvZ4OIXKNeGyPzOzteH2PzWzN0Q6ITWqOtrQzH7M4LDHQdx9QfjYPXpliYiMG58liHMaGDLcV7I8ytXXUNOm3AW83917zexc4B/d/aNm9knCKy8zawO6gdPDAOBbgIvMbDnw58Dh7u4lU6E8QBD862Z2AfA54H9GqLsmtQyV1zxeIiLxOZjgiqtUX7g8ikeBfzazq4GfADuAtwErw0T5DPDcEPu9FXjK3Z8If76ZINPwG0A/8C0z+0n4nAAzgVvN7EBgAsFtz9hVbV7u/vNGFJJGaYyWkvjFOYBhzZ2X0L3tXnpajM6i0zXz3Zz63qVV97v+vs+yfOs99OF0YCyc9R4uOv0ro1KTRPYUwa3C0g9rdhCxCZRPmwLcD2x093fW+Xz7zOwdwOnAOcAngdOArwP/291XmNl84Mooddeq5ve8zOwtZnZbeL9zy8BXnMU1s4GYpu27+pnanmX7rn6uWLGRVZu3R9p36U+f4Nr7nySXDz7onMsXuPb+J1n60yeqPq8kbyAuqTfXy5QJU+jN9bJk7RLWbIs+AcOaOy9hyW9X0mswpej0Giz57UrW3Fl5Bufr7/ssy7beTQ6nFcjhLNt6N9ff99nINcmo+AowkaBhET5ODJfXbYhpU04AppvZO8P12TCvFmAXMDn8/jfAbDM7NPx5IfBzM5sE7O/udwKfAY4J1+8PDHwqvWGpSyMZsPGvwPXAPuBPgVuA78RRVBqkMVpK4hdnXFL3tnvJutNOkOLUDmTd6d52b8X9lm+9ByO4zVL6uHzrPZFrkujCUYUXE9zCe134ePEojDY8CvgvM1sPfAm4guCK6Woz+zWwHjgp3LYb+Ga4rQEfAX5gZo8CRYLZRCYDPwmnVnkA+Otw3yvDbR8mmC6rIUYyk3K7u99nZubuzwBXhsVeEVNtTS2N0VISvzjjknpajClFHzQYu82D5ZX0hVdcpTLhcmkOYaMa1aHx7n4PMNS/UOYNse0PgR+WLLoPKP/o03PAO4bY9z+A/6i/0vqM5MprTzj52P81s0+a2Z8Dk2Kqq+mlMVpK4hdnXFJn0ekv61P9FiyvpAOj/J8+hXC5SFqNpHldSnAv9hLgeIL7oKlOlY8ijdFSEr8445K6Zr6bvBk5ghSnHJA3o2vmuyvut3DWe3CC+/2ljwtnvSdyTSJJGXHCRjivl7v7rnhKqq5ZEjbSGC0l8dNow3FNl7MNMpLJKOcSDNoYGJHyIvBRd384ptqG1SzNS0SkjJpXg4xkwMa3gU+4+xoAMzuFoJlVnOdLRERktI3kPa/CQOMCcPcHCG6fi4iINNRImtfPzWyZmc03sz8xs38BVpnZceGnuEVEZIwys8Vm9q469psfxkmNqpHcNhz4NPWXypa/nWDw0mnlO4QBj6sJPi3eCtzm7l8q22YiwQeejwf+AJzr7k+PoC4RkXS7cv/XTInClS82bEqUARaEHpp72XTggLs35DO9Ztbq7lXv6tXcvNz9T+uoYw9wmrvvNrMs8ICZ3eXuD5Vs81fADnc/1MzOA64Gzq3jWMOqlCMYJZ8wikojCqvVlFTNr05Y+EwwdfxIJn5cdTU8dB3s2Q0TJ8GJF8P8v4n8vJUmfqw26q/apJFx+fzqz3P303dT8AIZy3Dm7DO5at5Vwcoq5yLK661X1fOY0GSUkSbfbCZB47qO4O/LV6ZE4cr9L663gZnZVcBWd78u/PlKYDfBgJK/ILiguN3dv2Rmswk+zLyW4CLivWb2dwTTpDjwbXf/qpl1Az9x99vM7I8JpsHaL6z7dCBPkMI0l+Atpb9295+V1fU6gvEThxBkOV7o7hvC+v4oXP4s8JdVX+MIRhu+AVgCvMndzzKzI4F3uvu3aty/gyBS5CJ3X1uy/B7gSnf/TzNrBX4HTPcKhY1ktOFAjmA2Y7RnM+TyBfIFZ/GCINJruHVxNoOB/MIWC9Izih58XXraoRw9c2rFmiq9nlgb2BMr4a7LoGUCZNshnwumij/rmuqNZtXVsPrLgIFlwAuAw7zPwZuOq/t516xdypKNN5IF2qyFfi+SBy6f8zHofDtL1i4hm8nSlmmjv9BPvpDn8hMu59SZp1bcN84G9vnVn+eOp+54zfL3Hfw+rnrj6RXPRZTXW6+BrMZhz2OV9XGpdFwgkZpCIx9teOX+9zF0MO9zXPni6XUVYfZ24Gvu/ifhz48TXBicDCwK61wBfJmgWWwBTnL3h8zseOAqdz8j3Hequ+8caF7hfpsJ7pL9Mvz4VB/BZ4HnhFOsHA7cCxwGnEgw3crZZvZ14Hl3/zszO40g0PfYsHn9GXCKu+dqeY0jec+rm6A7vyn8+Qng09V2MrNMmJe1HVhZ2rhCncBWCFKLCYbgv34EdVVUKUcwSj5hFJXyC6vVlFTNPHht8JfqhI4gWG9CR/DzgzXMQfrQdYBBpjXo1pkwXe+h6yI9b/emW8gC7ZbBMNotQzZcXi1jsNK+cbr76eAf0lby3yvLq5yLKK+3XlXPY4xZjvXWlVRNERzM4MYFEadEcfdfATPM7E1mdgzBdChHAe8GfgU8AhwOvCXc5ZmSO2JbgEPM7OtmdibwUtnTvxV4zt1/GR7rpfDv7lMI827dfTPwDEHzKnUKsDzc5n7g9WHzA1hRa+OCkTWvA9z9+wQhjQONpmrgnrsX3P1Ygjlf3mFmbxvBMV9hZhea2TozW9fb21vzflt39NGeHRytNJAjWGldnF7eW6A8jm4gv7BaTUnVzM5ngquBUtl22Pls9X337A6uuEpZJlge4Xl7intos8H/C7dZCz3FPfTs7qEt0zZ4XUnGYKV941Twof/IFLxQ9VxEeb31qnoeYzpulLqSqimCp3g1UX5A5ClRgB8QBPGeC9xKcLX1T+5+bPh1aMmds5cHdnL3HQRjHFYBHwduilhHrV6uvsmrRtK8Xjaz1xPOqmxmJxJcJdXE3XcCPyOY3bNUDzArfM5Wgnj9Pwyx/w3uPtfd506fPr3moivlCEbJJ4yiUn5htZqSqpmpBwW3sUrlczD1zdX3nTgpvFVYwgvB8gjP29kykf6y95X7vUhny8SqGYOV9o1TpryJly6vci6ivN56VT2PMWY51ltXUjVFEMuUKAQN6zyCBvYDgjtnHw2nNsHMOs3sNe81mNkBQEsY1vtFoHw0+W+AA8P3vTCzyeHf3WuAD4bLDgPeHG5bqnSb+QS3EMuv7Goykub11wT3Ov/IzH5BMELwU5V2MLPpA1NFm1k7cAbBvdJSK3g1I/Ec4P5K73eNVKUcwSj5hFFUyi+sVlNSNXPSpcH7L3v7gmC9vX3BzyddWn3fEy8GHAr7gi5dCNP1Trw40vN2HXE+eSDnBRwn5wXy4fJqGYOV9o3TmbODf7t5yX+vLK9yLqK83npVPY8xZjnWW1dSNdUtGJTxmilRoo42dPeNBIlIPe7+nLvfC/wb8J/hVCe38WpiUqlOgo9BrSe4DfiFsufdS3A19/VwapWVQBvwL0BL+Ny3Al3uXn4r40rg+HBalauIkI87kgEb/4Ogc88CPkAwsdnfuvsjFfY5mmAK6QxBo/y+uy82s8XAunDmzTaCe6BvJxhpc567V3wDZ6TxUJVyBKPkE0ZRy2jD4WpKquZXR8I9G1wNjPpow5E/79gdbTj0udBow9rqSnC0oeKhGmQkzWuDux8dxkL9PXANcIW7nxBngUNRtqGINCk1rwYZUTxU+Pg+4EZ3vwOYMPoliYiIVDaS5tVjZssI7nXeGSZjjGR/ERGRUTGS5vMXBO95vSccOfg6gjgTERGRhhpJPFQf8KOSn58jGBWTaolFLY03lWKPosRO1XvMiDVVHBAQ5fXEdS4SkqqYJkmVEc+k3AxGa8BGYlFL402laCmoP3aq3mMedkakmipGIvX11/96okRwNaGkoqMSpgEbDTKu37NKLGppvKkUexQldqreY0asqWL8UJTXE9e5SEgKY5rGtDAq6rY69rtz4PO6Fbapa7qUKEYyJcqYs3VHH1Pbs4OWNSRqabzZ+Qy0TRu87JXYI6+wLq5jRqupZ3cPUyZMGbT6lfihnT31v55qNadMxfMkgxx181GvmRLl0Q8/OqpTorj7bwmCIAapNgWJu7+3huduyHQppcb1lVdiUUvjTaXYoyixU/UeM2JNFeOHoryeuM5FQlIY05SIsHFdR5As/8qUKOHyupjZVWZ2ccnPV5rZZWb2WPhzl5mtMLP7gfvMrMPMvm9mj5vZ7Wa21szmhts+bWYHmNlsM9tkZjea2UYzuzdMTsLMus3snPD7PzazB83s12b2X2F81GwzW2Nmj4RfJ9X72gaM6+aVWNTSeFMp9ihK7FS9x4xYU8X4oSivJ65zkZDUxTQl57MEc2IN3PLpC3+OMpr7VoIR4gP+gmC+rlLHAeeE06Z8gmBexSOBvyWY12sobwGuc/c5wE6CtKVXmNmE8NiXuvsxwLuAHMGsIme4+3EEH7daWv9LC4zr5jX/8BksXjCHGZPbeDGXZ8bkNg3WiMNhZwSDDia/Afp3Bo8DgxAqrYvrmBFrOnXmqVx+wuVMb5/OS3tfYnr79FcHIUR5PXGdi4RUPE9SqlFTomwt22ylu78Qfn8K8L1w38eADcM89VPuvj78/mFgdtn64aZLyQI3hrmHPwCOrPe1DRjX73lB0MDUrBpgoCmMdF1cx4xY06kzTx3+L+Eoryeuc5GQiudJBjzF0JNRjtaUKG8kuBoqN6IpSEKlQbsFoH24Dct8Bvg9wVQrLUB/5c2rG9dXXiIiTaBRU6JU8gvC24xmdiTBxJX1GG66lP0JrsiKwEKCsPZI1LxERBIUjip8zZQoUUcblk+JUmXzfwGmm9njwD8AGxnBfI0lx6w0XcqHw2WHU99V3yDj+kPKIiKjLJUfUjazDJB1934z+yPgp8Bbw2bUlMb9e14iIkIH8DMzyxI04E80c+MCNS9pkIoZdxUmqow1Gy/CcSutv3799SzftJy+fB8d2Q4WHrGQi4696NXjJpXl2ISSmkBTBnP3XcDcpOsYCd02lNhVzLh78kFY/WXAwDLgBcBh3udYc+hJ8WXjrbq67uNWej2PPf8YyzYsw8zIkKFAAXdn0dGLggYWV35hCnMR48o+TDhTMZW3DdNIAzYkdhUz7h66DjDItEJL+IjBQ9fFm40X4biV1i/ftBwzo9VaBz0u37Q8OG5SWY5NKK7frzIVxwc1L4ldz+4e2jJtg5a9knG3Z3dw5VPKMrBnd+X9oopw3Err+/J9ZMpGAWfI0JcPP8Kz85ngyqjUaGU5xvG8MYrr9xvr/zfSNNS8JHYVM+4mTgpv2ZXwAkycFG82XoTjVlrfke2gwODnLVCgIxt+hCepLMcmFNfvV5mK44Oal8SuYsbdiRcDDoV9UAwfcTjx4niz8SIct9L6hUcsxN3Z5/sGPS48YmFw3KSyHJtQXL9fZSqODxqwIQ2h0YZDjTZ8NrgyGvXRhqP8vDEag6MNNWCjQdS8RERGj5pXg+i2oYiIpI6al4iIpI6al4iIpI7iocabpCKE6jzu9bd/kOU71tPXYnQUnYXTjuWiP/9urMesRaUBAYomEomfBmyMJ0lFCNV53Otv/yDLXvw15k6GYOY7N2PR/sdUb2AxvtZK8UNAktFEkjwN2GgQ3TYcT5KKEKrzuMt3rMfcCYObgkd3lu9YH9sxa1EpfkjRRCKNoeY1niQVIVTncfta7DXTrWbC5XEdsxaV4ocUTSTSGGpe40lSEUJ1Hrej6GVBS8Gtw45iDbe6Y3ytleKHFE0k0hhqXuNJUhFCdR534bRjcTPC4Kbg0YyF046N7Zi1qBQ/pGgikcbQgI3xJqkIoTqPOzqjDUf/tWq0oQxDAzYaRM1LRGT0qHk1iG4biohI6qh5iYhI6qh5iYhI6qh5iYhI6ijbUEZHlBzBSvsmlE8o6abf7din0YYSXZQcwUr7QiL5hPpLLt0S/t1qtGGD6LahRBclR7DSvgnlE0q66Xc7Pqh5SXRRcgQr7ZtQPqGkm36344Oal0QXJUew0r4J5RNKuul3Oz6oeUl0UXIEK+2bUD6hpJt+t+ODBmzI6IiSI1hp34TyCSXdEvzdasBGg6h5iYiMHjWvBtFtQxERSR01LxERSR01LxERSZ1Y46HMbBZwC/AGgslwb3D3a8u2mQ/8B/BUuOhH7r44zrqaXaQ3m2OMU4p03Arr0zhwYs3apXRvuoWe4h46WybSdcT5nHrCJfEfN4XnSiQOsQ7YMLMDgQPd/REzmww8DPw3d3+8ZJv5wGXufnatzzuWB2xEiraJEtMURbXjVli/pqMtdTFNa9YuZcnGG8kCbdZCvxfJA5fP+VisDUyRVqmgARsNEuttQ3d/zt0fCb/fBWwC9EnBCiJF28QYpxTpuBXWpzHKp3vTLWSBdstgGO2WIRsuj/W4KTxXInFp2HteZjYbeDuwdojV7zSzX5vZXWY2Z5j9LzSzdWa2rre3N85SExUp2ibGOKVIx62wPo1RPj3FPbTZ4D86bdZCT3FPvMdN4bkSiUtDmpeZTQJ+CHza3V8qW/0IcJC7HwN8Hfg/Qz2Hu9/g7nPdfe706dNjrTdJkaJtYoxTinTcCuvTGOXT2TKRfi8OWtbvRTpbJsZ73BSeK5G4xN68zCxL0Li+6+4/Kl/v7i+5++7w+zuBrJkdEHddzSpStE2McUqRjlthfRqjfLqOOJ88kPMCjpPzAvlweazHTeG5EolL3AM2DLgZeMHdPz3MNm8Efu/ubmbvAG4juBIbtrCxPGADRmu04ejHKUU6boX1aRxBp9GGMgwN2GiQuJvXKcAa4FFg4D7L5cCbAdz9m2b2SeAiYB+QA/7a3R+s9LxjvXmJSGqpeTVIrJ/zcvcHqPLLdPdvAN+Isw4RERlblLAhIiKpo+YlIiKpo+YlIiKpE+t7XlKnpPIJo1h1NTx0HezZDRMnwYkXw/y/SboqERmjdOXVbAZyAHf9HtqmBY93XRYsb1arrobVXw4+v9WSDR5XfzlYLiISAzWvZpNUPmEUD10HGGRaoSV8xMLlIiKjT82r2SSVTxjFnt1gmcHLLBMsFxGJgZpXs0kqnzCKiZPAC4OXeSFYLiISAzWvZpNUPmEUJ14MOBT2QTF8xMPlIiKjT82r2Rx2RjCJ4+Q3QP/O4DHuySSjmv83MO9zwftzxXzwOO9zGm0oIrGJNdswLso2FJEmpWzDBtGVl4iIpI6al4iIpI6al4iIpI7ioSpYtXk7y1ZvYeuOPmZN62DRvEOYf/iMpMtqzvioZqwpTuPt9Yo0GV15DWPV5u1csWIj23f1M7U9y/Zd/VyxYiOrNm9PtrBmjI9qxpriNN5er0gTUvMaxrLVW8hmjI4JrZgFj9mMsWz1lmQLa8b4qGasKU7j7fWKNCE1r2Fs3dFHe3Zw5FF7NsO2HX0JVRRqxvioZqwpTuPt9Yo0ITWvYcya1kEuPzjyKJcvMHNaR0IVhZoxPqoZa4rTeHu9Ik1IzWsYi+YdQr7g9O3dh3vwmC84i+YdkmxhzRgf1Yw1xWm8vV6RJqTmNYz5h89g8YI5zJjcxou5PDMmt7F4wZzkRxs2Y3xUM9YUp/H2ekWakOKhRERGj+KhGkRXXiIikjpqXiIikjpqXiIikjpqXiIikjrKNpQxa83apXRvuoWe4h46WybSdcT5nHrCJbXtrOxCkaamKy8Zk9asXcqSjTfSW9zDFGuht7iHJRtvZM3apdV3VnahSNNT85IxqXvTLWSBdstgGO2WIRsur0rZhSJNT81LxqSe4h7abPD/3m3WQk9xT/WdlV0o0vTUvGRM6myZSL8XBy3r9yKdLROr76zsQpGmp+YlY1LXEeeTB3JewHFyXiAfLq9K2YUiTU/NS8akU0+4hMvnfIzpLRN5yYtMb5nI5XM+VttoQ2UXijQ9ZRuKiIweZRs2iK68REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkddS8REQkdWJtXmY2y8x+ZmaPm9lGM3vNVLQWWGpmT5rZBjM7Ls6axoQnVkL32fC1o4LHJ1YmXZGISEPFfeW1D/if7n4kcCJwsZkdWbbNWcBbwq8LgetjrindnlgJd10Gu34PbdOCx7suUwMTkXEl1ubl7s+5+yPh97uATUBn2WbvB27xwEPAVDM7MM66Uu3Ba6FlAkzoALPgsWVCsFxEZJxo2HteZjYbeDuwtmxVJ7C15OdtvLbBYWYXmtk6M1vX29sbW51Nb+czkG0fvCzbDjufTaYeEZEENKR5mdkk4IfAp939pXqew91vcPe57j53+vTpo1tgmkw9CPK5wcvyOZj65mTqERFJQOzNy8yyBI3ru+7+oyE26QFmlfw8M1wmQznpUijuhb194B48FvcGy0VExom4Rxsa8C1gk7v/72E2WwGcH446PBF40d2fi7OuVDvsDDjrGpj8BujfGTyedU2wXERknGiN+flPBhYCj5rZ+nDZ5cCbAdz9m8CdwHuBJ4E+4CMx15R+h52hZiUi41qszcvdHwCsyjYOXBxnHSIiMrYoYUNERFJHzUtERFJHzUtERFJHzUtERFJHzUtERFJHzUtERFJHzUtERFJHzUtERFLHgs8Ip4uZ9QLPNPCQBwDPN/B4tWrGulRT7ZqxLtVUu6Hqet7dz0yimPEmlc2r0cxsnbvPTbqOcs1Yl2qqXTPWpZpq16x1jRe6bSgiIqmj5iUiIqmj5lWbG5IuYBjNWJdqql0z1qWaatesdY0Les9LRERSR1deIiKSOmpeIiKSOmpeZcwsY2a/MrOfDLGuy8x6zWx9+HVBg2p62sweDY+5boj1ZmZLzexJM9tgZsc1QU3zzezFknN1RQNqmmpmt5nZZjPbZGbvLFufxHmqVlMS5+mtJcdbb2Yvmdmny7Zp6LmqsaYkztVnzGyjmT1mZv9uZm1l6yea2a3heVprZrPjrkkCsc6knFKXApuAKcOsv9XdP9nAegb8qbsP90HNs4C3hF8nANeHj0nWBLDG3c9uQB0DrgXudvdzzGwC0FG2PonzVK0maPB5cvffAMdC8I81oAe4vWyzhp6rGmuCBp4rM+sELgGOdPecmX0fOA/oLtnsr4Ad7n6omZ0HXA2c24j6xjtdeZUws5nA+4Cbkq5lhN4P3OKBh4CpZnZg0kU1kpntD8wDvgXg7nvdfWfZZg09TzXWlLTTgf/n7uWJNUn+PzVcTUloBdrNrJXgHx6/LVv/fuDm8PvbgNPNzBpY37il5jXY14DPAcUK23wgvI1ym5nNakxZOHCvmT1sZhcOsb4T2Fry87ZwWZI1AbzTzH5tZneZ2ZyY6zkY6AX+Nbzte5OZ7Ve2TaPPUy01QWPPU7nzgH8fYnkS/08NGK4maOC5cvce4BrgWeA54EV3v7dss1fOk7vvA14EXh9nXRJQ8wqZ2dnAdnd/uMJmPwZmu/vRwEpe/RdX3E5x9+MIbuVcbGbzGnTcSqrV9AhwkLsfA3wd+D8x19MKHAdc7+5vB14GPh/zMauppaZGn6dXhLcxFwA/aNQxq6lSU0PPlZlNI7iyOhh4E7CfmX0ozmNK7dS8XnUysMDMnga+B5xmZt8p3cDd/+Due8IfbwKOb0Rh4b8AcfftBO8DvKNskx6g9CpwZrgssZrc/SV33x1+fyeQNbMDYixpG7DN3deGP99G0DhKNfo8Va0pgfNU6izgEXf//RDrGv7/VLWaEjhX7wKecvded88DPwJOKtvmlfMU3lrcH/hDjDVJSM0r5O5fcPeZ7j6b4LbF/e4+6F9ZZff8FxAM7IiVme1nZpMHvgfeDTxWttkK4PxwhNiJBLc3nkuyJjN748C9fzN7B8H/a7H9oXb33wFbzeyt4aLTgcfLNmvoeaqlpkafpzJ/yfC35xp6rmqpKYFz9Sxwopl1hMc9ndf+mV8BfDj8/hyCvzeU/NAAGm1YhZktBta5+wrgEjNbAOwDXgC6GlDCG4Dbwz+zrcC/ufvdZvZxAHf/JnAn8F7gSaAP+EgT1HQOcJGZ7QNywHkN+EP9KeC74a2nLcBHEj5PtdSUxHka+EfHGcCikmWJnqsaamrouXL3tWZ2G8Htyn3Ar4Abyv5O+Baw3MyeJPg74by46pHBFA8lIiKpo9uGIiKSOmpeIiKSOmpeIiKSOmpeIiKSOmpeIiKSOmpeIiKSOmpeMuZZMJXGa6a4KVnfZWbfiOG4XWb2ppKfn25geobImKbmJRKfLoJMPBEZZWpe0hTCyKk7wsTwx8zsXDM73sx+HibX3zMQz2Vmq8zsWgsmJHwsjArCzN5hZv8ZJrg/WBLLNJI6ppvZD83sl+HXyeHyK83s2+Gxt5jZJSX7/K2Z/cbMHrBgwsLLzOwcYC5BusZ6M2sPN/+UmT1iwUSeh0c+cSLjlJqXNIszgd+6+zHu/jbgboLk8HPc/Xjg28A/lmzf4e7HAp8I1wFsBk4NE9yvAJbUUce1wFfd/Y+BDzB4brfDgfcQhBB/ycyyZjaw3TEEobJzAdz9NmAd8EF3P9bdc+FzPB+m8V8PXFZHfSKCsg2leTwK/LOZXQ38BNgBvA1YGWYoZgjmVBrw7wDuvtrMppjZVGAycLOZvYVgvrFsHXW8CzjSXp1PcIqZTQq/vyOcVWCPmW0nyHg8GfgPd+8H+s3sx1We/0fh48PAf6+jPhFBzUuahLs/YWbHEYTB/gNwP7DR3d853C5D/Pz3wM/c/c/NbDawqo5SWoATw2b0irCZ7SlZVKC+Pz8Dz1Hv/iKCbhtKkwhH5fW5+3eArwAnANPN7J3h+qwNnjn33HD5KQTTdbxIMJfSwJxTXXWWci9BEvxAXcdW2f4XwJ+ZWVt4hXZ2ybpdBFeDIjLK9C8/aRZHAV8xsyKQBy4imIZiqZntT/D/6teAjeH2/Wb2K4Jbgx8Nl32Z4LbhF4E76qzjEuA6M9sQHnM18PHhNnb3X5rZCmAD8HuC258vhqu7gW+aWQ4Y7gpSROqgKVEkdcxsFXCZu69LuhYAM5vk7rvNrIOg2V3o7o8kXZfIWKYrL5HobjCzI4E24GY1LpH46cpLxg0z+whwadniX7j7xUnUIyL1U/MSEZHU0WhDERFJHTUvERFJHTUvERFJHTUvERFJnf8Pi8XozePdy5AAAAAASUVORK5CYII=)
%% Cell type:markdown id: tags:
Since this is a multiclass classification problem, we need to convert the label column to binary columns. This is equivalent to one-hot encoding.
%% Cell type:code id: tags:
``` python
#one-hot encoding
categories = iris.species.unique() # create a vector with the category names,
lencat = len(categories) # store the length of the categories vector to use later
# one hot encoding, create three new columns which are true for that particular species, and false if not
for species in categories: # loop over all the labels in categories
# + concatenates two strings
iris['species_'+species] = pd.Series(iris['species']==species).astype(int)
iris.describe() # show the statistics of the nunmerical features of the dataframe
```
%%%% Output: execute_result
sepal_length sepal_width petal_length petal_width species_setosa \
count 150.000000 150.000000 150.000000 150.000000 150.000000
mean 5.843333 3.057333 3.758000 1.199333 0.333333
std 0.828066 0.435866 1.765298 0.762238 0.472984
min 4.300000 2.000000 1.000000 0.100000 0.000000
25% 5.100000 2.800000 1.600000 0.300000 0.000000
50% 5.800000 3.000000 4.350000 1.300000 0.000000
75% 6.400000 3.300000 5.100000 1.800000 1.000000
max 7.900000 4.400000 6.900000 2.500000 1.000000
species_versicolor species_virginica
count 150.000000 150.000000
mean 0.333333 0.333333
std 0.472984 0.472984
min 0.000000 0.000000
25% 0.000000 0.000000
50% 0.000000 0.000000
75% 1.000000 1.000000
max 1.000000 1.000000
%% Cell type:markdown id: tags:
Now we create a features table XX which comprises the first four columns of the iris data set, and a target series Y which is the binary column we have just created to identify whether or not the species is versicolor. Here we will used the normalized version of the features matrix X.
%% Cell type:code id: tags:
``` python
# the table XX will have the first 3 columns (0 to 3) of iris, and Y will have the binary classification column
XX, Y = iris[iris.columns[:4]], iris.species_versicolor
X = (XX-XX.mean())/XX.std() # create a new feature matrix for analysis which has mean 0 and standard deviation 1
```
%% Cell type:markdown id: tags:
Then we split these data sets into training and testing sets.
%% Cell type:code id: tags:
``` python
from sklearn.model_selection import train_test_split # import the splitting method from sklearn
# split the data into 80% training and 20% testing, random_state=0 ensures that the results are repeatable
X_train,X_test,y_train,y_test=train_test_split(X,Y,train_size=0.8,random_state=0)
```
%% Cell type:markdown id: tags:
## Logistic Regression
%% Cell type:markdown id: tags:
We can now introduce the sigmoid or logistic function. This is a defined as
$$ f(x) = \displaystyle{\frac{1}{1+e^{-x}}}, $$
and is a smooth, one-to-one (every x value gives a unique y value) function with the domain $(-\infty,\infty)$ and range $(0,1)$. If x is negative then the logistic function is less than 0.5 but greater than 0, while if x is positive the logistic function is greater than 0.5 but less than 1.
%% Cell type:code id: tags:
``` python
xa = np.linspace(-6,6) # create a linear array of x values between -6 and 6
plt.plot(xa,1/(1+np.exp(-xa))); # plot the logistic function
plt.xlabel("x") # xlabel
plt.ylabel("y") # ylabel
plt.title("Logistic function y = 1/(1+exp(-x))"); # title
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAp/0lEQVR4nO3deXxU9bnH8c+THQJhX8Mui4CKS8StVauooC14u8q1Wru59Gq12lqtXaz13trd1mrVutZq0aJVWnHDrdYWBRRBwmIIW1gDhCVkzzz3j3PQMQYIIcOZyXzfr9e8MnPmzDnPmZnM95zfWX7m7oiISPrKiLoAERGJloJARCTNKQhERNKcgkBEJM0pCERE0pyCQEQkzSkIBDO708x+0IrXDTKzSjPLbMVrbzazzWa2YX9feyBau6ypysz+YmbnRl1HS5nZT83sqla+9nEzmxT3uI+ZLTaz3DYrsJ0ynUeQWsxsJfA1d5+VqvM2s0HAUmCwu29qi9r2MJ+LCOr9WKLmcTCZ2SeAHwJHAxXuPqSZcU4AfuXuJ5rZEcA0YKy7u5n1A+4CioB+wFB3X3mw6t8XM+sFzAeGu3t1K14/HviDux8TN+wOYLG739ZmhbZD2iKQKAwCtiQyBNqpXcB9wHf2Ms45wMzw/iXAw/7B2l4MeBb4zL5mZGYXmdkDrS+1VS4CZrYmBADc/U2gwMyK4gY/TPA+yF4oCNoJM8s1s1vNbF14uzV+k9jMrjWz9eFzXzMzN7Ph4XMPmNnN4f2eZvYPM9tmZlvN7DUzyzCzhwh+wP8eNgdda2ZDwulkha/tbmb3h/OoMLMnm6lzAvAC0D+czgNmdqqZlTUZb2U4LmZ2o5k9ZmZ/MrOdZrYo/p/dzAaa2RNmVm5mW8zs92Y2GrgTOCGcz7amyxo+/rqZlYTLOsPM+sc952Z2qZm9F74ft5uZNbNMfc2sysx6xA07Oqwne38+x71x9zfd/SGgdC+jnc0HQTAJeDXu9Rvd/Q5gzoHUYWZdzOze8Pu0NmzmyzSzHDObb2ZXhONlmtnrZvbD8PGNZjbdzB4NP8e3zGxc3KQ/VG8z8z0k/JyODh/3D9/jU+NGe4UgDHd7AxhmZoMPZJnbOwVB+3EDcDxwJDAOGA98H8DMJgJXAxOA4cCpe5nONUAZ0AvoA3wPcHe/AFgNfMrdO7n7z5t57UNAR2As0Bv4TdMRwmalScC6cDoXtXD5JhM0c3QFZgC/D5ctE/gHsAoYAhQC09x9MXAp8J9wPl2bTtDMTgN+CnyeoKlkVTiPeJ8EjgWOCMc7q5ll2kDwA/T5uMEXhHXUNzPf/w6DZU+3QS16Rz463X4En9nbZpYPDCVogmtrDwANBN+lo4AzCZrg6oAvAjeFQXwdkAn8b9xrpwB/BboDjwBPxoXl4Xur192XA98F/mxmHYH7gQfd/ZW40RYTfP93v6YBKIkfJh+lIGg/zgducvdN7l4O/JjgxwiCH6j73X2Ru1cBN+5lOvUEP4qD3b3e3V+La1rYo/BHaBJwqbtXhK/d49pdK/zL3We6eyNB4Oz+xx4P9Ae+4+673L3G3f/VwmmeD9zn7m+5ey1wPcEWxJC4cW5x923uvhp4mSBom/MgwY/g7nCaGtb5Ee7+iLt33cttdQvrb+ps4Nnw8+oaDtvZymk1y8z6hPO5Kny/NxEE/nkA7v4ucDPwJPBt4ILwM9ttnrtPDwPy10AewQoMYc17rdfd/0jww/4Gwff0hiaj7OSDZd/bMImjIGg/+hOs0e62Khy2+7k1cc/F32/qFwT/aM+bWamZXdfC+Q8Etrp7RQvH31/xRxdVAXlhk9RAYFW45re/PvSeuXslsIVgq2JP8+20h2k9BYwxs6HAGcD2sM36YIpvFtoW/u3c0heb2R27t0qAO4D4LZcF4WiDgWxgfdy4dxFsAe72YDjeTHd/r8ls3v/uuXuMYOtz9/e0Ir5eM3smbNarNLPz46bxR+Aw4LYwwON1jlv2vQ2TOAqC9mMdwT/fboPCYQDrgQFxzw3c00Tcfae7X+PuwwiaY642s9N3P72X+a8BuptZ1/0tnGAnaMfdD8I16l4tfO0aYNDu/RRN7GtL5kPvWdic0gNY28J5fzAj9xrgMYKtggvYw9ZAOJ/z437gmrvtd9NQ2LxyCsH+F9x9F7AcGLkfy/CN3VslwDeA+C2XI8LR1gC1QM+45wrcfWzcpO4gaK47y8yaHrH1/nfPzDIIvpe7v6cL4ut190lhs14nd384fE0n4FbgXuBGM+veZPqjgXfi5pFF0IT1DrJHCoLUlG1meXG3LOAvwPfNrJeZ9SQ4zPDP4fiPAV82s9Fh2+oej6M3s0+a2fBwp+h2oJHgaBOAjcCw5l7n7uuBZ4A7zKybmWWb2cktXJ5lBGv454Q/aN8HWnrs95sEQXeLmeWH78dJcfUOMLOcPbz2LwTvy5EW7Fj/P+CNAzik8k8ER75MZi9B4O4Px/3ANXdrtmnIgp32eQRr5BYu6+5l+xiwwN13xL1kJkE4xE8jjw/e29zwcYuFn/PzwK/MrCCs6RAzOyWc/gXAMeH78E3gwfDHe7djzOzT4Xf2KoJQmb2nepvxW2Cuu38NeJrggIB4pxB8D3cbD6x091XIHikIUtNMoDrudiNBu+xcgrWqhcBb4TDc/RngdwRt3CV88I/XdLMaYAQwC6gE/gPc4e4vh8/9lCBstpnZt5t57QUE+xiWAJsI/tH3yd23E6yB3kOwNr6LoMmgJa9tBD5FsNa3OnzdF8KnXwIWARvMbHMzr51FEIqPE4TJIYRt3a3h7q8ThOZbCfrhOZng855JsMVXTfCjDB8+bHS3u4Hzw1DfrZrgs4Xgc2rNoZoXAjlAMUFzznSgX7glcytwobtXuvsjBN/J+IMGniL4fCoIvi+fjtuh/ifgbDPr0NxMzWwKMBG4LBx0NXD07mYjMzsWqGzSJHc+Hw0LaUInlKWh8IiOd4HcVratyx6Y2UsETSr3HOT5FgOfdffiJsMfAR5z9ycPZj3NMbMbCU4W++Jexvk/YJO739qK6T8O3OvuM8PHvQkORz0qbLqTPVAQpAkz+y+CNcaOBDvzYu5+bqRFtTPhGukLwEB3b9OjdfYx3xzgane/5WDNszVaEgQSDTUNpY9LCJprlhO0+1+299Flf5jZgwRNalcdzBAAcPe6ZA8BSW7aIhARSXPaIhARSXPNHXud1Hr27OlDhgyJugwRkZQyb968ze7e7Pk5KRcEQ4YMYe7cuVGXISKSUsxsj4c0q2lIRCTNKQhERNKcgkBEJM0pCERE0pyCQEQkzSUsCMzsPjPbZGbv7uF5M7PfWdBN4AILu58TEZGDK5FbBA8QXClwTyYRXOlyBHAx8IcE1iIiInuQsPMI3P2fTbr8a2oK8KewW73ZZtbVzPqF1zsXEWkX3J3ahhg19Y3U1Id/GxqprY9R1xijriFGbUNj+Dd4XN/oNMTi7jfGqG+McfroPowb2LXNa4zyhLJCPtxlYlk47CNBYGYXE2w1MGhQq/r1FhHZL+5OdX0jFVX1bKuqY1tVPduq6tlRU8/OmnoqaxrYUdPAzpoGdtbUU1XXyK66Bqrj/9Y2Ul3fuO+ZtVCfLnntLghazN3vJuhkg6KiIl0lT0RarTHmbNpZw9qKajbsqKF8Z+0Ht8rg7+bKWiqq6qlriO1xOmbQKSeLznlZdM7LJj83k/ycLHp1yiU/N4sOOZnk52TSITuT3Ozgb152JnnZGeRlZ5KblUFuViY5WRnBLTP4m5uVQXZmBtmZRlZmMDw708jMMD7cx1DbiTII1vLhvnMH0Iq+YkVEmtpRU8+K8l2s2LyL0vJK1lRUs3ZbNWsrqtm4o4aG2IfXJ7MyjJ6dcunVOZc+BXmM6VdA9/wcunbMoWvHbLp1zH7/fkFeNp3zssjPySIjIzE/zAdblEEwA7jczKYBxwHbtX9ARPZHVV0Di9fvpHj9DorX7WB5eSWl5bvYXPlBL6wZBv26dKCwaweOHdKNwm4d6N81uPXrkkfvznl07ZDdbn7UWyNhQWBmfwFOBXqaWRnwI4JOt3H3Owl6yzqboA/dKuDLiapFRFJfbUMjC8q2M29VBe+u3U7x+h2s2LyL3V2qdOmQzcg+nTjt0F4M69WJoT3zGdYzn0E9OpKblRlt8UkukUcNTd3H8w78T6LmLyKpbWdNPW+t3sabK7YwZ0UF88u2vd9mX9i1A2P7FzB5XH/G9u/CmP4F9O+Sl7A29PYuJXYWi0j75+4sL69k1uJNvLh4I/NWVRBzyMwwDutfwIXHD+bYod0pGtyNHp1yoy63XVEQiEhk6htjvLliK7MWb+TFxZtYvbUKgLH9C/jGqcM5flgPjhrUlfxc/VQlkt5dETnoitftYPq8Mp6av5Ytu+rIycrgpEN6cPHJwzh9dG/6dekQdYlpRUEgIgfFlspanpq/junzyihev4PsTGPC6D5MObKQk0f2pGOOfo6iondeRBJq/ppt3P3P5Ty/aCMNMefwwi78ePJYJo/rT7f8nKjLExQEIpIA7s6ry8q589XlzC7dSkFeFhedOITPFg3g0L4FUZcnTSgIRKTNNDTGeHrheu58tZTF63fQtyCPG84ezdTjBtFJO3yTlj4ZETlg7s6z727glmeXsGpLFYf0yufnnz2Cc48sJCdL/V8lOwWBiByQReu2c9Pfi3ljxVZG9enMXRccwxmj+6T1JRtSjYJARFplc2Utv3p+KdPmrKFrh2xuPvcwzjt2IFmZ2gJINQoCEdkv9Y0x7n99Bbe9WEJ1fSNfOWko3zx9BF06ZEddmrSSgkBEWqy0vJJvPTqfd8q2c9qhvbnhnNEc0qtT1GXJAVIQiMg+uTuPvLmam/+xmNzsDO44/2jOPrxf1GVJG1EQiMhele+s5brHF/Dikk18fERPfvm5cfQpyIu6LGlDCgIR2aNZxRv57uML2FnbwI8+NYYvnTBERwO1QwoCEfmIxphzyzOL+eNrKxjTr4C/nHckI/t0jrosSRAFgYh8yK7aBq6cNp9Zizdy4QmDueGc0erhq51TEIjI+9Zvr+arD8xlyYYd/HjyWL504pCoS5KDQEEgIgAsLNvO1/40h121jdx70bF8YlTvqEuSg0RBICI8t2gDV02bT/f8HB6/7DhG9dX+gHSiIBBJc/f9awU/ebqYcQO6cveFx9C7sw4NTTcKApE0ds9rpdz89GImju3LrecdSV62dgqnIwWBSJp64PUV3Pz0YiYd1pffTT2KbF0sLm3pkxdJQw/9ZyU3/r2Ys8b2UQiIgkAk3Tz8xip+8NQiJozuw21Tj1YIiIJAJJ1Me3M1N/ztXU47tDe3n3+Ueg8TQEEgkjb+OncN1/9tIaeM7MUd5x+ts4XlfQoCkTTw6rJyvvv4Aj42vCd3XXCMjg6SD1EQiLRz723cyeUPv8XIPp2584sKAfkoBYFIO7Z1Vx1ffXAuudkZ3HvRseTn6ohx+aiEBoGZTTSzpWZWYmbXNfP8IDN72czeNrMFZnZ2IusRSSe1DY1c+tA8Nuyo4e4Liyjs2iHqkiRJJSwIzCwTuB2YBIwBpprZmCajfR94zN2PAs4D7khUPSLpxN35/t/e5c2VW/nFZ4/g6EHdoi5JklgitwjGAyXuXurudcA0YEqTcRwoCO93AdYlsB6RtHH3P0v567wyvnnacKYcWRh1OZLkEhkEhcCauMdl4bB4NwJfNLMyYCZwRXMTMrOLzWyumc0tLy9PRK0i7cYLxRu55dklnHN4P66aMDLqciQFRL2zeCrwgLsPAM4GHjKzj9Tk7ne7e5G7F/Xq1eugFymSKko2VXLltLc5vLALv/zcOPUvLC2SyCBYCwyMezwgHBbvq8BjAO7+HyAP6JnAmkTarZr6Rq74y9vkZmVw9wVFdMjRYaLSMokMgjnACDMbamY5BDuDZzQZZzVwOoCZjSYIArX9iLTCLc8sYfH6Hfzyc+Po20V9CkjLJSwI3L0BuBx4DlhMcHTQIjO7ycwmh6NdA3zdzN4B/gJc5O6eqJpE2qtZxRt54N8r+fJJQzh9dJ+oy5EUk9CzS9x9JsFO4PhhP4y7XwyclMgaRNq7jTtq+M70dxjTr4DrJh0adTmSgqLeWSwiB6Ax5lw1bT419TF+N/UoXUhOWkXnm4uksDtfXc5/Srfw888cwfDenaIuR1KUtghEUtS8VRX8+oVlfGpcfz5XNCDqciSFKQhEUtCOmnqunPY2/brk8b//dRhmOl9AWk9NQyIp6P+eXsy6bdVMv+xECvKyoy5HUpy2CERSzL+Xb2banDV8/ePDdDE5aRMKApEUUlPfyPVPLGRwj466jpC0GTUNiaSQ38xaxqotVTzy9eN0CQlpM9oiEEkR767dzj2vreALRQM58RBdkkvajoJAJAXUN8a4dvoCuufn8L2zR0ddjrQzahoSSQF/fK2U4vU7uPOLR9Olo44SkralLQKRJFdaXsmts95j4ti+TDysX9TlSDukIBBJYrGYc/0TC8nNyuDHU8ZGXY60UwoCkST22Nw1vLFiKzecPZo+BepjQBJDQSCSpLZX1fOzZ5cwfkh3vnDswH2/QKSVFAQiSeo3s5axvbqeH00eo2sJSUIpCESS0NINO3lo9iqmjh/E2P5doi5H2jkFgUiScXdu+sci8nMyuebMUVGXI2lAQSCSZJ5btJHXS7ZwzZmj6J6fE3U5kgYUBCJJpKa+kf+dWcyoPp05/7hBUZcjaUJnFoskkXteK2XN1moe+dpxZGVqPU0ODn3TRJLE+u3V3P7yciYd1pcTh+uicnLwKAhEksRPZy4h5q6LyslBpyAQSQJzVm5lxjvruOSUQxjYvWPU5UiaURCIRCwWc276ezH9u+Rx2SmHRF2OpCEFgUjE/rFwPQvXbufbZ41Sr2MSCQWBSITqGmL88rmljO5XwLlHFkZdjqQpBYFIhB5+YxWrt1Zx3aRDycjQ9YQkGgoCkYjsrKnntpdKOGl4D04eocNFJToKApGI3PVqKVt31XHdxNG6uqhESkEgEoGNO2q451+lTB7Xn8MH6OqiEq2EBoGZTTSzpWZWYmbX7WGcz5tZsZktMrNHElmPSLK4ddYyGmPOd87S1UUlegm71pCZZQK3A2cAZcAcM5vh7sVx44wArgdOcvcKM+udqHpEkkXJpp08OmcNXzpxiE4ek6SQyC2C8UCJu5e6ex0wDZjSZJyvA7e7ewWAu29KYD0iSeFnzy6lY04WV5w2IupSRIDEBkEhsCbucVk4LN5IYKSZvW5ms81sYnMTMrOLzWyumc0tLy9PULkiiTd35VZeKN7IpacMU18DkjSi3lmcBYwATgWmAn80s65NR3L3u929yN2LevXqdXArFGkj7s5Pn1lC7865fOVjQ6MuR+R9iQyCtcDAuMcDwmHxyoAZ7l7v7iuAZQTBINLuvLRkE/NWVXDVhJF0zFFXIJI8EhkEc4ARZjbUzHKA84AZTcZ5kmBrADPrSdBUVJrAmkQiEYs5v3x+GYN7dORzRQOiLkfkQxIWBO7eAFwOPAcsBh5z90VmdpOZTQ5Hew7YYmbFwMvAd9x9S6JqEonKzHfXs3j9Dq6aMIJs9TwmScbcPeoa9ktRUZHPnTs36jJEWqyhMcaZt/6TTDOevepkMnVNIYmAmc1z96LmntOqiUiCPTl/HaXlu7jmzJEKAUlKCgKRBKpriHHrrGUcXtiFs8b2jbockWYpCEQS6NG5ayirqOaaM0fqwnKStBQEIglSU9/I7196j2OHdOOUkTr/RZKXgkAkQR76zyo27qjlmjNHaWtAkpqCQCQBKmsb+MOry/n4iJ4cP6xH1OWI7JWCQCQB7v/XCrbuquOaM3WZaUl+CgKRNra9qp67XyvljDF9OHJg16jLEdknBYFIG/vja6XsrGng6jNGRl2KSIsoCETa0NZdddz/+grOOaIfo/sVRF2OSIsoCETa0F2vLqe6vpFvTdBFdCV1KAhE2simnTU8+J+VTDmykOG9O0ddjkiLKQhE2sgfXllOfaNz5enaGpDUss8gMLMrzKzbwShGJFWt317Nw2+s5jNHFzKkZ37U5Yjsl5ZsEfQB5pjZY2Y20XSKpMhH3P5yCe6uDuklJe0zCNz9+wTdR94LXAS8Z2b/Z2aHJLg2kZRQVlHFo3PW8PmigQzs3jHqckT2W4v2EXjQe82G8NYAdAOmm9nPE1ibSEq47cUSzIzLTxsedSkirbLPHrTN7ErgQmAzcA9Bd5L1ZpYBvAdcm9gSRZLXys27mP5WGRccP5h+XTpEXY5Iq+wzCIDuwKfdfVX8QHePmdknE1OWSGr43YvvkZ1pfOMTaimV1LXPIHD3H+3lucVtW45I6ijZVMmT89fytY8Po3fnvKjLEWk1nUcg0kq/mbWMvOxMLjl5WNSliBwQBYFIKyxat52nF6znKycNpUen3KjLETkgCgKRVvjNC8soyMvi69oakHZAQSCyn95aXcGsxZu45JRD6NIhO+pyRA6YgkBkP/3q+aX07JTDRScOiboUkTahIBDZD/9evpnXS7Zw2anDyc9tydHXIslPQSDSQu7OL59bSt+CPM4/blDU5Yi0GQWBSAu9vHQTb63exhWnDycvOzPqckTajIJApAViMeeXzy1jUPeOfL5oYNTliLQpBYFICzy7aAPF63dw1YQRZGfq30baF32jRfahMeb8+oVljOjdiSlHFkZdjkibS2gQhB3ZLDWzEjO7bi/jfcbM3MyKElmPSGs8+fZaSjZVcvUZI8nMUL9M0v4kLAjMLBO4HZgEjAGmmtmYZsbrDFwJvJGoWkRaq6a+kV+/sIzDCguYeFjfqMsRSYhEbhGMB0rcvdTd64BpwJRmxvsJ8DOgJoG1iLTKn2evYu22aq6fNBr10irtVSKDoBBYE/e4LBz2PjM7Ghjo7k/vbUJmdrGZzTWzueXl5W1fqUgztlfX8/uXSzh5ZC9OGt4z6nJEEiayncVhD2e/Bq7Z17jufre7F7l7Ua9evRJfnAjwh1eWs726nu9OHBV1KSIJlcggWAvEH3A9IBy2W2fgMOAVM1sJHA/M0A5jSQbrtlVz/+srOPfIQsb27xJ1OSIJlcggmAOMMLOhZpYDnAfM2P2ku293957uPsTdhwCzgcnuPjeBNYm0yG9eWIY7XH3GyKhLEUm4hAWBuzcAlwPPAYuBx9x9kZndZGaTEzVfkQO1dMNOHn+rjAtPGMzA7h2jLkck4RJ6+UR3nwnMbDLsh3sY99RE1iLSUj97dgn5uVn8zyeGR12KyEGhM4tF4swu3cJLSzbxjVOH0y0/J+pyRA4KBYFIyN356TNL6FuQx5dPGhJ1OSIHjYJAJPTMuxt4Z802rj5jpC4zLWlFQSAC1DY08rNnlzCyTyc+c8yAqMsROagUBCLAvf9awaotVXz/nDG6sJykHQWBpL2NO2r4/UslTBjdh5NH6sx1ST8KAkl7P3tmCQ2Nzg8+OTrqUkQioSCQtDZvVQVPvL2Wr318KIN75EddjkgkFASStmIx58d/X0SfglydPCZpTUEgaWv6vDIWlG3n+kmjyc9N6En2IklNQSBpaUdNPT9/bgnHDO7GlCP7R12OSKS0GiRp6Xez3mPLrjruv2i8eh6TtKctAkk7JZsqeeDfK/lC0UAOH6C+BkQUBJJW3J2b/lFMh5xMvn2Weh4TAQWBpJm/L1jPP5eVc9WEkfTslBt1OSJJQUEgaaNiVx0/nrGIIwZ04UsnDI66HJGkoZ3FkjZ+8nQx26vr+fPXjiMrU+tAIrvpv0HSwqvLynnirbVcesohjO5XEHU5IklFQSDt3q7aBr73xEKG9crn8tN0BrFIU2oaknbvV88vY+22ah675AR1OCPSDG0RSLv29uoK7v/3Cr54/CDGD+0edTkiSUlBIO1WXUOM6x5fSN+CPL478dCoyxFJWmoaknbrzleXs3TjTu79UhGd87KjLkckaWmLQNqlJRt28PuXSvjUuP6cPrpP1OWIJDUFgbQ71XWNfPMvb1PQIZsffWpM1OWIJD01DUm785Oni1m2sZI/fWW8LiMh0gLaIpB25ZmF63nkjdVccsowdUQv0kIKAmk31m6r5ruPL2DcgC5cc4auLCrSUgoCaRcaGmNcNe1tYg6/m3oUOVn6aou0lPYRSLtw20slzFlZwa1fOJLBPfKjLkckpWi1SVLe7NIt3PbSe3z66ELOPaow6nJEUk5Cg8DMJprZUjMrMbPrmnn+ajMrNrMFZvaimeki8bJftlXV8a1H5zOoe0dumnJY1OWIpKSEBYGZZQK3A5OAMcBUM2t6UPfbQJG7HwFMB36eqHqk/alvjHH5I2+zubKW26YeTadctXSKtEYitwjGAyXuXurudcA0YEr8CO7+srtXhQ9nAwMSWI+0I+7Oj/++iH+VbOZ//+twdUIvcgASGQSFwJq4x2XhsD35KvBMc0+Y2cVmNtfM5paXl7dhiZKqHvz3Sv48ezWXnDyMzxcNjLockZSWFDuLzeyLQBHwi+aed/e73b3I3Yt69dJJQunulaWbuOkfxUwY3YdrdVVRkQOWyEbVtUD8qtqAcNiHmNkE4AbgFHevTWA90g68t3EnVzzyNqP6FvDb844kM8OiLkkk5SVyi2AOMMLMhppZDnAeMCN+BDM7CrgLmOzumxJYi7QDW3fV8dUH55Kbnck9XyoiXzuHRdpEwoLA3RuAy4HngMXAY+6+yMxuMrPJ4Wi/ADoBfzWz+WY2Yw+TkzRX29DIpQ/NY8OOGv544TEUdu0QdUki7UZCV6ncfSYws8mwH8bdn5DI+Uv70Bhzrp2+gDdXbuV3U4/iqEHdoi5JpF1Jip3FInvSGHO+M/0dnpq/jmsnjmLyuP5RlyTS7igIJGnFYs51jy/gibfWcvUZI/nGqcOjLkmkXVIQSFKKxZzv/W0hf51XxpWnj+Cbp4+IuiSRdktBIEknFnO+/9S7TJuzhitOG85VExQCIomkIJCk4u78aMYiHnljNZedeghXnzESM50rIJJICgJJGo0x54dPLeKh2au45ORhXHvWKIWAyEGgM3IkKeyqbeDKaW8za/EmLjl5GNdNOlQhIHKQKAgkcuu3V/PVB+ayZMMObpoylgtPGBJ1SSJpRUEgkVpYtp2vPjiHqrpG7rvoWE4d1TvqkkTSjoJAIvPsuxv41qPz6Z6fw+OXHceovp2jLkkkLSkI5KBzd+76Zyk/e3YJ4wZ05Y8XFtGrc27UZYmkLQWBHFSbdtZw7fQFvLK0nHOO6MevPjeOvOzMqMsSSWsKAjlonlu0geufWMiu2gZumjKWC44frCODRJKAgkASbldtAzf9vZhH567hsMICbv3CkQzvrf0BIslCQSAJNW9VBVc/Np/VW6v4xqmHcNWEkeRk6TxGkWSiIJCEqNhVx62zlvHQ7FX069KBRy8+gfFDu0ddlog0Q0Egbaq+McbDs1fxm1nvsbOmnvOPG8x3Jo6iIC876tJEZA8UBNJmXlm6iZufXkzJpko+NrwnP/jkGJ0bIJICFARywBav38HPn13Cy0vLGdKjI/dcWMTpo3vriCCRFKEgkFZxd95YsZU7X13OK0vL6ZybxffOPpSLThyqncEiKUZBIPslFnOeL97Ina8uZ/6abfTIz+HbZ47kguOH0KWj9gOIpCIFgbTI9qp6ZixYx/2vr6C0fBeDunfkJ+cexueOGaAzg0VSnIJA9qihMcZrJZuZPq+MF4o3UtcQY2z/Am6behSTDutLVqaagETaAwWBfIi7s2xjJU+8Xcbf3lrLpp21dO2YzX+PH8RnjxnA2P4F2gks0s4oCIS6hhhvrNjCi4s38eKSjazZWk1mhvGJUb357DGFfOLQ3uRmqflHpL1SEKSpsooqZpdu5aUlG/nnss1U1jaQm5XBx4b35NJTDuHMMX11aWiRNKEgSAOxmFNSXsmbK7YyZ+VW5qzYyrrtNQD07pzLp8b14/RD+3DS8J50yNGav0i6URC0M40xZ8XmShat20Hxuh0Ur9/BwrXb2VZVDwQ//McO7c4lQ7pz7JDuHNq3MxkZavMXSWcKghRV29DI6i1VlG7exYrNuygtr2TZxkqWbNhBTX0MgJzMDEb27cRZY/pyzJBuHDe0O4O6d9TOXhH5EAVBkqpriLFhew1rt1Wzdls167ZVs7aimnXbq1m1pYqyiipi/sH4vTrnckivfP57/GDG9i9gTP8ChvfuRLYO8RSRfUhoEJjZROC3QCZwj7vf0uT5XOBPwDHAFuAL7r4ykTVFpa4hxrbqOrZV1VOxq45t1fVsq6qjoqqezTtrKa+spXxneKusfb8pJ17PTrkUduvAEQO6cO5RhQzrmc+wXvkM6Zmvq3uKSKslLAjMLBO4HTgDKAPmmNkMdy+OG+2rQIW7Dzez84CfAV9IVE174+7UNcaoawhutXF/q+sbqXn/FqO2oZGquvBW20BVffB3V10ju2ob2FnTwM6aenbG3d/dXNOcDtmZ9OqcG67Vd+L4YT3o2SmXfl3zKOzagcKuHejbJU9n8IpIQiRyi2A8UOLupQBmNg2YAsQHwRTgxvD+dOD3Zmbu7rSxx+as4a5/Lqe+0WlojFHX6NQ3xmhojFHfGIRAa2VnGh1zsuiYk0l+bhad87Io6JDNgG4d6ZyXFd6y6dYxm64dc+jaMZtucX/zc9VCJyLRSeQvUCGwJu5xGXDcnsZx9wYz2w70ADbHj2RmFwMXAwwaNKhVxXTLz+HQfgVkZxjZmRlkZWaQkxl3PyuD3PCWk5VBTjgsLzuTvOwM8rIyyd19PzuTjjmZdMzOokNOpq62KSIpLSVWRd39buBugKKiolZtLZwxpg9njOnTpnWJiLQHiVyVXQsMjHs8IBzW7DhmlgV0IdhpLCIiB0kig2AOMMLMhppZDnAeMKPJODOAL4X3Pwu8lIj9AyIismcJaxoK2/wvB54jOHz0PndfZGY3AXPdfQZwL/CQmZUAWwnCQkREDqKE7iNw95nAzCbDfhh3vwb4XCJrEBGRvdPhLiIiaU5BICKS5hQEIiJpTkEgIpLmLNWO1jSzcmBVK1/ekyZnLacwLUvyaS/LAVqWZHUgyzLY3Xs190TKBcGBMLO57l4UdR1tQcuSfNrLcoCWJVklalnUNCQikuYUBCIiaS7dguDuqAtoQ1qW5NNelgO0LMkqIcuSVvsIRETko9Jti0BERJpQEIiIpLm0DAIzu8LMlpjZIjP7edT1HCgzu8bM3Mx6Rl1La5jZL8LPY4GZ/c3MukZd0/4ys4lmttTMSszsuqjraS0zG2hmL5tZcfj/cWXUNR0IM8s0s7fN7B9R13IgzKyrmU0P/08Wm9kJbTn9tAsCM/sEQV/J49x9LPDLiEs6IGY2EDgTWB11LQfgBeAwdz8CWAZcH3E9+8XMMoHbgUnAGGCqmY2JtqpWawCucfcxwPHA/6TwsgBcCSyOuog28FvgWXc/FBhHGy9T2gUBcBlwi7vXArj7pojrOVC/Aa4FUnavv7s/7+4N4cPZBL3ZpZLxQIm7l7p7HTCNYGUj5bj7end/K7y/k+AHpzDaqlrHzAYA5wD3RF3LgTCzLsDJBP234O517r6tLeeRjkEwEvi4mb1hZq+a2bFRF9RaZjYFWOvu70RdSxv6CvBM1EXsp0JgTdzjMlL0xzOemQ0BjgLeiLiU1rqVYCUpFnEdB2ooUA7cHzZz3WNm+W05g5TovH5/mdksoG8zT91AsMzdCTZ7jwUeM7NhydpF5j6W5XsEzUJJb2/L4e5PhePcQNA08fDBrE0+ysw6AY8DV7n7jqjr2V9m9klgk7vPM7NTIy7nQGUBRwNXuPsbZvZb4DrgB205g3bH3Sfs6Tkzuwx4Ivzhf9PMYgQXcio/WPXtjz0ti5kdTrCm8I6ZQdCc8paZjXf3DQexxBbZ22cCYGYXAZ8ETk/WUN6LtcDAuMcDwmEpycyyCULgYXd/Iup6WukkYLKZnQ3kAQVm9md3/2LEdbVGGVDm7ru3zKYTBEGbScemoSeBTwCY2UgghxS8MqG7L3T33u4+xN2HEHxZjk7GENgXM5tIsAk/2d2roq6nFeYAI8xsqJnlEPS9PSPimlrFgrWKe4HF7v7rqOtpLXe/3t0HhP8b5wEvpWgIEP5PrzGzUeGg04HitpxHu9wi2If7gPvM7F2gDvhSCq6Btje/B3KBF8Ktm9nufmm0JbWcuzeY2eXAc0AmcJ+7L4q4rNY6CbgAWGhm88Nh3wv7H5foXAE8HK5olAJfbsuJ6xITIiJpLh2bhkREJI6CQEQkzSkIRETSnIJARCTNKQhERNKcgkBEJM0pCERE0pyCQOQAmdmxYV8KeWaWH17H/7Co6xJpKZ1QJtIGzOxmgmvadCC4LsxPIy5JpMUUBCJtIDz1fw5QA5zo7o0RlyTSYmoaEmkbPYBOQGeCLQORlKEtApE2YGYzCHomGwr0c/fLIy5JpMXS8eqjIm3KzC4E6t39kbD/4n+b2Wnu/lLUtYm0hLYIRETSnPYRiIikOQWBiEiaUxCIiKQ5BYGISJpTEIiIpDkFgYhImlMQiIikuf8HcbgT2KnchIEAAAAASUVORK5CYII=)
%% Cell type:markdown id: tags:
Logistic regression assumes that the data is described by m instances or points, and at each of these there are n features. At each point from i=1,2,...,m, the probability of the data being classified as true or false is modelled by
$$ y_i = f(c_0 + c_1 x_{i,1} + c_2 x_{i,2} + \cdots + c_n x_{i,n}), $$
where $f$ is the **logistic** function. The algorithm then aims to calculate the coefficients $c_1$, $c_2$,..., $c_n$ and the intercept $c_0$. For linear regression the coefficients were calculated by minimizing the sum of the square of the errors. This is equivalent to maximimizing the likelihood of the observed data, assuming the data points are distributed with mean 0 and some standard deviation $\sigma$. Therefore, the observed data is the most likely data. For logistic regression these two statements are no longer equivalent, however, the coefficients can be calculated by again maximizing the likelihood of the observed data. This is done by using optimization algorithms.
The process for logistic regression using `sklearn` is the same as the other models we have so far considered, however the options change when instantiating the model. We will use the default optimisation solver `lbfgs` to fit the data and show the model partameters.
%% Cell type:code id: tags:
``` python
from sklearn.linear_model import LogisticRegression # import the LogisticRegression model
# instantiate the model (using the default parameters)
# penalty='none' implies no regularization (to be covered next week) and solver='lbfgs' is the default solver
# different solvers can be used, dependent on the type of penalties that are implemented
logreg = LogisticRegression(solver='lbfgs',penalty='none')
logreg.fit(X_train,y_train) # fit the training data to the model
print('Model coefficients are',np.round(logreg.coef_,3)) # print the model coefficients c1,...,c4
print('Model intercept is',np.round(logreg.intercept_,3)) # print the model intercept c0
```
%% Cell type:markdown id: tags:
Next we investigate the probabilities that are output. The second column, which is the probability of correctly picking the species, is the output of the logistic function. Since the problem is binary, the first column is just 1 minus the second column.
%% Cell type:code id: tags:
``` python
y_preda = logreg.predict_proba(X_test) # calculate the probabilities for the test features
# print out the probability table with a header
print('Probability table for testing set is:')
print(y_preda)
```
%% Cell type:markdown id: tags:
Using a classification threshold of 0.5, i.e., $P>0.5$ indicates true and $P<0.5$ indicates false, the probabilities are then used to predict the target test values, which can be compared against the actual target test values. We can see there are correct predictions at (1,1) (true positives) and (0,0) (true negative), but there are also false positives at (1,0) and false negatives at (0,1). False positives indicates that the model incorrectly predicts a positive value, similarly false negatives indicates that the model incorrectly predicts a negative value.
%% Cell type:code id: tags:
``` python
y_pred=logreg.predict(X_test) # calculate the predicted values of the model for the test features
plt.scatter(y_pred,y_test) # plot the predicted values against the actual test values
plt.xlabel('Predicted') # xlabel
plt.ylabel('Actual') # ylabel
plt.title('Logistic Regression'); # add a title
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYOklEQVR4nO3de5wlZX3n8c+XGS4qV5nGKAwOhsF1gkbcDrK6MbioC2ggJq6ZSVBxUdSsuiqbiIElSnS9rpcoXkDxgpGLusuOCmJUFHQdQqOCgppMAGGAyCAXQVRuv/xRNXpoumfOzHSdnu76vF+vflHnqTpVv6d76G/X85yqSlUhSeqvrWa7AEnS7DIIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCzWlJPpjkf27C+/ZMckeSBV3UtaVKcm6SF8x2HdqyxOsINCpJrgZeVFVfnqvHTnIk8BHgF8B9wFXAcVX1+c2tUZotnhFIG+9bVbU9sDPwfuCMJDvP9EH6drai2WMQaNYl2TbJu5Nc3369O8m2A+v/KskN7boXJakke7frPpbkje3yoiSfT3JrkpuTXJhkqySnAXsCn2uHg/4qyZJ2Pwvb9z40yUfbY9yS5OwN1V1V9wGnAQ8Blg705R1Jrknyk3bo6kEb0ZcPJDknyc+BpyZ5RJLPJlmb5KokrxzY1/5JJpL8rD3WO9v27ZJ8MslP2+/FxUke1q77WpIXtctbJTk+yY+T3JjkE0l2atet+/68oO3LTUmO2+QfsrZoBoG2BMcBBwCPB34X2B84HiDJwcBrgKcBewMHrmc/xwBrgDHgYcBfA1VVzwOuAf6wqravqrdN8d7TgAcDvwPsBrxrQ0W3f7G/ELgb+HHb/BZgn7YvewO7AydsRF/+DHgTsAPw/4HPAZe2+zkIeFWS/9xu+x7gPVW1I/DbwFlt+wuAnYDFwK7AS2mGsiY7sv16KvAoYHvgfZO2+Y/Ao9tjn5DkMdN/RzRXGQTaEvw5cGJV3VhVa4E3AM9r1z0X+GhVXV5VdwKvX89+7gYeDjyyqu6uqgtriEmwJA8HDgFeWlW3tO/9+nreckCSW4FfAu8AjqiqG5MEOBp4dVXdXFW3A/8LWL4Rffl/VfXN9mzjscBYVZ1YVXdV1ZXAKQP7uxvYO8miqrqjqlYNtO8K7F1V91bVJVX1symO9efAO6vqyqq6A3gdsHzdWVLrDVX1i6q6lCaQfnc93xfNUQaBtgSP4Dd/UdMuP2Jg3bUD6waXJ3s7sBr4UpIrkxw75PEXAzdX1S1Dbr+qqnYGdgFWAr/fto/RnFVc0g7J3Ap8sW2H4foy2PZI4BHr9tXu769pznYAjqI5+/hhO/zzrLb9NOA8mrmL65O8LcnWUxxrqu/7woH9A/zrwPKdNGcNmmcMAm0Jrqf5pbfOnm0bwA3AHgPrFk+3k6q6vaqOqapHAYcBr0ly0LrV6zn+tcBDN3bCt/0r+mXA85LsB9xEMwTzO1W1c/u1UzuxPGxfBuu8FrhqYF87V9UOVXVoe/x/rqoVNENZbwU+k+Qh7RnNG6pqGfAk4FnA86c41lTf93uAn2zM90Fzn0GgUdu6ncxc97UQOB04PslYkkU0Y+qfbLc/C3hhksckeTAw7TUDSZ6VZO92iOY24F6aj3hC88vtUVO9r6puAM4F3p9klyRbJ3nKMJ2pqpuBDwMntMM5pwDvSrJbW9PuA2P6Q/el9Y/A7Ulem+RBSRYk2TfJ77X7PiLJWHvcW9v33JfkqUke285h/IxmqOi+KfZ/OvDqJHsl2Z5mGOvMqrpnmL5r/jAINGrn0PzVvO7r9cAbgQngMuB7wLfbNqrqXODvgPNphn3WjYP/aop9LwW+DNwBfAt4f1Wd3657M03Y3Jrkf0zx3ufR/ML8IXAj8KqN6NO7gUOTPA547bo6k/ysrefRm9AXqupemr/mH09zvcJNNKGzU7vJwcDlSe6gmTheXlW/AH4L+AxNCPwA+DrNcNFkp7btF7T7/yXwio3ot+YJLyjTnNJ+auX7wLZz/S/X+dQXzW2eEWiLl+TZ7efzd6EZC//cXP3FOZ/6ovnDINBc8BKa4Zp/oRn3f9nslrNZ5lNfNE84NCRJPecZgST13MINb7JlWbRoUS1ZsmS2y5CkOeWSSy65qarGplo354JgyZIlTExMzHYZkjSnJPnxdOscGpKknjMIJKnnDAJJ6jmDQJJ6ziCQpJ7r7FNDSU6luWHWjVW17xTrQ3OjrENp7nN+ZFV9u4ta9n7dF7hn4Lq5hYHVb35mF4eSpBn3xDf9Az+5/a5fv37YDttw0XFPn7H9d3lG8DGauyNO5xCau0UupXmq0we6KGJyCADcU027JG3pJocAwE9uv4snvukfZuwYnQVBVV0A3LyeTQ4HPlGNVcDO7SMDZ9TkENhQuyRtSSaHwIbaN8VszhHszv0fy7embXuAJEcnmUgysXbt2pEUJ0l9MScmi6vq5Koar6rxsbEpr5CWJG2i2QyC67j/M1v3aNtm1MJsXLskbUketsM2G9W+KWYzCFYCz0/jAOC29tmxM2r1m5/5gF/6fmpI0lxx0XFPf8Av/Zn+1FBnzyNIcjpwILCI5sHhfwNsDVBVH2w/Pvo+mk8W3Qm8sKo2eDe58fHx8qZzkrRxklxSVeNTrevsOoKqWrGB9QX8t66OL0kazpyYLJYkdccgkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnOg2CJAcn+VGS1UmOnWL9nknOT/KdJJclObTLeiRJD9RZECRZAJwEHAIsA1YkWTZps+OBs6pqP2A58P6u6pEkTa3LM4L9gdVVdWVV3QWcARw+aZsCdmyXdwKu77AeSdIUugyC3YFrB16vadsGvR44Iska4BzgFVPtKMnRSSaSTKxdu7aLWiWpt2Z7sngF8LGq2gM4FDgtyQNqqqqTq2q8qsbHxsZGXqQkzWddBsF1wOKB13u0bYOOAs4CqKpvAdsBizqsSZI0SZdBcDGwNMleSbahmQxeOWmba4CDAJI8hiYIHPuRpBHqLAiq6h7g5cB5wA9oPh10eZITkxzWbnYM8OIklwKnA0dWVXVVkyTpgRZ2ufOqOodmEniw7YSB5SuAJ3dZgyRp/WZ7sliSNMsMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6rlOgyDJwUl+lGR1kmOn2ea5Sa5IcnmST3VZjyTpgRZ2teMkC4CTgKcDa4CLk6ysqisGtlkKvA54clXdkmS3ruqRJE2tyzOC/YHVVXVlVd0FnAEcPmmbFwMnVdUtAFV1Y4f1SJKm0GUQ7A5cO/B6Tds2aB9gnyTfTLIqycFT7SjJ0UkmkkysXbu2o3IlqZ9me7J4IbAUOBBYAZySZOfJG1XVyVU1XlXjY2Njo61Qkua5LoPgOmDxwOs92rZBa4CVVXV3VV0F/BNNMEiSRqTLILgYWJpkryTbAMuBlZO2OZvmbIAki2iGiq7ssCZJ0iSdBUFV3QO8HDgP+AFwVlVdnuTEJIe1m50H/DTJFcD5wF9W1U+7qkmS9ECpqtmuYaOMj4/XxMTEbJchSXNKkkuqanyqdbM9WSxJmmUGgST13LRXFid5LzDtuFFVvbKTiiRJI7W+W0w4EC9JPTBtEFTVx0dZiCRpdmzwpnNJxoDXAsuA7da1V9V/6rAuSdKIDDNZ/Pc01wHsBbwBuJrmYjFJ0jwwTBDsWlUfAe6uqq9X1X8FPBuQpHlimOcR3N3+94YkzwSuBx7aXUmSpFEaJgjemGQn4BjgvcCOwKs7rUqSNDIbDIKq+ny7eBvw1G7LkSSN2jCfGvooU1xY1s4VSJLmuGGGhj4/sLwd8GyaeQJJ0jwwzNDQZwdfJzkd+EZnFUmSRmpTbjq3FNhtpguRJM2OYeYIbuf+cwT/SnOlsSRpHhhmaGiHURQiSZodGxwaSvKVYdokSXPT+p5HsB3wYGBRkl2AtKt2BHYfQW2SpBFY39DQS4BXAY8ALuE3QfAz4H3dliVJGpX1PY/gPcB7kryiqt47wpokSSM0zMdH70uy87oXSXZJ8hfdlSRJGqVhguDFVXXruhdVdQvw4s4qkiSN1DBBsCDJuvkBkiwAtumuJEnSKA1zr6EvAmcm+VD7+iXAud2VJEkapWGC4LXA0cBL29eXAb/VWUWSpJHa4NBQVd0HXETzrOL9aR5T+YNuy5Ikjcr6LijbB1jRft0EnAlQVT6cRpLmkfUNDf0QuBB4VlWtBkjiIyolaZ5Z39DQHwM3AOcnOSXJQfzm6mJJ0jwxbRBU1dlVtRz4d8D5NLeb2C3JB5I8Y0T1SZI6Nsxk8c+r6lNV9YfAHsB3GPJ5BEkOTvKjJKuTHLue7f4kSSUZH7pySdKM2KgnlFXVLVV1clUdtKFt2wvPTgIOAZYBK5Ism2K7HYD/TvPJJEnSiG3KoyqHtT+wuqqurKq7gDOAw6fY7m+BtwK/7LAWSdI0ugyC3YFrB16vYdJzDJI8AVhcVV9Y346SHJ1kIsnE2rVrZ75SSeqxLoNgvZJsBbwTOGZD27bDUeNVNT42NtZ9cZLUI10GwXXA4oHXe7Rt6+wA7At8LcnVwAHASieMJWm0ugyCi4GlSfZKsg2wHFi5bmVV3VZVi6pqSVUtAVYBh1XVRIc1SZIm6SwIquoe4OXAeTT3Jjqrqi5PcmKSw7o6riRp4wxz99FNVlXnAOdMajthmm0P7LIWSdLUZm2yWJK0ZTAIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5zoNgiQHJ/lRktVJjp1i/WuSXJHksiRfSfLILuuRJD1QZ0GQZAFwEnAIsAxYkWTZpM2+A4xX1eOAzwBv66oeSdLUujwj2B9YXVVXVtVdwBnA4YMbVNX5VXVn+3IVsEeH9UiSptBlEOwOXDvwek3bNp2jgHOnWpHk6CQTSSbWrl07gyVKkraIyeIkRwDjwNunWl9VJ1fVeFWNj42NjbY4SZrnFna47+uAxQOv92jb7ifJ04DjgD+oql91WI8kaQpdnhFcDCxNsleSbYDlwMrBDZLsB3wIOKyqbuywFknSNDoLgqq6B3g5cB7wA+Csqro8yYlJDms3ezuwPfDpJN9NsnKa3UmSOtLl0BBVdQ5wzqS2EwaWn9bl8SVJG7ZFTBZLkmaPQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzy3scudJDgbeAywAPlxVb5m0flvgE8C/B34K/GlVXT3TdSw59gsPaLv6Lc+c6cNIUieOP/t7nH7RtdxbxYKEFU9czBv/6LEztv/OzgiSLABOAg4BlgErkiybtNlRwC1VtTfwLuCtM13HVCGwvnZJ2pIcf/b3+OSqa7i3CoB7q/jkqms4/uzvzdgxuhwa2h9YXVVXVtVdwBnA4ZO2ORz4eLv8GeCgJOmwJkmaU06/6NqNat8UXQbB7sBgpWvatim3qap7gNuAXSfvKMnRSSaSTKxdu7ajciVpy7PuTGDY9k0xJyaLq+rkqhqvqvGxsbHZLkeSRmbBNIMk07Vvii6D4Dpg8cDrPdq2KbdJshDYiWbSWJIErHji4o1q3xRdBsHFwNIkeyXZBlgOrJy0zUrgBe3yc4CvVs3g+Q7TfzrITw1Jmgve+EeP5YgD9vz1GcCChCMO2HNGPzWUGf69e/+dJ4cC76b5+OipVfWmJCcCE1W1Msl2wGnAfsDNwPKqunJ9+xwfH6+JiYnOapak+SjJJVU1PtW6Tq8jqKpzgHMmtZ0wsPxL4L90WYMkaf3mxGSxJKk7BoEk9ZxBIEk9ZxBIUs91+qmhLiRZC/x4E9++CLhpBsuZC+xzP9jnfticPj+yqqa8InfOBcHmSDIx3cen5iv73A/2uR+66rNDQ5LUcwaBJPVc34Lg5NkuYBbY536wz/3QSZ97NUcgSXqgvp0RSJImMQgkqefmZRAkOTjJj5KsTnLsFOu3TXJmu/6iJEtmocwZNUSfX5PkiiSXJflKkkfORp0zaUN9HtjuT5JUkjn/UcNh+pzkue3P+vIknxp1jTNtiH/beyY5P8l32n/fh85GnTMlyalJbkzy/WnWJ8nftd+Py5I8YbMPWlXz6ovmltf/AjwK2Aa4FFg2aZu/AD7YLi8HzpztukfQ56cCD26XX9aHPrfb7QBcAKwCxme77hH8nJcC3wF2aV/vNtt1j6DPJwMva5eXAVfPdt2b2eenAE8Avj/N+kOBc4EABwAXbe4x5+MZwf7A6qq6sqruAs4ADp+0zeHAx9vlzwAHJTP43LfR22Cfq+r8qrqzfbmK5olxc9kwP2eAvwXeCvxylMV1ZJg+vxg4qapuAaiqG0dc40wbps8F7Ngu7wRcP8L6ZlxVXUDzfJbpHA58ohqrgJ2TPHxzjjkfg2B34NqB12vatim3qap7gNuAXUdSXTeG6fOgo2j+opjLNtjn9pR5cVV9YZSFdWiYn/M+wD5JvplkVZKDR1ZdN4bp8+uBI5KsoXn+yStGU9qs2dj/3zeo0wfTaMuT5AhgHPiD2a6lS0m2At4JHDnLpYzaQprhoQNpzvouSPLYqrp1Novq2ArgY1X1v5P8B+C0JPtW1X2zXdhcMR/PCK4DBp/qvEfbNuU2SRbSnE7+dCTVdWOYPpPkacBxwGFV9asR1daVDfV5B2Bf4GtJrqYZS105xyeMh/k5rwFWVtXdVXUV8E80wTBXDdPno4CzAKrqW8B2NDdnm6+G+v99Y8zHILgYWJpkryTb0EwGr5y0zUrgBe3yc4CvVjsLM0dtsM9J9gM+RBMCc33cGDbQ56q6raoWVdWSqlpCMy9yWFXN5QdeD/Nv+2yaswGSLKIZKlrvc8C3cMP0+RrgIIAkj6EJgrUjrXK0VgLPbz89dABwW1XdsDk7nHdDQ1V1T5KXA+fRfOLg1Kq6PMmJwERVrQQ+QnP6uJpmUmb57FW8+Ybs89uB7YFPt/Pi11TVYbNW9GYass/zypB9Pg94RpIrgHuBv6yqOXu2O2SfjwFOSfJqmonjI+fyH3ZJTqcJ80XtvMffAFsDVNUHaeZBDgVWA3cCL9zsY87h75ckaQbMx6EhSdJGMAgkqecMAknqOYNAknrOIJCknjMI1DtJ7k3y3STfT/LpJA/ejH19LMlz2uUPJ1m2nm0PTPKkTTjG1e01AVInDAL10S+q6vFVtS9wF/DSwZXt1eYbrapeVFVXrGeTA4GNDgKpawaB+u5CYO/2r/ULk6wErkiyIMnbk1zc3vP9JfDre8G/r70//peB3dbtKMnX1t3Cor2H/reTXNo+/2EJTeC8uj0b+f0kY0k+2x7j4iRPbt+7a5Ivtc8T+DDN7Yalzsy7K4ulYbV/+R8CfLFtegKwb1VdleRomkv3fy/JtsA3k3wJ2A94NM197x8GXAGcOmm/Y8ApwFPafT20qm5O8kHgjqp6R7vdp4B3VdU3kuxJc/XsY2iuJP1GVZ2Y5Jk099KROmMQqI8elOS77fKFNLcceRLwj+2N2gCeATxu3fg/zY0Jl9I8NOT0qroXuD7JV6fY/wHABev2VVXT3Vv+acCygUdh7Jhk+/YYf9y+9wtJbtm0bkrDMQjUR7+oqscPNrS/jH8+2AS8oqrOm7TdTD4GcSvggKq630Nz5vYzkjQXOUcgTe084GVJtgZIsk+Sh9A89vJP2zmEh9M8AnSyVcBTkuzVvvehbfvtNLfHXudLDDxEJcnj28ULgD9r2w4BdpmpTklTMQikqX2YZvz/22keIv4hmjPo/wv8c7vuE8C3Jr+xqtYCRwP/J8mlwJntqs8Bz143WQy8EhhvJ6Ov4DefXnoDTZBcTjNEdE1HfZQA7z4qSb3nGYEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLP/Ruv00B9QxIrZwAAAABJRU5ErkJggg==)
%% Cell type:markdown id: tags:
To determine the number of correct predictions, false positives and false negatives, we can construct a confusion matrix. From the confusion matrix it can be seen that there are 15 true negatives, 5 true positives, 2 false positives where versicolor was incorrectly predicted, and 8 false negatives, where the other species was incorrectly predicted. The false positive and negatives will typically relate to the virginica species, rather than setosa.
%% Cell type:code id: tags:
``` python
from sklearn.metrics import confusion_matrix # import the confusion matrix function
cnf_matrix = confusion_matrix(y_test, y_pred) # create a confusion matrix for our actual and predicted values
# create a data frame from the confusion matrix with the column and row names being the class_names
class_names=['other', 'versicolor'] # names of the binary classes for plotting
cmatrix = pd.DataFrame(cnf_matrix,columns=class_names,index=class_names)
f, ax = plt.subplots(figsize=(7,6)) # initialise the plots and axes
sns.heatmap(cmatrix, annot=True, linewidths=.5) # plot the confusion matrix as a heatmap
plt.title('Confusion matrix') # add a title
plt.ylabel('Actual label') # add a ylabel
plt.xlabel('Predicted label') # add a xlabel
# adjust the bottom and top of the figure, so we can view all of it
bottom, top = ax.get_ylim() # get the y axis limits
ax.set_ylim(bottom + 0.5, top - 0.5); # adjust the y axis limits
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZ0AAAGDCAYAAADwA81JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjU0lEQVR4nO3deZhlVX3v//eHeR4bQQHHK3gVIyoag0oIqNcpgok3OF4DaBtnwSl6VcR7c6/5aQwEjdKAoGL4tQQ1DpFBDKMyNI0yOwSUyCDzKDJ0fe8fZxeUnabr1LB31a7zfj3Pfvqcffbe61vV9dS31lrfvXaqCkmSurDGXAcgSRodJh1JUmdMOpKkzph0JEmdMelIkjpj0pEkdcako3kryfpJvp3k9iTHz+A6r0ty8mzGNleSPD/JT+c6Dmm64n06mqkkrwUOBJ4E3An8GPibqjprhtd9A/BOYNeqemCmcc53SQp4YlX9Yq5jkdpiT0czkuRA4BDg/wBbA48G/hHYaxYu/xjgZ6OQcIaRZK25jkGaKZOOpi3JpsAngLdX1der6u6qur+qvl1V72+OWTfJIUmubbZDkqzbfLZ7kl8neW+SG5Jcl2Tf5rODgY8B+yS5K8n+ST6e5NgJ7T82SY3/Mk7yl0muTHJnkquSvG7C/rMmnLdrkvObYbvzk+w64bPTkvyvJGc31zk5yaKH+frH4//AhPj3TvLSJD9LckuSD084/tlJfpTktubYzyZZp/nsjOawnzRf7z4Trv/BJNcDR4/va855QtPGM5r3j0pyY5LdZ/L/KrXJpKOZ+CNgPeAbqznmfwLPAXYGngY8G/jIhM+3ATYFtgX2Bz6XZPOqOohB72lpVW1UVUetLpAkGwL/ALykqjYGdmUwzLfycVsA322O3RL4DPDdJFtOOOy1wL7AI4B1gPetpultGHwPtmWQJI8AXg88E3g+8NEkj2uOXQEcACxi8L3bE3gbQFXt1hzztObrXTrh+lsw6PUtnthwVf078EHg2CQbAEcDX6qq01YTrzSnTDqaiS2BmyYZ/nod8ImquqGqbgQOBt4w4fP7m8/vr6p/Be4CdpxmPGPATknWr6rrqurSVRzzMuDnVfWVqnqgqo4DrgD+dMIxR1fVz6rqHuBrDBLmw7mfwfzV/cD/zyChHFpVdzbtX8Yg2VJVF1TVOU27vwQOB/54iK/poKq6t4nn91TVEcAvgHOBRzJI8tK8ZdLRTNwMLJpkruFRwK8mvP9Vs+/Ba6yUtH4LbDTVQKrqbmAf4K+A65J8N8mThohnPKZtJ7y/fgrx3FxVK5rX40nhNxM+v2f8/CQ7JPlOkuuT3MGgJ7fKobsJbqyq301yzBHATsBhVXXvJMdKc8qko5n4EXAvsPdqjrmWwdDQuEc3+6bjbmCDCe+3mfhhVZ1UVS9k8Bf/FQx+GU8Wz3hM10wzpqn4PIO4nlhVmwAfBjLJOastL02yEYNCjqOAjzfDh9K8ZdLRtFXV7QzmMT7XTKBvkGTtJC9J8v81hx0HfCTJVs2E/MeAYx/umpP4MbBbkkc3RQwfGv8gydZJ9mrmdu5lMEw3topr/CuwQ5LXJlkryT7Ak4HvTDOmqdgYuAO4q+mFvXWlz38DPH6K1zwUWFZVb2IwV/WFGUcptcikoxmpqr9jcI/OR4Abgf8A3gF8sznkfwPLgIuAi4Hlzb7ptHUKsLS51gX8fqJYo4njWuAWBnMlK/9Sp6puBl4OvJfB8OAHgJdX1U3TiWmK3segSOFOBr2wpSt9/nHgS011219MdrEkewEv5qGv80DgGeNVe9J85M2hkqTO2NORJHXGpCNJ6oxJR5LUGZOOJKkzJh1JUmfm86q1ltVJWsgmuzF4aPffdOWMf1+uvejxsxbP6sznpCNJGsbYismPmSdMOpLUd7WqxTfmJ5OOJPXdWH+SjoUEkqTO2NORpJ4rh9ckSZ3p0fCaSUeS+q5HPR3ndCRJnbGnI0l95306kqTO9Gh4zaQjSX1nIYEkqSt9Kpm2kECS1Bl7OpLUdw6vSZI606PhNZOOJPWdJdOSpM70qKdjIYEkqTP2dCSp7ywkkCR1xuE1SVJnxsZmvk0iyReT3JDkklV89t4klWTRZNcx6UiShnEM8OKVdybZHngRcPUwFzHpSFLPVa2Y8TZ5G3UGcMsqPvp74ANADROrczqS1HdzNKeTZC/gmqr6SZKhzjHpSFLfzUL1WpLFwOIJu5ZU1ZLVHL8B8GEGQ2tDM+lIUt/NQk+nSTAPm2RW4QnA44DxXs52wPIkz66q6x/uJJOOJGnKqupi4BHj75P8Etilqm5a3XkWEkhS342tmPk2iSTHAT8Cdkzy6yT7TydUezqS1HcdFBJU1Wsm+fyxw1zHpCNJfecyOJKkzrgMjiRJ/5k9HUnqO4fXJEmdMelIkroyzNpp84VzOpKkztjTkaS+c3hNktSZHpVMm3Qkqe/s6UiSOtOjno6FBJKkztjTkaS+c3hNktSZHg2vmXQkqe961NNxTkeS1Bl7OpLUdz3q6Zh0JKnvnNORJHXGno4kqTM96ulYSCBJ6ow9HUnqO4fXJEmd6dHwmklHkvrOno4kqTMmndlx/01XznUIGhFrL3o8a62z7VyHoRHywH3XzHUIc2JeJx1J0hCq5jqCoZl0JKnvHF6TJHWmR0nHm0MlSZ2xpyNJfed9OpKkzvRoeM2kI0l9Z/WaJKkzPerpWEggSeqMPR1J6rse9XRMOpLUd1avSZK6UmMWEkiSutKj4TULCSRJnbGnI0l916M5HXs6ktR3YzXzbRJJvpjkhiSXTNj3qSRXJLkoyTeSbDbZdUw6ktR3Y2Mz3yZ3DPDilfadAuxUVX8A/Az40GQXMelIkiZVVWcAt6y07+SqeqB5ew6w3WTXcU5HkvpuflSv7QcsnewgezqS1HdVM96SLE6ybMK2eNjmk/xP4AHgq5Mda09HkvpuFno6VbUEWDLV85L8JfByYM+qyZe7NulIUt/N0YoESV4MfAD446r67TDnOLwmSZpUkuOAHwE7Jvl1kv2BzwIbA6ck+XGSL0x2HXs6ktR3HdwcWlWvWcXuo6Z6HZOOJPWdC35KkrpS86NkeigmHUnqux71dCwkkCR1xp6OJPVdj1aZNulIUt/1aHjNpCNJfdejQgLndCRJnbGnI0l95/CaJKkzFhJIkjpjT0eS1JU+rUhgIYEkqTP2dCSp7xxekyR1xqQjSeqM1WuSpM70qKfTWiFBkjWTHNDW9SVJ/dNa0qmqFcCqHm8qSZpFNVYz3rrS9vDa2Uk+CywF7h7fWVXLW25XkkZHj4bX2k46Ozf/fmLCvgL2aLldSRodPbo5tNWkU1V/0ub1JUn90uqKBEm2TnJUku8175+cZP8225SkkTNWM9860vYyOMcAJwGPat7/DHhPy21K0mjpUdJpe05nUVV9LcmHAKrqgSQrWm5zZHzk/3yGM84+jy0234xvHvsFAD531LGc8K0T2XyzTQF491veyG67Pnsuw9QCtN12j+KYLx7KI7ZeRFVx5JFf5bDPHjXXYY2sKgsJxt2dZEsGxQMkeQ5we8ttjoy9X/pCXvvnr+DD/+vTv7f/Dfvszb6vfdUcRaVR8MADD/D+DxzMhT++hI022pDzzj2R7596Bpdf/vO5Dm00Wb32oAOBbwFPSHI2sBXgb8NZssvOT+Wa634z12FoBF1//Q1cf/0NANx1191cccXP2fZR25h0NKm2q9eWJ/ljYEcgwE+r6v422xQcd8K3+daJp/KUJz2R97/jzWy6ycZzHZIWsMc8Zjt2ftpOnHvehXMdyujqUU+ni+fpPBt4GvAM4DVJ/sfDHZhkcZJlSZYtWbKkg9AWnn1e+TK+97UvcsIxn2OrLbfgU589Yq5D0gK24YYb8LWlR3Dg+w7izjvvmutwRpYrEjSSfAV4AvBjYLyAoIAvr+r4qloCjGebuv+mK9sMb0FatMXmD75+1Stewtvff9AcRqOFbK211uL4pUdw3HHf4Jvf/N5chzPaetTTaXtOZxfgydWn0oqeu/GmW9hq0RYAnHr6D/kvj3/MHEekheqIJX/H5Vf8gkMOdVRizvVnQYLWk84lwDbAdS23M5Lef9AnOf/Ci7jttjvYc+/X87b938D5F17ET39+JQS23WZrDvrAu+Y6TC1Az931Wbzh9a/ioosvY9n5JwPw0Y9+ku+d+IM5jkzzXdrohCT5NoNhtI0ZrL92HnDv+OdV9YohLuPwmjqz9qLHs9Y62851GBohD9x3TWbrWre9bo8Z/yLf7Ks/mLV4Vqetns6nJz9EkjQrRn1Op6pOB0jyt1X1wYmfJflb4PQ22pWkkdSjOZ22S6ZfuIp9L2m5TUnSPNVKTyfJW4G3MViJ4KIJH20M/LCNNiVpVHV5n81MtTWn80/A94D/C3wS2K3Zf1ZVeduyJM2mUR9eq6rbq+qXwDnAscAiBuuufSnJO9toU5JGlSsSPGR/4DlVdTc8WETwI+CwltuVpNEx6j2dCcJDy9/QvO6kFlySNP+03dM5Gjg3yTea93sDPulJkmZR2dMZqKrPAPsCtzTbvlV1SJttStLIGZuFbRJJvpjkhiSXTNi3RZJTkvy8+Xfz1V0DOni0QVUtr6p/aDYr1yRpltXYzLchHAO8eKV9fw2cWlVPBE5t3q9WF8/TkSS1qYOeTlWdwWDEaqK9gC81r7/EYApltUw6kqTfe4hmsy0e4rStq2r8KQLXA1tPdkLbhQSSpJbNRiHBSg/RnM75lWTSG35MOpLUc3NYvfabJI+squuSPBK4YbITHF6TpJ7rqJBgVb4FvLF5/UbgXyY7waQjSZpUkuMYrCizY5JfJ9mfwdqaL0zyc+AFzfvVcnhNkvqu2l/opape8zAf7TmV65h0JKnn+rQigUlHknquxvqzpKVJR5J6rk89HQsJJEmdsacjST1XHRQSzJaHTTpJ7gTG7y4d/4qqeV1VtUnLsUmShtCn4bWHTTpVtXGXgUiSpqdPhQRDzekkeV6SfZvXi5I8rt2wJEkL0aRzOkkOAnYBdmTwJNB1gGOB57YbmiRpGDXpMpvzxzCFBK8Eng4sB6iqa5M49CZJ80SfhteGSTr3TVyyOsmGLcckSZqChZZ0vpbkcGCzJG8G9gOOaDcsSdKwFtTwWlV9OskLgTuAHYCPVdUprUcmSVpwhr059GJgfQb36VzcXjiSpKnq0/DapCXTSd4EnAf8GfAq4Jwk+7UdmCRpOFWZ8daVYXo67weeXlU3AyTZEvgh8MU2A5MkDWdBrEgwwc3AnRPe39nskyTNA2MLZO21A5uXvwDOTfIvDOZ09gIu6iA2SdICs7qezvgNoP/ebOP+pb1wJElTtSBWma6qg7sMRJI0PX2qXhtm7bWtgA8ATwHWG99fVXu0GJckaUh9ujl0mFWmvwpcATwOOBj4JXB+izFJkhaoYarXtqyqo5K8u6pOB05PYtKRpHliQQ2vAfc3/16X5GXAtcAW7YUkSZqKBVEyPcH/TrIp8F7gMGAT4IBWo5IkDW1BVK+Nq6rvNC9vB/6k3XAkSVPVp0KC1d0cehiDm0FXqare1UpEkqQFa3U9nWWdRSFJmrYFMadTVV/qMhBJ0vQsqDkdSdL8tiDmdCRJ/dCn4bVhViSQJGlWzOvqtbUXPb7tJqQHPXDfNXMdgjQtC2VOZ86r107dep+5DkEjYs/fLOW1j3nlXIehEfJPv/rGrF2rT8NrVq9JUs/1qI5g6EcbfBB4Mj7aQJI0A8M+2uByfLSBJM1LY5UZb10ZJulsWVVHAfdX1elVtR9gL0eS5omqzHjrio82kKSeG5vrAKbARxtIUs8VC6B6bZyPNpAkzZZhqteOZhUVec3cjiRpjo11UDOd5ADgTQzywcXAvlX1u6leZ5jhte9MeL0e8EoG8zqSpHlgrOXhtSTbAu8CnlxV9yT5GvBq4JipXmuY4bUTVmr8OOCsqTYkSWpHR3M6awHrJ7kf2IBpdj6ms+DnE4FHTKcxSdLsG5uFLcniJMsmbIvHr19V1wCfBq4GrgNur6qTpxPrMHM6d/L7czrXM1ihQJK0QFTVEmDJqj5LsjmwF4NFAm4Djk/y+qo6dqrtDDO8tvFULypJ6k4Hw2svAK6qqhsBknwd2BWYctKZdHgtyanD7JMkzY3ZGF6bxNXAc5JskCTAngyWR5uy1T1PZz0Gk0WLmq7VeCrdBNh2Oo1JkmZf2ysSVNW5Sf4ZWA48AFzIwwzFTWZ1w2tvAd4DPAq4gIeSzh3AZ6fTmCSpn6rqIOCgmV5ndc/TORQ4NMk7q+qwmTYkSWpHn5bBGaZkeizJZuNvkmye5G3thSRJmoqxzHzryjBJ581Vddv4m6q6FXhzaxFJkqZkjMx468owy+CsmSRVVQBJ1gTWaTcsSdKwFtTjqoETgaVJDm/ev6XZJ0nSlAyTdD4ILAbe2rw/BTiitYgkSVOyoB7iVlVjwBeajSTPZ/Awt7e3G5okaRhj6U/12jA9HZI8HXgN8BfAVcDX2wxKkjS8BTGnk2QHBonmNcBNwFIgVeXTQyVpHlkow2tXAGcCL6+qX8CDT46TJGlaVnefzp8xeG7CvyU5Isme0KPbXiVpRCyIm0Or6ptV9WrgScC/MViH7RFJPp/kRR3FJ0maRJ9uDp10RYKquruq/qmq/hTYjsHqoj7ETZLmiZqFrStTelx1Vd1aVUuqas+2ApIkLVxTSjrDSrJmkq+2cW1J0u/r05zOUPfpTFVVrUjymCTrVNV9bbQhSRpYKCXTM3UlcHaSbwF3j++sqs+02KYkjZwFcXPoLPj3ZlsD2LjFdiRppHU5PDZTrSWdqjoYIMlGzfu72mpLktQPrRQSACTZKcmFwKXApUkuSPKUttqTpFE1NgtbV9ocXlsCHFhV/waQZHcGj0TYtcU2JWnkWEgwsOF4wgGoqtOSbNhie5I0kso5HQCuTPJR4CvN+9czqGiTJM2iPvV0WpvTAfYDtmLw7J2vN6/3a7E9SdI812b12q3Au9q6viRpoE89nVlPOkm+zWruVaqqV8x2m5I0ykb95tBPt3BNSdLDGOmbQ6vq9PHXSdYBdmje/rSq7p/t9iRJ/dHanE5zX86XgF8yeOLo9kneWFVntNWmJI2ikZ7TmeDvgBdV1U8BkuwAHAc8s8U2JWnkmHQG1h5POABV9bMka7fYniSNpFEvJBi3LMmRwLHN+9cBy1psT5JG0kgXEkzwVuDtPHSvzpnAP7bY3sjb/i0v5VGv3QOAuy6/msvf/XnG7rV2Q+049KzDuefuexhbMcbYihV85E/fP9chqQfaTDprAYeOP7QtyZrAui22N9LW3WZztn/TSzjn+Qcy9rv72WnJe9h67125bunpk58sTdPfvPqj3HnrnXMdxsjr05xOm8vgnAqsP+H9+sD3W2xv5GXNNVhjvXXImmuw5gbrcO/1t851SJI6ULOwdaXNns56Ex/cVlV3JdmgxfZG2r3X38rVn/8Oz13+j4zdcx+3nH4Rt5x+0VyHpQWsKP762IOg4NSvnsQPjjtlrkMaWWM9KiVoM+ncneQZVbUcIMkzgXtWd0KSxcBigMMPP5wntBjcQrPWphuy6MW78MNnvYMHbv8tTz3yALb58+dx/QlnzXVoWqAO/vMPc+tvbmGTLTflQ8cexLX/fg1XnHfZXIelea7N4bX3AMcnOTPJWcBS4B2rO6GqllTVLlW1y+LFi1sMbeHZYren8rurb+D+m++kHljBDd89j02fteNch6UF7Nbf3ALAHTffzrKTzuUJOz9xjiMaXT45FKiq85M8CRj/zecyOC363TU3sckznsga66/D2D33scXzd+KOn/j4IrVj3fXXJWuE3939O9Zdf12eutvOfP3Qr811WCOrP4Nr7awyvUdV/SDJn6300Q5JqKqvz3abgjuW/4IbvnMuzz7lk9SKMe68+Cqu+Yp1G2rHpos244AlHwRgzbXW5Ox/OZOLTr9wjqMaXX2qXmujp/PHwA+AP13FZ8XggW5qwVWfOp6rPnX8XIehEXDDf/yGD73kwLkOQ42Rvjm0qg5q/t13tq8tSZobSTYDjgR2YtCB2K+qfjTV67RWSJDk3Uk2ycCRSZYneVFb7UnSqBqjZrwN4VDgxKp6EvA04PLpxNpm9dp+VXUH8CJgS+ANwCdbbE+SRlLbN4cm2RTYDTgKoKruq6rbphNrm0lnfJTxpcCXq+rSCfskSbNkNkqmkyxOsmzCNvG+lccBNwJHJ7mwGb3acDqxtpl0LkhyMoOkc1KSjelXkYUk9cJsDK9NvE+y2ZZMaGIt4BnA56vq6cDdwF9PJ9ZW7tNJEuBjwFbAlVX12yRbAhYXSFL//Br4dVWd27z/Z+ZT0qmqSvKvVfXUCftuBm5uoz1JGmVt3xxaVdcn+Y8kOzYP59wTmNaaR22uvbY8ybOq6vwW25CkkdfRvMU7ga8mWQe4kmmOXLWZdP4QeH2SXzIY/wuDTtAftNimJI2cLlaZrqofA7vM9DptJp3/1uK1JUk91Fr1WlX9Ctge2KN5/ds225OkUeVD3IAkBzHoiu0IHA2sDRwLPLetNiVpFPXpXpQ2h9deCTwdWA5QVdc29+pIkmZR9ejhBm0mnfua0ukCmO7dq5Kk1etTT6fNOZavJTkc2CzJm4HvA0e02J4kaZ5rs6dTwFnAHcAOwMeq6pQW25OkkdRFyfRsaTPpbATsB9wCLAUuarEtSRpZ/Uk57ZZMH1xVTwHeDjwSOD2Jz0+WpFnW0fN0ZkWbPZ1xNwDXM1h37REdtCdJI8VCAiDJ25KcBpzK4CFub3YJHEkabW32dLYH3tOs1yNJaon36QBV9aG2ri1Jekifhte6mNORJLWoTz0dF+CUJHXGno4k9ZzDa5KkzoxVf4bXTDqS1HP9STkmHUnqvT6tvWYhgSSpM/Z0JKnn+lQybdKRpJ6zek2S1Jk+zemYdCSp5/o0vGYhgSSpM/Z0JKnnnNORJHWmXJFAktSVPhUSOKcjSeqMPR1J6jnndCRJnelTybRJR5J6rk9zOiYdSeq5PlWvWUggSeqMPR1J6jkLCSRJnbGQQJLUGQsJJEmdsZBAkqRVsKcjST3Xp+G1zONu2bwNTJJmQWbrQrtv94IZ/7487dffn7V4VseejiT13FhHnYckawLLgGuq6uXTuYZzOpKkYb0buHwmFzDpSFLP1Sxsk0myHfAy4MiZxOrwmiT1XEeFBIcAHwA2nslF7OlIUs+NUTPekixOsmzCtnj8+kleDtxQVRfMNFZ7OpLUc7NRhVxVS4AlD/Pxc4FXJHkpsB6wSZJjq+r1U23Hno4kabWq6kNVtV1VPRZ4NfCD6SQcsKcjSb3Xp5tDTTqS1HNdrjJdVacBp033fJOOJPXcPF5Z5j8x6UhSz/VpeM1CAklSZ+zpSFLPObwmSepMn4bXTDqS1HNdVq/NlHM6kqTO2NORpJ7r6nk6s8GkI0k916fhNZOOJPWcPR1JUmf61NOxkECS1Bl7OpLUcw6vSZI606fhNZOOJPWcPR1JUmf61NOxkECS1Bl7OpLUc1Vjcx3C0Ew6ktRzrjItSepMn56n45yOJKkz9nQkqeccXpMkdaZPw2smHUnqOW8OlSR1xptDJUlaBXs6ktRzzulIkjpj9ZokqTN96uk4pyNJ6ow9HUnqOUumJUmd6dPwmklHknrOQgJJUmf61NOxkECS1Bl7OpLUcxYSSJI606e110w6ktRz9nQkSZ2xkECSpFWwpyNJPdenOR17OpLUc1U14211kmyf5N+SXJbk0iTvnm6s9nQkqec6mNN5AHhvVS1PsjFwQZJTquqyqV7Ino4kabWq6rqqWt68vhO4HNh2OteypyNJPdfljE6SxwJPB86dzvnzOelkrgPooySLq2rJXMeh0eHP3Nx74L5rZvz7MsliYPGEXUtW/n9NshFwAvCeqrpjWu30qb5bk0uyrKp2mes4NDr8mRsNSdYGvgOcVFWfme51nNORJK1WkgBHAZfPJOGASUeSNLnnAm8A9kjy42Z76XQuNJ/ndDQ9jq2ra/7MLXBVdRazNM/unI4kqTMOr0mSOmPS6aEkmyV524T3uyf5zlzGpNGR5BNJXjCN8/w5lUmnpzYD3jbZQcNK4tyefk8GVvn7oao+VlXf7yAGfy4XIJNODyQ5MMklzfYe4JPAE5oKkk81h22U5J+TXJHkq02JI0memeT0JBckOSnJI5v9pyU5JMkyYNqL92l+S/LJJG+f8P7jSd6X5P1Jzk9yUZKDm88em+SnSb4MXAJsn+SY5ufu4iQHNMcdk+RVzetnJflhkp8kOS/JxknWS3J0c86FSf5kFXFtkeSbTfvnJPmDCfF9JcnZwFc6+BapY/4lMc8leSawL/CHDKpHzgVeD+xUVTs3x+zOYFmKpwDXAmcDz01yLnAYsFdV3ZhkH+BvgP2ay6/jTX0L3lLgEOBzzfu/AP6WQQnssxn8TH0ryW7A1cATgTdW1TnNz962VbUTDIZ1J144yTrN9fepqvOTbALcw+CPmKqqpyZ5EnBykh1Wiutg4MKq2jvJHsCXgZ2bz54MPK+q7pml74HmEZPO/Pc84BtVdTdAkq8Dz1/FcedV1a+bY34MPBa4DdgJOKXp+KwJXDfhnKVtBa35oaouTPKIJI8CtgJuBZ4KvAi4sDlsIwbJ5mrgV1V1TrP/SuDxSQ4DvgucvNLldwSuq6rzm7buAEjyPAZ/7FBVVyT5FbBy0nke8OfNMT9IsmWTtAC+ZcJZuEw6C8e9E16vYPB/G+DSqvqjhznn7taj0nxwPPAqYBsGf2g8Bvi/VXX4xIOahRwf/JmoqluTPA34b8BfMegl7Uf7/LlcwJzTmf/OBPZOskGSDYFXMhg+23iIc38KbJXkj2CwdlKSp7QXquappcCrGSSe44GTgP2axRtJsm2SR6x8UpJFwBpVdQLwEeAZKx3yU+CRSZ7VHL9xM/l/JvC6Zt8OwKObYyeaeMzuwE3TXUBS/WJPZ55rHpp0DHBes+vIqrogydlJLgG+x2DoY1Xn3tdM+P5Dkk0Z/H8fAlzafuSaL6rq0gwevHVNVV0HXJfkvwI/aoZd72IwT7hipVO3BY6eUMX2oZWue18zT3hYkvUZzOe8APhH4PNJLmbw8K+/rKp7m7bGfRz4YpKLgN8Cb5y1L1jzmisSSJI64/CaJKkzJh1JUmdMOpKkzph0JEmdMelIkjpj0tGcSLKiWTvukiTHJ9lgBteauBbYkUmevJpjd0+y6zTa+GVz38pQ+1c65q4ptvXxJO+baoxSH5h0NFfuqaqdm3W97mNwx/uDprvCcFW9qaouW80huwNTTjqSZodJR/PBmcB/aXohZyb5FnBZkjWTfGrCashvgQeX3f9ssyLy94EH76ZvVs/epXn94iTLmxWQT22Wefkr4ICml/X8JFslOaFp4/wkz23O3TLJyUkuTXIkQzyqt1k1+YLmnMUrffb3zf5Tk2zV7HtCkhObc85sFseUFjRXJNCcano0LwFObHY9g8EK2lc1v7hvr6pnJVkXODvJyQxW1N6RwWrEWwOXAV9c6bpbAUcAuzXX2qKqbknyBeCuqvp0c9w/AX9fVWcleTSDJWL+K3AQcFZVfSLJy4D9h/hy9mvaWB84P8kJVXUzsCGwrKoOSPKx5trvAJYAf1VVP0/yhwzu5N9jGt9GqTdMOpor6zerYcOgp3MUg2Gv86rqqmb/i4A/GJ+vATZlsBrybsBxVbUCuDbJD1Zx/ecAZ4xfq6pueZg4XgA8ecISLZs0a5LtBvxZc+53k9w6xNf0riSvbF5v38R6MzDGQyt6Hwt8vWljV+D4CW2vO0QbUq+ZdDRX7hl/HtC45pfvxBWGA7yzqk5a6biXzmIcawDPqarfrSKWoTWLVr4A+KOq+m2S04D1Hubwatq9beXvgbTQOaej+ewk4K1J1obBisXNSttnAPs0cz6PBP7TkymBc4DdkjyuOXeLZv+d/P4K3ScD7xx/k2Tn5uUZwGubfS8BNp8k1k2BW5uE8yQGPa1xazBY4Znmmmc1KypfleS/N20kg8cISAuaSUfz2ZEM5muWNytqH86gd/4N4OfNZ18GfrTyiVV1I7CYwVDWT3hoeOvbwCvHCwmAdwG7NIUKl/FQFd3BDJLWpQyG2a6eJNYTgbWSXM7gceLnTPjsbuDZzdewB/CJZv/rgP2b+C4F9hrieyL1mqtMS5I6Y09HktQZk44kqTMmHUlSZ0w6kqTOmHQkSZ0x6UiSOmPSkSR1xqQjSerM/wOZJpaBS2BFfwAAAABJRU5ErkJggg==)
%% Cell type:markdown id: tags:
We can calculate various scores to determine the accuracy of the classifier. The accuracy score is the number of correct scores divided by the number of samples. The precision score is the number of true positives divided by the true positives plus false positives. The recall score is the number of true positives divided by the true positives plus false negatives. The precision and recall scores give an indication of how well the algorithm is able to pick positive samples, whereas the accurary gives an indication of how well the algorithm can predict correct samples overall.
%% Cell type:code id: tags:
``` python
from sklearn.metrics import accuracy_score, precision_score, recall_score # import the score functions
print("Accuracy:",np.round(accuracy_score(y_test, y_pred),3)) # calculate and print the accuracy score
print("Precision:",np.round(precision_score(y_test, y_pred),3)) # calculate and print the precision score
print("Recall:",np.round(recall_score(y_test, y_pred),3)) # calculate and print the recall score
```
%% Cell type:markdown id: tags:
To further quantify the accurary of the classifier we can construct a Response Operating Characteristic (ROC) curve, which plots the false positive rate (FPR) against the true positive rate (TPR) as the classification threshold is changed. The TPR is the same as the precision score, while the FPR is the number of false positives divided by the false positives plus true negatives. For a perfect classifier we would expect that TPR=1 and FPR=0. We have considered up until this point a classification threshold of 0.5. If the classification threshold is very high (close to 1), then we are not going to flag any positive respones, so TPR=0 and FPR=0. If the classification threshold is very low (close to 0), then we will only flag positive responses and TPR=1 and FPR=1. Hence if our classifier is good, as the threshold changes the ROC curve will hug the left and top boundaries of the figure. If this occurs the area under the ROC curve (AUC) will approach 1. Therefore classifiers can be compared for a particular problem by measuring the AUC. Any classifier which is below the line y=x and has AUC less than 0.5 will be useless.
%% Cell type:code id: tags:
``` python
# import the functions to calculate the parameters for the ROC curve and the AUC
from sklearn.metrics import roc_curve, roc_auc_score
y_pred_proba = logreg.predict_proba(X_test)[::,1] # extract the second column of the model probabilities
# calculate the false positive and true positive rates as the threshold is varied, we don't use thresholds
fpr, tpr, thresholds = roc_curve(y_test, y_pred_proba)
auc = roc_auc_score(y_test, y_pred_proba) # calculate the area under the ROC curve (AUC)
# plot the FPR vs TPR and format label with AUC to 3 decimal places
plt.plot(fpr,tpr,label="Logistic Regression, auc = %0.3f " % auc)
plt.plot([0,1],[0,1],'k--') # plot x = y for comparison
plt.plot([0, 0, 1], [0, 1, 1], 'g--', label="'Perfect' Classifier")
plt.xlabel('False Positive Rate') # add xlabel
plt.ylabel('True Positive Rate') # add ylabel
plt.title('ROC curve') # add title
plt.legend(loc='best'); # add legend
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABBjUlEQVR4nO3dd3gU1dfA8e8hoSogAhakt0ACIUgQAWnSm6AoIKDgLwqIUqVYaAIKIkiTXowgCAiiNBsoomIhYggQiiGUhB6E0CHlvn/sJm9C2gaymSR7Ps+zDzs7d2fOJGHP3rkz54oxBqWUUq4rl9UBKKWUspYmAqWUcnGaCJRSysVpIlBKKReniUAppVycJgKllHJxmgiUUsrFaSJQOY6IHBWR6yJyRUROi4i/iNx7W5t6IvKjiFwWkUgR2SAinre1KSQi00XkuH1bh+3LxTL3iJRyLk0EKqdqb4y5F/ABagJvxa0QkbrA98DXQAmgHLAb+E1Eytvb5AG2Al5AK6AQUBc4DzzmrKBFxN1Z21YqJZoIVI5mjDkNfIctIcSZDCw1xswwxlw2xvxnjBkJ/AGMtbd5ESgNPG2MCTbGxBpjzhpjxhtjNie3LxHxEpEfROQ/ETkjIm/bX/cXkQkJ2jUWkfAEy0dFZISIBAFX7c/X3LbtGSIy0/68sIgsFpFTInJCRCaIiNvd/aSUK9NEoHI0ESkJtAZC7MsFgHrAF8k0Xw00tz9vBnxrjLni4H4KAluAb7H1Mipi61E46nmgLXAfsBJoY98m9g/5zsAKe1t/INq+j5pAC+DldOxLqUQ0Eaic6isRuQyEAWeBMfbX78f2d38qmfecAuLO/xdNoU1K2gGnjTFTjTE37D2NP9Px/pnGmDBjzHVjzDFgF/C0fd2TwDVjzB8i8iDQBhhkjLlqjDkLTAO6pmNfSiWiiUDlVB2NMQWBxkAV/v8D/gIQCzyczHseBiLsz8+n0CYlpYDDdxSpTdhtyyuw9RIAuvH/vYEyQG7glIhcFJGLwHzggbvYt3JxmghUjmaM+RnbqZQp9uWrwO/Ac8k078z/n87ZArQUkXsc3FUYUD6FdVeBAgmWH0ou1NuWvwAa209tPc3/J4Iw4CZQzBhzn/1RyBjj5WCcSiWhiUC5gulAcxGpYV9+E+gpIgNEpKCIFLEP5tYF3rW3WYbtQ3etiFQRkVwiUlRE3haRNsnsYyPwsIgMEpG89u3Wsa8LxHbO/34ReQgYlFbAxphzwDbgE+CIMWa//fVT2K54mmq/vDWXiFQQkUbp/aEoFUcTgcrx7B+qS4HR9uVfgZbAM9jGAY5hG3R9whjzr73NTWwDxgeAH4BLwF/YTjElOfdvjLmMbaC5PXAa+BdoYl+9DNvlqUexfYivcjD0FfYYVtz2+otAHiAY26muNaTvNJZSiYhOTKOUUq5NewRKKeXiNBEopZSL00SglFIuThOBUkq5uGxX4KpYsWKmbNmyVoehlFLZyt9//x1hjCme3LpslwjKli1LQECA1WEopVS2IiLHUlqnp4aUUsrFaSJQSikXp4lAKaVcnCYCpZRycZoIlFLKxTktEYjIEhE5KyJ7U1gvIjJTREJEJEhEHnVWLEoppVLmzB6BP7ZJv1PSGqhkf/QG5joxFqWUUilw2n0ExpjtIlI2lSYdsE0gboA/ROQ+EXnYXm/dKQJPBzLo20FJXn+/6fvUK1WPHWE7eHvr20nWT281HZ+HfNgSuoUJ2yckWT+/3Xw8inmw4eAGpv4+Ncn6ZU8vo1ThUqzau4q5AUnz3ZrOayhWoBj+gf74B/onWb+5+2YK5C7AnJ1zWL1vdZL123ptA2DKjilsPLQx0br8ufPzTfdvABj/83i2Hkk8jW7RAkVZ23ktAG9teYvfw39PtL5koZJ89sxnAAz6dhCBpwMTra9ctDIL2i8AoPeG3hw6fyjRep+HfJjeajoAPb7sQfil8ETr65asy8RmEwHotLoT56+dT7S+abmmjGo0CoDWy1tzPep6ovXtKrdjaL2hADT2b8ztyt3TjGv/NSXa3OCX/95Isr5s/jaUK9CWm7EX2XHhnSTrKxR4mtL5m3Et5gx/XhyXZL3HPc9TIt8TXIo+xt+Rk5Os97y3Fw/mrc2FqEMEXpqRZH31gn0plqc6Ebf2sOfyvCTrfQoNpEjuypy5uZPgK/5J1tcqPJxC7mU4eeNXDl79PMn6OveNpoDbgxy/voXD19YlWV+vyHvkzXUfR65t4uj1zUnWN7h/Ku6Sj5Brawm7/mOS9U2KzgbgwNUVnLrxW6J1bpKXhvd/BMC+K59w9mbi+3/y5CpM/SLvAxB0eS7nbyU+eZDfrTiP3zcWgH8uTedi1L+J1hd0L4Vv4TcBCIicxOXoxJO83Ze7EjULDQLgj4tjuR5zLtH6onmq4V3wVQB+u/A2t2IjE61/IK8vXve+BMD2/4YQY24mWv9wvvpUuacbAD+df43blcr/JBULdMqQv70/LowlNjqaAY+uYkz7jJ+DyMobyh4h8fR84fbXkiQCEemNrddA6dKlMyU4lTMEhV/kVuQlKj+Ux+pQlLoj50KC+O/KfiSXG7GxsU7Zh1PnI7D3CDYaY6ols24jMMk+SQgishUYYYxJ9bZhX19fcyd3Fm8J3QJAs/LN0v1elX11mW/r4azqU9fiSJRKn4sXLzJs2DAWLVpExYoVWbRoEY0a3flEdCLytzHGN7l1VvYITmCb8DtOSftrThF3SkcTgVIqq4uJiaFevXocPHiQ4cOHM3bsWPLnz++0/VmZCNYDr4vISqAOEOnM8QGllMrqzp8/z/3334+bmxvvvfcepUqVwtc32S/xGcqZl49+DvwOeIhIuIj4iUhfEelrb7IZCAVCgIVAP2fFopRSWZkxhs8++4zKlSuzaNEiAJ5++ulMSQLg3KuGnk9jvQGSDrUrpZQLCQsLo2/fvmzevJnHH3+c+vXrZ3oMemexUkpZ5PPPP8fLy4tt27Yxffp0fv31Vzw9PTM9jmw3H8Gdmt9uvtUhKKVUIkWKFKFOnTosWLCAcuXKWRaHyyQCj2IeVoeglHJx0dHRTJs2jVu3bvHOO+/QqlUrWrZsiYhYGpfLnBracHADGw5usDoMpZSL2r17N48//jjDhw8nKCiIuHu4rE4C4EKJYOrvU5Mt/6CUUs508+ZNRo0aha+vL2FhYXzxxResXLkySySAOC6TCJRSygr//vsvH3zwAd26dSM4OJhnn302SyUBcKExAqWUyixXrlzh66+/pnv37lSrVo0DBw5Qvnx5q8NKkfYIlFIqA/3www9Ur16dF154gf379wNk6SQA2iNQWdSKP4/zdeDdl54KPnUJz4cLZUBESqXuwoULDB06lCVLllC5cmV+/vlnqlatanVYDnGZRLDs6WVWh6DS4evAExnyIe75cCE6+DySQVEplbyYmBjq16/PoUOHeOuttxg9ejT58uWzOiyHuUwiKFW4VNqNVJbi+XAhLR+tsrSIiIj4InHvv/8+pUuX5tFHs9+suy4zRrBq7ypW7V1ldRhKqRzAGMPSpUsTFYnr2LFjtkwC4EKJYG7A3GSniVRKqfQ4duwYrVu3pmfPnlStWpWGDRtaHdJdc5lEoJRSd+uzzz6jWrVq/Prrr8yaNYtffvmFKlWqWB3WXXOZMQKllLpbxYsXp379+syfP58yZcpYHU6G0USglFIpiIqKYurUqURFRTFq1ChatmxJixYtstydwXdLTw0ppVQy/vnnH+rUqcNbb71FcHBwlioSl9FcpkewpvMaq0NQSmUDN27cYNy4cUyePJlixYqxdu1annnmGavDciqXSQTFChSzOgSlVDYQEhLClClTePHFF5k6dSpFihSxOiSnc5lE4B/oD0Avn16WxqGUynquXLnCunXreOGFF6hWrRoHDx60dMawzOYyYwT+gf7xyUAppeJ89913eHl50bNnz/gica6UBMCFEoFSSiV0/vx5evbsSatWrShQoAC//PJLtikSl9Fc5tSQUkrFiSsSFxISwjvvvMPIkSOzVZG4jKaJQCnlMs6dO0fRokVxc3Pjgw8+oEyZMvj4+FgdluU0EWRzGVW3P6vReQRURjLG4O/vz5AhQ5g0aRJ9+vShQ4cOVoeVZbhMItjcfbPVIThFRtXtz2p0HgGVUY4ePUrv3r354YcfaNCgAU2aNLE6pCzHZRJBgdwFrA7BabRuv1LJW7ZsGa+++ioiwpw5c+jTpw+5cuk1MrdzmUQwZ+ccAPrV7mdxJEqpzPLggw/SsGFD5s2bR+nSpa0OJ8tymUSwet9qQBOBUjlZVFQUkydPJiYmhtGjR9OiRQtatGhhdVhZnvaRlFI5wq5du6hduzYjR47k4MGD8UXiVNo0ESilsrXr16/z5ptv8thjj3HmzBnWrVvH8uXLc2SVUGdxaiIQkVYiclBEQkTkzWTWlxaRn0TkHxEJEpE2zoxHKZXzhIaG8tFHH9GrVy+Cg4Pp2LGj1SFlO05LBCLiBswGWgOewPMi4nlbs5HAamNMTaArMMdZ8Silco5Lly7h7+8PgJeXF//++y+LFi1yiUqhzuDMweLHgBBjTCiAiKwEOgDBCdoYIO4C+MLASWcFs63XNmdtWimViTZv3kzfvn05ceIEderUoWrVqjlq2kgrOPPU0CNAWILlcPtrCY0FeohIOLAZ6J/chkSkt4gEiEjAuXPnnBGrUiqLi4iI4IUXXqBt27YULFiQ3377zWWLxGU0qweLnwf8jTElgTbAMhFJEpMxZoExxtcY41u8ePE72tGUHVOYsmPK3UWrlLJEXJG4lStXMnr0aHbt2sXjjz9udVg5hjNPDZ0ASiVYLml/LSE/oBWAMeZ3EckHFAPOZnQwGw9tBGBovaEZvWmllJOcOXOG4sWL4+bmxpQpUyhTpgze3t5Wh5XjOLNHsBOoJCLlRCQPtsHg9be1OQ40BRCRqkA+QM/9KOXijDEsXrwYDw8PFixYAED79u01CTiJ0xKBMSYaeB34DtiP7eqgfSIyTkSesjd7A3hFRHYDnwO9jN4FopRLCw0NpVmzZrz88sv4+PjQrFkzq0PK8ZxaYsIYsxnbIHDC10YneB4M1HdmDEqp7OPTTz+lX79+uLm5MW/ePF555RUtEpcJXKbWUP7c+a0OQSmVhhIlSvDkk08yd+5cSpYsaXU4LsNlEsE33b+xOgSl1G1u3brFpEmTiI2NZezYsTRv3pzmzZtbHZbL0T6XUsoSO3fupFatWowZM4bQ0FAtEmchl0kE438ez/ifx1sdhlIu79q1awwdOpTHH3+cCxcusH79epYuXapF4izkMolg65GtbD2y1eowlHJ5R44cYdasWbzyyivs27eP9u3bWx2Sy3OZMQKllHUiIyP58ssveemll/Dy8iIkJIRSpUql/UaVKVymR6CUssamTZvw8vLi5Zdf5sCBAwCaBLIYTQRKKac4d+4c3bt3p127dhQpUoTff/+dKlWqWB2WSobLnBoqWqCo1SEo5TJiYmJ44oknOHLkCO+++y5vvvkmefLksToslQKXSQRrO6+1OgSlcrzTp0/zwAMP4ObmxtSpUylbtizVqlWzOiyVBodPDYlIAWcGopTKvmJjY5k/fz6VK1dm/vz5ALRr106TQDaRZiIQkXoiEgwcsC/XEJFsN6XkW1ve4q0tb1kdhlI5TkhICE2bNqVv377Url2bli1bWh2SSidHTg1NA1piLyFtjNktIg2dGpUT/B7+u9UhKJXjfPLJJ/Tr1488efKwcOFC/Pz89MawbMihMQJjTNhtv9wY54SjlMpOSpcuTcuWLZk9ezaPPHL7TLQqu3AkEYSJSD3AiEhuYCC2+QWUUi7m5s2bTJw4kdjYWMaNG0fTpk1p2rSp1WGpu+TIYHFf4DVsE8+fAHyAfk6MSSmVBf3555/UqlWLd999l+PHj2uRuBzEkUTgYYzpbox50BjzgDGmB1DV2YFltJKFSlKykNY3Vyq9rl69ypAhQ6hbty6RkZFs3LgRf39/HQvIQRw5NTQLeNSB17K0z575zOoQlMqWjh07xpw5c+jbty+TJk2iUKFCVoekMliKiUBE6gL1gOIiMiTBqkKAm7MDU0pZ5+LFi6xZs4aXX34ZT09PQkJCdMawHCy1U0N5gHuxJYuCCR6XgGedH1rGGvTtIAZ9O8jqMJTK8r7++ms8PT3p27dvfJE4TQI5W4o9AmPMz8DPIuJvjDmWiTE5ReDpQKtDUCpLO3v2LAMGDGDVqlV4e3uzfv16LRLnIhwZI7gmIh8CXkC+uBeNMU86LSqlVKaKiYmhfv36HD9+nAkTJjB8+HBy585tdVgqkziSCJYDq4B22C4l7Qmcc2ZQSqnMcfLkSR566CHc3NyYMWMGZcuWxdPT0+qwVCZz5PLRosaYxUCUMeZnY8z/AO0NKJWNxcbGMnfuXKpUqcK8efMAaNOmjSYBF+VIjyDK/u8pEWkLnATud15IzlG5aGWrQ1AqSzh06BCvvPIK27dvp1mzZrRu3drqkJTFHEkEE0SkMPAGtvsHCgGDnBmUMyxov8DqEJSy3OLFi3n99dfJly8fS5YsoVevXnpjmEo7ERhjNtqfRgJNAESkvjODUko5R9myZWndujWzZ8/m4YcftjoclUWkdkOZG9AZW42hb40xe0WkHfA2kB+omTkhZozeG3oD2jNQruXmzZuMHz8egAkTJmiROJWs1HoEi4FSwF/ATBE5CfgCbxpjvsqE2DLUofOHrA5BqUy1Y8cO/Pz8OHDgAP/73/8wxuhpIJWs1BKBL+BtjIkVkXzAaaCCMeZ85oSmlLoTV65c4Z133mHWrFmUKlWKb7/9VmcNU6lK7fLRW8aYWABjzA0gNL1JQERaichBEQkRkTdTaNNZRIJFZJ+IrEjP9pVSSR0/fpz58+fz2muvsXfvXk0CKk2p9QiqiEiQ/bkAFezLAhhjjHdqG7aPMcwGmgPhwE4RWW+MCU7QphLwFlDfGHNBRB64i2NRymVduHCBL774gt69e+Pp6UloaCglSpSwOiyVTaSWCO52zoHHgBBjTCiAiKwEOgDBCdq8Asw2xlwAMMacvct9psjnIR9nbVopS61bt45+/fpx7tw5GjVqhIeHhyYBlS6pFZ2720JzjwBhCZbDgTq3takMICK/YSttPdYY8+3tGxKR3kBvsM2Reiemt5p+R+9TKqs6ffo0/fv3Z82aNfj4+LBp0yY8PDysDktlQw5NXu/k/VcCGgMlge0iUt0YczFhI2PMAmABgK+vr86Pp1xeTEwMDRo0ICwsjPfff5+hQ4dqkTh1x5yZCE5gu/w0Tkn7awmFA38aY6KAIyJyCFti2JnRwfT4sgegM5Wp7C08PJwSJUrg5ubGzJkzKVeunJaKVnfNkaJziEh+EUlvn3MnUElEyolIHqArsP62Nl9h6w0gIsWwnSoKTed+HBJ+KZzwS+HO2LRSThcbG8usWbOoUqUKc+fOBaB169aaBFSGSDMRiEh7IBD41r7sIyK3f6AnYYyJBl4HvgP2A6uNMftEZJyIPGVv9h1wXkSCgZ+AYXqfglKJHThwgIYNGzJgwACeeOIJ2rVrZ3VIKodx5NTQWGxXAG0DMMYEikg5RzZujNkMbL7ttdEJnhtgiP2hlLrNokWLeP311ylQoACffvopL7zwgt4drDKcQ2WojTGRt/3x6YCtUpmgQoUKtG/fno8//pgHH3zQ6nBUDuVIItgnIt0AN/sNYAOAHc4NK+PVLVnX6hCUStONGzcYN24cAO+//z5NmjShSZMmFkelcjpHEkF/4B3gJrAC23n9Cc4MyhkmNptodQhKpeq3337Dz8+PgwcP8vLLL2uROJVpHLlqqIox5h1jTG37Y6S99pBSKgNcvnyZ/v3706BBA27evMl3333HwoULNQmoTONIIpgqIvtFZLyIVHN6RE7SaXUnOq3uZHUYSiURHh7OokWL6N+/P3v27KFFixZWh6RcjCMzlDURkYewTVIzX0QKAauMMdnq9ND5a1nrqtQVfx7n68Db769Lv+BTl/B8uFAGRKQy0/nz51m9ejWvvvoqVatWJTQ0VGcMU5Zx6IYyY8xpY8xMoC+2ewpGp/4OlZavA08QfOrSXW/H8+FCdPB5JAMiUpnBGMOaNWvw9PRkwIABHDx4EECTgLJUmj0CEakKdAE6AeeBVdgmsld3yfPhQqzqo1czuYpTp07x2muvsW7dOmrVqsX333+vReJUluDIVUNLsH34tzTGnHRyPErlSHFF4k6cOMHkyZMZPHgw7u5W13xUysaRMYIc8ZW1aTmdsFtlvrCwMB555BHc3NyYPXs25cqVo3LlylaHpVQiKY4RiMhq+797RCQowWNPgpnLso1RjUYxqtEoq8NQLiImJoaZM2cmKhLXsmVLTQIqS0qtRzDQ/q9WuFIqHfbv34+fnx+///47rVu3pn379laHpFSqUuwRGGNO2Z/2M8YcS/gA+mVOeBmn9fLWtF7e2uowVA63YMECfHx8OHToEMuWLWPTpk13PKueUpnFkctHmyfzWrb7RL0edZ3rUdetDkPlcJUqVeLpp58mODiYHj166N3BKltI8dSQiLyK7Zt/+dvGBAoCvzk7MKWyg+vXrzN27FhEhEmTJmmROJUtpdYjWAG0xzarWPsEj1rGmB6ZEJtSWdr27dupUaMGkydPJjIyEtv0GkplP6klAmOMOQq8BlxO8EBE7nd+aEplTZcuXaJfv340atSImJgYtm7dyty5c/U0kMq2UrtqaAW2K4b+xjYRTcK/cgOUd2JcGa5dZb34SWWMkydP4u/vz5AhQxg3bhz33HOP1SEpdVdSTATGmHb2fx2aljKrG1pvqNUhqGwsIiKC1atX069fP6pUqcKRI0d0xjCVYzgyeX19EbnH/ryHiHwkIno9nHIJxhhWrVqFp6cngwYN4tChQwCaBFSO4sjlo3OBayJSA1uxucPAMqdG5QSN/RvT2L+x1WGobOTkyZN07NiRrl27UqZMGf7++2+9M1jlSI5UvYo2xhgR6QB8bIxZLCJ+zg5MKSvFxMTQsGFDTpw4wZQpUxg4cKAWiVM5liN/2ZdF5C3gBaCBiOQCcjs3LKWscezYMUqWLImbmxtz5syhfPnyVKxY0eqwlHIqR04NdcE2cf3/jDGngZLAh06NSqlMFhMTw0cffUTVqlXji8S1aNFCk4ByCWkmAvuH/3KgsIi0A24YY5Y6PTKlMsnevXupV68eb7zxBk2bNqVjx45Wh6RUpnJkhrLO2HoA27DdSzBLRIYZY9Y4ObYM1dmrs9UhqCxo3rx5DBgwgMKFC7NixQq6du2qN4Ypl+PIGME7QG1jzFkAESkObAGyVSLoVzvbFUxVTmSMQUSoWrUqzz33HNOnT6d48eJWh6WUJRxJBLnikoDdeRyc9D4ruRZ1DYACuQtYHImy0rVr1xg9ejRubm588MEHNGrUiEaNGlkdllKWcuQD/VsR+U5EeolIL2ATsNm5YWW8Nsvb0GZ5G6vDUBbatm0b3t7eTJ06lStXrmiROKXsHBksHgbMB7ztjwXGmBHODkypjBIZGUmfPn3iy0P/+OOPzJ49W8cClLJLbT6CSsAUoAKwBxhqjDmRWYEplVFOnTrFZ599xtChQ3n33XcpUEBPDyqVUGo9giXARqATtgqks9K7cRFpJSIHRSRERN5MpV0nETEi4pvefSiVnHPnzjFrlu1PtkqVKhw9epQPP/xQk4BSyUgtERQ0xiw0xhw0xkwByqZnwyLiBszGNq2lJ/C8iHgm064gMBD4Mz3bVyo5xhhWrFhB1apVeeONN+KLxOkVQUqlLLWrhvKJSE3+fx6C/AmXjTG70tj2Y0CIMSYUQERWAh2A4NvajQc+AIalM/Z06eXTy5mbV1lAWFgYr776Kps2baJOnTosXrxYi8Qp5YDUEsEp4KMEy6cTLBvgyTS2/QgQlmA5HKiTsIGIPAqUMsZsEpEUE4GI9AZ6A5QufWcVsDUR5GzR0dE0btyY06dPM23aNPr374+bm5vVYSmVLaQ2MY1TZ+C2F6/7COiVVltjzAJgAYCvr+8dXfMXcS0CgGIFit3J21UWdfToUUqVKoW7uzvz58+nfPnylC+frSbPU8pyzrwx7ARQKsFySftrcQoC1YBtInIUeBxY76wB42dXP8uzq591xqaVBaKjo5kyZQpVq1Zlzpw5ADRr1kyTgFJ3wJkF1ncClUSkHLYE0BXoFrfSGBMJxH89F5Ft2C5RDXBiTCoHCAoKws/Pj4CAADp06ECnTp2sDkmpbM1pPQJjTDTwOvAdsB9YbYzZJyLjROQpZ+1X5Wxz5syhVq1aHDt2jFWrVrFu3TpKlChhdVhKZWuOVB8VoDtQ3hgzzj5f8UPGmL/Seq8xZjO3laMwxoxOoW1jhyJWLimuSFy1atXo2rUr06ZNo1gxHe9RKiM4cmpoDhCL7SqhccBlYC1Q24lxKQXA1atXGTlyJO7u7nz44Yc0bNiQhg0bWh2WUjmKI4mgjjHmURH5B8AYc0FE8jg5rgz3qu+rVoeg0mnr1q288sorHDlyhP79+8f3CpRSGcuRRBBlv0vYQPx8BLFOjcoJulTrYnUIykEXL15k6NChLF68mEqVKrF9+3YaNGhgdVhK5ViODBbPBNYBD4jIe8CvwPtOjcoJwiLDCIsMS7uhstyZM2dYuXIlI0aMYPfu3ZoElHKyNHsExpjlIvI30BRbeYmOxpj9To8sg72w7gUAtvXaZm0gKllxH/4DBw7Ew8ODo0eP6mCwUpkkzR6B/Sqha8AGYD1w1f6aUnfNGMNnn32Gp6cnw4cP599//wXQJKBUJnLk1NAmbOWoNwFbgVDgG2cGpVzD8ePHadu2LS+88AIeHh4EBgZSqVIlq8NSyuU4cmqoesJle6E4nQle3ZW4InFnz55l5syZ9OvXT4vEKWWRdJeYMMbsEpE6abdUKqnQ0FDKlCmDu7s7CxcupEKFCpQtW9bqsJRyaY7cWTwkwWIu4FHgpNMicpI36r5hdQguLTo6mqlTpzJmzBgmT57MgAEDaNq0qdVhKaVwrEdQMMHzaGxjBWudE47ztPdob3UILiswMBA/Pz927drF008/zXPPPWd1SEqpBFJNBPYbyQoaY4ZmUjxOczDiIAAexTzuajsr/jzO14En0m6YhuBTl/B8uNBdbyer+/jjjxk8eDBFixZlzZo1WilUqSwoxauGRMTdGBMD1M/EeJymz8Y+9NnY566383XgCYJPXbrr7Xg+XIgOPo/c9XayKmNs8wd5e3vTvXt3goODNQkolUWl1iP4C9t4QKCIrAe+AK7GrTTGfOnk2LIsz4cLsapPXavDyJKuXLnCO++8Q+7cuZkyZYoWiVMqG3DkPoJ8wHls1UfbAe3t/yqVyPfff0+1atWYNWsWUVFR8b0CpVTWllqP4AH7FUN7sRWcS1j2Uf+Hq3gXLlxgyJAh+Pv74+Hhwfbt23niiSesDksp5aDUegRuwL32R8EEz+MeSgFw9uxZ1qxZw1tvvUVgYKAmAaWymdR6BKeMMeMyLRInG9lwpNUh5CinT5/m888/Z/DgwfFF4ooWLWp1WEqpO5BaIshRM4A0K9/M6hByBGMMS5cuZfDgwVy7do127dpRqVIlTQJKZWOpnRrKUbd9Bp4OJPB0oNVhZGtHjx6lVatW9OrVC09PTy0Sp1QOkWKPwBjzX2YG4myDvh0E6HwEdyo6OpomTZoQERHB7Nmz6du3L7lyOXLRmVIqq0t30TnlWkJCQihXrhzu7u4sWbKE8uXLU6ZMGavDUkplIP1Kp5IVFRXF+++/j5eXF7NnzwagSZMmmgSUyoG0R6CS2LVrF35+fgQGBvLcc8/RpUsXq0NSSjmR9ghUIjNnzuSxxx7j9OnTfPnll6xevZoHH3zQ6rCUUk7kMj2C95u+b3UIWZoxBhGhZs2avPjii0ydOpUiRYpYHZZSKhO4TCKoV6qe1SFkSZcvX+att94ib968TJ06lQYNGtCgQQOrw1JKZSKXOTW0I2wHO8J2WB1GlvLtt99SrVo15syZgzFGi8Qp5aJcpkfw9ta3Ab2PAOD8+fMMGTKEpUuXUrVqVX777Tfq1tWy2kq5KpfpEaj/d/78edatW8eoUaP4559/NAko5eKcmghEpJWIHBSREBF5M5n1Q0QkWESCRGSriOhF6k5y6tQppkyZgjGGypUrc+zYMcaNG0fevHmtDk0pZTGnJQL7fMezgdaAJ/C8iHje1uwfwNcY4w2sASY7Kx5XZYxhyZIlVK1alVGjRhESEgKgVwQppeI5s0fwGBBijAk1xtwCVgIdEjYwxvxkjLlmX/wDKOnEeFzOkSNHaNGiBX5+ftSoUYPdu3drkTilVBLOHCx+BAhLsBwO1EmlvR/wTXIrRKQ30BugdOnSdxTM9FbT7+h92VV0dDRPPvkk58+fZ+7cufTu3VuLxCmlkpUlrhoSkR6AL9AoufXGmAXAAgBfX987usbR5yGfOw0vW/n3338pX7487u7ufPLJJ1SoUIFSpUpZHZZSKgtz5lfEE0DCT6CS9tcSEZFmwDvAU8aYm84KZkvoFraEbnHW5i0XFRXFhAkTqFatGh9//DEAjRs31iSglEqTM3sEO4FKIlIOWwLoCnRL2EBEagLzgVbGmLNOjIUJ2ycAOXOmsoCAAPz8/AgKCqJr1648//zzVoeklMpGnNYjMMZEA68D3wH7gdXGmH0iMk5EnrI3+xC4F/hCRAJFZL2z4smpZsyYQZ06dYiIiODrr7/m888/54EHHrA6LKVUNuLUMQJjzGZg822vjU7wPOd9Pc8kcUXifH198fPzY/Lkydx3331Wh6WUyoayxGCxctylS5cYMWIE+fLlY9q0adSvX5/69etbHZZSKhvT6wmzkc2bN+Pl5cWCBQtwd3fXInFKqQzhMj2C+e3mWx3CHYuIiGDQoEEsX74cLy8v1qxZQ506qd2SoZRSjnOZROBRzMPqEO7YhQsX2LBhA2PGjOHtt98mT548VoeklMpBXCYRbDi4AYD2Hu0tjsQxJ06cYPny5QwbNoxKlSpx7NgxHQxWSjmFy4wRTP19KlN/n2p1GGkyxrBw4UI8PT0ZO3Yshw8fBtAkoJRyGpfpEWQHhw8f5pVXXuGnn36icePGLFy4kIoVK1odVpYSFRVFeHg4N27csDoUpbKkfPnyUbJkSXLnzu3wezQRZBHR0dE0bdqU//77j/nz5/Pyyy9rkbhkhIeHU7BgQcqWLYuIWB2OUlmKMYbz588THh5OuXLlHH6fJgKLHTx4kAoVKuDu7s6nn35KhQoVKFlSq3Gn5MaNG5oElEqBiFC0aFHOnTuXrvfpV06L3Lp1i3fffZfq1asze/ZsABo1aqRJwAGaBJRK2Z38/3CZHsGyp5dZHUK8v/76Cz8/P/bu3Uu3bt3o3r271SEppVyYy/QIShUuRanC1pdknj59OnXr1o2/N2D58uUUK1bM6rBUOtx77713vY2AgAAGDBiQ4vqjR4+yYsUKh9vfrnHjxnh4eFCjRg1q165NYGDg3YSbodavX8+kSZOsDsMpJk6cSMWKFfHw8OC7775Lts3WrVt59NFH8fHx4YknnoifPjbO2rVrERECAgIA+OGHH6hVqxbVq1enVq1a/PjjjxkfuDEmWz1q1apl7sTKPSvNyj0r7+i9CXWet8N0nrcj3e+LjY01xhjz22+/mT59+piLFy/edSyuKDg42OoQzD333OP0ffz000+mbdu2d/z+Ro0amZ07dxpjjFmyZIlp1qxZhsQVHR2dIdvJifbt22e8vb3NjRs3TGhoqClfvnyyP69KlSrF/x3Pnj3b9OzZM37dpUuXTIMGDUydOnXif3+7du0yJ06cMMYYs2fPHlOiRIk0Y0nu/wkQYFL4XHWZU0NzA+YC0KVal0zdb2RkJMOHDyd//vxMnz6devXqUa9evUyNIad6d8M+gk9eytBtepYoxJj2Xul+X2BgIH379uXatWtUqFCBJUuWUKRIEXbu3Imfnx+5cuWiefPmfPPNN+zdu5dt27YxZcoUNm7cyM8//8zAgQMB2/nd7du38+abb7J//358fHzo2bMnNWvWjG9/5coV+vfvT0BAACLCmDFj6NSpU4qx1a1blw8//BCAq1ev0r9/f/bu3UtUVBRjx46lQ4cOXLt2jV69erF37148PDw4efIks2fPxtfXl3vvvZc+ffqwZcsWZs+ezdGjR5k5cya3bt2iTp06zJkzBwA/P7/4mP73v/8xePBgZs6cybx583B3d8fT05OVK1fi7+9PQEAAH3/8MUePHuV///sfERERFC9enE8++YTSpUvTq1cvChUqREBAAKdPn2by5Mk8++yzqf4OXn31VXbu3Mn169d59tlneffddwEoW7YsAQEBFCtWjICAAIYOHcq2bdvS/XNMy9dff03Xrl3Jmzcv5cqVo2LFivz111/UrVs3UTsR4dIl299tZGQkJUqUiF83atQoRowYEf/7AqhZs2b8cy8vL65fv87NmzfJmzfvHcd6O5c5NWSFDRs24OnpyaJFi8ibN68WicvBXnzxRT744AOCgoKoXr16/IfQSy+9xPz58wkMDMTNzS3Z906ZMoXZs2cTGBjIL7/8Qv78+Zk0aRINGjQgMDCQwYMHJ2o/fvx4ChcuzJ49ewgKCuLJJ59MNbZvv/2Wjh07AvDee+/x5JNP8tdff/HTTz8xbNgwrl69ypw5cyhSpAjBwcGMHz+ev//+O/79V69epU6dOuzevZuiRYuyatUqfvvtt/hjWr58OYGBgZw4cYK9e/eyZ88eXnrpJQAmTZrEP//8Q1BQEPPmzUsSW//+/enZsydBQUF079490emvU6dO8euvv7Jx40befPPNNH8H7733HgEBAQQFBfHzzz8TFBSUantHfo6DBw/Gx8cnySO5U1snTpxINCNgyZIlOXEiyaSMLFq0iDZt2lCyZEmWLVsWf2y7du0iLCyMtm3bphjz2rVrefTRRzM0CYALDRZnpnPnzjFw4EA+//xzqlevzldffUXt2rWtDivHuZNv7s4QGRnJxYsXadTINuV2z549ee6557h48SKXL1+O/0bYrVs3Nm7cmOT99evXZ8iQIXTv3p1nnnkmzSvHtmzZwsqVK+OXixQpkmy77t27c+vWLa5cuRI/RvD999+zfv16pkyZAtguxz1+/Di//vprfK+kWrVqeHt7x2/Hzc0t/pvy1q1b+fvvv+P/nq9fv84DDzxA+/btCQ0NpX///rRt25YWLVoA4O3tTffu3enYsWN8Mkro999/58svvwTghRdeYPjw4fHrOnbsSK5cufD09OTMmTOp/kwAVq9ezYIFC4iOjubUqVMEBwcnOo7bOfJznDZtWpr7Ta9p06axefNm6tSpw4cffsiQIUNYsGABQ4YMwd/fP8X37du3jxEjRvD9999neEzaI3CCyMhINm/ezLvvvktAQIAmAZWqN998k0WLFnH9+nXq16/PgQMHMmS7y5cvJzQ0lJ49e9K/f3/ANia4du1aAgMDCQwM5Pjx41StWjXV7eTLly++N2OMoWfPnvHvP3jwIGPHjqVIkSLs3r2bxo0bM2/ePF5++WUANm3axGuvvcauXbuoXbs20dHRDsef8FtvWr3pI0eOMGXKFLZu3UpQUBBt27aNv/vc3d2d2NhYgHTfkZ6eHsEjjzxCWFhY/HJ4eDiPPPJIojbnzp1j9+7d8dWDu3Tpwo4dO7h8+TJ79+6lcePGlC1blj/++IOnnnoqfsA4PDycp59+mqVLl1KhQoV0HYMjNBFkkLCwMCZOnIgxhooVK3Ls2DFGjx6tlUJdQOHChSlSpAi//PILAMuWLaNRo0bcd999FCxYkD///BMg0bfPhA4fPkz16tUZMWIEtWvX5sCBAxQsWJDLly8n27558+bx956ArTptSkSE8ePH88cff3DgwAFatmzJrFmz4j9Y//nnH8DWK1m9ejUAwcHB7NmzJ9ntNW3alDVr1nD2rG2K8f/++49jx44RERFBbGwsnTp1YsKECezatYvY2FjCwsJo0qQJH3zwAZGRkVy5ciXR9urVqxf/c1m+fDkNGjRI8VjiVKlSJclrly5d4p577qFw4cKcOXOGb775Jn5d2bJl4091rV27Nv51R36O06ZNi096CR/Jnap66qmnWLlyJTdv3uTIkSP8+++/PPbYY4naFClShMjISA4dOgTYrgiqWrUqhQsXJiIigqNHj3L06FEef/xx1q9fj6+vLxcvXqRt27ZMmjTJaZNQuUwiWNN5DWs6r8nw7cbGxjJv3jy8vLyYMGFCfJG4woULZ/i+VNZw7do1SpYsGf/46KOP+PTTTxk2bBje3t4EBgYyerRtRtbFixfzyiuv4OPjw9WrV5P9u5g+fXr86ZjcuXPTunVrvL29cXNzo0aNGklOT4wcOZILFy5QrVo1atSowU8//ZRqvPnz5+eNN97gww8/ZNSoUURFReHt7Y2XlxejRo0CoF+/fpw7dw5PT09GjhyJl5dXsrF6enoyYcIEWrRogbe3N82bN+fUqVOcOHGCxo0b4+PjQ48ePZg4cSIxMTH06NGD6tWrU7NmTQYMGJCkeOKsWbP45JNP8Pb2ZtmyZcyYMSPVY4mIiEi2d1CjRg1q1qxJlSpV6NatW6IPzDFjxjBw4EB8fX0TjdOk9+eYFi8vLzp37oynpyetWrVi9uzZ8ftr06YNJ0+exN3dnYULF9KpUydq1KjBsmXLEg0MJ+fjjz8mJCSEcePGxfdI4hJxhknpcqKs+rjTy0czSsLLRw8dOmQaNWpkANO0aVNz+PBhS2NzBVnh8tH0uHz5cvzziRMnmgEDBlgYTcqio6PN9evXjTHGhISEmLJly5qbN29aHFVSGzZsMDNmzLA6jCxPLx9NgX+gPwC9fHplyPaio6Np3rw5Fy9eZPHixbz00kta+kAlsWnTJiZOnEh0dDRlypRJdTDQSteuXaNJkyZERUVhjGHOnDlZ8rRmu3btrA4hR9JEkE6XTh3l3gdK4u7uzrJly6hQoUKi64CVSqhLly506ZK5967ciYIFC8YPTCrX4zJjBHfr5s2bjBkzhu/Gv0DINtuAU4MGDTQJKKWyPZfpEdyNP/74Az8/P4KDgylTpxVl6rSyOiSllMow2iNIw9SpU6lXrx6XL1+23QTy0mjy3qtXBCmlcg5NBCmIuwGlbt269O3bl71799K6dWuLo1JKqYznMolgc/fNbO6+Oc12Fy9exM/PL/52+3r16jFnzhwKFSrk7BBVNlG2bFnAVio6f/78+Pj44OnpSd++feO/QDjq+eefx9vbO92lDC5evBhf7C0ulsaNGyfb9tChQ7Rp04ZKlSrx6KOP0rlzZ86cOcO2bdsy9Cqcl19+meDgYAC++OILqlatSpMmTdJdQltlPpcZIyiQu0Cabb766iv69evH2bNnGT58OMYYvSRUpapChQoEBgYSHR3Nk08+yVdffcUzzzyT5vuio6OJiIhg586dSerROyIuEfTr1y/Vdjdu3KBt27Z89NFHtG/fHoBt27aleypDRyxatCj++eLFi1m4cCFPPPEEAL6+vg5vJzo6Gnd3l/loyhJc5qc9Z6ft21O/2kn/45w9e5bXX3+dL774Ah8fHzZu3Mijjz6a2SGqO9DYv3GS1zp7daZf7X5ci7pGm+Vtkqzv5dOLXj69iLgWwbOrE5c23tZrW5r7LF68eJLX3N3dqVevHiEhIZw7d46+ffty/PhxwHbncP369Rk7diyHDx8mNDSU0qVLs2/fPk6cOIGPjw+zZs2iRIkSvPbaa5w7d44CBQqwcOFCqlSpwpkzZ+jbty+hoaEAzJ07l5kzZ3L48GF8fHxo3rw5AwYM4P77708S14oVK6hbt258EgDiew7btv3/sf71118MHDiQGzdukD9/fj755BM8PDzYt28fL730Erdu3SI2Npa1a9dSokQJOnfuTHh4ODExMYwaNYouXbrQuHFjpkyZwubNm/n111/x8/Pjqaeeom3btvEltFMqg+3v78+XX37JlStXiImJ4eeff07z96AyjsskgtX7bHVUkksEly5d4ocffuC9995j2LBh5M6dO7PDU9nIzp07k7x27do1tm7dyrhx4xg4cCCDBw/miSee4Pjx47Rs2ZL9+/cDtjo+v/76K/nz5+fo0aO0a9cuvjJo06ZNmTdvHpUqVeLPP/+kX79+/PjjjwwYMIBGjRqxbt06YmJiuHLlCpMmTWLv3r2JZh6Lq+KZ0N69e6lVq1aax1SlShV++eUX3N3d2bJlC2+//TZr165l3rx5DBw4ML6SaUxMDJs3b6ZEiRJs2rQJsBVZTGj06NH8+OOPTJkyBV9f30QJJ64M9pIlS7h48SKPPfYYzZo1A2xlmIOCgpJNaMq5XCYR3O748eMsW7aMt99+m4oVK3L8+HEKFixodVgqnVL7Bl8gd4FU1xcrUMyhHkBq4r6ViwgdOnSgdevW9OzZM/5cOdi+aMQVW3vqqafInz9/ku1cuXKFHTt28Nxzz8W/dvPmTQB+/PFHli5dCthKQhcuXDjVQnN3IjIykp49e/Lvv/8iIkRFRQG2iyXee+89wsPDeeaZZ6hUqRLVq1fnjTfeYMSIEbRr186hQnFxUiqDDbYicJoErOHURCAirYAZgBuwyBgz6bb1eYGlQC3gPNDFGHPUmTHFFYkbMWIEsbGxdOnShYoVK2oSUHckbowgodjYWP744w/y5cuXpP0999yT7HZiY2O57777MnxuYS8vL4dOs4waNYomTZqwbt26RAPP3bp1o06dOmzatIk2bdowf/58nnzySXbt2sXmzZsZOXIkTZs2jS+ylxZjL4Pt4eGR6PU///wzxZ+Ncj6nXTUkIm7AbKA14Ak8LyKetzXzAy4YYyoC04APnBUP2LrvjRs35rXXXqNu3brs27ePihUrOnOXygW1aNGCWbNmxS878uFeqFAhypUrxxdffAHYPjB3794N2E4ZzZ1rm2o1JiaGyMjIVMtUJ9StWzd27NgRfxoHYPv27ezduzdRu8jIyPja+QnrIYWGhlK+fHkGDBhAhw4dCAoK4uTJkxQoUIAePXowbNgwdu3alWYccVIqg62s5cwewWNAiDEmFEBEVgIdgOAEbToAY+3P1wAfi4iYuL+SDHQ04gphB4Nw+zsftV98h/vqtmHEd6eAU+naTvCpS3g+rJeSqpTNnDmT1157DW9vb6Kjo2nYsGGy0zTebvny5bz66qtMmDCBqKgounbtSo0aNZgxYwa9e/dm8eLFuLm5MXfuXOrWrUv9+vWpVq0arVu3TrGUcf78+dm4cSODBg1i0KBB5M6dG29vb2bMmEFERER8u+HDh9OzZ08mTJiQaKrE1atXs2zZMnLnzs1DDz3E22+/zc6dOxk2bBi5cuUid+7c8UnKEaNGjWLQoEF4e3sTGxtLuXLlkp21TWUuccJnrm3DIs8CrYwxL9uXXwDqGGNeT9Bmr71NuH35sL1NxG3b6g30BihdunStY8eOpTuedzfs4+ftv3Bv8UfIX7jYnR4WAB18HqFbndJ3tQ11Z/bv35/mjFpKubrk/p+IyN/GmGSv480Wg8XGmAXAAgBfX987ylxj2ntBFpnjVimlshJn3ll8AiiVYLmk/bVk24iIO1AY26CxUkqpTOLMRLATqCQi5UQkD9AVWH9bm/VAT/vzZ4EfnTE+oHIW/RNRKmV38v/DaYnAGBMNvA58B+wHVhtj9onIOBF5yt5sMVBUREKAIUDSGaGVSiBfvnycP39ek4FSyTDGcP78+WQvXU6N0waLncXX19foTEquKyoqivDwcG7cuGF1KEplSfny5aNkyZJJKiRk+8FipeLkzp2bcuXKWR2GUjmKy5ShVkoplTxNBEop5eI0ESillIvLdoPFInIOSP+txTbFgIg0W+UsesyuQY/ZNdzNMZcxxiSdTINsmAjuhogEpDRqnlPpMbsGPWbX4Kxj1lNDSinl4jQRKKWUi3O1RLDA6gAsoMfsGvSYXYNTjtmlxgiUUkol5Wo9AqWUUrfRRKCUUi4uRyYCEWklIgdFJEREklQ0FZG8IrLKvv5PESlrQZgZyoFjHiIiwSISJCJbRaSMFXFmpLSOOUG7TiJiRCTbX2royDGLSGf773qfiKzI7BgzmgN/26VF5CcR+cf+993GijgziogsEZGz9hkck1svIjLT/vMIEpFH73qnxpgc9QDcgMNAeSAPsBvwvK1NP2Ce/XlXYJXVcWfCMTcBCtifv+oKx2xvVxDYDvwB+Foddyb8nisB/wBF7MsPWB13JhzzAuBV+3NP4KjVcd/lMTcEHgX2prC+DfANIMDjwJ93u8+c2CN4DAgxxoQaY24BK4EOt7XpAHxqf74GaCoikokxZrQ0j9kY85Mx5pp98Q9sM8ZlZ478ngHGAx8AOaFutSPH/Aow2xhzAcAYczaTY8xojhyzAQrZnxcGTmZifBnOGLMd+C+VJh2ApcbmD+A+EXn4bvaZExPBI0BYguVw+2vJtjG2CXQigaKZEp1zOHLMCflh+0aRnaV5zPYucyljzKbMDMyJHPk9VwYqi8hvIvKHiLTKtOicw5FjHgv0EJFwYDPQP3NCs0x6/7+nSecjcDEi0gPwBRpZHYsziUgu4COgl8WhZDZ3bKeHGmPr9W0XkerGmItWBuVkzwP+xpipIlIXWCYi1YwxsVYHll3kxB7BCaBUguWS9teSbSMi7ti6k+czJTrncOSYEZFmwDvAU8aYm5kUm7OkdcwFgWrANhE5iu1c6vpsPmDsyO85HFhvjIkyxhwBDmFLDNmVI8fsB6wGMMb8DuTDVpwtp3Lo/3t65MREsBOoJCLlRCQPtsHg9be1WQ/0tD9/FvjR2Edhsqk0j1lEagLzsSWB7H7eGNI4ZmNMpDGmmDGmrDGmLLZxkaeMMdl5nlNH/ra/wtYbQESKYTtVFJqJMWY0R475ONAUQESqYksE5zI1ysy1HnjRfvXQ40CkMebU3Wwwx50aMsZEi8jrwHfYrjhYYozZJyLjgABjzHpgMbbuYwi2QZmu1kV89xw85g+Be4Ev7OPix40xT1kW9F1y8JhzFAeP+TughYgEAzHAMGNMtu3tOnjMbwALRWQwtoHjXtn5i52IfI4tmRezj3uMAXIDGGPmYRsHaQOEANeAl+56n9n456WUUioD5MRTQ0oppdJBE4FSSrk4TQRKKeXiNBEopZSL00SglFIuThOBypJEJEZEAhM8yqbS9koG7M9fRI7Y97XLfodqerexSEQ87c/fvm3djruN0b6duJ/LXhHZICL3pdHeJ7tX41TOp5ePqixJRK4YY+7N6LapbMMf2GiMWSMiLYApxhjvu9jeXceU1nZF5FPgkDHmvVTa98JWdfX1jI5F5RzaI1DZgojca59HYZeI7BGRJJVGReRhEdme4BtzA/vrLUTkd/t7vxCRtD6gtwMV7e8dYt/WXhEZZH/tHhHZJCK77a93sb++TUR8RWQSkN8ex3L7uiv2f1eKSNsEMfuLyLMi4iYiH4rITnuN+T4O/Fh+x15sTEQesx/jPyKyQ0Q87HfijgO62GPpYo99iYj8ZW+bXMVW5Wqsrr2tD30k98B2V2yg/bEO213whezrimG7qzKuR3vF/u8bwDv2527Y6g0Vw/bBfo/99RHA6GT25w88a3/+HPAnUAvYA9yD7a7sfUBNoBOwMMF7C9v/3YZ9zoO4mBK0iYvxaeBT+/M82KpI5gd6AyPtr+cFAoByycR5JcHxfQG0si8XAtztz5sBa+3PewEfJ3j/+0AP+/P7sNUiusfq37c+rH3kuBITKse4bozxiVsQkdzA+yLSEIjF9k34QeB0gvfsBJbY235ljAkUkUbYJiv5zV5aIw+2b9LJ+VBERmKrU+OHrX7NOmPMVXsMXwINgG+BqSLyAbbTSb+k47i+AWaISF6gFbDdGHPdfjrKW0SetbcrjK1Y3JHb3p9fRALtx78f+CFB+09FpBK2Mgu5U9h/C+ApERlqX84HlLZvS7koTQQqu+gOFAdqGWOixFZRNF/CBsaY7fZE0RbwF5GPgAvAD8aY5x3YxzBjzJq4BRFpmlwjY8whsc110AaYICJbjTHjHDkIY8wNEdkGtAS6YJtoBWyzTfU3xnyXxiauG2N8RKQAtvo7rwEzsU3A85Mx5mn7wPq2FN4vQCdjzEFH4lWuQccIVHZRGDhrTwJNgCRzLottHuYzxpiFwCJs0/39AdQXkbhz/veISGUH9/kL0FFECojIPdhO6/wiIiWAa8aYz7AV80tuztgoe88kOauwFQqL612A7UP91bj3iEhl+z6TZWyzzQ0A3pD/L6UeV4q4V4Kml7GdIovzHdBf7N0jsVWlVS5OE4HKLpYDviKyB3gROJBMm8bAbhH5B9u37RnGmHPYPhg/F5EgbKeFqjiyQ2PMLmxjB39hGzNYZIz5B6gO/GU/RTMGmJDM2xcAQXGDxbf5HtvEQFuMbfpFsCWuYGCX2CYtn08aPXZ7LEHYJmaZDEy0H3vC9/0EeMYNFmPrOeS2x7bPvqxcnF4+qpRSLk57BEop5eI0ESillIvTRKCUUi5OE4FSSrk4TQRKKeXiNBEopZSL00SglFIu7v8A4P8QUT/TYaMAAAAASUVORK5CYII=)
%% Cell type:markdown id: tags:
## Exercise
%% Cell type:markdown id: tags:
In this exercise you are going to construct a simple one vs rest multiclass classifier based on the probabilities for each of the binary problems for the three iris species. A one-vs-rest classification scheme constructs multiple models (one for each class) where each model predicts a binary classification and combines the results to come up with a single prediction.
For this problem, only use the first two features, i.e., sepal width and sepal length.
The steps for this multiclass classifier are:
1. Create a new features dataframe only using the first two columns of the iris features.
2. For each of the iris categories, model the data using Logistic Regression and calculate the probabilities for the testing set of the category being correct. Store these probabilites in an array. For each category you will need to re-do the train-test split, but make sure you use the same random state. Alternatively, split the data set initially with all three labels, and then work on each category. The necessary one-hot encoding has already been performed.
3. For each instance in feature space, choose the category with the highest probability. The function `numpy.argmax()` will be useful.
4. Plot the final classification for the testing set in feature space, with colours based on the predicted category.
Q: After completing the code below and inspecting the results, comment on which class(es) appear most difficult to predict accurately. Can you suggest any improvements to the model which might result in better accuracy?
A:
%% Cell type:code id: tags:
``` python
# Part 1
```
%% Cell type:code id: tags:
``` python
# Part 2
```
%% Cell type:code id: tags:
``` python
# Part 3
```
%% Cell type:code id: tags:
``` python
# Part 4
```
%% Cell type:code id: tags:
``` python
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment