{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Classification of Spam\n", "\n", "We'll now be exploring the field of Natural Language Processing (NLP), which concerns itself with interpreting, predicting, and classifying the written word. As a first foray into this field, we'll construct a simple spam classifier.\n", "\n", "Our goal in this project is to classify text messages as either spam or not spam (also playfully known as \"ham\" messages). We'll be using a collection of English SMS messages from [Kaggle](https://www.kaggle.com/uciml/sms-spam-collection-dataset/home) as our dataset. First thing to notice is that this dataset is NOT homogeneous! First off, the number of spam/ham messages is unbalanced: only (425+322)=747 spam messages, compared to (3375+450+1002)=4827 ham messages. The other thing is that this data doesn't all come from the same source:\n", "\n", "* [Grumbletext](www.grumbletext.co.uk): UK forum for people complaining about spam messages. (425 spam messages)\n", "* [NUS SMS Corpus](www.comp.nus.edu.sg/~rpnlpir/downloads/corpora/smsCorpus/): Dataset of legitimate messages collected by the Dept. of Computer Science at the National University of Singapore. Vast majority of messages are from students attending the university. (3375 ham messages)\n", "* [Caroline Tag's PhD thesis](http://etheses.bham.ac.uk/253/1/Tagg09PhD.pdf): Dataset collected by Caroline Tag during her doctoral studies on the linguistic aspects of texting. (450 ham messages)\n", "* [SMS Spam Corpus v.0.1 Big](http://www.esp.uem.es/jmgomez/smsspamcorpus/): Dataset (itself composed of other datasets) collected by researchers at the Universidad Europea de Madrid. (1002 ham and 322 spam messages)\n", "\n", "We'll need to take special care to make sure that the unbalanced nature of the dataset, and the various sources and nationalities of the senders of these messages does not affect the classification. As such, a simple classification accuracy measure will not suffice, as our classifier could still get 4827/(4827+747) = 87% accuracy by just specifying everything as ham and misclassifying every spam sample." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## EDA and data preprocessing\n", "Ok, let's dive in. First, we'll download our dataset and extract the CSV file therein (this requires installing the [Kaggle API](https://github.com/Kaggle/kaggle-api)):" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Warning: Your Kaggle API key is readable by otherusers on this system! To fix this, you can run'chmod 600 /home/ecotner/.kaggle/kaggle.json'\n", "Archive: sms-spam-collection-dataset.zip\n", " inflating: Data/spam.csv \n" ] } ], "source": [ "%%bash\n", "kaggle datasets download uciml/sms-spam-collection-dataset --quiet\n", "mkdir -p Data\n", "unzip sms-spam-collection-dataset.zip -d Data\n", "rm sms-spam-collection-dataset.zip" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's upload this using pandas and take a look at the data." ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
v1v2Unnamed: 2Unnamed: 3Unnamed: 4
0hamGo until jurong point, crazy.. Available only ...NaNNaNNaN
1hamOk lar... Joking wif u oni...NaNNaNNaN
2spamFree entry in 2 a wkly comp to win FA Cup fina...NaNNaNNaN
3hamU dun say so early hor... U c already then say...NaNNaNNaN
4hamNah I don't think he goes to usf, he lives aro...NaNNaNNaN
5spamFreeMsg Hey there darling it's been 3 week's n...NaNNaNNaN
6hamEven my brother is not like to speak with me. ...NaNNaNNaN
7hamAs per your request 'Melle Melle (Oru Minnamin...NaNNaNNaN
8spamWINNER!! As a valued network customer you have...NaNNaNNaN
9spamHad your mobile 11 months or more? U R entitle...NaNNaNNaN
\n", "
" ], "text/plain": [ " v1 v2 Unnamed: 2 \\\n", "0 ham Go until jurong point, crazy.. Available only ... NaN \n", "1 ham Ok lar... Joking wif u oni... NaN \n", "2 spam Free entry in 2 a wkly comp to win FA Cup fina... NaN \n", "3 ham U dun say so early hor... U c already then say... NaN \n", "4 ham Nah I don't think he goes to usf, he lives aro... NaN \n", "5 spam FreeMsg Hey there darling it's been 3 week's n... NaN \n", "6 ham Even my brother is not like to speak with me. ... NaN \n", "7 ham As per your request 'Melle Melle (Oru Minnamin... NaN \n", "8 spam WINNER!! As a valued network customer you have... NaN \n", "9 spam Had your mobile 11 months or more? U R entitle... NaN \n", "\n", " Unnamed: 3 Unnamed: 4 \n", "0 NaN NaN \n", "1 NaN NaN \n", "2 NaN NaN \n", "3 NaN NaN \n", "4 NaN NaN \n", "5 NaN NaN \n", "6 NaN NaN \n", "7 NaN NaN \n", "8 NaN NaN \n", "9 NaN NaN " ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "raw_data = pd.read_csv(\"./Data/spam.csv\", encoding=\"latin\") # Need to use latin encoding since UTF throws error\n", "raw_data.head(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, looks like the class (spam/ham) is under column \"v1\", and the actual text of the message is in \"v2\". What's in the other 3 unnamed columns? Is there even anything in them? Also, you can already tell from this snippet that the use of texting slang is pervasive throughout this corpus. Let's look take a look at the \"unnamed\" columns." ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
v1v2Unnamed: 2Unnamed: 3Unnamed: 4
95spamYour free ringtone is waiting to be collected....PO Box 5249MK17 92H. 450Ppw 16\"NaN
281ham\\Wen u miss someonethe person is definitely special for u..... B...why to miss themjust Keep-in-touch\\\" gdeve..\"
444ham\\HEY HEY WERETHE MONKEESPEOPLE SAY WE MONKEYAR...HOWU DOIN? FOUNDURSELF A JOBYET SAUSAGE?LOVE ...NaNNaN
671spamSMS. ac sun0819 posts HELLO:\\You seem coolwanted to say hi. HI!!!\\\" Stop? Send STOP to ...NaNNaN
710hamHeight of Confidence: All the Aeronautics prof...this wont even start........ Datz confidence..\"NaNNaN
\n", "
" ], "text/plain": [ " v1 v2 \\\n", "95 spam Your free ringtone is waiting to be collected.... \n", "281 ham \\Wen u miss someone \n", "444 ham \\HEY HEY WERETHE MONKEESPEOPLE SAY WE MONKEYAR... \n", "671 spam SMS. ac sun0819 posts HELLO:\\You seem cool \n", "710 ham Height of Confidence: All the Aeronautics prof... \n", "\n", " Unnamed: 2 Unnamed: 3 \\\n", "95 PO Box 5249 MK17 92H. 450Ppw 16\" \n", "281 the person is definitely special for u..... B... why to miss them \n", "444 HOWU DOIN? FOUNDURSELF A JOBYET SAUSAGE?LOVE ... NaN \n", "671 wanted to say hi. HI!!!\\\" Stop? Send STOP to ... NaN \n", "710 this wont even start........ Datz confidence..\" NaN \n", "\n", " Unnamed: 4 \n", "95 NaN \n", "281 just Keep-in-touch\\\" gdeve..\" \n", "444 NaN \n", "671 NaN \n", "710 NaN " ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raw_data[pd.notna(raw_data[\"Unnamed: 2\"]) | pd.notna(raw_data[\"Unnamed: 3\"]) | pd.notna(raw_data[\"Unnamed: 3\"])].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looks like the extra columns contain other messages or metadata? Some of it may be relevant to the classification (e.g. it looks like spam messages have some kind of PO box), but it seems like only a very small fraction of the dataset (50/5572=0.9%) has these, so I'm just going to drop these extra columns." ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [], "source": [ "raw_data = raw_data[[\"v1\", \"v2\"]].rename(columns={\"v1\": \"y\", \"v2\": \"msg\"})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's take a look at some summary statistics for the dataset - what is the class balance, number of messages, etc. We'll make use of the python package NLTK (Natural Language ToolKit) to simplify the analysis. First, we'll have to do some tokenization. I'm going to convert all messages to lowercase, then split tokens based on non-alphanumeric characters, discarding the punctuation itself. Also, I'll convert the \"spam\"/\"ham\" labels to binary identifiers (where \"spam\" --> 1)." ] }, { "cell_type": "code", "execution_count": 198, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ymsg
00[go, until, jurong, point, crazy, available, o...
10[ok, lar, joking, wif, u, oni]
21[free, entry, in, 2, a, wkly, comp, to, win, f...
30[u, dun, say, so, early, hor, u, c, already, t...
40[nah, i, don, t, think, he, goes, to, usf, he,...
\n", "
" ], "text/plain": [ " y msg\n", "0 0 [go, until, jurong, point, crazy, available, o...\n", "1 0 [ok, lar, joking, wif, u, oni]\n", "2 1 [free, entry, in, 2, a, wkly, comp, to, win, f...\n", "3 0 [u, dun, say, so, early, hor, u, c, already, t...\n", "4 0 [nah, i, don, t, think, he, goes, to, usf, he,..." ] }, "execution_count": 198, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import nltk\n", "from nltk.tokenize import RegexpTokenizer\n", "tokenizer = RegexpTokenizer(r\"\\w+\")\n", "\n", "data = []\n", "for _, row in raw_data.iterrows():\n", " y, msg = row\n", " y = 0 if (y==\"ham\") else 1\n", " msg = msg.lower()\n", " data.append([y, tokenizer.tokenize(msg)])\n", "data = pd.DataFrame(data, columns=[\"y\", \"msg\"])\n", "data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, before we go any further and start to make any more assumptions about how to process the messages, we should split this into training, validation and test sets. Then we'll blind the test set. Even though we haven't got to the typical training procedure yet, any more decisions regarding how to process the data will cross-contaminate the test set so we need to keep it separate. Since we only have a couple thousand data points, I'll use a 80/10/10 train/val/test split to get good statistics when evaluating the validation/test sets." ] }, { "cell_type": "code", "execution_count": 310, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Training set:\n", "Size: 4458, # spam examples: 601 (13.5%), # ham examples: 3857 (86.5%)\n", "\n", "Validation set:\n", "Size: 557, # spam examples: 66 (11.8%), # ham examples: 491 (88.2%)\n", "\n", "Test set:\n", "Size: 557, # spam examples: 80 (14.4%), # ham examples: 477 (85.6%)\n", "\n" ] } ], "source": [ "import numpy as np\n", "np.random.seed(0)\n", "shuffled_data = data.sample(len(data), random_state=0)\n", "\n", "# Determine the split\n", "M = len(shuffled_data)\n", "train_idx = int(M*0.80)+1\n", "val_idx = int(M*(0.80+0.10))+1\n", "X_train = shuffled_data[\"msg\"].values[:train_idx]\n", "y_train = shuffled_data[\"y\"].values[:train_idx]\n", "X_val = shuffled_data[\"msg\"].values[train_idx:val_idx]\n", "y_val = shuffled_data[\"y\"].values[train_idx:val_idx]\n", "X_test = shuffled_data[\"msg\"].values[val_idx:]\n", "y_test = shuffled_data[\"y\"].values[val_idx:]\n", "\n", "# Check split stats\n", "for name, X, y in [(\"Training\", X_train, y_train), (\"Validation\", X_val, y_val), (\"Test\", X_test, y_test)]:\n", " print(name + \" set:\")\n", " M = len(X); nspam = np.sum(y); nham = M - nspam\n", " print(\"Size: {}, # spam examples: {} ({:.1f}%), # ham examples: {} ({:.1f}%)\\n\".format(M, nspam, 100*nspam/M, nham, 100*nham/M))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's compute those summary stats. We'll gather frequency data on all the tokens in the training set." ] }, { "cell_type": "code", "execution_count": 311, "metadata": {}, "outputs": [], "source": [ "# Concatenate all tokens together\n", "all_tokens = []\n", "all_spam_tokens = []\n", "all_ham_tokens = []\n", "for x, y in zip(X_train, y_train):\n", " if (y==1): all_spam_tokens.extend(x)\n", " elif (y==0): all_ham_tokens.extend(x)\n", " all_tokens.extend(x)\n", "\n", "# Get frequency distribution of words\n", "from nltk.probability import FreqDist\n", "freq_dist = FreqDist(token for token in all_tokens)\n", "spam_freq_dist = FreqDist(token for token in all_spam_tokens)\n", "ham_freq_dist = FreqDist(token for token in all_ham_tokens)" ] }, { "cell_type": "code", "execution_count": 312, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total # of tokens with frequency of at least...\n", "1: 7790, 2: 3793, 3: 2535, 4: 1969, 5: 1645, 10: 914, 50: 220, 100: 118, \n", "\n", "Most frequent tokens (both spam/ham):\n", " i to you a the u and in is me my it for your of call s that have on \n", "2388 1800 1761 1161 1067 942 783 716 712 648 614 590 561 555 482 473 461 459 444 441 \n", "\n", "Most freqent spam tokens:\n", " to a call å you your the 2 free for now or txt u is on ur 4 have and \n", " 552 288 274 244 238 220 176 169 168 168 161 150 133 132 128 119 117 115 103 103 \n", "\n", "Most frequent ham tokens:\n", " i you to the a u and in me my is it that of for s so but can have \n", "2337 1523 1248 891 873 810 680 652 624 607 584 561 442 405 393 385 365 359 357 341 \n", "\n", "LEAST frequent tokens:\n", " mega asda counts toyota camry olayiwola mileage landing kane shud \n", " 1 1 1 1 1 1 1 1 1 1 \n" ] } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "print(\"Total # of tokens with frequency of at least...\")\n", "freq_count = {}\n", "for token in freq_dist:\n", " f = freq_dist[token]\n", " freq_count[f] = freq_count.get(f, 0) + 1\n", "freq_counts = [sum([freq_count[f] for f in freq_count if f >= n]) for n in [1,2,3,4,5,10,50,100]]\n", "print(\"\".join([str(n)+\": {}, \" for n in [1,2,3,4,5,10,50,100]]).format(*freq_counts), end=\"\\n\\n\")\n", "print(\"Most frequent tokens (both spam/ham):\")\n", "freq_dist.tabulate(20)\n", "print(\"\\nMost freqent spam tokens:\")\n", "spam_freq_dist.tabulate(20)\n", "print(\"\\nMost frequent ham tokens:\")\n", "ham_freq_dist.tabulate(20)\n", "print(\"\\nLEAST frequent tokens:\")\n", "FreqDist(dict(freq_dist.most_common()[-10:])).tabulate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looking at the number of tokens with a given count, we can see that the vast majority of distinct tokens only appear a handful of times throughout the entire corpus. This means that we can safely ignore tokens that have counts less than a small number (say, 5?) because they are not distinctive enough to discriminate against in classification.\n", "\n", "Looking at the most frequent tokens, we can see that \"stop words\" such as \"i\", \"to\", \"you\", etc. appear quite often; however, the most frequent spam tokens are not exactly the same as the most frequent ham tokens! For example, the token \"i\" is the most frequent ham token, yet it doesn't even appear on the top 20 list of spam tokens. Likewise, the token \"å\" is one of the top 5 spam tokens, yet doesn't appear in the top 20 ham tokens.\n", "\n", "But just because a token is frequent throughout the corpus doesn't mean it appear in most messages (i.e. it has a high \"document frequency\"). Let's take a look at the document frequency and token frequency of the union of the top 10 tokens in either set:" ] }, { "cell_type": "code", "execution_count": 313, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "'me' - spam document frequency: 3.7%, ham document frequency: 13.7%\n", "\tspam token frequency: 0.2%, ham token frequency: 1.1%\n", "\n", "'the' - spam document frequency: 24.0%, ham document frequency: 17.8%\n", "\tspam token frequency: 1.1%, ham token frequency: 1.6%\n", "\n", "'for' - spam document frequency: 24.5%, ham document frequency: 9.1%\n", "\tspam token frequency: 1.1%, ham token frequency: 0.7%\n", "\n", "'i' - spam document frequency: 6.0%, ham document frequency: 41.9%\n", "\tspam token frequency: 0.3%, ham token frequency: 4.1%\n", "\n", "'in' - spam document frequency: 9.3%, ham document frequency: 15.3%\n", "\tspam token frequency: 0.4%, ham token frequency: 1.1%\n", "\n", "'u' - spam document frequency: 16.3%, ham document frequency: 14.6%\n", "\tspam token frequency: 0.9%, ham token frequency: 1.4%\n", "\n", "'your' - spam document frequency: 31.1%, ham document frequency: 7.6%\n", "\tspam token frequency: 1.4%, ham token frequency: 0.6%\n", "\n", "'free' - spam document frequency: 21.8%, ham document frequency: 1.3%\n", "\tspam token frequency: 1.1%, ham token frequency: 0.1%\n", "\n", "'my' - spam document frequency: 1.2%, ham document frequency: 12.7%\n", "\tspam token frequency: 0.0%, ham token frequency: 1.1%\n", "\n", "'call' - spam document frequency: 42.3%, ham document frequency: 4.8%\n", "\tspam token frequency: 1.8%, ham token frequency: 0.4%\n", "\n", "'2' - spam document frequency: 21.6%, ham document frequency: 5.7%\n", "\tspam token frequency: 1.1%, ham token frequency: 0.5%\n", "\n", "'you' - spam document frequency: 32.4%, ham document frequency: 27.6%\n", "\tspam token frequency: 1.5%, ham token frequency: 2.7%\n", "\n", "'a' - spam document frequency: 36.6%, ham document frequency: 18.6%\n", "\tspam token frequency: 1.9%, ham token frequency: 1.5%\n", "\n", "'å' - spam document frequency: 32.3%, ham document frequency: 0.1%\n", "\tspam token frequency: 1.6%, ham token frequency: 0.0%\n", "\n", "'and' - spam document frequency: 15.0%, ham document frequency: 14.1%\n", "\tspam token frequency: 0.7%, ham token frequency: 1.2%\n", "\n", "'to' - spam document frequency: 61.9%, ham document frequency: 25.5%\n", "\tspam token frequency: 3.6%, ham token frequency: 2.2%\n", "\n" ] } ], "source": [ "_ = [dist.most_common(10) for dist in [spam_freq_dist, ham_freq_dist]]\n", "for t in {item[0] for sublist in _ for item in sublist}:\n", " spam_doc_freq = sum([(t in x) for x in X_train[y_train==1]])/len(X_train[y_train==1])\n", " ham_doc_freq = sum([(t in x) for x in X_train[y_train!=1]])/len(X_train[y_train!=1])\n", " print(\"'{}' - spam document frequency: {:.1f}%, ham document frequency: {:.1f}%\".format(t, 100*spam_doc_freq, 100*ham_doc_freq))\n", " print(\"\\tspam token frequency: {:.1f}%, ham token frequency: {:.1f}%\".format(100*spam_freq_dist.freq(t), 100*ham_freq_dist.freq(t)), end=\"\\n\\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that just because a token appears frequently within a class doesn't mean that it is helpful in discriminating classes. Take a look at \"u\": even though it's in the top most frequent tokens in both classes, its document frequency in both is about 15%; it will not be useful to discriminate between spam and ham. However, there are a number of tokens that appear primarily in one class and not the other, such as \"me\", \"i\", \"free\", \"my\", \"call\", and \"å\". We can find the most discrimnative tokens by looking at their \"discriminative ratio\" (I just made that up), defined by\n", "\n", "$$\\text{dr}(t,D) = \\ln\\left(\\frac{\\text{df}(t,D_\\text{spam})}{\\text{df}(t,D_\\text{ham})}\\right),$$\n", "\n", "where $t$ is the token, $D_i$ is a set of documents, and $\\text{df}(t,D) = \\{d \\in D: t \\in d\\}/|D|$ is the document frequency. Tokens with large $|\\text{dr}(t,D)|$ should have good discriminative power since the document frequency within one class is significantly larger, and the sign of this metric should specify which class the token is biased towards." ] }, { "cell_type": "code", "execution_count": 603, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in log\n", " \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Top spam features:\n", "('å', 6.028295894068239)\n", "('uk', 5.810293742254448)\n", "('www', 5.496636183399406)\n", "('nokia', 4.474009801709218)\n", "('txt', 4.392746837630453)\n", "('mobile', 3.7367519227018002)\n", "('cash', 3.521597761421069)\n", "('won', 3.2630439618923117)\n", "('stop', 3.007672732915791)\n", "('free', 2.802421714149846)\n", "('reply', 2.5941614033319977)\n", "('1', 2.3460649987452924)\n", "('text', 2.2468155546817834)\n", "('call', 2.181448533082571)\n", "('our', 2.12992497780842)\n", "('new', 1.8590500236730203)\n", "('or', 1.6686030594351484)\n", "('from', 1.625435172491515)\n", "('4', 1.6067073176665483)\n", "('your', 1.4031832841877871)\n", "('now', 1.4022916211773053)\n", "('only', 1.3890463944272846)\n", "('ur', 1.352488798693487)\n", "('send', 1.3414676174441496)\n", "('2', 1.3329569277762412)\n", "('for', 0.9886963869858907)\n", "('with', 0.9196517167887156)\n", "('to', 0.8863179808939132)\n", "('this', 0.8620479114156929)\n", "('on', 0.7266507374275272)\n", "('have', 0.6851157284663115)\n", "('a', 0.6762080009771575)\n", "('get', 0.6758803469768475)\n", "('t', 0.5706811997234378)\n", "('just', 0.48577103066519994)\n", "('is', 0.43262356112681216)\n", "('no', 0.4076163603821683)\n", "('of', 0.2988023554296916)\n", "('the', 0.29798569552313564)\n", "('are', 0.2835136629146011)\n", "('s', 0.20226599064172146)\n", "('you', 0.16225891233517728)\n", "('be', 0.15826233265069203)\n", "('u', 0.10896325084553317)\n", "('and', 0.06359369256363905)\n", "('we', 0.01464682095954321)\n", "\n", "Top ham features:\n", "('lt', -inf)\n", "('gt', -inf)\n", "('my', -2.3853552331248133)\n", "('but', -2.189832164472323)\n", "('ll', -2.062923312608294)\n", "('how', -1.9782494355591895)\n", "('ok', -1.9695913728160748)\n", "('i', -1.9451402769519104)\n", "('me', -1.317108071510026)\n", "('when', -1.2812634786920782)\n", "('that', -1.2583791848584907)\n", "('it', -1.0919463560395792)\n", "('so', -1.0225111045409898)\n", "('do', -0.9518575628688977)\n", "('m', -0.9409376727548754)\n", "('what', -0.903876787655024)\n", "('not', -0.8882208905824712)\n", "('at', -0.7619888004395601)\n", "('can', -0.7217979198563149)\n", "('if', -0.5776097925755947)\n", "('know', -0.544178595100713)\n", "('up', -0.5251150563134485)\n", "('in', -0.49402446924341725)\n", "('will', -0.1181126688863976)\n" ] } ], "source": [ "dr = []\n", "_ = [dist.most_common(50) for dist in [spam_freq_dist, ham_freq_dist]]\n", "for t in {item[0] for sublist in _ for item in sublist}:\n", " spam_doc_freq = sum([(t in x) for x in X_train[y_train==1]])/len(X_train[y_train==1])\n", " ham_doc_freq = sum([(t in x) for x in X_train[y_train!=1]])/len(X_train[y_train!=1])\n", " #idf = -np.log(sum([(t in x) for x in X_train])/len(X_train))\n", " try:\n", " dr.append((t, np.log(spam_doc_freq/ham_doc_freq)))\n", " except ZeroDivisionError:\n", " pass\n", "dr.sort(key=lambda e: abs(e[1]), reverse=True)\n", "print(\"Top spam features:\", *filter(lambda e: e[1]>0, dr), sep=\"\\n\", end=\"\\n\\n\")\n", "print(\"Top ham features:\", *filter(lambda e: e[1]<0, dr), sep=\"\\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The presence of the above tokens should then be the top 20 most discriminative features when trying to classify spam messages. We can see that 'å' is clearly the top feature, with a heavy positive bias, as predicted. It appears that most of the features with strong negative bias (\"my\", \"i\", \"me\") are first-person pronouns - obviously whoever is sending spam messages doesn't like to talk about themselves that much! It also appears as though there are two tokens (\"lt\" and \"gt\") which only appear in ham messages, and are thus completely correlated with the ham class. If you look at the messages which contain these tokens, it's obvious that they're some kind of transcription error; I'm pretty sure that they are corrupted versions of '<' and '>' (i.e. the less/greater than symbols '<' = '<' and '>' = '>'). The uncorrupted versions are entirely within the spam class, so I don't think it would be fair to use these tokens as classification features. On that note, it also appears that 'å' is also a corrupted version of the symbol for GBP (£), but since it's essentially a stand-in, I think it'll be fine to keep it." ] }, { "cell_type": "code", "execution_count": 604, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Occurrences of lt and gt:\n", "[['ham'\n", " 'Great! I hope you like your man well endowed. I am <#> inches...']\n", " ['ham'\n", " 'A gram usually runs like <#> , a half eighth is smarter though and gets you almost a whole second gram for <#>']\n", " ['ham'\n", " 'Do you know what Mallika Sherawat did yesterday? Find out now @ <URL>']\n", " ['ham' 'Does not operate after <#> or what']\n", " ['ham'\n", " \"Turns out my friends are staying for the whole show and won't be back til ~ <#> , so feel free to go ahead and smoke that $ <#> worth\"]]\n", "\n", "Occurrences of '<' and '>':\n", "[['spam'\n", " 'SIX chances to win CASH! From 100 to 20,000 pounds txt> CSH11 and send to 87575. Cost 150p/day, 6days, 16+ TsandCs apply Reply HL 4 info']\n", " ['spam'\n", " 'XXXMobileMovieClub: To use your credit, click the WAP link in the next txt message or click here>> http://wap. xxxmobilemovieclub.com?n=QJKGIGHJJGCBL']\n", " ['spam'\n", " 'TheMob> Check out our newest selection of content, Games, Tones, Gossip, babes and sport, Keep your mobile fit and funky text WAP to 82468']\n", " ['spam'\n", " 'Please CALL 08712404000 immediately as there is an urgent message waiting for you.']\n", " ['spam'\n", " 'RT-KIng Pro Video Club>> Need help? info@ringtoneking.co.uk or call 08701237397 You must be 16+ Club credits redeemable at www.ringtoneking.co.uk! Enjoy!']]\n", "\n", "Occurrences of 'å':\n", "[['spam'\n", " \"FreeMsg Hey there darling it's been 3 week's now and no word back! I'd like some fun you up for it still? Tb ok! XxX std chgs to send, å£1.50 to rcv\"]\n", " ['spam'\n", " 'WINNER!! As a valued network customer you have been selected to receivea å£900 prize reward! To claim call 09061701461. Claim code KL341. Valid 12 hours only.']\n", " ['spam'\n", " 'URGENT! You have won a 1 week FREE membership in our å£100,000 Prize Jackpot! Txt the word: CLAIM to No: 81010 T&C www.dbuk.net LCCLTD POBOX 4403LDNW1A7RW18']\n", " ['ham' 'Fine if thatåÕs the way u feel. ThatåÕs the way its gota b']\n", " ['spam'\n", " 'Thanks for your subscription to Ringtone UK your mobile will be charged å£5/month Please confirm by replying YES or NO. If you reply NO you will not be charged']]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:2: UserWarning: This pattern has match groups. To actually get the groups, use str.extract.\n", " \n" ] } ], "source": [ "print(\"Occurrences of lt and gt:\")\n", "print(raw_data[raw_data.msg.str.contains(r\"\\W(gt|lt)\\W\")].values[:5], end=\"\\n\\n\")\n", "print(\"Occurrences of '<' and '>':\")\n", "print(raw_data[raw_data.msg.str.contains(r\"<|>\")].values[:5], end=\"\\n\\n\")\n", "print(\"Occurrences of 'å':\")\n", "print(raw_data[raw_data.msg.str.contains(r\"å\")].values[:5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also need to make sure that there is some subset of these discriminative tokens $T_\\text{disc}$ such that every document $d \\in D$ contains at least one element of $T_\\text{disc} = \\{t | t \\in d \\, \\forall t \\in T_\\text{disc}, d \\in D \\}$." ] }, { "cell_type": "code", "execution_count": 861, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of tokens which completely covers entire dataset: |T_disc| = 68\n", "Top 20 tokens: å, uk, www, nokia, txt, mobile, cash, won, stop, free, reply, my, 1, text, but, call, our, ll, how, ok\n" ] } ], "source": [ "# Iterate over all tokens\n", "T_disc = []\n", "X_all = list(X_train) + list(X_val)\n", "for t, ratio in dr:\n", " # Remove 'lt' and 'gt' from consideration because that's probably cheating\n", " if np.isinf(ratio): continue\n", " # Test to see which document the token is in, and remove it from consideration\n", " for i in reversed(range(len(X_all))):\n", " if t in X_all[i]:\n", " del X_all[i]\n", " T_disc.append(t)\n", " if len(X_all) == 0:\n", " break\n", "print(\"Number of tokens which completely covers entire dataset: |T_disc| =\", len(T_disc))\n", "print(\"Top 20 tokens: \", end=\"\")\n", "print(*T_disc[:20], sep=\", \")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic regression and assessment of performance metrics\n", "As a simple first attempt at classification, let's just take a subset of the top discriminative tokens (except for 'lt' and 'gt') and use those as features in a linear classifier. First, we'll need a function to convert our documents to (sparse) feature vectors." ] }, { "cell_type": "code", "execution_count": 862, "metadata": {}, "outputs": [], "source": [ "import scipy as sp\n", "from scipy.sparse import dok_matrix\n", "\n", "token_to_idx_dict = {T_disc[i]: i+1 for i in range(len(T_disc))} # Reserve 0 index for unknown token\n", "X_train_sparse = dok_matrix((len(X_train), len(token_to_idx_dict)+1), dtype=np.float32)\n", "X_val_sparse = dok_matrix((len(X_val), len(token_to_idx_dict)+1), dtype=np.float32)\n", "X_test_sparse = dok_matrix((len(X_test), len(token_to_idx_dict)+1), dtype=np.float32)\n", "for X, S in [(X_train, X_train_sparse), (X_val, X_val_sparse), (X_test, X_test_sparse)]:\n", " for i, x in enumerate(X):\n", " for t in x:\n", " S[i, token_to_idx_dict.get(t, 0)] = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll import a logistic regression classifier from scikit-learn and fit it to that. Since there is a class imbalance in the dataset, we'll rebalance it by assigning extra weight to the spam class (the ratio of spam:ham is about 6:1)." ] }, { "cell_type": "code", "execution_count": 863, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LogisticRegression(C=100.0, class_weight={0: 1, 1: 6}, dual=False,\n", " fit_intercept=True, intercept_scaling=1, max_iter=100,\n", " multi_class='ovr', n_jobs=1, penalty='l2', random_state=None,\n", " solver='liblinear', tol=0.0001, verbose=0, warm_start=False)" ] }, "execution_count": 863, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.linear_model import LogisticRegression\n", "\n", "model = {}\n", "model[1] = LogisticRegression(penalty=\"l2\", C=100.0, solver=\"liblinear\", class_weight={0:1, 1:6})\n", "model[1].fit(X_train_sparse, y_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's see how our classifier performs. First, let's look at some standard performance metrics across all the classes such as ROC, precision/recall, and $F_1$ score." ] }, { "cell_type": "code", "execution_count": 864, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAEXCAYAAAA+1KxEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd81fX1+PHXuTeLJCQhJGGFsKeIgIgbQcFdZ93WWf3aqq221tGqVWvV+mvraLW1Ks4ijlZFiwMVVARlKEOQvQkjg5Cd3HF+f3w+wCUkJIR7c2/CeT4e93Hv/ax7btD3/Zz3FFXFGGOMMcYYY0xs8kQ7AGOMMcYYY4wxDbOkzRhjjDHGGGNimCVtxhhjjDHGGBPDLGkzxhhjjDHGmBhmSZsxxhhjjDHGxDBL2owxxhhjjDEmhlnSZowxxhhjjDExLC7aARhjjDHGGNOWiEgHoBioCNm8A8hVWyTZNIO1tJmoE5G1IlIlIuUiskVEXhSR1JD9x4jIZyJSJiI7ROQ9ERlc5xppIvK4iKx3r7PKfZ/V8t/IGBNLRGSxiIxp5Jg8t+zwtlBYYSci94nIq+7rniKiImKVs8ZExzCgWFVTQx7dYiFhs3KhdbKkzcSKH6lqKk4hNxy4C0BEjgY+Bt4FugK9gAXAVyLS2z0mAfgUOAQ4FUgDjgaKgFEt+zWMMfujTqXN1rqVNuGgqoeo6vRGjlnv3lQFwvGZItJVRDaG41rGmFZpGLCkuSeLyB0issmtsF4mIie527uLyH9FpEBEikTk7yHnDBKR6SJS4lZWnRWyb617zYVAhYjEueXUf9xrrRGRXxzIFzaRZUmbiSmqugX4CKewA3gUeFlVn1DVMlUtVtW7ga+B+9xjrgDygHNVdYmqBlV1m6r+QVWntPR3MMbst52VNiOAkcDdoTvF0dp+r04HPqy70Wq4jTloDKeZSZuIDABuAo5Q1fbAKcBatyfA+8A6oCfQDZjknhMPvIdT0Z0D3Az8273WTpcAZwAZQNA9foF7nZOAW0TklObEbCKvtf0ImjZORHKB04CVIpIMHAO8Wc+hbwDj3dfjgA9VtbxlojTGRIKqbgI+AIa4tcV/FJGvgEqgt4iki8jzIrLZrYF+MLQ7o4hcJyI/uDXTS0RkhLt9rYiMc1+PEpG5IlLqtuz91d2+R3dCtwZ6sogUi8hKEbku5HPuE5E3RORl97MWi8jIOl/ndGBKyOc3uYZbRLwi8lu3m3eZiMwTke7uvidEZIMb/zwROT7s/xDGmHAYBvzEbfUqEZFFbhk22+1ZMGQf5waARGCwiMSr6lpVXYXTe6gr8BtVrVDValWd4Z5zFJAKPKKqtar6GU6Cd0nIdZ9U1Q2qWgUcAWSr6gPu8auBZ4GLw/pXMGFjSZuJFe+ISBmwAdgG/B7IxPlvdHM9x28Gdo5X69jAMcaYVsRNTE4HvnM3/QS4HmiPU7P8IuAH+uLUYp8M/NQ99wKc1vcrcLpIn4XTRbquJ4AnVDUN6INTAVSfScBGnBukHwMPiciJIfvPco/JACYDoV2U4oHRwNSQ4/enhvtX7vGnu9/lGpzEFWAOzs1gJjAReFNEkhr4DsaYKBCRRGAQMFpVM9zHoTj/H58BvLWv81V1JXALTpm2TUQmiUhXoDuwTlX99ZzWFdigqsGQbetwypidNoS87gF0DUkqS4DfAp3257ualmNJm4kV57hdAMYAA3ESsu04Nzdd6jm+C1Dovi5q4BhjTOvwjnvDMAP4HHjI3f6iqi52b1AycZKYW9wa5m3AY+yuFf4p8KiqzlHHSlVdV89n+YC+IpKlquWq+nXdA9zk8VjgDrcmez7wHE5CuNMMVZ3ijoF7BTgsZN9oYIGqloVs258a7p8Cd6vqMve7LFDVIgBVfVVVi1TVr6p/wamND+3+ZIyJviE49y8LQzeqqk9VC5pyAVWdqKrH4SRXCvwJJ+nKa6CbdT7QvU5X8jxgU+hlQ15vANaEJJUZqtpeVU9vSnym5VnSZmKKqn6OU5v+Z1WtAGYBF9Rz6IU4k48AfAKcIiIpLRKkMSbcznFvGHqo6s/dxAb2rhWOBzaH1Ao/gzN2A5wa6FVN+Kxrgf7AUhGZIyJn1nNMV5xZ30KTrro11ltCXlcCSSE3Uru6RobYnxruBr+LiNzmdgHd4Z6Xzu5eB8aY2DAc+F5Va5tzsogMEJET3Ra7aqAKJwmcjdOz6BERSRGRJBE51j3tG5yy6HYRiRdnxtwf4Y55q8dsoMztut3O7ZY9RESOaE7MJvIsaTOx6HFgvIgcBtwJXCkivxCR9iLSQUQexJkd8n73+Fdwboj+IyIDRcQjIh3dMSFWY2RM61W3VrgGyAqpFU5T1UNC9vdp9IKqK1T1Epxk70/AW/VU+OQDmSLSPmRb3RrrfakvadufGu56v4s7fu12nEqrDqqagbPukzQxLmNMyxgGzD2A8xOBR3B6FG3BKa/uclv2f4TTRXw9ThfuiwDcBPFHOPMCFAJPA1eo6tL6PsC91plurGvcc57DqQgyMchmsTIxR1ULRORl4F5VPd8d5/EgTpepIPAlcJyqrnCPr3EnGbgfZwxJB2ArzjIB30TjOxhjwktVN4vIx8BfROQeoBxnCZBct4X+OeCvIjID+BYn6fHV7SIpIpcDH7nlTIm7OXQMCKq6QURmAg+LyG04LXPXApc1FqeI9AISVfWHfRy2q4YbeBKoxRn/0k5V57jf5Q8isgRYCRyKkzC2xxnTVwDEicidOGPejDExRFVvOsDzF9LAkkWquh44p4F9i4ETGtjXs55t+ew5UYmJYZa0mahroCD5WcjrGThj3fZ1jR04g3ZvCXN4xpjYcQVO7fMSnARmNU5rGar6poh0xJmcoxuwFmcik7rj2k7FSe6S3X0Xq2qVyF6NVZcA/8RpddsO/F5VP2lCjGewdyvbHlQ14HbL/AtODXcisIzdSx381d32MU7Xx6XAuTjLoXwILAcqcMb0bcAY02qIyBSc1q0BIvKMqr4Y5ZBMKyEa/YXZjTHGmDbBvSH7u60RaYwxJpwiOqZNRCaIyDYR+b6B/SIiT7pr4CwUd00dd9+VIrLCfVwZsv1wd62Lle651pffGLNfIlE2GeOaDkyLdhDGGGPalkhPRPIiTleUhpwG9HMf1wP/ABCRTJx1uo7E6dP7exHp4J7zD+C6kPP2dX1jjKnPi4S/bDIGVX00ZPZLY4wxJiwimrSp6hdA8T4OORt42V2H5msgQ0S6AKcAU1W1WFW340wucaq7L01Vv1anX+fLNDAY0xhjGhLusinyERtjjDHmYBbtiUi6secg6o3utn1t31jP9r2IyPU4NeSkpKQcPnDgwPBFbUy01JRB0cpoR9GoioRsUrJyGz1u3rx5haqa3QIh7a/9LZv20qwyqHgN+GsoaNeLLaXVZKcm4vE0rwd4YpyH9HbxzTrXmINJDJdDYZOVlaU9e/aMdhjGmHo0tQyKdtIWMar6L+BfACNHjtS5cw9kuQzTWgWDyuOfrmB7RbPWt4w5/XfM5CdrfsMrvR5lS7t+TT7vnfmbOOuwrlx1bM/IBRciuX0H2qdnNnqciNSd2a/NaFYZ9NJZ4K/m3cNf4JeT5h/Q5/uAU4/KQ9wlvNoleLlxbF9L5Iypoy2XQzv17NkTuw8yJjY1tQyKdtK2Cege8j7X3baJPad4z8UZ3L3JfV33eGPqta64kic/XUFKgpfEeG+0wzlgO7SYnwAfrQuyRAJNPs+T3JnBAwfRKbdr5IJrW/a3bAqPmjJIzuTsYd0449Auzb7MO/PzuX/yYqYs2gJAsVtpMfGb9RzXN6vZ101JjOOGE3rTr1P7xg82xhhjTNhEO2mbDNwkIpNwBvbvcBdQ/Qh4KGSA/8k4K8EXi0ipiByFs2jyFcDfohK5iZp352/ipZlrm3Rstc9ZM/eh8w7l7GH19mJrXZYHYSK8eu0o6HZ4tKNpy/arbArbp9aUQYceAMR5mz/k+MeH5/Ljw3fXb20rreaGV+dRWu1nTWFFs6+bX1LF5AWb+NmYvtw4tg+Jca2/IsQYY4xpDSKatInIazi10lkishFn1rV4AFX9J84CpKcDK4FK4Gp3X7GI/AGY417qAVXdOWnAz3FmfmsHfOA+zEFk6pKtfLu+hOP7Nd5ikJII4wblMCLPJvgzu0WobDpwNWWQGP5WrJy0JP7782MP+DqF5TU8+P4Snvx0Be8vzOeR84Yyqlfj3WCNMcYYc2AimrSp6iWN7Ffgxgb2TQAm1LN9LjAkLAGa/aaqPPHpCp79YjX+YHQWZvcFgvTNSeWVa4+Myueb1i8SZVNY1JRBYlpELh0OWamJPH7xcM4dkcvv3l7Ehc/M4pJRedwyrh/tEqLb6pYY57GWP2OMMW1WtLtHmlYkEFTuefd7Jn6znnGDOtEnJyVqsYzqabX7po0JBsBXEZGWtnA7oX82H986mic+WcFzM9bw2uz10Q4JgMk3HUtKYv0/a13T20U9sTTGGGOay5I20yQ1/gC3vj6fKYu28LMxfbj9lAGING8qcmNMPWrKnOeyLVC4omU/u10HSNm/CUqSE+K46/RBnDO8G1+tLIxQYE3z3283sWRzKWf9/asGj0mK9/DCVaMASE7wMjQ33cowY4wxrYYlbaZR5TV+rn95LjNXFXH3GYP46fG9ox2SMW1P0J0NdN4LzqMleRPg9jWQmLrfpw7qksagLtHt0nn5UT2YtnQbtYFgvft/OWk+1b4glzz79a5tvx7fnx5ZTm+BeI9wwoBskhPsJ9EYY0xssl8os0+F5TVc/cIclmwu5a8XHsZ5IxpfMNkY0wwpHeGKyVBR0LKfu+Rd+GEyeFrvz0FSvJfT9rFEwvH9slm2xWnJVJRfTprPX6Yu3+OYh887lEtG5UU0TmOMMS2gtgJ8VXv1IPl+0w6WbSnjvBHdWmVPi9b7K20ibkNxJVdMmM3mHVU8e8XhnDiwU7RDMqZt631Cy3/mhtmQmA7xSS3/2S0kMyWBo/t03PV+6q2jKSx31q7bWlrNZc99w5RFm1lb5CyHEO/xcM1xvchMSYhKvMYYY/YWDCpriirIL6kC4Ng+WYhAaXkV6dsXwerpsOZz53ct6KM8aygleSezLnssr6xI4sMlW7k9bhLff7Ca7R2G0vmQE+g/8iRIzY7uF2siS9pMvZZtKeOKCd9QVRvg1WuPZKRN/GFM21SxrdX8YIVLRnICGclOQpaTlkiX9CRmrylm9hpn9YYaf5AvVxbSJ6vhyZay2ydy+6kD8XpaX22tMcbEElVlXVElH3y/heVbyxgzIJv5G0rwilBR6+e79SWsLqggKzWB/B3Vu84bml7FzZVPc7RnMUg1QYRV3j7Mjz+LteVexm37luGFfyaXP5NLF67oPJJjSt5joy+LAdt+IKFgIkyHHUm5fO8dyNp2h3DRuecT12UIeGJv4ipL2mLc299t5L/fbmrxz52/oYTkBC9v3nAMAzrH/mx2xphmKi+AlJxoRxE1aUnxzLrrpD22PTzlB6Z8v5niipp6zwkElPwd1cxeW0xKQhzD8zL49ckDWiJcY4xp9Wr9Qd5fmM+iTTtYurmMpVtK2V7p27X/7e923/fGeQR/UMlMSWBglzRuPqkfvbNSmLJoMzXz32S8dx7vBo+jqNs4Jmzqjs+bQUZSAh2zEhhy1AN8WrWVnoWf07NwGj3WfoCmdaPrz7+h3Kf89+OPWP3tp4wIrORwzxyOrfgEnn2C7XE5LD7zXY4bNjgaf54GWdIW4/4zbxPfrd/e4onTqJ6Z3HfWIXTPTG7RzzXGtLCKAsi2hCPUXacP4q7TBzW4PxhUfvPWQtYUljNjZSEzVhYyd+12Xrj6CJLiY6921hhjokFVmbW6iKLyWk4/tAsCTPl+MzdN/A6AdvFeBnRuz6lDOnNI13SG5qaTlhRPjT9I5/Qk0pLiqA0E612D88jeHaHbHHgfzr71KejQk2vqjaILMAz4JVTvQDSIJLUnLQkuPv8Cas8+n62l1SQnx3PdM++QtW0WD/M80958ilv/dy6F5TV0SUviwXOH1D9MKBhwumRqEPqOC98frx6WtMWQ/JIqLvrXLCpqAru2lVb5GNUrk4nXHRXFyIwxbVbFNuh1fLSjaFU8HuEvFx4GwPRl27jqhTnMWl3EwHs+pENyPP6A8v8uGMr4wZ3xCK1ywLsxxoQqq/axrqiSzTuqWVNYTmVtgPR28YzI60BRRQ0/bC6jqLwWrwe2lNawrbSajdur2OSOP7v5te/2uN4dpw7k+tG9G+1iXl/CRjAIG76G+f92eoqkN3ESqaT0vTYlxHl2NVD87efnUVL5IwKvzeWmspnsyLuWH7ZWUJK/igmf+MntcAx9slOdmItWwfyJsOA1KHVbBs9+GoZf1rRYmsGSthhRXuNn/oYSNhRXMW5QJ7qk754U4ORDbAIQY0wEBHxQtf2g7h55oMYMyGHGHWN5fsYa/AGlsjbAf77dyA2vfgvA2AHZPHTeoWSnJhLn9UQ5WmOMaVwgqKwqKGfu2u28vzCfytoAy7aUUeUL7PM8j0Cc10PX9CRy2icxrHsGVx7Tg4qaANOXF5DgFS4/qgc/GtoVT91kraYMilZC4UooXA6+Shh8NuQeATsrvjZ9CwsmOTMel20GbyKMvx884Slbk+K9dE73wvG30uHNq/hz34XQF/zv30Z5QQJPPnkerwTG88WAt+iybjKIB/qcBCf/Ab59BSbfBPHtYMh5YYmnLlHViFw4lowcOVLnzp0b7TAaVO0LMOqPn1Ba7QfgzRuO5ojWNPHHC2dAwQ/RjuLg4K+F2jK4bhp0GxHtaA6YiMxT1ZHRjiPSYrYMKtsKf+kP8cnOD40Ji2p/EA0qlSE3OF4RUpPi8HqE/W538ybAjydAj2PCGqdxHAzlUMyWQSamLM7fwZtzN/LG3A1U1jrlV7t4L72yUjikaxonDsyha0Y7unVoR0Kch8KyGpZvLadjagIDOrcnMc6DV6ThCqpgAHZscBKzohVOcla4wknWyjbvPk48zjI0gVrI7A1DL4buR8C/LwDxQr/xcMi50P8USIzA8CFVeOE0WD8LgJLOx7JyWzkjgwtYF8yhh2cbs1LHc/T1f4M0d7mZ2kqYcLJzn3bT7P36uKaWQdbSFmVz1xazfGs5pdV+LhmVx9F9OjK8e0a0w9o/676CLodBbpv+zYsdiWnQaUi0ozBtQWoOjPmt00XShM3OfhJVFbUUldcwa7UzKyXl0CsrheP6ZtHkHpOVRbD4bdixMRKhGmPaqGBQ8QXrHw+2k6qyfGs5s1YVMvWHrXy1sgiAUw/pzPjBnRiam07fnNQGu3inJcXTOzt17x3VO0ISsxUhz6sgEDLBU1I6dOwHvcdCVl/ndVY/J1Hz1zgtagsmwfSHdp9z09fOMZEkAuc/Bx/cAZ2GkHHC7YwUD6yYSu5Hv4OibTxYNJYjphdz31lu0paQDDmHwPqZEQvLkrYoWldUwY//OWvX+0tH5XFo7t79bVuFfifDib+LdhTGmP0hAmPuiHYUbVam+0gqruTb9dv55aT5sAWu7tOT/jntEWDswBw6pe1jjbw1XzpJW6p1kzfGNKyy1k9SnJfaQJAPvt/MY1NXsL64kqzUBLJSE0mI81Be40eA0f2zSU2M42+frdx1fm6Hdlx3fC+uPa43ndObuW7nV0/ArKegfOvubeKFDj2dRKvPic5zVn8nQUvJosEarLhEGH658yjZAIvecBbN7ti3ebHtr/RcuPjfe27rfzLePidSW7CC3tOqeXHmWm4c25fs9onOfm+8M+wgQixpi6IdVc4/7D1nDmbMgGz61FdbYYwxplXrnplM98xkanxBbv/PQl74au0e+28d17/Bcwdsm8+pwMvf17J99YpGPyszNYHLj8yzyU+MaaOqfQG+WlnItrIaKmr8rCmsYO7a7SzbWlbv8QM7p1FUUUtmSgLdM5NZsKFkjzLonjMHc/LgTgc+W7i/1pmYIz4Zxt23u9WsQy+ISziwa2d0h+N/fWDXCBdvHAmdB3HJqELeW5DPEX/8hLWPnOHsa9fBGSeu2nAyegAsaYuiWn8QgP6dUi1hM8aYNu7CI7ozbnCnXWX/GU9+SVFFLY99srzBc67zLuHUeHh0ZgnlNHxcqF4dUziuX1ZYYjbGtKzyGj8bt1eyuqCCpZtLWVNUSbxHWF9cybItZZTVOPMfJFJLNymkf0Ixl2aW06NbAVKyjj4JReQEtxHXsReeIefBIYc5rUYuX8BZI+2YHqnklC9FNn0AW7pBh7P2P9FQhdJ8+O5VmPu808J24t1w3K3h/JPEpKN7d9z1ekNxpZP0pmSBv9qZVCUpLeyfaUlbFPgCQe5993tWbC0HGpjO1BhjTJuTmbK7xnnu3eNobC4w+WgG+m0KC+85v9EbqhkrC7liwmwuf/4bzh3ejXvPHEyHlAOs4TbGRFQwqDw/Yw2TF+SztrBiV1IGTmLWS7YwIqWQ8xK3kuZZR5/2hfSKKySpKmQs8nacmRQz8qBDD0g/GjYvgI/vdh7dj3JmNGzXgfiNczh341x4bxEEQ7ry9T8NfvQEtO/kJGPVO6BsizNBSPlW53nn+9DnQK1zft/xcNQN0PvElvnDRZmIMOUXx3Pxv2ZxzlNf8fxVRzAsw116oGQddD407J9pSVsULNtSxmuzN9A9sx1H9OxA3xxrZTPGmIONiDResV2xFdp3QpqwXMBxfbMYN6gTS/J38PZ3m7j4iO7OArTGmJgRdKfT/25DCfM3lPDlD5vwlm1kRGoJP+teQXfdRCffBjIq1xNfvglBwY/zSOvmTNKRMcJJzjJ67H5O7bT31PdFq5wxsd//Fz643dkWn+LMPn3MTdBtpPN68dvw6QPw9yMgOdNJxvxVewefmAbtOzuPvKPd112cRaUjPTlIDBrcNY23bzyWS5/9mj99sJTXzuzl7CheY0lba3bPO9/z2VKnVqTG70yj+vSlh7feiUeMMcZEXtkW56aoCTwe4bkrR/LOd5u45fX5dExNjHBwxpimCNZWs2rJHGbMnsuO/OV08m+mh2zlWE8BD0oRnsQg+ICNQEKqM9lGr6PdcWHurIod+0BCyv59cMc+MPo251G4wpmRMWcQeOr08Dr6RqelbPrDTot++y67E7Kdz6mdINEaGerqk53KiLwOLN1SBonZzkZfPQlvGEQ0aRORU4EnAC/wnKo+Umd/D2ACkA0UA5er6kYRGQs8FnLoQOBiVX1HRF4ETgB2uPuuUtX5kfweB6rWH+SNuRvIy0xmaK4znX+H5HgGdYnA2hLGmEY1t2xy9wWARe6h61X1rBYL3Bx8yjY7M6ZNuqz+/UkZcPqje9zMlVQ63ZXG/fVzfnFSP341vuGJTowx4bVsSxmbSipZsbWcL1cUkhwPv1h7I0N0BTvboqqTM6FDTxKzxyCZvSCzlzPDYodeTpIUiYmEGmsJy+4PF7wQ/s89CPgCQdYUVrAwv5yhAEF/Y6c0S8SSNhHxAk8B43HqDuaIyGRVXRJy2J+Bl1X1JRE5EXgY+ImqTgOGudfJBFYCH4ec9xtVfStSsYfb2D9Pp8Yf5MeH5/J/J/SJdjjGHNQOpGxy91Wp6rAWDdocvPqOh7UznO42dVVth7J8p5a80+Bdm08Z0pnC8lr+Pm0lT366wpI2YyLs0x+2MnlBPu/OzyeNCvrJRgZ4NjJe1jM8bi1DWMH3A2+h+5Fnk961H0mRWBDaRM1vThnIJz9s45GPVjARQAMR+ZxItrSNAlaq6moAEZkEnA2E3hgNBn7lvp4GvFPPdX4MfKCqlRGMNWJq/UE2lVTRMSWBi47oHu1wjDHhK5uMibzTH21439f/gA/v3GsNty7p7bjtlAE888UqfAFlyqLNnH5o07pYGmMaV1nr58WZa3n0w2UM7NyetK2z+VncZG5P3Eg3Kdp1nCa0d7oj9rmTIWPujEwLmom6AZ3b89hFh/HI69MgCQhGJmlrfGRz83UDNoS83+huC7UAOM99fS7QXkTqjpq+GHitzrY/ishCEXlMRGK60/6/v1kHwDXH9SIj2WbxMiYGHGjZlCQic0XkaxE5J7KhGrMPZVvAE++sDVSPB84eAsDP//0tReU1LRmZMW1SrT/IU9NWcsSDn/Doh8sA6B1fzEspTzKqXT6dDz3RWaPs0jfglkXIXRuQn06FsXdZwtbGnTOsG/Hx8c6bVtjS1hS3AX8XkauAL4BNwK5vKiJdgEOBj0LOuQvYAiQA/wLuAB6oe2ERuR64HiAvLy8y0TdiSX4pny8vAODKY3pGJQZjTLPsq2zqoaqbRKQ38JmILFLVVXUvEAtlkGnjyrdBas7eM8a5LhrZnY8Wb2H6sgL+8P4S/nrhMDweu3E0Zn9tK61m8eZS/vi/H1i5rZyjemdy24k9OGz9S8TPehIEuG4qZA+IdqgmSkSEntntoRg0GCASJW0kk7ZNQGh/wFx32y6qmo9bmy0iqcD5qloScsiFwNuq6gs5Z7P7skZEXsC5udqLqv4LJ6lj5MiRjayEExl/+nApny8voHNaEsnxthabMTHigMomVd3kPq8WkenAcGCvpC0WyiDTxpVv2atrZCiPR7jj1IFMX1bAO/PzKa8J0Du7/tnnuqQncfWxvSIVqTGtSrUvwPwNJawtrOCjxVuYtsypgM9Ijue1MxI4quM25KuHYc3nMPgcGH+/M5GIOaiNG9qTwDShetNi9nOezyaJZNI2B+gnIr1wboguBi4NPUBEsoBiVQ3itKBNqHONS9ztoed0UdXNIiLAOcD3EYr/gCzauIPPlxdwWPcM3rrhaKvdNCZ2NLtsEpEOQKWq1rjHHAvsY9CRMRFUttVZo2kfBnVJ49Vrj+Tal+Ywbdk2vlq5d6tclc9pRE6I83DZkfu+njFt2bSl23jysxUsyS+lxh8EnETtlhN7M7qLn4GJxSRPDJkw+OynYXgDM7uag06HjAxeD4zlku//DaNvDvvadRFL2lTVLyI34XRt9AITVHWxiDwAzFXVycAY4GERUZwuSDfuPF9EeuLUhn9e59L/FpFsnMbo+cANkfoOB+Ld+U7F/Zj+2cQ3YVFUY0zLOMCyaRDwjIgEccYEP1Jn1kljWk75Fug+qtHDjuuXxbIHT2tw/4xu+WNAAAAgAElEQVQVhVz+/De8MWeDJW0xoAlLkuQBLwEZ7jF3quqUFg+0jdlaWs39r01jdMJSfta9gsFJ28nybyGxYiPy9QYI+vY84aifW8Jm9pDg9fCA/wIuajcb7/RH4MfPh/X6ER3T5hYiU+psuzfk9VtAvVP3q+pa9p4cAFU9MbxRHoDlHzmLFdZRVuNHZy3n/+I93JqyBmZ+XM/JbYn1/DKtS3PLJlWdiTPO1pjo8tdCZREUr4aZfz+gSx0H3NNxDVu21nL1E4Xcf/nJ5HVMDk+cZr80cUmSu4E3VPUfIjIYpyzr2eLBthE7Kmv5+vlf8WZhT17lH+T6CmEzkJzltGR3GQaDz4aMHpCRB+KBdV/B6NujHbqJMfFeD0WkU97pCNKL9s4PDlS0JyJp3d64AvzVe21uD9zjTiBDW8/Xdsqw5QyMMablKLTLdMbUrKnbIWX/XQvghScKyhn9/+I5eXAn/nXFyAO+rtlvTVmSRIE093U6kN+iEbYhX64o4J3/vsZfql7hlJ2jWC54CfqOg8TUhk/sM7ZF4jOtSw+3sqtU25FeXRr261vSdiCCfjj6Jjjhjj02n/LEF2zaXsXL14xiRF79UzG3KeLZd+FmjDEmvOIS4dfL6q04bJaaMnhsMEN7d4Pl8PGSrZRV+2ifFN/4uSac6luS5Mg6x9wHfCwiNwMpwLj6LmQz2O7bzLlzWPLO4/zJO4VdU/31OwUOsZVcTPP06OhMPzI730f35IqwX9+StgMVlwhJaXtsWl8eRxXJHNqnO9h4NmOMMZEQl+A8wqHSWRB47IjBXNmxBy/NWkdFTcCStth0CfCiqv5FRI4GXhGRIe7ESbvYDLb1KN1M2bevs23maxxTu5Rj4sB36CVw2kOwcQ50s9Zl03wJcc49f5lPIFAb9utb0hZmheU1VPkCnDQwxyYgMcYY0zpUFjvPyR05ZUhnXpq1jjfnbmBglzRy2idyWPeM6MZ38Gh0SRKc3qynAqjqLBFJArKAbS0SYWsTDMD8f8PCN9C1M2iPsi7Yky963cyoM68jKdudfKf/KdGN07QJndOS8FXGQcDX+MH7yZK2MAoElfcXOF3Lj+mbFeVojDHGmCZyW9pI7sjwnA6kJHj5y9TlAIjA9/edQkqi3TK0gEaXJAHWAycBL4rIICAJKGjRKFuTjXNg8s1Up/ViAj/mrZoj6dTrUF67+qhoR2baoMuOzMM3LQ4N+sK+wLaVwGG0JL+U+95zxgr372RjvIwxxrQSu5K2TNolePn012MoLK9h4uz1TPxm/a5uPyaymrgkya+BZ0XkVpxJSa5SVev+WJcqlObjm/E34oEzCn5OVVof7jp3EKcc0jna0Zk2ql2ClwrikECt89+ghC91s6Stmd5bkM/pQeWLZduYWroIgIKyGgD+cdkIju+XHc3wjDHGmKarLHSekzsC0Dk9ic7pSSS6ydrnywoYN7hTtKI7qDRhSZIlwLEtHVfMU4W1X8KGb2DTt86jfAvxwAO+n5DUZRCvXzOKrNTEaEdq2rB2CV5K1a3k0iCIN2zXtqStmR6bupxTFdYUVPBx8dZd2/MykxncNW0fZxpjjDExprLImQm4cCV446HTEPB4GNTZ+T179svVlrSZ2LZuJrz0I+d1Vn/Kux3HW1tzmLQ1j1FHHc97PzoEjyfcHdaM2VNKQlzEVi+2pK2Zqn0BPALXHNeLa06qd7ZdY4wxpnUI+Jxa4edOdN6f/TQMv4wLj+jOczNW882aYooraslMCdNslcaE05ov4IM7UYQPx7zL77/ysW2j0/vpzKFduPfMwZawmRYRybG/lrQ1U40/2PhBxhhjTGtwwu3Q6wQoy4f3frnHrk5pSSzfWs6G4kpL2kxs2b4O/4e/JW7Z+2yWbP7suY3/fFhOgtdDu3gvU355PL2yUqIdpTmIpFrSFlu+WV1EUUWtM1+TMcYY09olpUP/k2HjPOd9yu4ZkM8e1o0vVxTSIdkSNhM7tHQzpU+NJc5XzmP+C3k97ixyczL565k9OHd4NySME0AY01SWtMWYtUXOKudWIBhjjGlTQqb+36mk0lkkNj3ZFto2saHaF2Dpy7cxwFfGubUP8McbLuK2vA52X2aiLiUxZOKRME/qanP4HgArG4wxxrQpIVP/77R5RzUAb87dgC9gQwNMdG0ormTkPf+ld8GnfJU0mrfvv57De2RawmZiQmpiHBr2Fdoc1tJmjDHGGEc9LW1jBmTz/Iw1PPi/H/h8eQHnj8jd67SkeC8nDcoh3mt1wSZyluYX88k/b2NG4oekSRUnXfIrJCF8U6obc6BsIhJjjDHGRF5lEXjiIHH30jXH9c1i4k+P5NLnvuHLFYV8uaKw3lPHDMjmlEM6c8HhucRZ8mbCaEelj2VbSvnw+d9zb/x/KM0bB+PvQLqPinZoxuwhOcEbsZ54lrQ1w9eri6MdgjHGGBN+lUVOK1vIXYeIcEzfLObdPY7Sav9ep6wpLOeaF+cyfVkB05cV0CE5nlOHdGnJqE0b9tXKQv750kv8UiZxb/xyClIHkH31WzZGxcQkESEhQpVWlrQ1w1Hr/sEjiW8hQb+zGKkxxhjTFlTvgPKt8IfsvXZ1dB919QLWpMC23udx5KKzWLK5zJI2ExYrt5aR/8r1vOL9lPLEzlSP/hPZo660hM3EtIQ4DwTCf11L2pohp3Il1d5UEo+9FoZdFu1wjDHGmPA45mbI7LXfp8mCSWRWrwfgyU9XcNZhXeib0z7c0ZmDyDvfbeK+N2cyP+FTSvqcQ8bF/4T4dtEOy5hGxXt3Jm3hnT3SkrZmSPB6KCST9JPujXYoxhhjTPjkjnQe+2vFVOKSMxiRl8G360sY99cv6JeTynkjcvnZmD7hj9O0aarKwx/8wEnt10ENZIy62BI202okxHmgNvzXjWjfPhE5VUSWichKEbmznv09RORTEVkoItNFJDdkX0BE5ruPySHbe4nIN+41XxeRFl/tM6hKYrx1izSmNTvA8ulKEVnhPq5s2ciNiUHVO5CkDN74v6O5/Kg8Tj+0M1tLq5m5qv5JS4xpSCCo3PbmQpLK1vG7pP9CQnvoPTbaYRnTZPFxkckRIpZ5iIgXeAo4DRgMXCIig+sc9mfgZVUdCjwAPByyr0pVh7mPs0K2/wl4TFX7AtuBayP1HRriDyieCK3BYIyJvAMpn0QkE/g9cCQwCvi9iHRoqdiNiUnVpZCUTpzXw4PnHMrTlx1Op7QkUiM4/bVpe4JB5YZX5zH928V82O4eOlRvgLOegPikaIdmTJMleiOzDEUkm4tGAStVdbWq1gKTgLPrHDMY+Mx9Pa2e/XsQZ+XEE4G33E0vAeeELeImKCirwR/UMPdSNca0sAMpn04BpqpqsapuB6YCp7ZAzMbEpmAQapykLdT2Sh+bd1Sjar+YpmlemrWWqUu28kzP6SRpNXLtVBhyfrTDMma/xHsj07ATyaStG7Ah5P1Gd1uoBcB57utzgfYisnNyqiQRmSsiX4vIzsSsI1CiqjvnHK7vmgCIyPXu+XMLCgoO9Lvssr3S6aSaZN0jjWnNDqR8asq5ESuDjIk5NaWAQlLaHpu3V9Yyf0MJi/NLoxOXaTVUlb9/toJn3vuSmzJnc3jB28iwSyFnYLRDM2a/SYRmN4125nEbcIKIfAecAGxi9ySZPVR1JHAp8LiI7NdIZlX9l6qOVNWR2dl7T13cXIGgU2MYbwuHGtPW7at8alSkyiBjYk71Due5Tkvbb08fBMB7C/J59et1fLHcKi/M3jQYZPKz93PG9DP4Oulmbqt8HEnJgTF7DTU2plXw7MzZwtzLIJKdzTcB3UPe57rbdlHVfNyabBFJBc5X1RJ33yb3ebWITAeGA/8BMkQkzm1t2+uakba1tBrARrQZ07o1u3wSkU3AmDrnTo9ksMbEtBq3JW3m32HhG7s2X1DtZ0BCCczafei6zORdr9vFe8lpn7h755H/BwPPiHS0JsbMfuuvnJ3/GCtThqLH/wLpNRo6DQGPVY6b1kk8kckSIpm0zQH6iUgvnJuhi3FazXYRkSygWFWDwF3ABHd7B6BSVWvcY44FHlVVFZFpwI9xxqBcCbwbwe+wl1UFFeQB3gj9gxhjWkSzyyfgI+ChkMlHTnb3G3NwysiDfqc4LW7+ml2b0+LgqO6pBFUpqqhl4/ZKtm3fseepCRkkeD2waR6kdbOk7SCzvqiS0kVT2BTfhV63TUciNIGDMS2rlSVtquoXkZtwbnC8wARVXSwiDwBzVXUyTm31wyKiwBfAje7pg4BnRCSI04XzEVVd4u67A5gkIg8C3wHPR+o7hAoElbvfWcS360q4DXcNBmNMq3Qg5ZOqFovIH3ASP4AHVLW4xb+EMbEiKR0ue6PeXTtvMrq4j52en7GGP7y/hM6FScy660Tk//WxdbgOJqqsnPwo/nmvMN67geK8M/BawmbaiEi160R0Ll5VnQJMqbPt3pDXb7F7JsjQY2YChzZwzdU4M7+1qK2l1bw2ewNd05PIap+I1xOBVfOMMS2mueWTu28Cu1vejDH76eIjuvOH95ewpbSac56eyTu+asSStoNG+fQn6PvdQyz09OPzvndwwvk/j3ZIxoSNp41ORNJqbHHHst18Uj+Gd8+wMW3GGGNMM6UkxvHhLcczNDedBRu2g7/KWtoOIlWznmVWYDDBa6ZywuW/hXYZ0Q7JmLCJVI5gSVsTlVU7qwxktIuPciTGGGNM6zewcxpXH9uTBPyIBnl2Vj6+QDDaYZkIK1m/mOzajWzpMpZheR0aP8GYVmZ3Q1t4Z4+0pG0/5aQlRTsEY4wxpk0Y0z+HW8fkAbClUixpa8uqStCJF9N+wnHUaBxDx14Y7YiMiQiRyKRXlrQ1UUWNv/GDjDHGGNNkHVIS+NkxzhQlVSSSFGeTUbRVOu8lZPkHvOo/iacHvUyfQcOiHZIxERGhIW2RnYikzSheQ3D649wSV0rPRXOgcLn1vTfGGGPCwV8FwBjPfHTaw87Ua9kDYMh5UQ7MhE0wiP+Lv/BtcCCvZ/2C/110fLQjMiZiJEJZmyVtTTF3AmcWvej8tea62w69IIoBGWOMMW1EUgaV3jROZh58Oc/ZlpBqSVsbsnDx9wyt3cEHchFv/fyYiN3UGhMLojYRiYi0E5G7ROSf7vu+InJahOKJTRqkXJPoVTMR7tvhPM5/LtpRGWOMMa1fciY7frGcntUT6Vk9ke86X4h6rE65rQgGlU1v34MfD9dddhnJCfZva9q2nXUSqi0/EckEnKTxOPd9PvBQWKNoJX56XK9oh2CMMca0OV3S23HOsK54PcIPGwvwSUK0QzJh8sU3szktOJ1lPX9Ct/42js20feK2tUUjaeunqg8BPjeASiLX8hfTEm2AtDHGGBMRj188nBevPoJE8aHexGiHY8KgtNpH6dRHAeg39vIoR2NMC3Gb2oJRSNpqRSQJd7EBEekF1IY1ihi3848e77XJNo0xxphIqfEFSaKWzZXhvdkx0fHP6asY7l9AYd7pJPQYFe1wjGkR0Vxc+w/Ah0CuiLwETAN+G6F4YtKupC3uoGxgNMYYY1rEsLwMEvFR5o9j0D0f8sznq6j2Baj2BQgELZFrTbZX1PL6jO/p6ikiq8fgaIdjTIsLd4nV6GhQVf1AROYCx+Akj79R1W1hjiOmBd0figRraTPGGGMiJis1kWN6pFK03U9VUYCHP1jKwx8sBaB3dgqf/XpMdAM0TfbiVyt5TB5HPF4YeHq0wzGm1Ws0aRORj1X1ZODderYdFGoDQQDiPNbSZowxxkRSsvhIzunAy2eP4vv8HQBMX1bA/A0lUY7MNNWawgrKv3yG0d5FcObfoNvh0Q7JmBYX5iFtDSdtIpIAJAGdRKQ9u7topgF54Q0jti3csIOhQFq7+GiHYowxxrRt/mpISmd0/2xG988GoKi8liX5pVEOzDSFqvL3/37KnZ63qe0ykoThP4l2SMa0qF3LELbgRCQ3AouBge7zzsdHwD/DGkWMC7q9Un90WNcoR2KMMca0cf4aiNtz9kh/IEic9+Do7SIip4rIMhFZKSJ3NnDMhSKyREQWi8jElo5xXx79aBmnbfgr6XEBEs55MuQO1piDg+6c8j/Mo9oabGlT1ceAx0TkFlV9PKyf2gqJ2OyRxhhjTMT5qyAuaY9NtQGlpNJHeY2f1MS2uziziHiBp4DxwEZgjohMVtUlIcf0A+4CjlXV7SKSE51o97a5pILCGS8wLu47dMx90OmQaIdkTJvRlIlIHheRgcBgnO6SO7fHVM2OMcYYY9oAfw3E75m0JbitbF+vKmLc4E7RiKqljAJWqupqABGZBJwNLAk55jrgKVXdDhArk8NpMMD6p87h/8XNpibrEBKPuiHaIRkTVeEe09Zo05GI3A38C6dL5GnA48CPwxuGMcYYYwzOmLY6LW3njcgFwNP2O7x0AzaEvN/obgvVH+gvIl+JyNcicmp9FxKR60VkrojMLSgoiFC4u22b+jhH+mbzZY+fk/izLyC+XcQ/05hYFKkewU0p/i4CxgKbVfUnwGFASlMu3li/bBHpISKfishCEZkuIrnu9mEiMsvtq71QRC4KOedFEVkjIvPdx7AmfdMDsLW0JtIfYYxpYU0on/JEZJqIfOeWQ6e723uKSFVIGXRQjfE1JuJ81XuNafO5szjbMAXA6SXVDxgDXAI8KyIZdQ9S1X+p6khVHZmdnR3RgFQV/7xXmRscQL9z7wVv2+3Caky0NKX0q1LVAOB3Z5HcAvRo7KSQftmn4XStvERE6q6u+GfgZVUdCjwAPOxurwSuUNVDgFOBx+sUSL9R1WHuY34TvsMBSYzzhL2J0xgTPU0sn+4G3lDV4cDFwNMh+1aFlEHWB8iYcFF1W9r2bKXxBZwf4bi239S2Cege8j7X3RZqIzBZVX2qugZYjpPERc3KzdvJrNkIXYfROcNa2MzBbtf0kWG9alNKv+/chGkCMBeY7T4as6tftqrWAjv7ZYcaDHzmvp62c7+qLlfVFe7rfGAbENlqon0QlDiPZW3GtCFNKZ8UZ4kTgHQgvwXjM+bgFPABuldLW9CtOZ00Z30UgmpRc4B+ItLLXXrpYmBynWPewWllQ0SycLpLrm7JIOv6cNZ3tJNa8gbaemzGRMo+kzYREeA+VS1R1aeAM4D/U9UrmnDtpvTLXgCc574+F2gvIh3rxDAKSABWhWz+o9td6TER2bNk331e2Ppyn1HxNgnqO6BrGGNiSlPKp/uAy0VkIzAFuDlkXy+32+TnInJ8fR/Q0uNJjGkT/FXOc50xbSPyOgCwo8rH9oralo6qxaiqH7gJZ3mlH3Ba+xeLyAMicpZ72EdAkYgswanw/o2qFkUnYggGFV0wCYCcfiOjFYYxbd4+kzZVVWBqyPuVqvptGD//NuAEEfkOOAGnC0Bg504R6QK8AlytqkF38104a8cdAWQCdzQQe1j7cq/y9j7gaxhjWpVLgBdVNRc4HXhFRDzAZiDP7Tb5K2CiiKTVPbklx5MY02b43THkdWaPbJfgpW9OKtOXFXDaE19GIbCWo6pTVLW/qvZR1T+62+5V1cnua1XVX6nqYFU9VFUnRTPessoaLpUPWZ91PHQdEc1QjGnTmtI9cr6IDG/GtRvtl62q+ap6nnvz8zt3WwmAexP0P+B3qvp1yDmb3QKrBngBp5tTRPmI44dkqz0ypg1pyriRa4E3AFR1Fs6SJ1mqWrOzVltV5+H0Augf8YiNORj4q53nOi1tAP+8/HDGDsimosbfwkGZfVm9ejlZUkpxtxNtIW1jIqgpSdtwnMUdl4nIt26XoKa0tjXaL1tEstyaa3Ba0Ca42xOAt3EmKXmrzjld3GcBzgG+b0IsByw5wdsSH2OMaRlNGTeyHjgJQEQG4SRtBSKS7U5kgoj0xpkAIKrjSYxpM3a2tNWTtPXNSSW3QzJeryUGseSlmc7olUHdrUeBMZHUlDlZz2r8kL2pql9EdvbL9gITdvbLBua6zfxjgIdFRIEvgBvd0y8ERgMdReQqd9tV7kyR/xaRbJypWeYDLTJzm1jtkTFtRhPLp1/jTKV9K86kJFepqorIaOABEfEBQeAGVS2O0lcxpm3x1T+mbSd/UA+GGSRbjYoaPwvWFUEiJCbUO8WAMQevME8932jSpqqrGjtmH+dOwRnAH7rt3pDXbwFv1XPeq8CrDVzzxObGY4wxOzWhfFoCHFvPef8B/hPxAI05GO2jpQ3AHwgSby1tMWNtUQWpuIm2WDJtjCMyZZT9H9ZE9hNhjDHGRNiuMW31t9r4g0qcJW0xY+I36/mRdxaKQNfmTH9gjGkqS9qayHpHGmOMMRG2j4lIwEnaNhRX8fT0lS0YlKnP0i2lTJqznjPTViHZA6Fjn2iHZEyb1qSkTURyRWSs+zpRRFIiG1YssqzNGGOMiaidSZs3vt7dFxyeC8CMFYUtFZFpwKc/bCMnWETXyqVw2MXRDseYNq/RpE1ErsGZVe05d1MP4N1IBhVzVK2lzRhjjIm0rAHO84qP6909un82w7pn4PXYj3K0TV+2jaNz3IXOcwZFNxhjYlJ4JyJpSkvbL4CjgFIAVV0O5IQ1ihi2vaIWxRn8bIwxxpgIyu4PfcfD7H+Br7reQwJBJd5rozuiqbTax7frS7g0ZY6zIbljdAMyJoZoFCciqVbV2p1v3PWJDooqLl8gyNQftgLQMTUhytEYY4wxB4FjboaKAlj0Rr27fYEgX64oYHtFbb37TeRNW7qNbrqFEVv/AwPOgG6HRzskY9q8piRtX4nI7UCSO67tdeD9yIYVG75YXsDtby0EIKOdJW3GGGNMxPUaDZ0PhZl/h+DevVzSkuLxBZS/fWaTkUTLwo07ODZuGZ6gD8beZbO1GdMCmpK03Q6UAUuBXwKfAr+LZFCx4tv12wGI83romXUQzr1ijDHGtDQROPpmKFwGKz/Za/fzV40EwF9PQmdaxterizgstcR507FvdIMx5iDRlKTtDOA5VT1XVc9R1X+o6kFRUq7YWg44fUGtDskYY4xpIUPOg/ZdYdbf9trVPime9klxNhlJlJRW+1icX8o4zzzIHQXx7aIdkjEHhbgmHHMB8DcR+Qyna+RUVQ1ENqzY4BGhf6dUpCzakRhj6hKRX+1rv6r+taViMcaEmTcejroBpt4LmxdAl8P22B0IKnGWtEXF1MVbyWE7WRUr4Mh7oh2OMbFLwzt7ZKNJm6r+REQScVrcrgaeEZEPVPWGsEYSa1QZv30iUl0AQX+0ozHG7K19tAMwxkTQiCvh80fh/Vuh+1G7t2f2wh/MxeuxGSSj4anpK/lVxhdQDQw8M9rhGBN7IjTGsyktbahqjYi8C1QBXuBCoG0nbTs2cn7xc9QSD4lp0HlItCMyxoRQ1fujHYMxJoLaZcBxt8CMJ6BgubMtUAOBWiT4irW0RcG8ddtZXVDB8D5+KMmBnIHRDsmYg0ajSZuIjAcuAsYBM4CXgUsjHFcMcJo0n065kVt+Y/eGxsQaEXlyX/tV9RctFYsxJkJG/8Z57DTtIfj8T9QEPTamLQomfrOetCQP/XxLbSybMS2sKS1t1+OMZbtZVasiHE/MCQTD2x/VGBM286IdgDGmhQV8qLtcrLW0tbyvVxdxe6e5eLYuhpPujXY4xhxUmjKm7YKWCCTWqDrrmSfFW595Y2KRqr4U7RiMMS0s6HMmKQG8XkvaWlq1L8CQmgXQrgMcY50ZjGlJDWYkIvK5+7xdRIpDHttFpLjlQoyOxZtLAUhrFx/lSIwx+yIi2SLyZxGZIiKf7XxEOy5jTAQE/OBx6psXbypl8oJ8qn0HxYTWMcEXCJJTsx4y++xKno0xDQlvb719NSONdZ+zgOyQx873bVpFjfMjcFTvjlGOxBjTiH8DPwC9gPuBtcCcaAZkjIkQt6UtNTGO/y3azC9e+47py7ZFO6qDRiAYoHPVSuh5bLRDMSZmaYRWd24waQtZQPt5VQ2EPoDnIxJNDAm4ayvEW/cLY2JdR1V9HvCp6ueqeg1wYrSDMsZEQNCPeOKZedeJvHzNKAD8Nva8xRwVnI+HAHQeGu1QjIl5YV6mbZ8tbTvt8X+miHiBI5pycRE5VUSWichKEbmznv09RORTEVkoItNFJDdk35UissJ9XBmy/XARWeRe80mRCC2G4P6lPRG6vDEmbHzu82YROUNEhgOZjZ3UhPIpT0Smich3bhl1esi+u9zzlonIKeH7KsaYfQr4Iegjbc2H9C78jFM8c8if9Sb88N6ej5L10Y60QSLSKufJV1VO4Wuq4tJh0FnRDseYmBWpzKHBiUhE5A7gTqB9yBg2wemg2WhLm5vcPQWMBzYCc0RksqouCTnsz8DLqvqSiJwIPAz8REQygd8DI93Pm+eeux34B3Ad8A0wBTgV+GA/vnOTBNx2RkvajIl5D4pIOvBr4G9AGnDrvk5oYvl0N/CGqv5DRAbjlDc93dcXA4cAXYFPRKS/2wvBGBNJie2haju8fjm5wDMJwGacOa5D9ToBrpzc8vE1zcdAXrSD2F+l1X4O9ayhKLU/uXEJ0Q7HmIPOvmaPfBT4C04itasWej9uTEYBK1V1NYCITALOBkJvigYDv3JfTwPecV+fAkxV1WL33KnAqSIyHUhT1a/d7S8D5xCBpC3otrRZzmZMbFPV992XO9g9FrcxTSmfFCcBBEgH8t3XZwOTVLUGWCMiK93rzWr2lzDGNM34B2D4Zbve3v3O95RV+3ji4uG7j3nvFvBX///27jxOy7re//jrzbDvCLixCBYeRFGQySW09KeZWsnRLCXN9FgePWZp9XucFjO0Osc2j1lWai7lz4200nPCzAzTcgNkk3FDxRzgCKIiOwzz+f1xXQM3wz0z9zD3Pu/n43HDdV/L9/rc99zzmft7fZerBMFt18p9JAUMLGYs+fLAwuUcq3fQnh8odShmnVJrlbb3RsRLkm4juaIMQFNvxIhY0EbZw/HRo64AACAASURBVIDXM57XA4c122c+cCrwY+AUkla9wS0cOyx91GdZvxNJ55PcY46RI9t/Qasx7SPvOptZeZP0K+CLEfFO+nwQ8KN0bFtLcslP04A/SboY6AMcl3Hsk82O3SkPdTQHmVkWXbvDnuO3PX2t23rWNjbssI4efWFLyW8rey5J6/+mLNumFjmWvLj3mXpOVNC/n2+qbVYKrVXavgqcR9KFqLkA8nGp5SvATyWdAzwKLAXy0sUoIm4AbgCora1t91DApgO6+OadZuXuoKYKG0BEvJ2Oa+uoqcCtEfEjSUcAt0k6MNeDO5qDzKxtWxuDmvLsEjMLeDYiHm++QdK04ofTMavWbmLOa2/TozeoS02pwzHrlFqstEXEeen/R+1i2UuBERnPh6frMs+xjKSlDUl9gY9HxDuSlgJHNzv2kfT44c3W71BmvjTNHlmWfwrMLFMXSYPSMa+kY2JbuyAFOeQnkotWJwBExBOSepLc8iSXY82sCLY2RrleXD0NyNpHMyJGFzmWDnv4uRU0BnTrEiBX2sxKoc3ZIyWdKqlfuvxVSdMlHZxD2bOAMZJGS+pOMnB/h1HBkoZIaorha8DN6fKDwPGSBqVdnY4HHoyI5cC7kg5PZ408G7gvh1jaLbaNaSvLPwZmtt2PgCckfVvSt4HHScbktqbN/AT8AzgWQNL+QE9gZbrfGZJ6SBoNjAGezturMbOcNUbZtrT1jYj1pQ4iX25/+h/st0dfusRW6JLLxONmlm+5/OZNi4g1kt4PnERyI9vr2zooIhqAz5NUwJ4jmYVtkaQrJTXNFXs08IKkF4E9gO+mx74FfJvki9Us4MqmSUmAfwN+CSwGXqYAk5AA1C17F8jtDTKz0omIX5O02L+RPk6NiNvaOCaX/PRl4HOS5gN3AudEYhEwnWTSkj8CF3nmSLPS2NoY1JRnS1vTxGpIureUgXTUijUbmf/6O5x6yHAUW93SZlYibXUhgu1jzD4KXB8R9+XaHzsiZpBMk5257vKM5XuAe1o49ma2t7xlrp8N5DyuZJdsWsPzryb3eOndI5e3yMxKbDdgXUTcImmopNER8WprB+SQn+qAyS0c+13Si0xmVjpbo2zHnmcGtW/JosiD5e8kvTzfM7QvNG4Fj2kzK4lcaiTLJV0HnAhMSrsSVW8D1IrniJ9P5hfphfPePXuWOCAza42kpns6/hNwC9AN+H+0UOEys+rR2BjUlGWdjWhhueLULU96Hr13977gljazksml0vZJkm6RP0lnZdubjPu2VZ21b6DYyk0NJ7LXmImcNPYjpY7IzFp3CjAReAaSCY6axuGaWXUr4+6RB0t6l6TFrVe6TPo8IqJ/y4eWlxkLl7P3gJ6M2i29iO2WNrOSaLPSFhFrJS0CjpZ0NPBYRBRkHFk5+ePW93HdqRdBD7e0mZW5zRERkgJAUp9SB2RmxdEYQZcynIgkIqqiZrO5oZHZS97m9PeNQO+8lqzs3re0QZl1UrnMHvl54DfAyPQxXdK/FTqwcjCgV7dSh2BmbZsu6XpgoKTPAX8mmazIzKpcGbe0VYVZS95iw5atfGiPtXBtevvL9x5X2qDMOqlcxqadDxwaEV+PiK8DhwEXFDas0pOge031Dt0zqxYR8UOSCY3uJRnXdnlEXFvaqMysGLZG2d6nbZdJOkHSC5IWS2pxOIqkj0sKSbWFiuWvL66kW4143xt3bV+5+9hCnc7MWpHLmDYBmzOeb6ET3HO6pot8jzazChERDwEPAUjqIunMiLi9xGGZWYE1bA1WvJv1HtYVSVINcB3wIaAemCXp/nQ228z9+gFfBJ4qZDx/rnuDY/ftRffn09tYfn52IU9nZq3IpSnpNuApSZdJ+ibJjWt/VdiwSm9rY0VP9mRW9ST1l/Q1ST+VdLwSnwdeIZlAycyq3JqNW5i15G3eXLup1KHky6HA4oh4JSI2A3cBU7Ls923ge0DBaqwbt2zl1VXr+MDgd2H9m3DKDTBkTKFOZ2ZtaLPSFhHfB/4VWA+sBS5IuyNVtSF9e5Q6BDNr3W0k3SEXAp8FZgKfAP45IrJ9yTGzKvOpw0YCsGFz1dzffhjwesbz+nTdNpIOAUZExB9aK0jS+ZJmS5q9cuXKdgeyeMVaImDf3mlnq357trsMM8ufXAdtbQQ2Zfxf9Qb36V7qEMysdftGxDkRcT0wFRgHfDgi5pU4LjMrkpG79QbgN3Pqiaj+HjKSugBXA19ua9+IuCEiaiOidujQoe0+17NLVwOwX5elyYo9Dmx3GWaWP7nMHvkN4E5gL2A4cIekrxU6sFKr/tRvVvG2NC1ExFagPiKqZ3CLmbVp+KCk0nbtwy+xbHVV/PovBUZkPB+ermvSDzgQeETSEuBw4P5CTEayYOlq+vXsyqANS6CmO/QalO9TmFk75NLSdjbwvoi4LCK+QdLf+pyCRlVCWxuT/9dvbihtIGbWloMlvZs+1gAHNS1n3MjWzKrY5PcO4T9OGQ9QLS1ts4AxkkZL6g6cAdzftDEiVkfEkIgYFRGjgCeBkyMi7zOEPLt0NYfv1QXNux0O/Dh08YzaZqWUy+yRy5vt1zVdV5UibWOb/N4hJY7EzFpTLTevNbOOaZrxvxxvst1eEdGQTqj0IFAD3BwRiyRdCcyOiPtbLyE/Njc08vzyNXx63OuwfDMcdHoxTmtmrcil0vYWsEjSgyS9Bo8nmYL2aoCI+FIB4yu6xnTWSFX/XQ3MzMwqXtNkz9Vyk+2ImAHMaLbu8hb2PboQMbz4xho2b21kgl6ELl1hxGGFOI2ZtUMulbY/pI8mTxYolrLQ1LvCvQDMzMzKX2P6h7sKGtrKxsJ0EpIRaxfAXgdD994ljsjM2qy0RcRNxQikXDSm3SOr5IKdmZlZVWuqtFVD98hysXDpagb1hB4r5kPteaUOx8zIfcr/TqNpIhJ3jzQzMyt/TcMaalxpy5tnl67mY0NXoIaNMNJdI83KgSttzXR97VHAV+zMzMwqQdOYNv/dzo+tjcGLb6zhyB4vJytGHF7agMwqVOT5BmI5V9ok9cjrmctUlxV1AKzrM6KNPc3MzKzUto1p82XovHj9rfVs3NLIuC11MGgU9Nuj1CGZGbndXPtQSQuBl9LnB0v6SS6FSzpB0guSFkv6apbtIyXNlDRX0gJJJ6Xrz5Q0L+PRKGlCuu2RtMymbbu36xW3IdSFRY37sKm3k5RZtcohN/1XRo55UdI7Gdu2ZmwryvTbZtaypkqbu0fmxwtvrAGCPVbPcyubWRnJZfbIa4GPAr8HiIj5ko5p6yBJNcB1wIeAepLbBNwfEXUZu10GTI+In0saRzLF7aiIuB24PS1nPPD7iJiXcdyZhbiRJGxvyvREJGbVKZfcFBGXZux/MTAxo4gNETGhWPGaWevcPTK//rFqPUN5h24bV8GwSaUOx8xSuXQm6BIRrzVbtzWH4w4FFkfEKxGxGbgLmNJsnwD6p8sDgGVZypmaHlsUG7ckL81X7MyqVi65KdNU4M6iRGZm7dbU0jZryVsljqQ6LH1nA2N6vJ08GbRPaYMxs21yqbS9LulQICTVSLoEeDGH44YBr2c8r0/XZZoGnCWpnqSV7eIs5ZzOzl+Ybkm7Jn1Tyl67knS+pNmSZq9cuTKHcBMbNifTR+7Wp3vOx5hZRcklNwEgaR9gNPCXjNU909zypKR/bukku5qDzKx9Dhs9GIC/vujfs3xY9s4GJvVekTzZ7T2lDcbMtsml0nYh8CVgJPAGcHi6Lh+mArdGxHDgJOA2aftQYkmHAesj4tmMY86MiPHAUenj09kKjogbIqI2ImqHDh3ajpCSK3Zj9ujXvldiZtXoDOCeiMjsXbBPRNQCnwKukZT1W82u5yAza49J+wyiV7caajyuIS/q395Abc1L0K0P7LZvqcMxs1SblbaIWBERZ0TEkPRxRkS8mUPZS4HMKRiHp+synQdMT8/zBNATGJKx/QyatbJFxNL0/zXAHSRdnfIm7WVBVyd/s2qVS25q0loOegV4hB3Hu5lZCTRG4FEN+bH0nQ0ctHkejDwcunhKTrNy0eZEJJJuhJ1vNBAR57dx6CxgjKTRJF+IziC5Mp3pH8CxwK2S9ieptK1Mz9sF+CRJa1pTLF2BgRHxpqRuJBOk/Lmt19AeTROR+IqdWdXKJTchaSwwCHgiY90gktb/TZKGAJOB7xclajNrUYQnIsmHTQ1bWbdhAwPiDRiWtSOTmZVILrNHZlaKegKnsON4kKwiokHS54EHgRrg5ohYJOlKYHZE3A98GbhR0qUkFcNzIprauvgA8Hp6NbtJD+DBtMJWk8Z2Yw6vIWeNbmkzq2o55iZIKnN3ZeQkgP2B6yU1kvRUuKrZjLhmVgKNEZ71OQ/eWreZvbSKLjR6EhKzMtNmpS0i7s58Luk24G+5FB4RM0gmGMlcd3nGch3Jlepsxz5CMn4uc906oLDzz6Zfz9zSZla92spN6fNpWY57HBhf0ODMrN2SSpv/bnfUqrWb2Vf/mzwZNKqksZjZjnals/JooGrvPL12UwPgSpuZmVmlaAxoYTJpa4dV6zYzWsuTJ0PHljYYM9tBm5U2SW9Leit9vAM8BHyt8KGVxqaGZMr/3t1z6TlqZmZmpdTUg/me2W2O3LA2vLVuEyO1gsZufaD34FKHY2YZWq2ZpPdAO5jtM6s1NhvfUVUaG4OVazYxqlsXunf1jElmZmaVYtnqjWxq2EqPUgdSwVav38JeWkVj/+HubmpWZlqtmaQVtBkRsTV9VG2FDWDl2k2Au1iYmZlVCkl85fj9AM8g2VHvbmxgqFbTpd/upQ7FzJrJpTlpnqROcR+ixrROOrhP9xJHYmZmZrlqmvnZlbaOeXfDFkZqJV0Gjix1KGbWTIvdIyV1jYgGkhvHzpL0MrAOEEkj3CFFirFoqrsd0czMrDo1XXR1la1j1m7YwO56GwaMKHUoZtZMa2PangYOAU4uUiwl5zqbmZlZ5Wm66OqGto5pXLsyWeg7tLSBmNlOWqu0CSAiXi5SLCXX2Ohqm5mZWaWJCCSPSe+omvVvJgt9PKbNrNy0VmkbKulLLW2MiKsLEI+ZmZlZuzSGx7PlQ836FclCvz1LG4iZ7aS1SlsN0JdO1EXcY9rMzMwqT2NE5/myUkBr3l4J3YFeu5U6FDNrprVK2/KIuLJokZSBwAOZzczMKk3glrZ8+OL7h8BsoNegUodiZs20NuV/p8t+24a0dbpXbmZmVrka0zFt1jH79t6QLPQcUNpAzGwnrVXaji1aFGWiyu8dbmZmVpXCY9ry49EfJP/XtNYRy8xKocVKW0S8VcxAyoGrbGZmZpWnsTHo4jpbx+0xvtQRmFkLfCklQ2y7Oaczv5mZWaVoDE/3nxef/TNs3VzqKMwsC1faMrh3pJmZWeUJPKYtL7r1TB5m1mH5rle0Nqat03GdzczMrPJ4TJuZlYtCpSJX2jK4pc3MzKzyNEawesMWNm7ZWupQzMwKoqCVNkknSHpB0mJJX82yfaSkmZLmSlog6aR0/ShJGyTNSx+/yDhmkqSFaZnXKo+d2BvD92kz6yxyyE//lZGDXpT0Tsa2z0h6KX18priRm1lzvbrXAHD7U/8ocSRmZoVRsEqbpBrgOuBEYBwwVdK4ZrtdBkyPiInAGcDPMra9HBET0scFGet/DnwOGJM+TshXzOH7tJl1Crnkp4i4tCkHAT8BfpseuxvwLeAw4FDgW5J8J1qzErr4/4wBcEubmVWtQra0HQosjohXImIzcBcwpdk+AfRPlwcAy1orUNJeQP+IeDKSqR5/DfxzvgIOj2oz6yxyyU+ZpgJ3pssfBh6KiLci4m3gIfJ48cjM2q9rOt9/NYxry6EXwJck1aU9lB6WtE8p4jSz4ipkpW0Y8HrG8/p0XaZpwFmS6oEZwMUZ20an3Sb/KumojDLr2ygTAEnnS5otafbKlStzCthj2sw6jVzyEwDpF6LRwF/ac+yu5CAz2zXbhjdUeJ0tx15Kc4HaiDgIuAf4fnGjNLNSKPVEJFOBWyNiOHAScJukLsByYGTabfJLwB2S+rdSzk4i4oaIqI2I2qFDh+Z4TPuCN7NO4QzgnohoV7+rXclBZrZrmv5+V8ENttvsBRARMyNiffr0SWB4kWM0sxIoZKVtKTAi4/nwdF2m84DpABHxBNATGBIRmyJiVbp+DvAysF96fGZyylbmLmvqHln5Od/M2pBLfmpyBtu7Rrb3WDMrgu0TiVX8X/CcewGkzgMeKGhEZlYWCllpmwWMkTRaUneSLz73N9vnH8CxAJL2J6m0rZQ0NO0igKR9SSYceSUilgPvSjo8nTXybOC+fAX88sq1gO/XZtYJ5JKfkDQWGAQ8kbH6QeB4SYPSCUiOT9eZWYlsm0es4utsuZN0FlAL/KCF7e6ibVZFClZpi4gG4PMkX2aeI5klcpGkKyWdnO72ZeBzkuaTXMk+J51g5APAAknzSPprXxARb6XH/BvwS2AxSQtc3q4wNTYm/3evKXWvUTMrpBzzEySVubvSvNR07FvAt0kqfrOAKzPyk5mVQKR/v6tgIpKcWvIlHQd8Azg5IjZlK8hdtM2qS9dCFh4RM0gmGMlcd3nGch0wOctx9wL3tlDmbODA/Ea6oypI+mbWhrbyU/p8WgvH3gzcXLDgzKxdqmUiEjJ6AZBU1s4APpW5g6SJwPXACRGxovghmlkpuEnJzMzMKlpTU3ilX3TNsRfAD4C+wG8kzZO0U9duM6s+BW1pMzMzMyu0Kmppy6WX0nFFD8rMSs4tbWZmZlbRtlfaqqDWZmaWhSttZmZmVtmq5z5tZmZZudJmZmZmFa0xrbRVwX3azMyycqXNzMzMKlqkTW1uaTOzauVKm5mZmVW0xm3dI11rM7Pq5EqbmZmZVbTG7f0jzcyqkittZmZmVhXc0mZm1cqVtgwbtmwtdQhmZmbWTtum/C9xHGZmheJKW4bX314PVMfNOc3MzDqLaBrT5m81ZlalupY6gHLSt3vydnT19FOd2pYtW6ivr2fjxo2lDqVq9OzZk+HDh9OtW7dSh1I2/DnLP3/OOq+mljZ3j8ydc1BhOA9ZobjSliEy/rXOq76+nn79+jFq1CjkLwAdFhGsWrWK+vp6Ro8eXepwyoY/Z/nlz1nn1jQPyT1z6pni78o5cQ7KP+chKyR3JMgQAd1pgJrupQ7FSmjjxo0MHjzYf8TyRBKDBw/21dxm/DnLL3/OOrc9+vcA4KlX3ipxJJXDOSj/nIcsU76bgVxpyxBEUmnr6kpbZ+c/Yvnl9zM7vy/55fez8+rXsxvnTh5Fj27+WtMe/p3JP7+nVqhPgLNbhgjori3QtUepQzEzM7N2iPDskWZWvVxpyxBANxqQu0daCa1atYoJEyYwYcIE9txzT4YNG7bt+ebNm3Mq49xzz+WFF15odZ/rrruO22+/PR8hWwXy58yqTUSwcUtjqcOwHDkHmbWPJyLJFOExbVZygwcPZt68eQBMmzaNvn378pWvfGWHfSKCiKBLC/Nb33LLLW2e56KLLup4sFax/DmzatO1pgubtzayfPVG9upd6misLc5BZu1T0EqbpBOAHwM1wC8j4qpm20cCvwIGpvt8NSJmSPoQcBXQHdgM/N+I+Et6zCPAXsCGtJjjI2JFPuINoAfuHmnbXfHfi6hb9m5eyxy3d3++9bED2n3c4sWLOfnkk5k4cSJz587loYce4oorruCZZ55hw4YNnH766Vx++eUAHHnkkfz0pz/lwAMPZMiQIVxwwQU88MAD9O7dm/vuu4/dd9+dyy67jCFDhnDJJZdw5JFHcuSRR/KXv/yF1atXc8stt/D+97+fdevWcfbZZ/Pcc88xbtw4lixZwi9/+UsmTJiQ1/ek2NrKTek+nwSmkaSG+RHxqXT9VmBhuts/IuLkjsbjz1l1fs6suM47cjQ3/e3VtLXNHSXbwznIOcjKX8G6R0qqAa4DTgTGAVMljWu222XA9IiYCJwB/Cxd/ybwsYgYD3wGuK3ZcWdGxIT0kZcKGyT94bvJLW1Wvp5//nkuvfRS6urqGDZsGFdddRWzZ89m/vz5PPTQQ9TV1e10zOrVq/ngBz/I/PnzOeKII7j55puzlh0RPP300/zgBz/gyiuvBOAnP/kJe+65J3V1dXzzm99k7ty5BX19xZBLbpI0BvgaMDkiDgAuydi8ISP/dLjCVo78ObNKtPfAXnSrEZ4HovI5B5ntrJAtbYcCiyPiFQBJdwFTgMzftAD6p8sDgGUAEZH527II6CWpR0RsKmC822ePdKXNUrtylbCQ3vOe91BbW7vt+Z133slNN91EQ0MDy5Yto66ujnHjdrw20qtXL0488UQAJk2axGOPPZa17FNPPXXbPkuWLAHgb3/7G//+7/8OwMEHH8wBB5TX+7GLcslNnwOui4i3AfJ5cSgbf86q8nNmJRC+1eoucQ5yDrLyV8hK2zDg9Yzn9cBhzfaZBvxJ0sVAH+C4LOV8HHimWYXtlrSL0r3AdyLyk6aT+7S5e6SVrz59+mxbfumll/jxj3/M008/zcCBAznrrLOy3hume/ftFyFqampoaGjIWnaPHj3a3KdK5JKb9gOQ9HeSLpTTIuKP6baekmYDDcBVEfH7AsdbdP6cWaVyna06OAeZ7azUs0dOBW6NiOHAScBtkrbFJOkA4HvAv2Ycc2babfKo9PHpbAVLOl/SbEmzV65cmVMwTbNHuqXNKsG7775Lv3796N+/P8uXL+fBBx/M+zkmT57M9OnTAVi4cGHWLilVqiswBjiaJE/dKGlgum2fiKgFPgVcI+k92QrYlRxUjvw5s0qSp2u4Vkacg8wShWxpWwqMyHg+PF2X6TzgBICIeEJST2AIsELScOB3wNkR8XLTARGxNP1/jaQ7SLo6/br5ySPiBuAGgNra2pyyeNLS5kqbVYZDDjmEcePGMXbsWPbZZx8mT56c93NcfPHFnH322YwbN27bY8CAAXk/T5HlkpvqgaciYgvwqqQXSSpxszJy0CvpxEgTgZebHb9LOagc+XNmlSTwFCTVxjnILNU0nWq+HyQVwleA0SSzQM4HDmi2zwPAOeny/iRj2kQym+R84NQsZQ5Jl7sB9wAXtBXLpEmTIhf/8YdnI77VP2Lmf+a0v1Wnurq6UodQNrZs2RIbNmyIiIgXX3wxRo0aFVu2bNmlsrK9r8DsKFAOaumRY246AfhVujyEpDvlYGAQ0CNj/UvAuLbOmS0H+XO2XaE/Z9Z57PPv/xNLrj4u4pcfyvmYUuShYj+cg1qXzxwU4fe2s3vitmkR3+ofa1avymn/XHNQwVraIqJB0ueBB0nGhNwcEYskXZkGdz/wZZJuR5eSXCA7JyIiPe69wOWSLk+LPB5YBzwoqVta5p+BG/MVc03jlnShW76KNKtoa9eu5dhjj6WhoYGI4Prrr6dr18q+vWOOuelB4HhJdcBWktuOrJL0fuB6SY0k3cuvigj3o+mgavycmVnlcA6ySlDQT2REzABmNFt3ecZyHbBTO3dEfAf4TgvFTspnjJm2V9o8EYkZwMCBA5kzZ06pw8i7HHJTAF9KH5n7PA6ML0aMnUm1fs6suJJfW7P2cw6ySlDqiUjKSk1sThY8e6SZmVlF2VZn86A2M6tCrrRl2Gf17GTB3SPNzMwqiutsZlbNXGnLsOe655KFEYeXNhAzMzNrF3ePNLNq5kpbM+ujB+w+ttRhmJmZWTu4ymZm1cyVNrMyc8wxx+x089BrrrmGCy+8sMVj+vbtC8CyZcs47bTTsu5z9NFHM3v27FbPfc0117B+/fptz0866STeeeedXEO3CuLPmVUbN7RVFucgs/ZxpS2D872Vg6lTp3LXXXftsO6uu+5i6tSpbR679957c8899+zyuZv/IZsxYwYDBw7c5fKsfPlzZtUm0r/i8qi2iuAcZNY+vglFpsAjmG1HD3wV/ndhfsvcczyceFWLm0877TQuu+wyNm/eTPfu3VmyZAnLli1j4sSJHHvssbz99tts2bKF73znO0yZMmWHY5csWcJHP/pRnn32WTZs2MC5557L/PnzGTt2LBs2bNi234UXXsisWbPYsGEDp512GldccQXXXnsty5Yt45hjjmHIkCHMnDmTUaNGMXv2bIYMGcLVV1/NzTffDMBnP/tZLrnkEpYsWcKJJ57IkUceyeOPP86wYcO477776NWrV37fs2rnz5k/Z9ZhbmnrAOcg5yAre25pMyszu+22G4ceeigPPPAAkFx5/OQnP0mvXr343e9+xzPPPMPMmTP58pe/3OrA+5///Of07t2b5557jiuuuGKHe9B897vfZfbs2SxYsIC//vWvLFiwgC984QvsvffezJw5k5kzZ+5Q1pw5c7jlllt46qmnePLJJ7nxxhuZO3cuAC+99BIXXXQRixYtYuDAgdx7770FeFcs3/w5s6rli68VwTnIqtWqvvvx64YPQZf8zkbvlrYMvkhnO2nlKmEhNXUbmTJlCnfddRc33XQTEcHXv/51Hn30Ubp06cLSpUt544032HPPPbOW8eijj/KFL3wBgIMOOoiDDjpo27bp06dzww030NDQwPLly6mrq9the3N/+9vfOOWUU+jTpw8Ap556Ko899hgnn3wyo0ePZsKECQBMmjSJJUuW5Old6ET8OQP8ObOOcUtbBzgHAc5Blh/LdjuM7zb059Ru+W2JdUubWRmaMmUKDz/8MM888wzr169n0qRJ3H777axcuZI5c+Ywb9489thjDzZu3Njusl999VV++MMf8vDDD7NgwQI+8pGP7FI5TXr02H4z+pqaGhoaGna5LCsuf86smmwf02aVwjnILHeutGXwVTorF3379uWYY47hX/7lX7YNyl69ejW777473bp1Y+bMmbz22mutlvGBD3yAO+64A4Bnn32WBQsWAPDuu+/Sp08fBgwYwBtvvLGtawpAv379WLNmzU5lHXXUUfz+979n/fr1rFu3jt/97nccdsZfBgAAEWpJREFUddRR+Xq5ViL+nFk1qZa/4ZJOkPSCpMWSvpplew9Jd6fbn5I0qvhR5odzkFnu3D0yw7sbt5Q6BLNtpk6dyimnnLJtdq0zzzyTj33sY4wfP57a2lrGjm39foIXXngh5557Lvvvvz/7778/kyZNAuDggw9m4sSJjB07lhEjRjB58uRtx5x//vmccMIJ2/r7NznkkEM455xzOPTQQ4FkcPbEiRPdPaQK+HNm1aIa6mySaoDrgA8B9cAsSfdHRF3GbucBb0fEeyWdAXwPOL340eaHc5BZbtTa4M5qUVtbG23dswPg7u9+ho9ufoA+V6woQlRWrp577jn233//UodRdbK9r5LmRERtiUIqmmw5yJ+zwvD72nmt2biF8dP+xN/2+jHD+wLn/Smn48opD0k6ApgWER9On38NICL+M2OfB9N9npDUFfhfYGi08oXOOai4/N52bn+ue4Pfzq3nR5+YQK/uNW3un2sOcktbhq5dRE0X94Y3MzOrNE01FlX2n/FhwOsZz+uBw1raJyIaJK0GBgNvZu4k6XzgfICRI0cWKl4za+a4cXtw3Lg98l6ux7SZmZlZxesEHYfaJSJuiIjaiKgdOnRoqcMxsw5ypc0si87QbbiY/H5m5/clv/x+dnLV8eNfCozIeD48XZd1n7R75ABg1a6czL8z+ef31ArFlTazZnr27MmqVaucePMkIli1ahU9e/YsdShlxZ+z/PLnzKI6am2zgDGSRkvqDpwB3N9sn/uBz6TLpwF/aW08W0ucg/LPecgKyWPammxtYGzD83QN33ejsxs+fDj19fWsXLmy1KFUjZ49ezJ8+PBSh1FW/DnLP3/OOrdqqHukY9Q+DzwI1AA3R8QiSVcCsyPifuAm4DZJi4G3SCp27eYcVBjOQ1YorrQ1WfRbDtj6HBvU229KJ9etWzdGjx5d6jCsyvlzZpZf2yYiKWkUHRcRM4AZzdZdnrG8EfhER8/jHGRWWQraPTKHG0SOlDRT0lxJCySdlLHta+lxL0j6cK5l7rJNyU0Wf7zXVXkr0szKVy65RNInJdVJWiTpjoz1n5H0Uvr4TLZjzay4tnfzq/Rqm5nZzgrWqJTjDSIvA6ZHxM8ljSO5sjQqXT4DOADYG/izpP3SY9oqs0NWddsrX0WZWZnKJT9JGgN8DZgcEW9L2j1dvxvwLaCW5OL+nPTYt4v9Osxsu229I11nM7MqVMiWtkOBxRHxSkRsBu4CpjTbJ4D+6fIAYFm6PAW4KyI2RcSrwOK0vFzKNDNrSy655HPAdU2VsYhYka7/MPBQRLyVbnsIOKFIcZtZG1xnM7NqVMjhW7ncIHIa8CdJFwN9gOMyjn2y2bHD0uW2ygR2vKkksFbSCznEPAQ+/OYP/yWHPUtrCM1uolnGKiVWx5lfuca5T6EDaUEu+Wk/AEl/J5kQYFpE/LGFY4c1O3ZXc1BHVMpnozV+DeWhol/D3sl/Q/iscn0NpcpDRTNnzpw3Jb2W4+6V8vN3nPnlOPOrPXHmlINKPefGVODWiPiRpCNIZkM6MB8FR8QNwA3tOUbS7Iiozcf5C6lS4oTKidVx5lelxNmGrsAY4GiSeyU9Kml8rgfvSg7qiGp4z/0ayoNfQ/WJiJzvrl0p753jzC/HmV+FiLOQ3SNzuUHkecB0gIh4AuhJUjNt6dhcyjQza0suuaQeuD8itqTdtF8kqcQ5D5mZmVlRFbLSlssNIv8BHAsgaX+SStvKdL8zJPWQNJrki9LTOZZpZtaWXHLJ70la2ZA0hKS75Csk9086XtIgSYOA49N1ZmZmZgVRsO6ROd4g8svAjZIuJZmU5JxI5uxdJGk6UAc0ABdFxFaAbGXmMeyidWXqoEqJEyonVseZX2UdZ475qalyVgdsBf5vRKwCkPRtkoofwJUR8VbxX8VOyvo9z5FfQ3nwa+jcKuW9c5z55TjzK+9xavt9TczMzMzMzKzcFPTm2mZmZmZmZtYxrrSZmZmZmZmVsU5ZaZN0gqQXJC2W9NUs23tIujvd/pSkUcWPMqc4vySpTtICSQ9LKsm9ZtqKM2O/j0sKSSWZqjWXOCV9Mn1PF0m6o9gxpjG09XMfKWmmpLnpz/6kEsV5s6QVkp5tYbskXZu+jgWSDil2jNWoUvJCayolZ7SmUvJJayol17TEOahj/F2ouHFm7OfvQjmohPxU9BwUEZ3qQTLpwMvAvkB3YD4wrtk+/wb8Il0+A7i7TOM8BuidLl9YrnGm+/UDHiW5aXptOcZJMkvpXGBQ+nz3Mo3zBuDCdHkcsKTYcabn/gBwCPBsC9tPAh4ABBwOPFWKOKvpUSl5oaOvId2vpDkjDz+HkueTPLyGssg1rbwG56DC/vz9XSiPcab7+btQ/uIseX4qdg7qjC1thwKLI+KViNgM3AVMabbPFOBX6fI9wLGSVMQYIYc4I2JmRKxPnz5Jcr+oYsvl/QT4NvA9YGMxg8uQS5yfA66LiLcBImJFkWOE3OIMoH+6PABYVsT4tgcR8SjQ2qyJU4BfR+JJYKCkvYoTXdWqlLzQmkrJGa2plHzSmorJNS1xDuoQfxfKr0rJa5WSuyoiPxU7B3XGStsw4PWM5/Xpuqz7REQDsBoYXJTossSQyhZnpvNIavPF1macaXPwiIj4QzEDayaX93M/YD9Jf5f0pKQTihbddrnEOQ04S1I9MAO4uDihtVt7P8PWtkrJC62plJzRmkrJJ62pplzTEueglvm7UH5VSl6rlNxVLfkprzmoYPdps+KRdBZQC3yw1LE0J6kLcDVwTolDyUVXkm4BR5NcqXtU0viIeKekUe1sKnBrRPxI0hHAbZIOjIjGUgdm5aOc80JrKixntKZS8klrnGusYpRzzquwvFYpuavT5afO2NK2FBiR8Xx4ui7rPpK6kjS7ripKdFliSGWLE0nHAd8ATo6ITUWKLVNbcfYDDgQekbSEpE/v/SUYgJvL+1kP3B8RWyLiVeBFksRVTLnEeR4wHSAingB6AkOKEl375PQZtnaplLzQmkrJGa2plHzSmmrKNS1xDmqZvwvlV6XktUrJXdWSn/KbgwoxMK+cHyRXEF4BRrN9cOMBzfa5iB0H304v0zgnkgzUHFPO72ez/R+hNINvc3k/TwB+lS4PIWnSHlyGcT4AnJMu70/Sj1sl+vmPouUBuB9hxwG4T5cixmp6VEpe6OhraLZ/SXJGHn4OJc8neXgNZZNrWnkdzkGF+/n7u1Ae42y2f0nyWqXkrkrKT8XMQUV9YeXyIJnN5cX0l/wb6borSa7QQFJb/w2wGHga2LdM4/wz8AYwL33cX45xNtu3JIkqx/dTJN0X6oCFwBllGuc44O9pEpsHHF+iOO8ElgNbSK7MnQdcAFyQ8X5el76OhaX6uVfbo1LyQkdeQ7N9S5YzOvhzKIt80sHXUBa5ppX4nYMK+/P3d6E8xtls35LltUrJXZWQn4qdg5QWamZmZmZmZmWoM45pMzMzMzMzqxiutJmZmZmZmZUxV9rMzMzMzMzKmCttZmZmZmZmZcyVNjMzMzMzszLmSlsnJ2mrpHkZj1Gt7DtK0rPFi65lkmolXZsuHy3p/RnbLpB0dhFjmSDppGKdz8xyk5HfnpX0G0m981DmttzTwva9Jd3T0fOYWWWSNDjjO9X/SlqaLr8jqa4A5zta0v+085hHst3YW9I5kn6av+gsn7qWOgAruQ0RMaHUQbRXRMwGZqdPjwbWAo+n236R7/NJ6hoRDS1sngDUAjPyfV4z65Bt+U3S7ST3z7m6aaMkkdyMtTHXApvlnmzblwGn7XLEZlbRImIVyfcCJE0D1kbED9OL4m1Wrtr4vmGdmFvabCdpi9pjkp5JH+/Pss8Bkp5Orx4tkDQmXX9WxvrrJdVkOXaJpO9LWpju+96M8/4lLe9hSSPT9Z9Ir5TPl/Rouu5oSf+TJsELgEvTcx4laZqkr0gaK+npZq9rYbo8SdJfJc2R9KCkvbLEeaukX0h6Cvi+pEMlPSFprqTHJf2TpO4kN3s8PT3/6ZL6SLo5fW1zJU3p8A/FzDrqMeC9aR54QdKvgWeBEZKOT3+3n0lb5PoCSHpf+rs+P/197pd5VVvSBzOuqM9Nt2/rkSCpp6Rb0lw3V9Ix6fpzJP1W0h8lvSTp+yV6T8ysuGok3ShpkaQ/SeoF21q+rpE0G/iipKGS7pU0K31MTvfbKeek5faVdI+k5yXdnl6QQtKx6X4L0+8lPZoHJOlcSS+m35cmF+l9sF3gSpv1ykgAv0vXrQA+FBGHAKcD2boCXQD8OL2KXQvUS9o/3X9yun4rcGYL510dEeOBnwLXpOt+AvwqIg4Cbs847+XAhyPiYODkzEIiYgnwC+C/ImJCRDyWse15oLuk0emq04G7JXVLz3VaREwCbga+20Kcw4H3R8SXgOeBoyJiYhrTf0TE5nT57vT8dwPfAP4SEYcCxwA/kNSnhfLNrMAkdQVOBBamq8YAP4uIA4B1wGXAcWnOmw18Kb0gczfwxTT3HAdsaFb0V4CL0nx3VJbtFwGR5rqpwK8k9Uy3TSDJSeNJLvqMyNsLNrNyNQa4Ls097wAfz9jWPSJqI+JHwI9Jvte8L93nl+k+LeWcicAlwDhgX2BymmtuBU5Pc1BX4MLMYNIL1leQVNaOTI+3MuXukZate2Q34KeSmipe+2U57gngG5KGA7+NiJckHQtMAmalF3l6kVQAs7kz4///SpePAE5Nl28Dmq4+/x24VdJ04LfteXHAdJIvRlel/58O/BNwIPBQGmcNsLyF438TEVvT5QEkX7rGAEHyPmVzPHCypK+kz3sCI4Hn2hm7mXVML0nz0uXHgJuAvYHXIuLJdP3hJF9U/p7mg+4k+e2fgOURMQsgIt4FSPdp8nfgaiVdL38bEfXNth9JcoGIiHhe0mtsz6cPR8TqtMw6YB/g9Ty9bjMrT69GRFNOmgOMyth2d8byccC4jHzSP+0B0FLOeToi6gHSnDcKWJOe78W0jF+RXEi6hu0OAx6JiJXpsXeT/TuflQFX2iybS4E3gINJWmM3Nt8hIu5Iuw1+BJgh6V8BkbSUfS2Hc0QLyzvvGHGBpMPSc82RNCm3lwEkSfA3kn6bFBUvSRoPLIqII3I4fl3G8reBmRFxipJumY+0cIyAj0fEC+2I08zyb6eLUukXnMzfawEPRcTUZvuNb6vwiLhK0h+Ak0gqfR8mS75swaaM5a3477FZZ9D8975XxvPMvNQFODwimueTbDknW7nOJ1XI3SMtmwEkV5gbgU+TtETtQNK+wCsRcS1wH3AQ8DBwmqTd0312k7RPC+c4PeP/J9Llx4Ez0uUzSa6MI+k9EfFURFwOrASadyNaA/Qji4h4mSSBfZPtV7FeAIZKOiItv5ukA1qIM9MAYGm6fE4r538QuDijT/nEHMo2s9J4kqQrUdPY2j6S9iPJE3tJel+6vl/azXKbNDctjIjvAbOAsc3Kfoy0i3ha5si0XDOz1vwJuLjpSdrzKZeck+kFYFRTbiP5PvfXZvs8BXxQyYyX3YBP5OsFWP650mbZ/Az4jKT5JAlhXZZ9Pgk8mzbDHwj8OiLqSMaG/EnSAuAhYKcJPlKD0n2+SNKyB0mCOjdd/+l0GyRjwhamg/sfB+Y3K+u/gVPScXlHZTnX3cBZJF0lScehnQZ8L32N84CdJlvJ4vvAf0qay45XsWaSdGOYJ+l0kha5bsACSYvS52ZWhtJuQecAd6a55wlgbJonTgd+kuaJh0i6Ome6RMkkSQuALcADzbb/DOiiZAKku4FzImITZmat+wJQq2RitjqSeQSg7ZyzTdpKdy5Jb6OFQCPJHACZ+ywHppHkvb/jYRxlTRGt9kwzyztJS4DaiHiz1LGYmZmZmZU7t7SZmZmZmZmVMbe0mZmZmZmZlTG3tJmZmZmZmZUxV9rMzMzMzMzKmCttZmZmZmZmZcyVNjMzMzMzszLmSpuZmZmZmVkZ+/+gB/DLq2mxUwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Training statistics:\n", " precision recall f1-score support\n", "\n", " ham 0.99 0.96 0.98 3857\n", " spam 0.80 0.94 0.86 601\n", "\n", "avg / total 0.96 0.96 0.96 4458\n", "\n", "Training classification accuracy: 0.960\n", "\n", "Validation statistics:\n", " precision recall f1-score support\n", "\n", " ham 0.99 0.97 0.98 491\n", " spam 0.79 0.94 0.86 66\n", "\n", "avg / total 0.97 0.96 0.97 557\n", "\n", "Validation classification accuracy: 0.964\n" ] } ], "source": [ "from sklearn.metrics import classification_report, roc_curve, precision_recall_curve, accuracy_score\n", "\n", "# Get training/validation predictions\n", "y_train_pred = model[1].predict_proba(X_train_sparse)\n", "y_val_pred = model[1].predict_proba(X_val_sparse)\n", "# Plot ROC curves\n", "def roc_plot(y_train, y_train_pred, y_val, y_val_pred):\n", " fpr, tpr, _ = roc_curve(y_train, y_train_pred[:,1])\n", " plt.plot(fpr, tpr, label=\"Training\")\n", " fpr, tpr, _ = roc_curve(y_val, y_val_pred[:,1])\n", " plt.plot(fpr, tpr, label=\"Validation\")\n", " plt.ylim(ymin=0.8, ymax=1)\n", " plt.xlabel(\"False positive rate\")\n", " plt.ylabel(\"True positive rate\")\n", " plt.title(\"ROC\")\n", " plt.legend(loc=4)\n", "# Plot precision/recall curves\n", "def pr_plot(y_train, y_train_pred, y_val, y_val_pred):\n", " precision, recall, _ = precision_recall_curve(y_train, y_train_pred[:,1])\n", " plt.plot(precision, recall, label=\"Training\")\n", " precision, recall, _ = precision_recall_curve(y_val, y_val_pred[:,1])\n", " plt.plot(precision, recall, label=\"Validation\")\n", " plt.ylim(ymin=0.6, ymax=1)\n", " plt.xlabel(\"Precision\")\n", " plt.ylabel(\"Recall\")\n", " plt.title(\"Precision/recall\")\n", " plt.legend(loc=3)\n", "def f1_plot(y_train, y_train_pred, y_val, y_val_pred):\n", " p, r, thresh = precision_recall_curve(y_train, y_train_pred[:,1]); p=p[:-1]; r=r[:-1]\n", " plt.plot(thresh, 2*p*r/(p+r), label=\"Training\")\n", " p, r, thresh = precision_recall_curve(y_val, y_val_pred[:,1]); p=p[:-1]; r=r[:-1]\n", " plt.plot(thresh, 2*p*r/(p+r), label=\"Validation\")\n", " #plt.ylim(ymin=0.6, ymax=1)\n", " plt.xlabel(\"Threshold\")\n", " plt.ylabel(\"F1\")\n", " plt.title(\"$F_1$ score\")\n", " plt.legend(loc=3)\n", "\n", "plt.subplot(131)\n", "roc_plot(y_train, y_train_pred, y_val, y_val_pred)\n", "plt.subplot(132)\n", "pr_plot(y_train, y_train_pred, y_val, y_val_pred)\n", "plt.subplot(133)\n", "f1_plot(y_train, y_train_pred, y_val, y_val_pred)\n", "plt.subplots_adjust(right=2, wspace=0.25)\n", "plt.show()\n", "\n", "# Show classification reports\n", "print(\"Training statistics:\")\n", "print(classification_report(y_train, np.argmax(y_train_pred, axis=1), target_names=[\"ham\", \"spam\"]))\n", "print(\"Training classification accuracy: {:.3f}\".format(accuracy_score(y_train, np.argmax(y_train_pred, axis=1))))\n", "print(\"\\nValidation statistics:\")\n", "print(classification_report(y_val, np.argmax(y_val_pred, axis=1), target_names=[\"ham\", \"spam\"]))\n", "print(\"Validation classification accuracy: {:.3f}\".format(accuracy_score(y_val, np.argmax(y_val_pred, axis=1))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks pretty good at a first glance; 96% classification accuracy doesn't look too bad. But since there is a class imbalance, accuracy isn't everything. Also, we have to consider the risk of making misclassifications. Specifically, if we want to make sure we're catching all the spam (true positives), we want to make sure we have a high recall (true positive rate). However, we REALLY don't want to misclassify a ham message as spam (false positive), since if a legitimate message is accidentally deleted, it could have grave consequences, whereas if a spam message gets through the filter, it is simply a minor nusiance. As it is now, with a threshold of 0.5, our classifier will catch 91% of the spam, and incorrectly identifies ham messages 2% of the time, which doesn't sound bad. But the fraction of messages which it classifies as spam that are actually ham is as high as 20%!" ] }, { "cell_type": "code", "execution_count": 865, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fraction of ham messages incorrectly identified as spam (FPR): 0.0326\n", "Fraction of spam messages correctly identified (TPR): 0.9394\n", "Fraction of messages identified as spam which are actually ham (FDR): 0.2051\n" ] } ], "source": [ "thresh = 0.5\n", "print(\"Fraction of ham messages incorrectly identified as spam (FPR): {:.4f}\".format(np.mean(y_val_pred[y_val==0][:,1] > thresh)))\n", "print(\"Fraction of spam messages correctly identified (TPR): {:.4f}\".format(np.mean(y_val_pred[y_val==1][:,1] > thresh)))\n", "print(\"Fraction of messages identified as spam which are actually ham (FDR): {:.4f}\".format(np.mean(y_val[y_val_pred[:,1]>thresh] == 0)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our primary goal is to minimize the number of false positives. Since the ratio of ham:spam is already relatively large (about 6:1), we don't expect to be making a lot of positive identifications anyway. Therefore, the false positive rate ($FPR=FP/N$) is always going to be relatively small, regardless of performance. The false discovery rate ($FDR=FP/(TP+FP)$ or 1-precision), on the other hand, will be relatively large since the denominator $TP+FP$ is smaller than the total number of negative examples $N=TN+FP$. So minimizing the FDR should be our primary goal. This is equivalent to maximizing the precision ($TP/(TP+FP)$), but large relative changes in small FDR may not amount to significant changes in precision (like changing FDR from 0.01 -> 0.03 only changes precision from 0.99 -> 0.97). However, in the end, we really don't want to end up throwing a large fraction of legitimate messages away, so FPR will be useful to us as a constraint (e.g. don't want to discard any more than 1/1000 ham messages).\n", "\n", "Our secondary goal is to positively identify as much actual spam as possible. This is simply measured by the recall or true positive rate $TPR=TP/T=TP/(TP+FN)$.\n", "\n", "Ideally, we would be able to combine these two metrics (1/FDR and TPR) into a single one so that we may measure performance at a glance, but since $1/FDR \\in [1, \\infty)$ and $TPR \\in [0,1]$, this does not seem feasible. However, what we can do is try to maximize recall, while maintaining a constraint on the false positive rate. To evaluate our model using this metric, we can plot FPR vs TPR vs different thresholds (the typical ROC curve) on a semilog plot, decide on a maximum allowable cutoff for the FPR (which determines the classification threshold), and then check the TPR at this threshold. Let's see how our current model is evaluated using this metric:" ] }, { "cell_type": "code", "execution_count": 792, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEaCAYAAADg2nttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FfW9//HXJyEQdgigKGHTYtlkjVh3qYpob0UtKrhrlatVe6v13trWq5bW6k+9Vq3UihZRq1LUqlhRxIo7yiKCAiKLKAFZZIewZPn8/phJPIQk5yScyTkk7+fjcR6Z+c53Zj4ZQj75znfm+zV3R0REpCoZqQ5ARETSn5KFiIjEpWQhIiJxKVmIiEhcShYiIhKXkoWIiMSlZCEiInEpWYjUgJktN7MdZrbNzFab2Xgzaxaz/Wgze9PMtprZZjN72cx6ljtGCzO7z8y+Do+zNFxvW/vfkUjVlCxEau7H7t4M6Af0B34NYGZHAa8DLwEHA12BucD7ZnZIWKch8G+gFzAUaAEcBawHBtXutyESn+kNbpHqM7PlwBXu/ka4fhfQy91/ZGbvAp+6+8/K7fMqsM7dLzazK4DbgUPdfVsthy9SbWpZiOwjM8sFTgOWmFkT4Gjg2QqqTgROCZdPBl5TopD9hZKFSM29aGZbgRXAWuBWIIfg/9U3FdT/Bijtj2hTSR2RtKRkIVJzZ7p7c+BEoDtBItgIlAAHVVD/IODbcHl9JXVE0pKShcg+cve3gfHAPe6+HZgOnFNB1XMJOrUB3gBONbOmtRKkyD5SshBJjvuAU8ysL3ATcImZ/dzMmptZazP7A8HTTr8L6z9JcPvqeTPrbmYZZtbGzH5jZqen5lsQqZyShUgSuPs64AngFnd/DzgVOJugX+Irgkdrj3X3xWH9XQSd3J8DU4EtwAyCW1kf1fo3IBKHHp0VEZG41LIQEZG4IksWZjbOzNaa2WeVbDcze8DMlpjZPDMbELPtEjNbHH4uiSpGERFJTJQti/EEwxhU5jSgW/gZBTwEYGY5BM+rH0kw7MGtZtY6wjhFRCSOyJKFu78DbKiiyjDgCQ98CLQys4MIOganuvsGd99I0PlXVdIREZGIpbLPogPBo4Ol8sOyyspFRCRFGqQ6gH1hZqMIbmHRtGnTgd27d09xRCIisG1XEYVFJWXr+Zt21Og4zdhBV/uG1Z7DLrIAaNk4i4YNMveol5GZSXbTljU6x+zZs79193bx6qUyWawEOsas54ZlKwmGT4gtf6uiA7j7WGAsQF5ens+aNSuKOEWkHtq2q4iJM1ewo7C4WvsVlzj3Tv1ij7KDgOO6teWGUw6r1rGar3yX7025iC9+9DTbD8yjaaMGdDugGWZWreNUxcy+SqReKpPFJOBaM5tA0Jm92d2/MbMpwB9jOrWHEM4TICJSG778djtXPjGLJWtrPijwr4Z254x+BwNgwEEts6v/S74wmE/rsAObQ6fUPucTWbIws2cIWghtzSyf4AmnLAB3/yswGTgdWAIUAJeF2zaY2e+BmeGhRrt7VR3lIiJJM23RWn7+zBwaZBhPXXEkR3TJqfYxzCArs269xhZZsnD3kXG2O3BNJdvGAeOiiEtEpCLuzkNvL+XuKYvo3r4FYy8aSMecJqkOK23s1x3cIlK/3PjsXGZ8uYEk3rIvU1hUwqrNO/mPPgdx9/C+NG6YGX+nekTJQkTSys7CYnYVlnDGmPdYtWnHHrdzCnYHnc1nhn0ByTawSw4XHtkpqR3IdYWShYik3Pptu9iys4hFq7fws6c+piQc37Rpw0zOP7JTWb0MM849oiOHtmtWe8EtfROevRSKi2rvnKVKwnNa6vs/lCxEJKXWb9vFD+74N4XFQYZo2jCT/zq5Gw0yMjirfwdaN22Y2gDXfQE7N8MRV0CD7No/f6MWcFDf2j9vOUoWIpISn63czPL128nfuIPCYmfU8YfQ86AWHNKuKX1yW6U6vL0N/i00qf6TUXWFkoWIRGrO1xv5ZMWmPcrc4Z7XF5X1QQAMH5gbvE8gaUnJQkQi8c4X65i+bD0TZnzNxoLCvbZnGIy9aCBd2zalaaMGHNyqcQqilEQpWYhI0r39xTouGTcDgGaNGvDE5YPok7vn2EVZmRk0baRfQfsL/UuJSFKUlDi/eeFTVm3eyTtfrAPg5h/14IrjDklxZJIMShYiss+e+ugr/jJtKSvD0VX7dmzF0F7tlSjqECULkfpq+fvw3p/AS+LXrURRifPx1xvJ3V3MH4G2OY3o2rYpTbIyg1lpnkxatKmz9N+pjiAtKFmI1FeLJsOSN6DDwBrtvrFgN8vXF5AFNM40OuY0oVVjh+JtUL1RvdNbh4HQsiNk12y+iLpCyUKkPmvYFK5M/C/n4hLn89VbKCx2zhzzPgA//+H3uPL4Q2ienRVVlJIGlCxE6pGdhcV89OUGiktK+N6G7Rxc4rzz+ZqE93/knS+Zvmx92foPux/ADUO+H0WokmaULETqkdtfWciTHwYTo/22wRpGZhZz+fjqzzD5yMV5NMgwjjq0TbJDlDSlZCFSD2wuKGTCzK/LEsWzVx1F51lv0nhRJi9dcUy1jtW+ZTYHtkjBGEmSUkoWIvXAlAWruePVzwG45KjOwexvi7LBjL4d03AcJkk7ShYi9UBJOOb3u/8zWLO/SY1EmizMbChwP5AJPOrud5bb3plg+tR2wAbgQnfPD7cVA5+GVb929zOijFWkrvr7lPfYPusZrsncTatZ8yA7/G+fP7PqHUViRJYszCwTGAOcAuQDM81skrsviKl2D/CEuz9uZj8E7gAuCrftcPd+UcUnUh+4OxvefYSfN3gRsoAPylVIg3kSZP8QZctiELDE3ZcBmNkEYBgQmyx6AjeEy9OAFyOMR6TemTJ/DZmUUEQmDf63gkdkTfNMS2KinKuvA8EL/6Xyw7JYc4Gzw+WzgOZmVvosXraZzTKzD83szAjjFKmTFn6zhav+PhuAzIwMyMza+5OR+uk6Zf+Q6p+UG4ETzGwOcAKwku8GCujs7nnA+cB9ZnZo+Z3NbFSYUGatW7eu1oIWSXfFJc6tL80HYECn1liK45H9X5S3oVYCHWPWc8OyMu6+irBlYWbNgJ+4+6Zw28rw6zIzewvoDywtt/9YYCxAXl6eR/JdiKSpJWu3sfCbLRVue33BGmYs3wDAD7q2gdW1GZnURVEmi5lANzPrSpAkRhC0EsqYWVtgg7uXAL8meDIKM2sNFLj7rrDOMcBdEcYqsl/5av12Tr737bj1nr/6aGxx9d/QFikvsmTh7kVmdi0wheDR2XHuPt/MRgOz3H0ScCJwh5k58A5wTbh7D+BhMyshuFV2Z7mnqETqtTHTlgDB2Ey/Ob17hXVaNM7igObZsLg2I5O6KtL3LNx9MjC5XNktMcvPAc9VsN8HwOFRxiaS9op2w8rZ4N+N971tZxETZn3NjhWbOLlxMY+c0AIrmFvx/gXAemDzioq3i1SD3uAWSVezxsFrv9qjqBlwRWzB4wkeq3HrJAUl9ZWShUi62r01+HrRC6zeWsT1Ez8BoE9uK6774fdo1qga/31b5kYQoNQnShYi6a7LcXyy8FumlxRyzzl9GT5Qv/il9qX6PQsRqUQ49h+rN+/k7S++pWGDDIb2bp/aoKTeUstCJE1NXbCGU4Fj73qTIhowpOeB1bv1JJJE+skTSVPbdxUB8Icze+MZWZz4/XYpjkjqMyULkTS0YkMBS9dthywYcUSnYBwnkRRSn4VIunHn7tcWkkFJqiMRKaOWhUg6KdxJyf19eWDbasgCxzANAyhpQMlCJJ3s3k7GttW8WdyPTTl9OPuk4yFT/00l9fRTKJKG3irpy8UX/BEOaJbqUEQA9VmIiEgClCxE0siu4u8GDcxQV4WkEd2GEqklxSXO8vXb8Uqm6SosLuH8+//NnOxgvUubprUXnEgcShYiEXJ35uVvpmB3Mfe8vojZX22ssn7p2LBXHncIGWpaSBpRshCJ0AtzVnLDxD3nm3hgZP9K6zcp3AT/go6tm0Qdmki1KFmIJNPaz+GVX0LxbgB6rt7C8w2L6ZTThOysTBo1yKDhzCq6CksKaylQkepRshBJphUfwVfv4Z2OYsXWEjYXNyIzy2iXU43Jh7oNga7HRxejSA0oWYgkyb8XruHbmSs4D7h0y3/y9uqGAFx2TBfyftwrtcGJ7KNIH501s6FmtsjMlpjZTRVs72xm/zazeWb2lpnlxmy7xMwWh59LooxTJBken/4V8/I3A7B5RyHdDmjGGzecwK1KFFIHRNayMLNMYAxwCpAPzDSzSe6+IKbaPcAT7v64mf0QuAO4yMxygFuBPMCB2eG+VT9KIpJiB7XKhm3w4s+OgZYdUh2OSNJE2bIYBCxx92XuvhuYAAwrV6cn8Ga4PC1m+6nAVHffECaIqcDQCGMV2Sfujlf2AoVIHRBlsugArIhZzw/LYs0Fzg6XzwKam1mbBPfFzEaZ2Swzm7Vu3bqkBS5SXdf/4xPeXfwtGaZ3I6RuSnUH943Ag2Z2KfAOsBIornKPGO4+FhgLkJeXpz/rpFaUlDifr95KYfF38018tmoLXds25fTD28MHKQxOJCJRJouVQMeY9dywrIy7ryJsWZhZM+An7r7JzFYCJ5bb960IY5W6avd2+OBBKNxerd227ypm3dZdFW5bum4bi9du26NsOHBIu6Z0WaMWrtRNUSaLmUA3M+tKkCRGAOfHVjCztsAGdy8Bfg2MCzdNAf5oZqUPpw8Jt4tUz9cfwlt/hMyGYInddXUgs7CY9pVsbw8ckwlZmRnEzkuUsdVgK9CqEzSuxnsVIvuByJKFuxeZ2bUEv/gzgXHuPt/MRgOz3H0SQevhDjNzgttQ14T7bjCz3xMkHIDR7r4hqlilDivtdL50MnQ8Im71l+eu4tZJ89mwazed2zThwZEDKqyX06whHVo1TmakImkt0j4Ld58MTC5XdkvM8nPAc5XsO47vWhoikVm7dScPvbWU3UUlPPXR1wC0bJzF1OtPoGEDjeIvAqnv4BZJubcXreOx95fTukkWbZo25PazenNqr/aYnmwSKaNkIfXOhBlfM2Hmd09mr98edGS/fN2x5Gq0V5EKKVlIvfPa/NUsXbuN/p2DTugWjbM4oksO7VtkpzgykfSlZCH10iHtmvLE5YNSHYbIfkPJQvZfJSXw1XvBuxSVWfXJHquFxSUsW7edVk2yIg5OpG5RspD918rZ8PiPE6ub3YKC3UWc89fpfL2hgPYtc6KNTaSOiZsszKwx8Augs7tfZWbfA7q5+6uRRydSlcKC4OuP74eD+u61eWNBIXPzN1HUoCnbVzbj7nHvkL9xBwB3/aRPbUYqst9LpGUxDvgUODZcXwU8CyhZSHpo0w0O3nte6z+99BlPTF8HbAK+ux310W9O4kB1ZotUSyLJopu7jzSzcwDcvcD0ALrsB3YVltC2WUMm/udRZWXtmjeiebb6K0SqK5FksdvMsgmGzCEc62l3pFGJ7INJc1fx5sI1fPz1JhpkZHBIu2apDklkv5dIsvg98BqQa2aPAycAV0QalUg1/OGVBSxslFm2/v6S9QB0btOEwd3bpSoskTolbrJw91fNbBZwNMEYm//t7msjj0zqpd1FJVw+fiZrtuyMW7df0VzuBublb6ak03dzS+R1bs05ebmcd0SnCCMVqV8SeRrqdXcfArxUQZlIUm0q2M17S76l18Et6Nym6qE3OhQ0hgI494iODD/76FqKUKR+qjRZmFlDIBs40Mya893I/S0A/ckmkTr/yE5ccGTnqist2wZPwPABubUTlEg9VlXL4hrgBuAAYD7fJYstwF8jjktERNJIpcnC3f8E/MnMfuHu99ViTFJPrdhQwPRl61MdhohUIJEO7vvMrDvQk+C2VGn501EGJvXPf02Yw8dfbwKghd6FEEkriXRw30wwB3Z3gilSTwXeA5QsJKl2FJZwZNccfn9mb7odoHcjRNJJInNGngcMBr5x94uAvkDTRA5uZkPNbJGZLTGzmyrY3snMppnZHDObZ2anh+VdzGyHmX0SftRHUk+0aJzFYQc21yx1ImkmkZfydrh7sZkVhU9FrQbiPKYCZpYJjAFOAfKBmWY2yd0XxFS7GZjo7g+ZWU+C+bq7hNuWunu/anwvIiISkURaFnPMrBXBgIKzgBnhJ55BwBJ3X+buu4EJwLBydZzgUVyAlgSDFEo989pnqzlzzPt8+e22VIciIpWosmURDhh4m7tvAsaY2RSghbt/nMCxOwArYtbzgSPL1bkNeN3MriO4tXVyzLauZjaH4FHdm9393QriGwWMAujUSa9+7K/e/mItC1Zt4ahD23BW/w6pDkdEKlBlsnB3N7OpQO9wfUmSzz8SGO/u/2dmRwFPmllv4Bugk7uvN7OBwItm1svdt5SLbywwFiAvL8+THJsk2ctzV3Hjs3Mp8T3/qYpKnAObZ/O4pjkVSVuJ9Fl8Ymb93X1ONY+9EugYs54blsX6KTAUwN2nh6Pbtg3HntoVls82s6XAYQS3wWQ/tWTtNnYVlfCzEw/da1uf3JYpiEhEEpVIsuhP0Dm9FNhO8Ca3u/uAOPvNBLqFQ5qvBEYA55er8zVwEjDezHoQvMexzszaARvCjvVDgG7AskS/KUlv/zO0e6pDEJFqSiRZnFGTA7t7kZldS/BuRiYwzt3nm9loYJa7TwJ+CTxiZtcTdHZfGt76Oh4YbWaFQAlwlbtvqEkcIiKy7xJ5g3tpTQ/u7pMJHoeNLbslZnkBcEwF+z0PPF/T84qISHIl8uisiIjUc4nchhKpkb+8tYRFq7eWrX/+zdYqaotIOksoWZhZLtDN3aeZWSOggbtvjzY02d/dN3Ux2VkZ5DRtWFZ2co8Dk3eCwh3JO5aIVCmRgQQvB64leMP6UIKhPv7Cni/QiVTogh905ldRPP306XPwrxugUQtos/ejuCKSXIn0Wfwc+AHBm9S4+xcEEyKJ1L6dm+Gfo+D5n8IB3eGqd6F5+1RHJVLnJXIbaqe77y4dBTQcIFBDgsoetu0qYlPB7j3KnCS/VP/1h/DPK2HzSjjxN3DcLyFT3W4itSGR/2nvm9n/ANlmNphgutV/RRuW7E9KSpzj75rGhu2799rWMDMJD9wVF8E7d8E7d0PLjnD5a9BRQ4OI1KZEksX/EAzW9znwXwQv2T0cZVCyf3Fgw/bdDOl5ICf3/K4DO8OMk7rv4x3LDV8GrYn8mdB3JJx2F2S3iL+fiCRVIsniR8Cj7v5Q1MHI/q13h5acm9cxfsVEuMPcCTD5RrBM+Mnf4PDhyTm2iFRbIvcIzgGWmNlj4cx3mVEHJfXcjo3w3OXw4lXQvg9c/b4ShUiKJTLcx0XhuxU/Ai4DHjazV939qsijk/3PkjegYB+G8SosgLfvhm2r4Yf/C8deDxn6+0Qk1RJ6lMTdd5nZS8AOgkEBzwWULGQPzXauhr//ZN8PlHMIXP465A7c92OJSFIk8lLeKcB5BC/hvQc8wd5DjYvQoGRXsDDkdvj+aTU/UMuO0KBh/HoiUmsSaVmMAv4BXOfuGl9B4mt2oN6qFqljEumzOKc2AhERkfRVabIws7fd/QQz2wh7vIpbOlNeTuTRSUqUlDjLvt1GcUli9YtLNP25SF1XVcticPi1bW0EIunj7x99xS0vza/2fo2yND2KSF1VabJw99K/K//m7pfGbjOz8cClSJ20qaAQgD+P7E9mRmLDgGWYcXzOJvgwyshEJFUS6eDuE7sSvpR3RCIHN7OhwP0Ej9s+6u53ltveCXgcaBXWuSmcihUz+zXwU6AY+Lm7T6nyZMW7YdOKRMKScnYUFvOvuavYFd53WrZ8IwezgdM7FZFp1RgzcuO6iCIUkVSrqs/iV8BNQHMzK33Lygj6L/4W78BhUhkDnALkAzPNbFI473apm4GJ7v6QmfUkmK+7S7g8AugFHAy8YWaHuXtxpSdcMx/u6x0vLKlAY4LX9EtdCJBNkOZrQo+9itQ5VbUs7gL+D7iDIGkAUOUv7D0NApa4+zIAM5sADANik4UDpaPCtQRWhcvDgAnuvgv40syWhMebXunZWnWCM0YnGJrE+nTlZp748Cuu++H3aN8iG4DMDKteq6JUg2zodmqSIxSRVKsqWXzP3Reb2ZMEf+EDUDqvhbvPi3PsDkDsfaF84MhydW4DXjez64CmfDf7Xgf2vPudH5btwcxGEbwHQqdOnWDARXFCkorkZ33Ds+9/zE97HUfD9hrRVUT2VlWyuImgz2BMBdscOD4J5x8JjHf3/zOzo4AnzSzhe0nuPhYYC5CXl6fnN0VEIlLV01A/Db8eV8NjrwRix6vODcti/RQYGp5nupllEzyqm8i+IiJSS+I+GG9mZ5tZ83D5JjObaGZ9Ezj2TKCbmXU1s4YEHdaTytX5GjgpPHYPgm7VdWG9EWbWyMy6At2AGYl+UyIiklyJvEV1m7tvNbOjgdOBp0hgpjx3LwKuJZhZbyHBU0/zzWy0mZ0RVvslcKWZzQWeAS71wHxgIkFn+GvANdXoWBcRkSRL5D2L0l/S/wE87O4vmdltiRw8fGdicrmyW2KWFwDHVLLv7cDtiZxHRESilUiy+MbMxgCnAQPDW0oa10FEpB5J5Jf+ucDbwOnuvpGgA/qmqncREZG6JG6ycPdtwHzgRDO7Cmjt7q9GHpmIiKSNRJ6GuhZ4FugUfiaa2c+iDkxERNJHojPlDQpbGJjZH4EPgL9EGZiIiKSPRPosDNgds14YlomISD2RSMviSeAjM3ueIEmcSTCsuOynngonN3IPRkgpHSelRgMHiki9kMgc3HeZ2VvAsQS/V65y95lRBybRWbJ2G5lm/OeJh5aVtWycxaHtmqUwKhFJZ4m0LAB2AruAkvCrpKFNBbuZl785br2VG3fQKCuDXw75fi1EJSJ1QdxkYWa/Bc4HXiC4DfW0mT3l7ndEHZxUz+iXF/DPOYmNt3hwy+yIoxGRuiSRlsXFQH93LwAws9uBOQSTIkkaKdhdTMecxtx3Xr+4dTu0alILEYlIXZHQcB/l6jUIyyQNNclqwMDOOakOQ0TqmESSxQZgvplNIejgHkIwn/a9AO5+Q4TxiYhIGkgkWbwSfkp9WFlFERGpmxJ5dPZvtRGIiIikLw01LiIicSlZiIhIXAknCzNrFGUgIiKSvhIZonyQmX0KLA7X+5rZnxM5uJkNNbNFZrbEzPaaMMnM/mRmn4SfL8xsU8y24phtk6rxPYmISJIl8jTUAwTzb78I4O5zzWxwvJ3MLBMYA5wC5BM8bjspnHeb8FjXx9S/Dugfc4gd7h7/7TIREYlcIskiw92/sj1HJC1OYL9BwBJ3XwZgZhOAYcCCSuqPBG5N4Lj13tQFa5i7YtNe5V+s3UpWhrqhRCT5EkkWK8xsEOBha+E64IsE9usArIhZzweOrKiimXUGugJvxhRnm9ksoAi4091frGC/UQSTM9GpU6cEQqobfvfyfPI37iAzY+8hxX90+EEpiEhE6rpEksXVBLeiOgFrgDfCsmQaATzn7rEtls7uvtLMDgHeNLNP3X1p7E7uPhYYC5CXl+fUE+4wfGAu95zTN9WhiEg9kchLeWsJfplX10qgY8x6blhWkRHANeXOuzL8uiycT6M/sHTvXUVEJGqJDFH+CN9NplbG3UfF2XUm0M3MuhIkiREEQ52XP353oDUwPaasNVDg7rvMrC1wDHBXvFhFRCQaidyGeiNmORs4iz37Iirk7kVmdi0wBcgExrn7fDMbDcxy99LHYUcAE7x0js9AD+BhMysheLz3ztinqEREpHYlchvqH7HrZvYk8F4iB3f3ycDkcmW3lFu/rYL9PgAOT+QcIiISvZo8Z9kVODDZgYiISPpKpM9iI9/1WWQQzG+x19vYkjwlJc7rC9awdWdhhdu37y6q5YhEpL6rMllY8CZeX757iqmkXN+CRGDh6i1c9ffZVdZp06xhLUUjIhInWbi7m9lkd+9dWwEJ7C4qAeDu4X34wSFtKqzToVXj2gxJROq5RJ6G+sTM+rv7nMijkT20bd6IjjlNUh2GiEjlycLMGrh7EcHLcDPNbCmwHTCCRseAWopRRERSrKqWxQxgAHBGLcUiIiJpqqpkYQDlx2MSEZH6p6pk0c7Mbqhso7vfG0E8IiKShqpKFplAM8IWhoiI1F9VJYtv3H10rUUiIiJpq6rhPtSiEBERoOpkcVKtRSEiImmt0mTh7htqMxAREUlfNRl1VkRE6hklCxERiUvJQkRE4kpkIEFJsq07C9lRWFzp9k0FFc9jISKSKpEmCzMbCtxP8ILfo+5+Z7ntfwIGh6tNgAPcvVW47RLg5nDbH9z98ShjrS2rNu3g+LumUVQSf1qQhplq+IlIeogsWZhZJjAGOAXIJxi5dpK7Lyit4+7Xx9S/jmCEW8wsB7gVyCOYpW92uO/GqOKtLRu276aoxLngyE70OKhFpfWaNMzkiC45tRiZiEjlomxZDAKWuPsyADObAAwDFlRSfyRBggA4FZha+viumU0FhgLPRBhvrTrhsHYM6dU+1WGIiCQkyvscHYAVMev5YdlezKwz0BV4s7r7iohI9NLlpvgI4Dl3r7zXtwJmNsrMZpnZrHXr1kUUmoiIRJksVgIdY9Zzw7KKjGDPW0wJ7evuY909z93z2rVrt4/hiohIZaJMFjOBbmbW1cwaEiSESeUrmVl3oDUwPaZ4CjDEzFqbWWtgSFgmIiIpEFkHt7sXmdm1BL/kM4Fx7j7fzEYDs9y9NHGMACa4u8fsu8HMfk+QcABGa6wqEZHUifQ9C3efDEwuV3ZLufXbKtl3HDAusuBERCRh6dLBLSIiaUzJQkRE4lKyEBGRuJQsREQkLiULERGJS8lCRETiUrIQEZG4lCxERCQuJQsREYlLyUJEROJSshARkbgiHRtKRKQmCgsLyc/PZ+fOnakOpc7Izs4mNzeXrKysGu2vZJEES9ZuZfZXiU0PvnLjjoijEdn/5efn07x5c7p06YKZpTqc/Z67s379evLz8+natWuNjqFkkQQ3v/gZHy6r3gjqbZs3iigakf3fzp07lSiSyMxo06YN+zKjqJJFEhQWO0d0ac39I/onVL9RgwzaNFOyEKmKEkVy7ev1VLJIkkZCvlK2AAAQVElEQVQNMjm4VeNUhyEiEgk9DSUiUs769evp168f/fr1o3379nTo0KFsfffu3Qkd47LLLmPRokVV1hkzZgxPPfVUMkKOXL1vWRQWl3DBox+xenPNn7pYvWUng7rkJDEqEUmlNm3a8MknnwBw22230axZM2688cY96rg77k5GRsV/cz/22GNxz3PNNdfse7C1JNJkYWZDgfsJ5uB+1N3vrKDOucBtgANz3f38sLwY+DSs9rW7nxFFjFt2FDLjyw30zW3JIe2a1fg4Pzr8oCRGJSLpaMmSJZxxxhn079+fOXPmMHXqVH73u9/x8ccfs2PHDs477zxuuSWYOfrYY4/lwQcfpHfv3rRt25arrrqKV199lSZNmvDSSy9xwAEHcPPNN9O2bVt+8YtfcOyxx3Lsscfy5ptvsnnzZh577DGOPvpotm/fzsUXX8zChQvp2bMny5cv59FHH6Vfv361+r1HlizMLBMYA5wC5AMzzWySuy+IqdMN+DVwjLtvNLMDYg6xw91r7Wr8ZGAuFx/VpbZOJyIJ+t3L81mwaktSj9nz4Bbc+uNeNdr3888/54knniAvLw+AO++8k5ycHIqKihg8eDDDhw+nZ8+ee+yzefNmTjjhBO68805uuOEGxo0bx0033bTXsd2dGTNmMGnSJEaPHs1rr73Gn//8Z9q3b8/zzz/P3LlzGTBgQI3i3ldR9lkMApa4+zJ33w1MAIaVq3MlMMbdNwK4+9oI4xER2WeHHnpoWaIAeOaZZxgwYAADBgxg4cKFLFiwYK99GjduzGmnnQbAwIEDWb58eYXHPvvss/eq89577zFixAgA+vbtS69eNUty+yrK21AdgBUx6/nAkeXqHAZgZu8T3Kq6zd1fC7dlm9ksoAi4091fjDBWEUlTNW0BRKVp06Zly4sXL+b+++9nxowZtGrVigsvvLDCt84bNmxYtpyZmUlRUVGFx27UqFHcOqmS6qehGgDdgBOBkcAjZtYq3NbZ3fOA84H7zOzQ8jub2Sgzm2Vms/blZRMRkZrYsmULzZs3p0WLFnzzzTdMmTIl6ec45phjmDhxIgCffvpphS2X2hBly2Il0DFmPTcsi5UPfOTuhcCXZvYFQfKY6e4rAdx9mZm9BfQHlsbu7O5jgbEAeXl5HsU3ISJSmQEDBtCzZ0+6d+9O586dOeaYY5J+juuuu46LL76Ynj17ln1atmyZ9PPEY+7R/I41swbAF8BJBEliJnC+u8+PqTMUGOnul5hZW2AO0A8oAQrcfVdYPh0YFts5Xl5eXp7PmjWr2nGu37aLgX94g9HDeqmDWyRNLFy4kB49eqQ6jLRQVFREUVER2dnZLF68mCFDhrB48WIaNKj+3/oVXVczmx3exalSZC0Ldy8ys2uBKQT9EePcfb6ZjQZmufukcNsQM1sAFAP/7e7rzexo4GEzKyG4VXZnVYlCRKSu2rZtGyeddBJFRUW4Ow8//HCNEsW+ivSM7j4ZmFyu7JaYZQduCD+xdT4ADo8yNhGR/UGrVq2YPXt2qsNIeQe3iIjsB5QsREQkLiULERGJS8lCRETiUrIQESln8ODBe71gd99993H11VdXuk+zZsFApKtWrWL48OEV1jnxxBOJ94j/fffdR0FBQdn66aefzqZNmxINPTJKFiIi5YwcOZIJEybsUTZhwgRGjhwZd9+DDz6Y5557rsbnLp8sJk+eTKtWrarYo3YoWYiIlDN8+HBeeeWVsomOli9fzqpVq+jfvz8nnXQSAwYM4PDDD+ell17aa9/ly5fTu3dvAHbs2MGIESPo0aMHZ511Fjt27Cird/XVV5OXl0evXr249dZbAXjggQdYtWoVgwcPZvDgwQB06dKFb7/9FoB7772X3r1707t3b+67776y8/Xo0YMrr7ySXr16MWTIkD3Okyz1fvIjEUlzr94Eqz+NX6862h8Op+01vU6ZnJwcBg0axKuvvsqwYcOYMGEC5557Lo0bN+aFF16gRYsWfPvtt/zgBz/gjDPOqHR+64ceeogmTZqwcOFC5s2bt8fw4rfffjs5OTkUFxdz0kknMW/ePH7+859z7733Mm3aNNq2bbvHsWbPns1jjz3GRx99hLtz5JFHcsIJJ9C6dWsWL17MM888wyOPPMK5557L888/z4UXXpicaxVSy0JEpAKxt6JKb0G5O7/5zW/o06cPJ598MitXrmTNmjWVHuOdd94p+6Xdp08f+vTpU7Zt4sSJDBgwgP79+zN//vy4AwS+9957nHXWWTRt2pRmzZpx9tln8+677wLQtWvXssmQqhoCfV+oZSEi6a2KFkCUhg0bxvXXX8/HH39MQUEBAwcOZPz48axbt47Zs2eTlZVFly5dKhySPJ4vv/ySe+65h5kzZ9K6dWsuvfTSGh2nVOnQ5hAMbx7FbSi1LEREKtCsWTMGDx7M5ZdfXtaxvXnzZg444ACysrKYNm0aX331VZXHOP7443n66acB+Oyzz5g3bx4QDG3etGlTWrZsyZo1a3j11VfL9mnevDlbt27d61jHHXccL774IgUFBWzfvp0XXniB4447LlnfblxqWYiIVGLkyJGcddZZZbejLrjgAn784x9z+OGHk5eXR/fu3avc/+qrr+ayyy6jR48e9OjRg4EDBwLBjHf9+/ene/fudOzYcY+hzUeNGsXQoUM5+OCDmTZtWln5gAEDuPTSSxk0aBAAV1xxBf3794/kllNFIhuivLZpiHKRukNDlEdjX4Yo120oERGJS8lCRETiUrIQkbRUV26Rp4t9vZ5KFiKSdrKzs1m/fr0SRpK4O+vXryc7O7vGx9DTUCKSdnJzc8nPz2fdunWpDqXOyM7OJjc3t8b7K1mISNrJysqia9euqQ5DYkR6G8rMhprZIjNbYmY3VVLnXDNbYGbzzezpmPJLzGxx+LkkyjhFRKRqkbUszCwTGAOcAuQDM81skrsviKnTDfg1cIy7bzSzA8LyHOBWIA9wYHa478ao4hURkcpF2bIYBCxx92XuvhuYAAwrV+dKYExpEnD3tWH5qcBUd98QbpsKDI0wVhERqUKUfRYdgBUx6/nAkeXqHAZgZu8DmcBt7v5aJft2KH8CMxsFjApXt5nZonJVWgKb48TZEth8yf+DSu51VXaMisoTKWsLfBsnpmRK5Bokc/+Er3k1tyVaXlG9+n7Na7q9vvyMV/cYUf2MV7Yt6mveOaFa7h7JBxgOPBqzfhHwYLk6/wJeALKArgQJohVwI3BzTL3/BW6sQQxj97VOZdsrKk+kDJgV1TWv6TVI5v77es2rc70rub4V/RvU62uun/HkHiOqn/F0v+ZR3oZaCXSMWc8Ny2LlA5PcvdDdvwS+ALoluG8iXk5Cncq2V1SeaFlt2tfzV3f/fb3m1bneFZWn+npD+l1z/Ywn9xhR/YxXti0trnlkAwmaWQOCX/4nEfyinwmc7+7zY+oMBUa6+yVm1haYA/Qj7NQGSqeV+hgY6O4bIgm2FpnZLE9g0C5JHl3z2qXrXftq45pH1mfh7kVmdi0whaA/Ypy7zzez0QRNpknhtiFmtgAoBv7b3dcDmNnvCRIMwOi6kChCY1MdQD2ka167dL1rX+TXvM4MUS4iItHR2FAiIhKXkoWIiMSlZCEiInEpWaQRM+thZn81s+fM7OpUx1MfmNmZZvaImf3DzIakOp66zswOMbO/mdlzqY6lLjOzpmb2ePizfUEyjqlkkSRmNs7M1prZZ+XK4w6mWMrdF7r7VcC5wDFV1ZWkXfMX3f1K4CrgvCjj3d8l6Xovc/efRhtp3VTN63828Fz4s31GMs6vZJE84yk3flXMYIqnAT2BkWbW08wON7N/lfuUDqJ4BvAKMLl2w98vjScJ1zx0c7ifVG48ybveUn3jSfD6E7zIXDpkUnEyTq75LJLE3d8xsy7lissGUwQwswnAMHe/A/iPSo4zCZhkZq8AT1dURwLJuOZmZsCdwKvu/nG0Ee/fkvUzLjVTnetPMDpGLvAJSWoUqGURrYQGRCxlZiea2QNm9jBqWdRUta45cB1wMjDczK6KMrA6qro/423M7K9AfzP7ddTB1QOVXf9/Aj8xs4dI0tAgalmkEXd/C3grxWHUK+7+APBAquOoL8IRGpSUI+bu24HLknlMtSyilawBESVxuua1S9c7tWrt+itZRGsm0M3MuppZQ2AEMCnFMdV1uua1S9c7tWrt+itZJImZPQNMB75vZvlm9lN3LwJKB1NcCEyMHXVX9o2uee3S9U6tVF9/DSQoIiJxqWUhIiJxKVmIiEhcShYiIhKXkoWIiMSlZCEiInEpWYiISFxKFrLfMbNiM/sk5tOlirpdyg/pnCpmlmdmD4TLJ5rZ0THbrjKzi2sxln5mdnptnU/2fxobSvZHO9y9X6qDqC53nwXMCldPBLYBH4Tb/prs85lZg/ClrYr0A/LQgJWSILUspE4IWxDvmtnH4efoCur0MrMZYWtknpl1C8svjCl/OJwjoPy+y83sLjP7NKz7vZjzvhke799m1iksP8fMPjOzuWb2Tlh2YjivQxeCwfSuD895nJndZmY3mll3M5tR7vv6NFweaGZvm9lsM5tiZgdVEOd4C2Zb/Ai4y8wGmdl0M5tjZh+Y2ffDYSFGA+eF5z/PgpnVxoXf2xwzG7bP/yhSt7i7PvrsVx+CyVw+CT8vhGVNgOxwuRswK1zuAnwWLv8ZuCBcbgg0BnoQDOGcFZb/Bbi4gnMuB34bLl8M/Ctcfhm4JFy+HHgxXP4U6BAutwq/nhiz323AjTHHL1sPv6+u4fKvCCZmyiJohbQLy88DxlUQ53jgX0BmuN4CaBAunww8Hy5fCjwYs98fgQtL4wW+AJqm+t9an/T56DaU7I8qug2VBTxoZv0IkslhFew3HfitmeUC/3T3xWZ2EjAQmBnMg0RjYG0l530m5uufwuWjCKawBHgSuCtcfh8Yb2YTCeYWqI6JBMngzvDrecD3gd7A1DDOTOCbSvZ/1t1LZ0drCTwetqKc4DpVZAhwhpndGK5nA50IxhsSUbKQOuN6YA3Ql+D26s7yFdz96fD2zI+AyWb2n4ABj7t7IhPxeCXLe1d0v8rMjgzPNdvMBib2bQDwD+BZM/tncChfbGaHA/Pd/agE9t8es/x7YJq7nxXe/nqrkn0M+Im7L6pGnFKPqM9C6oqWwDfuXgJcRPCX9x7M7BBgmQcTHr0E9AH+TTBLXukc6Dlm1rmSc5wX83V6uPwBwbDQABcA74bHOdTdP3L3W4B17DnnAMBWoHlFJ3H3pQSto/8lSBwAi4B2ZnZUePwsM+tVSZyxWvLd/AaXVnH+KcB1FjZbzKx/AseWekTJQuqKvwCXmNlcoDt7/nVd6lzgMzP7hOCWzhPuvoCgT+B1M5sHTAX26jgOtQ7r/BdBSwaCaVkvC8svCrcB3B12hn9GkFDmljvWy8BZpR3cFZzrH8CFBLekcPfdwHDg/4Xf4yfAXp34FbgLuMPM5rDnnYRpQM/SDm6CFkgWMM/M5ofrImU0RLlIAsxsOZDn7t+mOhaRVFDLQkRE4lLLQkRE4lLLQkRE4lKyEBGRuJQsREQkLiULERGJS8lCRETiUrIQEZG4/j8XeH9D8eGEjAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fpr, tpr, _ = roc_curve(y_train, y_train_pred[:,1])\n", "plt.semilogx(fpr, tpr, label=\"Training\")\n", "fpr, tpr, _ = roc_curve(y_val, y_val_pred[:,1])\n", "plt.semilogx(fpr, tpr, label=\"Validation\")\n", "plt.ylim(ymin=0.6, ymax=1)\n", "plt.xlabel(\"False positive rate\")\n", "plt.ylabel(\"True positive rate\")\n", "plt.title(\"ROC\")\n", "plt.legend(loc=4)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I think a false positive rate of 1 in 100 (0.01) is probably a pretty good constraint (also our validation set only has ~500 samples in it so it's impossible to measure an FPR less than ~1/500). We'll use that from now on. Looks like this corresponds to a recall of roughly 0.80-0.85. Let's write a function to do the evaluation for us:" ] }, { "cell_type": "code", "execution_count": 793, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Training - FPR=9.852e-03, TPR=0.875, thresh=0.824\n", "Validation - FPR=8.147e-03, TPR=0.879, thresh=0.880\n" ] } ], "source": [ "def evaluate_model(model, X, y, FPR_max=0.01):\n", " # Find y_pred\n", " y_pred = model.predict_proba(X)\n", " \n", " # Get FPR, TPR, and thresholds from sklearn's roc_curve function\n", " fpr, tpr, thresh = roc_curve(y, y_pred[:,1])\n", " \n", " # Find the threshold with the largest FPR below the max\n", " for i in reversed(range(len(fpr))):\n", " if fpr[i] <= FPR_max: break\n", " return fpr[i], tpr[i], thresh[i]\n", "print(\"Training - FPR={:.3e}, TPR={:.3f}, thresh={:.3f}\".format(*evaluate_model(model[1], X_train_sparse, y_train)))\n", "print(\"Validation - FPR={:.3e}, TPR={:.3f}, thresh={:.3f}\".format(*evaluate_model(model[1], X_val_sparse, y_val)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have our metric down we should start tuning hyperparameters. I'll just roll the dice on a random set of hyperparameters a number of times, evaluate them using k-fold cross-validation, then pick the best model." ] }, { "cell_type": "code", "execution_count": 830, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Model 0: FPR=9.174e-03, TPR=0.849±0.069, thresh=0.921\n", "Model 2: FPR=9.335e-03, TPR=0.859±0.070, thresh=0.532\n", "Model 17: FPR=9.259e-03, TPR=0.860±0.041, thresh=0.774\n", "Model 87: FPR=5.807e-03, TPR=0.866±0.052, thresh=0.633\n", "Model 94: FPR=8.226e-03, TPR=0.873±0.056, thresh=0.679\n" ] } ], "source": [ "from sklearn.model_selection import KFold\n", "kf = KFold(n_splits=5, shuffle=True, random_state=None)\n", "\n", "best_hyperparameters = None\n", "best_tpr = 0\n", "X = sp.sparse.vstack([X_train_sparse, X_val_sparse]).todok()\n", "y = np.concatenate([y_train, y_val], axis=0)\n", "\n", "for _ in range(100):\n", " # Pick hyperparameters\n", " C = 10**np.random.uniform(-2,2)\n", " spam_weight = np.random.uniform(0.5, 10)\n", " penalty = np.random.choice([\"l1\", \"l2\"])\n", " \n", " # Instantiate, train, and evaluate model using k-fold CV\n", " tpr_list = []\n", " for train_idx, val_idx in kf.split(X):\n", " temp_model = LogisticRegression(penalty=penalty, C=C, solver=\"liblinear\", class_weight={0:1, 1:spam_weight})\n", " temp_model.fit(X[train_idx], y[train_idx])\n", " fpr, tpr, thresh = evaluate_model(temp_model, X[val_idx], y[val_idx])\n", " tpr_list.append(tpr)\n", " if np.mean(tpr_list) > best_tpr:\n", " best_tpr = np.mean(tpr_list)\n", " best_hyperparameters = (C, spam_weight, penalty)\n", " print(\"Model {}: FPR={:.3e}, TPR={:.3f}\\u00B1{:.3f}, thresh={:.3f}\".format(_, fpr, best_tpr, 2*np.std(tpr_list), thresh), end=\"\\n\")" ] }, { "cell_type": "code", "execution_count": 870, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "C=9.16e-01, spam_weight=3.00, penalty=l1\n", "Test set: FPR=0.4864, TPR=0.988, threshold=0.013\n", "Test set: FPR=0.0964, TPR=0.938, threshold=0.104\n", "Test set: FPR=0.0042, TPR=0.825, threshold=0.857\n", "Test set: FPR=0.0000, TPR=0.750, threshold=0.964\n" ] } ], "source": [ "C, spam_weight, penalty = best_hyperparameters\n", "model[2] = LogisticRegression(penalty=penalty, C=C, solver=\"liblinear\", class_weight={0:1, 1:spam_weight})\n", "model[2].fit(X_train_sparse, y_train)\n", "print(\"C={:.2e}, spam_weight={:.2f}, penalty={}\".format(*best_hyperparameters))\n", "for FPR_max in [1/2, 1/10, 1/100, 1/1000]:\n", " print(\"Test set: FPR={:.4f}, TPR={:.3f}, threshold={:.3f}\".format(*evaluate_model(model[2], X_test_sparse, y_test, FPR_max=FPR_max)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have a recall of 82.5% and a false positive rate of <1% on the test set (which has been completely blinded up until this point). Looking around for other studies/implementations to benchmark against, I found [this study](https://arxiv.org/abs/cs/0009009), which compared several spam filtering techniques. For comparison, our constraint of max FPR of 1/100 would be roughly comparable to one of their models with parameter $\\lambda = 100$ (this comparison falls apart for $\\lambda\\sim 1$). The value of $\\lambda$ can be understood as the relative cost of misclassifying ham as spam vs misclassifying spam as ham - it is $\\lambda$ times worse to make a false positive identification than a false negative. As we can see from their Table 1, our logistic regression model outperforms all the models they considered with regards to spam recall (TPR) for equivalent values of $\\lambda$." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }