lectures/lec21/cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CS 429: Information Retrieval\n",
    "<br>\n",
    "\n",
    "## Clustering Words\n",
    "\n",
    "<br>\n",
    "\n",
    "### Dr. Aron Culotta\n",
    "### Illinois Institute of Technology\n",
    "### Spring 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Motivation\n",
    "\n",
    "Often, we want to know which features appear together.\n",
    "\n",
    "- If you liked *Twilight* you might like *Nosferatu*.\n",
    "- \"happy\" is a synonym of \"glad.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use k-means to cluster together related words from Twitter.\n",
    "\n",
    "**Caution:** This uses live Twitter data, which often contains profanity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n",
      "700\n",
      "800\n",
      "900\n",
      "1000\n",
      "1100\n",
      "1200\n",
      "1300\n",
      "1400\n",
      "1500\n",
      "1600\n",
      "1700\n",
      "1800\n",
      "1900\n",
      "2000\n",
      "2100\n",
      "2200\n",
      "2300\n",
      "2400\n",
      "2500\n",
      "2600\n",
      "2700\n",
      "2800\n",
      "2900\n",
      "3000\n",
      "3100\n",
      "3200\n",
      "3300\n",
      "3400\n",
      "3500\n",
      "3600\n",
      "3700\n",
      "3800\n",
      "3900\n",
      "4000\n",
      "4100\n",
      "4200\n",
      "4300\n",
      "4400\n",
      "4500\n",
      "4600\n",
      "4700\n",
      "4800\n",
      "4900\n",
      "5000\n"
     ]
    }
   ],
   "source": [
    "# Get some tweets containing the word 'i'.\n",
    "\n",
    "import os\n",
    "from TwitterAPI import TwitterAPI\n",
    "\n",
    "# Read Twitter credentials from environmental variables.\n",
    "api = TwitterAPI(os.environ.get('TW_CONSUMER_KEY'),\n",
    "                 os.environ.get('TW_CONSUMER_SECRET'),\n",
    "                 os.environ.get('TW_ACCESS_TOKEN'),\n",
    "                 os.environ.get('TW_ACCESS_TOKEN_SECRET'))\n",
    "\n",
    "# Collect 5000 tweets.\n",
    "tweets = []\n",
    "while True: \n",
    "    r = api.request('statuses/filter', {'track':'i',\n",
    "                                      'language':'en'})\n",
    "    if r.status_code != 200: # error\n",
    "        break\n",
    "    else:\n",
    "        for item in r.get_iterator():\n",
    "            tweets.append(item)\n",
    "            if len(tweets) > 5000:\n",
    "                break\n",
    "            elif len(tweets) % 100 == 0:\n",
    "                print(len(tweets))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5003\n"
     ]
    }
   ],
   "source": [
    "print(len(tweets))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "text @Maddy_Valone thanks madelynn, I love you💖\n",
      "description: actions speak louder than words// 10.23.15❤️ 03.07.16 AB❤️❤️\n",
      "name: Tay\n",
      "location: None\n"
     ]
    }
   ],
   "source": [
    "# Each tweet is a Python dict.\n",
    "print('text', tweets[0]['text'])\n",
    "print('description:', tweets[0]['user']['description'])\n",
    "print('name:', tweets[0]['user']['name'])\n",
    "print('location:', tweets[0]['user']['location'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tweets = [t for t in tweets if 'text' in t]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4623"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(tweets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['thanks', 'madelynn', 'i', 'love', 'you']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Tokenize each tweet text.\n",
    "import re\n",
    "tokens = []\n",
    "for tweet in tweets:\n",
    "    text = tweet['text'].lower()\n",
    "    text = re.sub('@\\S+', ' ', text)  # Remove mentions.\n",
    "    text = re.sub('http\\S+', ' ', text)  # Remove urls.\n",
    "    tokens.append(re.findall('[A-Za-z]+', text)) # Retain words.\n",
    "tokens[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Count words.\n",
    "from collections import Counter\n",
    "\n",
    "word_counts = Counter()\n",
    "for tweet in tokens:\n",
    "    word_counts.update(tweet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7502 unique terms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[('i', 5702),\n",
       " ('rt', 1845),\n",
       " ('to', 1465),\n",
       " ('the', 1384),\n",
       " ('a', 1220),\n",
       " ('you', 1085),\n",
       " ('m', 964),\n",
       " ('and', 941),\n",
       " ('my', 884),\n",
       " ('t', 780)]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Inspect word counts.\n",
    "import math\n",
    "\n",
    "print(len(word_counts), 'unique terms')\n",
    "word_counts.most_common(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3145 words occur at least two times.\n"
     ]
    }
   ],
   "source": [
    "# Retain in vocabulary words occurring more than twice.\n",
    "vocab = set([w for w, c in word_counts.items() if c > 1])\n",
    "print('%d words occur at least two times.' % len(vocab))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Prune tokens.\n",
    "newtoks = []\n",
    "for i, tweet in enumerate(tokens):\n",
    "    newtok = [token for token in tweet if token in vocab]\n",
    "    if len(newtok) > 0:\n",
    "        newtoks.append(newtok)\n",
    "tokens = newtoks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['thanks', 'i', 'love', 'you']"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A sample pruned tweet.\n",
    "tokens[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# For each term, create a context vector, indicating how often\n",
    "# each word occurs to the left or right of it.\n",
    "from collections import defaultdict\n",
    "import numpy as np\n",
    "\n",
    "# dict from term to context vector.\n",
    "contexts = defaultdict(lambda: Counter())\n",
    "window = 2\n",
    "for tweet in tokens:\n",
    "    for i, token in enumerate(tweet):\n",
    "        features = []\n",
    "        for j in range(np.amax([0, i-window]), i):\n",
    "            features.append(tweet[j] + \"@\" + str(j-i))\n",
    "        for j in range(i+1, min(i + window, len(tweet))):\n",
    "            features.append(tweet[j] + \"@\" + str(j-i))\n",
    "        contexts[token].update(features)\n",
    "        # Optionally: ignore word order\n",
    "        # contexts[token].update(tweet[:i] + tweet[i+1:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('m@1', 949),\n",
       " ('rt@-1', 791),\n",
       " ('and@-1', 233),\n",
       " ('love@1', 219),\n",
       " ('ve@1', 209),\n",
       " ('can@1', 208),\n",
       " ('rt@-2', 201),\n",
       " ('don@1', 172),\n",
       " ('have@1', 169),\n",
       " ('just@1', 166),\n",
       " ('but@-1', 164),\n",
       " ('ll@1', 162),\n",
       " ('when@-1', 139),\n",
       " ('the@-2', 136),\n",
       " ('am@1', 134),\n",
       " ('was@1', 125),\n",
       " ('me@-2', 118),\n",
       " ('want@1', 110),\n",
       " ('i@-2', 109),\n",
       " ('hate@1', 106)]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contexts['i'].most_common(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('i@-1', 5675),\n",
       " ('i@-2', 5464),\n",
       " ('i@1', 4350),\n",
       " ('rt@-1', 1844),\n",
       " ('rt@-2', 1838)]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute the number of different contexts each term appears in.\n",
    "tweet_freq = Counter()\n",
    "for context in contexts.values():\n",
    "    tweet_freq.update(context)\n",
    "tweet_freq.most_common(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('m@1', 0.6074600339659124),\n",
       " ('rt@-1', 0.46777716298830435),\n",
       " ('love@1', 0.16287334308387164),\n",
       " ('ve@1', 0.16207223522779549),\n",
       " ('and@-1', 0.14995192581405317)]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Transform each context vector to be term freq / tweet frequency. \n",
    "# Also then normalize by length.\n",
    "for term, context in contexts.items():\n",
    "    for term2, frequency in context.items():\n",
    "        context[term2] = frequency / (1. + math.log(tweet_freq[term2]))\n",
    "    length = math.sqrt(sum([v*v for v in context.values()]))\n",
    "    for term2, frequency in context.items():\n",
    "        context[term2] = 1. * frequency / length\n",
    "    \n",
    "contexts['i'].most_common(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('fake@-2', 0.5975239422063642),\n",
       " ('friends@-1', 0.5372074585396509),\n",
       " ('being@1', 0.4509488564161085),\n",
       " ('is@1', 0.15957649683887706),\n",
       " ('go@-2', 0.11810454155883642),\n",
       " ('rt@-1', 0.10816530593283122),\n",
       " ('to@-1', 0.09752316082466882),\n",
       " ('middle@-1', 0.08828866275699869),\n",
       " ('project@1', 0.0720574292682831),\n",
       " ('incident@-2', 0.06803418701026041)]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contexts['school'].most_common(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we have a list of dictionaries, one per term, indicating the terms that co-occur (weighted by inverse tweet frequency).\n",
    "\n",
    "Next, we have to cluster these vectors. To do this, we'll need to be able to compute the euclidean distance between two vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.4142135623730951\n",
      "2.23606797749979\n"
     ]
    }
   ],
   "source": [
    "# This is not efficient!\n",
    "def distance(c1, c2):\n",
    "    if len(c1.keys()) == 0 or len(c2.keys()) == 0:\n",
    "        return 1e9\n",
    "    keys = set(c1.keys()) | set(c2.keys())\n",
    "    distance = 0.\n",
    "    for k in keys:\n",
    "        distance += (c1[k] - c2[k]) ** 2\n",
    "    return math.sqrt(distance)\n",
    "\n",
    "print(distance({'hi':10, 'bye': 5}, {'hi': 9, 'bye': 4}))\n",
    "print(distance({'hi':10, 'bye': 5}, {'hi': 8, 'bye': 4}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['love', 'hope', 'miss', 'm', 'swear', 'am', 'll', 'bet', 've',\n",
       "       'hate'], \n",
       "      dtype='<U22')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def find_closest(term, n=5):\n",
    "    terms = np.array(list(contexts.keys()))\n",
    "    context = contexts[term]\n",
    "    distances = []\n",
    "    for term2, context2 in contexts.items():\n",
    "        distances.append(distance(context, context2))\n",
    "    return terms[np.argsort(distances)][:n]\n",
    "\n",
    "find_closest('love', n=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3143 nonzero contexts\n"
     ]
    }
   ],
   "source": [
    "nz_contexts = [t for t, context in contexts.items()\n",
    "               if len(context) > 1]\n",
    "contexts = dict([(term, contexts[term]) for term in nz_contexts])\n",
    "print(len(nz_contexts), 'nonzero contexts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['a@-1' 'a@-2' 'a@1' 'abgehakt@-1' 'abgehakt@-2' 'abgehakt@1' 'abi@1'\n",
      " 'ability@-1' 'ability@-2' 'ability@1']\n",
      "  (0, 2730)\t0.356687862264\n",
      "  (0, 5155)\t0.660034013128\n",
      "  (0, 6855)\t0.238859146112\n",
      "  (0, 7255)\t0.425839590179\n",
      "  (0, 8208)\t0.207871150335\n",
      "  (0, 8237)\t0.394367096803\n"
     ]
    }
   ],
   "source": [
    "# Transform context dicts to a sparse vector\n",
    "# for sklearn.\n",
    "from sklearn.feature_extraction import DictVectorizer\n",
    "\n",
    "vec = DictVectorizer()\n",
    "X = vec.fit_transform(contexts.values())\n",
    "names = np.array(vec.get_feature_names())\n",
    "print(names[:10])\n",
    "print(X[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "828\n",
      "  (0, 1)\t0.00477244645424\n",
      "  (0, 2)\t0.00940532605072\n",
      "  (0, 15)\t0.00614343176802\n",
      "  (0, 74)\t0.00835941470593\n",
      "  (0, 79)\t0.0181453447983\n",
      "  (0, 138)\t0.00729571555916\n",
      "  (0, 177)\t0.00841294973494\n",
      "  (0, 201)\t0.0422557059747\n",
      "  (0, 230)\t0.00886440634587\n",
      "  (0, 244)\t0.00679373548211\n",
      "  (0, 246)\t0.00681831071018\n",
      "  (0, 251)\t0.0201877524318\n",
      "  (0, 268)\t0.00629437188756\n",
      "  (0, 269)\t0.0190704309034\n",
      "  (0, 274)\t0.00972830351336\n",
      "  (0, 275)\t0.029380834053\n",
      "  (0, 276)\t0.00974457462836\n",
      "  (0, 400)\t0.0104642604187\n",
      "  (0, 465)\t0.00650715082774\n",
      "  (0, 548)\t0.0181453447983\n",
      "  (0, 590)\t0.0180522229886\n",
      "  (0, 591)\t0.00920809982738\n",
      "  (0, 631)\t0.022490687173\n",
      "  (0, 731)\t0.00704309616134\n",
      "  (0, 762)\t0.0129629837028\n",
      "  :\t:\n",
      "  (0, 8846)\t0.022490687173\n",
      "  (0, 8853)\t0.0292866721661\n",
      "  (0, 8864)\t0.00993423536806\n",
      "  (0, 8876)\t0.006409009642\n",
      "  (0, 8897)\t0.00679373548211\n",
      "  (0, 8911)\t0.0186546052707\n",
      "  (0, 8912)\t0.012498065172\n",
      "  (0, 8913)\t0.00620620836335\n",
      "  (0, 8963)\t0.00566090463241\n",
      "  (0, 8964)\t0.0387255296187\n",
      "  (0, 8970)\t0.00941521382615\n",
      "  (0, 9057)\t0.0471457683681\n",
      "  (0, 9067)\t0.00965410893082\n",
      "  (0, 9098)\t0.0859912121048\n",
      "  (0, 9107)\t0.0247317853321\n",
      "  (0, 9113)\t0.00886440634587\n",
      "  (0, 9131)\t0.00846880452335\n",
      "  (0, 9137)\t0.00894297328297\n",
      "  (0, 9161)\t0.0102695598793\n",
      "  (0, 9169)\t0.0193139622343\n",
      "  (0, 9170)\t0.0195729138021\n",
      "  (0, 9171)\t0.316065755412\n",
      "  (0, 9181)\t0.00610605274655\n",
      "  (0, 9182)\t0.00625534964522\n",
      "  (0, 9183)\t0.0366051781502\n",
      "you@1\n"
     ]
    }
   ],
   "source": [
    "# Which row of X is the word \"love\"?\n",
    "love_idx = list(contexts.keys()).index('love')\n",
    "print(love_idx)\n",
    "# What are the context feature values for love?\n",
    "print(X[love_idx])\n",
    "# Print a highly ranking feature.\n",
    "print(names[9171])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KMeans(copy_x=True, init='k-means++', max_iter=300, n_clusters=20, n_init=10,\n",
       "    n_jobs=1, precompute_distances='auto', random_state=None, tol=0.0001,\n",
       "    verbose=0)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's cluster!\n",
    "# http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html\n",
    "from sklearn.cluster import KMeans\n",
    "num_clusters = 20\n",
    "kmeans = KMeans(num_clusters)\n",
    "kmeans.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 to@-1 t@-1 i@-1 i@-2 me@1\n",
      "1 and@1 rt@-1 i@1 of@-1 to@-2\n",
      "2 the@-1 i@1 and@1 a@-1 in@-2\n",
      "3 and@-1 i@1 to@-1 in@1 the@1\n",
      "4 m@-2 so@-1 i@-2 not@-1 m@-1\n",
      "5 the@-2 i@1 and@1 a@-1 most@-1\n",
      "6 it@-1 like@1 that@-1 to@1 the@1\n",
      "7 the@-1 in@-2 i@1 of@1 for@-2\n",
      "8 i@-1 rt@-2 i@-2 t@1 a@1\n",
      "9 my@-1 and@1 to@-2 your@-1 i@1\n",
      "10 it@-2 do@-2 s@-1 but@-1 art@1\n",
      "11 an@-1 for@1 in@-2 i@1 of@-2\n",
      "12 a@-2 i@1 i@-2 in@-1 the@-1\n",
      "13 of@1 a@-1 the@-1 i@-2 the@-2\n",
      "14 to@1 i@-2 i@-1 to@-1 m@-1\n",
      "15 and@-2 i@-1 the@1 it@1 for@-1\n",
      "16 i@-2 i@-1 ve@-1 just@-1 in@1\n",
      "17 i@1 rt@-1 i@-2 to@-2 me@-1\n",
      "18 i@-2 m@-1 i@1 be@-1 the@1\n",
      "19 a@-1 i@1 and@1 have@-2 in@-2\n"
     ]
    }
   ],
   "source": [
    "# Let's print out the top features for each mean vector.\n",
    "# This is swamped by common terms\n",
    "for i in range(num_clusters):\n",
    "    print(i, ' '.join(names[np.argsort(kmeans.cluster_centers_[i])[::-1][:5]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.07242557,  0.99899023,  1.02986095,  1.02875689,  1.06510437,\n",
       "        1.02811637,  1.0538714 ,  1.14194475,  1.12543407,  1.06321653,\n",
       "        1.02268792,  1.07671936,  1.01856784,  1.08200541,  1.01522164,\n",
       "        1.04871631,  1.04493059,  1.04957583,  1.10209854,  1.05484867])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# .transform will compute the distance from each context to each cluster.\n",
    "distances = kmeans.transform(X)\n",
    "distances[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 see do make say sleep hear be tell know \n",
      "\n",
      "1 me like and you that work this a on \n",
      "\n",
      "2 same sky fact finale world freezer bathroom ground couch \n",
      "\n",
      "3 families bea sliding hug dang yup taylor cristal protein \n",
      "\n",
      "4 thankful excited grateful glad bummed pissed bored tired much \n",
      "\n",
      "5 person thing ensure bell vlogs time record cl of \n",
      "\n",
      "6 tooo feels seems bby takes looks bringing doc upstairs \n",
      "\n",
      "7 best sky finale world fact freezer bathroom ground res \n",
      "\n",
      "8 love swear am miss ll ve have just was \n",
      "\n",
      "9 life ability boyfriend countrymen son ear sister parents fault \n",
      "\n",
      "10 whatever wrong bigger cold bea within prolly magic fails \n",
      "\n",
      "11 argument hour interpreter attitude athlete audition era article outspoken \n",
      "\n",
      "12 from photo cs of sending mother bowl town boy \n",
      "\n",
      "13 type bunch ton member perk part plenty bottle piece \n",
      "\n",
      "14 trying listening supposed able listen tend want embarrassed tryin \n",
      "\n",
      "15 freaked petty votes keeping arguing cant making among loads \n",
      "\n",
      "16 a never crossing you finished curious received witnessed still \n",
      "\n",
      "17 but sometimes when ooh and things yeah if goodbye \n",
      "\n",
      "18 hungry buying tryin loyal leaving uploading sorry crossing shocked \n",
      "\n",
      "19 headache few little hustle burrito blessing sweetheart dickhead blogger \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Finally, we'll print the words that are closest\n",
    "# to the mean of each cluster.\n",
    "terms = np.array(list(contexts.keys()))\n",
    "for i in range(num_clusters):\n",
    "    print(i, ' '.join(terms[np.argsort(distances[:,i])[1:10]]), '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, interpreting these results requires a bit of investigation.\n",
    "\n",
    "Some patterns do emerge:\n",
    " - Cluster 7 has Polish\n",
    " - Cluster 8 are about communication acts\n",
    " - Cluster 13 is family\n",
    " - Cluster 14 seems to be abbreviations/[disfluencies](http://en.wikipedia.org/wiki/Speech_disfluency)\n",
    "\n",
    "As the number of tweets increases, we expect these clusters to become more coherent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**How does error decrease with number of cluster?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2857.5917288226187"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kmeans.score(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=5 score=2970.76\n",
      "k=10 score=2907.65\n",
      "k=20 score=2859.23\n",
      "k=50 score=2774.11\n",
      "k=100 score=2686.25\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEPCAYAAACHuClZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVXW9//HXe0A0RMVLyjmgIkdH0UQG1FQUxntqXvLC\noWNhSf08oidIy9RzVLqaZh000+yCtzRCvIAEaAqjJ0tEGIQAHQpFQcFLyEVSpPn8/lhr2JtxkGH2\n7MvMfj8fj/1gz3evtddnVsi7tb6XpYjAzMwsFxXFLsDMzNo+h4mZmeXMYWJmZjlzmJiZWc4cJmZm\nljOHiZmZ5SyvYSJpW0kzJNVKmi/pB2n7zpIel/SSpMck7ZS1z1WSFklaKOmkrPZ+kuZKqpM0Op91\nm5nZ1slrmETEB8CxEVEF9AGOkzQAuBJ4IiL2B6YBVwFIOhAYDPQGTgFuk6T0624HhkVEJVAp6eR8\n1m5mZs2X99tcEbEufbtteryVwJnA3Wn73cBZ6fszgLERsSEiXgEWAYdL6gbsEBEz0+3uydrHzMyK\nLO9hIqlCUi2wHKiJiAXAHhGxAiAilgO7p5t3B17L2n1Z2tYdWJrVvjRtMzOzEtAx3weIiHqgStKO\nwGOSqoHGa7h4TRczszYs72HSICJWS5oMHAqskLRHRKxIb2G9mW62DNgza7ceadvm2j9CkoPJzKwF\nIkJb3qpp+R7NtVvDSC1JnwBOBGqBicCX0s0uACak7ycCQyR1krQPsC/wXHorbJWkw9MO+aFZ+3xE\nRPgVwXXXXVf0Gkrl5XPhc+Fz8fGvXOX7yuRfgLvTAKgA7o2IJ9M+lHGSLgSWkIzgIiIWSBoHLAA+\nBIZH5re8BLgL2A6YHBFT81y7mZk1U17DJCLmAf2aaP87cMJm9rkeuL6J9lnAwa1do5mZ5c4z4Nux\n6urqYpdQMnwuMnwuMnwuWo9a415ZKZEU7e13MjPLN0lEqXbAm5lZeXCYmJlZzhwmZmaWM4eJmZnl\nzGFiZmY5c5iYmVnOHCZmZpYzh4mZmeWsYKsGl7r6+npqa2sBqKqqoqLCOWtm1lz+FxOorZ1P//4j\nGThwCQMHLqF//5HU1s4vdllmZm1G2S+nUl9fT//+I5kzZzSZbK2nb9+RzJo12lcoZlYWvJxKjmpr\na6mrq2bTU1FBXd2gjbe9zMzs45V9mJiZWe7KPkyqqqqorKwB6rNa66msfIqqqqriFGVm1saUfZ8J\nJB3wF154B3V1g6ivhw0bavjTn/6Tww47KE9VmpmVllz7TBwmqeyhwV//ehWXXlrB4MGtXZ2ZWWly\nmDTSGg/HGjcObrsNampapyYzs1LnMGmkNcJk/Xro2RMefxw+9anWqcvMrJR5aHAedOoEX/0q3H57\nsSsxM2sbfGWyGcuWwcEHw5IlsMMOrVCYmVkJ85VJnnTvDscdB/feW+xKzMxKX17DRFIPSdMkzZc0\nT9LX0vZDJP1J0guSJkjqkrbvLWmdpNnp67as7+onaa6kOkmj81l3g+HDk474dnbxZmbW6vJ9ZbIB\nuCwiDgKOBIZL6g38ErgiIg4BHgauyNrnrxHRL30Nz2q/HRgWEZVApaST81w7xx4L//wnPP10vo9k\nZta25TVMImJ5RMxJ368FXgS6A/tFxB/TzZ4Azsna7SP37CR1A3aIiJlp0z3AWXkrfONxM1cnZma2\neQXrM5HUE+gLPAvMl3RG+tFgoEfWpj3TW1zTJR2dtnUHlmZtszRty7uhQ5Mhwm+8UYijmZm1TQV5\nOFbaJzIeGBERayUNA26RdA0wEVifbvoGsFdErJTUD3hE0oFbe7xRo0ZtfF9dXU11dXWLa99pJ/j3\nf4df/hKuvbbFX2NmVlJqamqoacWZ2XkfGiypIzAJmBIRNzfx+X7AvRFxRBOfTQcuB14HpkdE77R9\nCDAoIi5uYp9WGRqcbe5cOOUUeOUV2GabVv1qM7OS0BaGBo8BFmQHiaRPpn9WAP8D/Dz9ebe0DUm9\ngH2BxRGxHFgl6XBJAoYCEwpQOwB9+kCvXjBxYqGOaGbWtuR7aPAA4HzgOEm1aV/IZ4DPS3oJWAAs\ni4i70l0GAnMlzQbGARdFxLvpZ5cAvwbqgEURMTWftTd2ySXws58V8ohmZm2HZ8A30/r1sPfeMG0a\n9O7d6l9vZlZUbeE2V7vQqRN85SseJmxm1hRfmWyF116DQw6BV1+FLl3ycggzs6LwlUkB7bknVFfD\nffcVuxIzs9LiMNlKw4cnHfHt7ILOzCwnDpOtdPzxSWf8M88UuxIzs9LhMNlKElx8sYcJm5llcwd8\nC7z7LuyzDyxcCN265fVQZmYF4Q74IujaFc47D371q2JXYmZWGnxl0kJz5sDpp8PLL0PHgiyXaWaW\nP74yKZK+fWGvveDRR4tdiZlZ8TlMcnDJJZ4Rb2YGvs2Vkw8+SNbreuop2H//ghzSzCwvfJuriLbd\nFoYNg9tvL3YlZmbF5SuTHL36KlRVJX9uv33BDmtm1qp8ZVJke+0FxxwD999f7ErMzIrHYdIKvF6X\nmZU7h0krOOEEWLcO/vznYldiZlYcDpNWUFGRrNflYcJmVq7cAd9KVq6Effap54EHatllF6iqqqKi\nwlltZm2DO+BLxCuvzKdDh5GcdtoSBg5cQv/+I6mtnV/ssszMCsJXJq2gvr6e/v1HMmfOaDL5XE/f\nviOZNWu0r1DMrOT5yqQE1NbWUldXzaans4K6ukHU1tYWqSozs8JxmJiZWc7yGiaSekiaJmm+pHmS\nvpa2HyLpT5JekDRBUpesfa6StEjSQkknZbX3kzRXUp2k0fmse2tVVVVRWVkD1Ge11rP33k9RVVVV\nnKLMzAoo31cmG4DLIuIg4EhguKTewC+BKyLiEOBh4AoASQcCg4HewCnAbZIa7uHdDgyLiEqgUtLJ\nea692SoqKhgz5iL69h1J584P0rnzg/ToMYK33rqIuXN98Wdm7V9BO+AlPQLcCjwQETunbT2AxyLi\nIElXAhERN6SfTQFGAUuAaRFxYNo+BBgUERc3cYyiDA2GpCO+oY+kqqqKhx6qYPhwGD8eBg4sSklm\nZs2Sawd8wZ4RKKkn0Bd4Fpgv6YyImEhyJdIj3aw7kD2PfFnatgFYmtW+NG0vKRUVFfTv33/jz+ee\nmzzi99xz4de/Tp7MaGbWHhUkTNI+kfHAiIhYK2kYcIuka4CJwPrWPN6oUaM2vq+urqa6uro1v36r\nnHACTJoEZ5wBN94IQ4cWrRQzs41qamqoqalpte/L+20uSR2BScCUiLi5ic/3A+6NiCOauM01FbiO\n5DbX9IjonbaX5G2uj7NwIZx8Mlx2GYwcWexqzMw21RbmmYwBFmQHiaRPpn9WAP8D/Dz9aCIwRFIn\nSfsA+wLPRcRyYJWkw9MO+aHAhALU3mp694Y//hF+/nO45hqvMGxm7Uter0wkDQCeBuYBkb6uBiqB\nS9KfH4qIq7P2uQoYBnxIclvs8bS9P3AXsB0wOSJGbOaYJXll0uCtt+CUU+Cww+DWW6FDh2JXZGaW\n+5WJl1MpgtWr4ayz4JOfhHvvhU6dil2RmZW7tnCbyxrZcUeYPBnWr09GeK1dW+yKzMxy4zApku22\ngwcegB49khFf77xT7IrMzFrOYVJEHTvCr36VTGgcOBCWLSt2RWZmLVOwSYvWNCmZf7LbbnD00fDY\nY1BZWeyqzMy2jsOkRFxxBey6K1RXw+9/D14f0szaEodJCRk2DHbeOZncOG5cEixmZm2B+0xKzNln\nw9ixMHgwTGhT0zLNrJz5yqQEHXdcMnT49NNh5Ur40peKXZGZ2cdzmJSoQw+F6dOTW15//3uyppeZ\nWanyDPgS99prcNJJ8LnPwfe/n4z+MjNrbV5OpZH2FiYAb78Np54KffvC7bd7PS8za31eTqUM7LYb\nPPkkLF4MQ4bABx8UuyIzs005TNqIHXZI5p/U18Npp8GaNcWuyMwsw2HShmy7bTL/pFcvOP745PaX\nmVkpcJi0MR06wB13JGFyzDFJB72ZWbF5aHAbJMH11yd9Kccck6zntf/+xa7KzMqZw6QNu/xy2GWX\nZNmVSZOgf/9iV2Rm5cph0sZ9+cvJel6nnAK/+x0ceyzU19dTW1sLQFVVFRUVvptpZvnleSbtRE1N\nsp7Xt741n9/85g7q6qoBqKysYcyYi6iqOqiY5ZlZifOkxUbKNUwAnn++niOPHMmGDaPJjK2op2/f\nkcyaNdpXKGa2WZ60aBtJtXTsWM2m/7NWUFc3aONtLzOzfHCYtDO++DCzYvA/Pe1IVVUVlZU1QH1W\naz3bbPMUXbr40Y1mlj95DRNJPSRNkzRf0jxJX0vbD5f0nKTa9M9D0/a9Ja2TNDt93Zb1Xf0kzZVU\nJ2l0PutuqyoqKhgz5iL69h1J584P0rnzgxx88Ai++MWLGDCggksugRUril2lmbVHee2Al9QN6BYR\ncyR1AZ4HPgfcBlwfEY9LOgW4IiKOlbQ38GhE9Gniu2YAl0bETEmTgZsj4rEmtivbDvgGTQ0Nfued\nZAn7u++GSy+Fb3wjWe/LzAxKvAM+IpZHxJz0/VrgReBfgTeArulmXYFlWbt95JdJQ2mHiJiZNt0D\nnJWvutu6iooK+vfvT//+/TeO4Np1V/jJT2DWLHj5ZdhvP7j1Vli/vsjFmlm7ULA+E0k9gb7ADOBK\n4MeSXgVuBK7K2rRneotruqSj07buwNKsbZambbaVevaEe+5JlmD5/e/hwAOTyY719Vvc1cxsswoy\nAz69xTUeGBERayU9DPxXRDwi6VxgDHAiyRXLXhGxUlI/4BFJB27t8UaNGrXxfXV1NdXV1a3wW7Qv\nhxwCU6bAtGlwxRXwox/BDTckC0iaWftXU1NDTU1Nq31f3ictSuoITAKmRMTNadvqiNgxa5tVEbFT\nE/tOBy4HXgemR0TvtH0IMCgiLm5in7LvM9la9fUwfjxcfTXsuy/88IfJUx3NrHyUdJ9JagywoCFI\nUoskDQKQdDxQl77fTVJF+r4XsC+wOCKWA6vSUWAChgITClB7WaioSJZiWbAATj8dPvMZ+MIX4JVX\nil2ZmbUV+R7NNQB4GpgHRPq6GniLZERXJ+B9YHhE1Eo6G/gOsJ5kssS1ETE5/a7+wF3AdsDkiBix\nmWP6yiRHa9bAj38MP/0pDB0K//3fyXL3ZtZ+eW2uRhwmrWfFCvjud2HsWLjsMhgxArbfvthVmVk+\ntIXbXNZG7bFHMnz4z3+GF16Aykr45S9hw4ZiV2ZmpWaLYSKpg6SvF6IYK0377ZcMH37kEfjtb+FT\nn4KHHwZfAJpZg2bd5pL0XEQcXoB6cubbXPkVkcxR+da3klteN94IRx+95f3MrLQVpM9E0v8C2wC/\nA95raI+I2S09cL44TAqjvh7uuw+uuQb69EmeSX+Qn79l1mYVKkymN9EcEXFcSw+cLw6Twnr/fbjt\ntmRuyumnw7e/DT16FLsqM9taHs3ViMOkON59N7nldccd8NWvJrfBdt652FWZWXMVZDSXpJ0k/UTS\n8+nrx5I+MmPdylfXrvCDH8DcufDOO7D//nDTTcmVi5m1f80dGjwGWAMMTl+rgTvzVZS1Xd27J8OH\nn3oK/vjHJFTuvhv++c9iV2Zm+dTcPpM5EdF3S22lwLe5SsszzyQLSa5Zk/SrnHIKqMUX0maWL4Wa\ntPiPrOXgG5ZJ+UdLD2rlY8CA5Arlu9+Fyy+HY4+FGTOKXZWZtbbmXpkcQvJAqoZ+kpXABRExN4+1\ntYivTErXhg1w110wahQceWTy5MfKymJXZWZQgCuTdBXf/SPiEKAP0CciqkoxSKy0dewIX/kK1NVB\n//7JVcvw4bB8ebErM7NcbTFMIqIeuCJ9vzoiVue9KmvXOneGK6+EF1+ET3wimex47bWw2n+zzNqs\n5vaZPCHpG5L2lLRLwyuvlVm7t+uuyVL3s2cnz06prEyWvfdz6c3anub2mbzcRHNERK/WLyk37jNp\nu154Aa66Cl56KelPGTw4eXCXmeVf3mfAp30mR0bEMy09SCE5TNq+6dOT4cT19clz6U84odgVmbV/\nhVqbqzYiqlp6kEJymLQPEfDAA8lz6f/t35I5KlVt4m+gWdtUqHkmT0o6J33+ulneScltroUL4cwz\n4dRT4fzz4eWmbriaWdE1N0wuAsYBH0haLWmNJI+9sbzbZptk+PCiRUkH/aGHwsiR8Pbbxa7MzLI1\nN0x2Ar4EfC8idgQOAk7MV1FmjXXpAtddBwsWJJMfDzgg6aR/770t72tm+dfcMPkZcATw+fTnNcCt\neanI7GM0PJf+2Wdh3rzkauWOO/xcerNia26YfDoiLgHeB4iIlUCnvFVltgX77gtjxybPpf/d75Ln\n0j/0kJ9Lb1YszQ2TDyV1AAJA0ieB+rxVZdZMhx0GTz4JN98M3/kOHHUU/N//Fbsqs/LT3DC5BXgY\n2F3S94E/Aj/Y0k6SekiaJmm+pHmSvpa2Hy7pOUm16Z+HZu1zlaRFkhZKOimrvZ+kuZLqJI3eqt/S\n2jUJTj45mUl/ySXwxS8mjxD+y1+KXZlZ+Wj2Y3slHQAcDwh4MiIWNmOfbkC3iJgjqQvwPPA54Dbg\n+oh4XNIpwBURcaykA4H7gMOAHsATwH4REZJmAJdGxExJk4GbI+KxJo7peSZl7oMPkufSX389fPaz\nyXPp99yz2FWZlbZCzTMhIl6MiJ9FxK3NCZJ0n+URMSd9vxZ4EfhX4A2ga7pZV2BZ+v4MYGxEbIiI\nV4BFwOFpKO0QETPT7e4Bzmpu7VZett0Wvv71ZDhxt27Qt2/yTPqVK4tdmVn7VbCVjyT1BPoCM4Ar\ngR9LehW4Ebgq3aw78FrWbsvStu7A0qz2pWmb2WbttFPmufQrVyYjv370Iz+X3iwfOhbiIOktrvHA\niIhYK+lh4L8i4hFJ55I8Y77V5q2MGjVq4/vq6mqqq6tb66utDereHX7xi+Rq5eqrk5WJv/OdpG+l\nQ4diV2dWHDU1NdTU1LTa9zW7z6TFB5A6ApOAKRFxc9q2Op382LDNuxHRVdKVJKsR35C2TwWuA5YA\n0yOid9o+BBgUERc3cTz3mdjH+tOfkoUkV61K1vw69VQ/l96sYH0mORgDLGgIktQiSYMAJB1P0jcC\nMBEYIqmTpH2AfYHnImI5sCodBSZgKDChALVbO9QwfPj734dvfhOqq5NJkGbWcnm9MpE0AHgamEcy\nRyWAq4G3SEZ0dSKZCDk8ImrTfa4ChgEfktwWezxt7w/cBWwHTI6IEZs5pq9MrNk2bIC7706Wajni\niCRg9t+/2FWZFV5BlqBvSxwm1hLr1sEtt8BNN8F55yWPEf6Xfyl2VWaF0xZuc5mVvIbn0r/0Emy/\nfbI8yzXX+Ln0Zs3lMDHLsuuuydXJ7Nnw6quw337JFcsHHxS7MrPS5jAxa8Leeyd9KX/4A0ydCr17\nw/33J48SNrOPcp+JWTNMn57Mot+wIXku/Yl+mo+1M+6Ab8RhYvkSAePHJxMf99knmaPSr1+xqzJr\nHe6ANysQKRnptWABfO5zcNpp8B//AYsXF7sys+JzmJhtpW22gYsvThaSPOCA5JkqI0bAW28VuzKz\n4nGYmLVQly7JfJSFC5OO+d694Xvf83PprTw5TMxytPvuyeKRzz6bPJBrv/3g5z+HDz9MPq+vr2fW\nrFnMmjWLeg8Hs3bKHfBmrez555ORX0uXwle+Mp/77ruDRYuqAaisrGHMmIuoqjqoqDWaNebRXI04\nTKwURMDUqfWcffZI3n9/NJmbAPX07TuSWbNGU1HhGwNWOjyay6wESbD77rVUVFSz6X9mFbz00iBq\na2uLU5hZnjhMzArsH/+A//zPpF/lzTeLXY1Z63CYmOVJVVUVlZU1QHanez19+jzFlVdW8fTTyaOE\nTzgB7rjDQ4utbXOfiVke1dbO58IL76CubhAA++1Xw513/ufGDvh//AOmTIEHHkj+POwwGDw4mRS5\n227FrNzKjTvgG3GYWKmpr6/f2EdSVVW12Y73deuSQBk3Lllc8tOfzgTLrrsWsmIrRw6TRhwm1h6s\nWweTJyfB8thjyVMgBw+Gs85ysFh+OEwacZhYe/Pee5lgefxxOPLITLDsskuxq7P2wmHSiMPE2rP3\n3oPf/z4Jlj/8AY46KhMsO+9c7OqsLXOYNOIwsXKxdm0mWJ54AgYMSILlzDMdLLb1HCaNOEysHK1Z\nkwmWJ5+Eo4/OBEvXrsWuztoCh0kjDhMrd2vWwKRJSbBMmwbHHJMEyxlnOFhs8xwmjThMzDJWr84E\ny/TpMHBgJlh22qnY1VkpKekwkdQDuAfYg2Qa8C8i4qeSxgKV6WY7Aysjop+kvYGFwIvpZ89GxPD0\nu/oBdwHbAZMjYuRmjukwMWvC6tXw6KNJsNTUwKBBmWDZccdiV2fFVuph0g3oFhFzJHUBZgFnRsSL\nWdvcBLwbEd9Lw+TRiOjTxHfNAC6NiJmSJgM3R8RjTWznMDHbglWrMsHy1FNQXZ0Ey+mnO1jKVUmv\nGhwRyyNiTvp+LclVR/dGmw0Gfpv180d+mTSUdoiImWnTPcBZrV+xWXnYaSf4whdg4kRYsgTOOQfG\njoU990xm3N9/f9L3YtZcBVvoUVJPoC8wI6vtGGB5RPwta9OekmZLmi7p6LStO7A0a5ulfDSUzKwF\nunaFoUOTK5UlS5I5K/ffDz16wNlnw29/62CxLetYiIOkt7jGAyPSK5QGn2fTq5LXgb0iYmXaR/KI\npAO39nijRo3a+L66uprq6uqWlG1Wdrp2hQsuSF4rV8KECfCb3yRL5p9wApx3Hnz2s9ClS7ErtVzV\n1NRQU1PTat+X99FckjoCk4ApEXFzVnsHYBnQLyJe38y+04HLSUJmekT0TtuHAIMi4uIm9nGfiVkr\n+/vfk2B54AF45hk48cQkWE47zcHSXpR0n0lqDLAgO0hSJwILs4NE0m6SKtL3vYB9gcURsRxYJelw\nSQKGAhMKULuZkawB9uUvJ2uEvfwynHoq3HUXdO+ehMq4cclSL1a+8j2aawDwNDAPiPR1dURMlXQn\n8OeI+EXW9mcD3wHWkwwlvjYiJqef9WfTocEjNnNMX5mYFcg778AjjyRh8uyzcPLJyaiwU0+Fzp2L\nXZ1tjZIeGlwMDhOz4nj77UywzJgBn/lMEiynnOJgaQscJo04TMyK7+234eGHk2CZOXPTYPnEJ4pd\nnTXFYdKIw8SstLz1ViZYnn8+CZTBg5OAcbCUDodJIw4Ts9L15puZYJk1K+lbaQiW7bYrdnXlzWHS\niMPErG1480146KEkWGbPToYZDx6cdOI7WArPYdKIw8Ss7VmxIhMsc+ZkguWkkxwsheIwacRhYta2\nLV+eCZYXXkhm3DcEy7bbFru69sth0ojDxKz9eOONTLDMm5cJlhNPdLC0NodJIw4Ts/bp9dczwfKX\nvyTL5TcES6dOxa6u7XOYNOIwMWv/Xn8dHnwwCZb585Nn3Z93XrIYpYOlZRwmjThMzMrLsmUwfnyy\nCOXChZlgOf54B8vWcJg04jAxK19Ll2aC5cUXk2ezNATLNtsUu7rS5jBpxGFiZgCvvZYJlrq6TLAc\nd5yDpSkOk0YcJmbW2KuvZoJl0aLk0cTnnQfHHutgaeAwacRhYmYfZ8mSTLD87W+bBkvHgjx7tjQ5\nTBpxmJhZc73yShIs48YlD/06++xkuPGgQeUXLA6TRhwmZtYSL7+cCZYlSzLBMnBgeQSLw6QRh4mZ\n5Wrx4kywvPbapsHSoUOxq8sPh0kjDhMza02LFyf9K+PGJUOPzzknCZZjjmlfweIwacRhYmb58re/\nZYLl9dczwXL00W0/WBwmjThMzKwQ/vrXTLAsX54JlgED2mawOEwacZiYWaEtWpQJljff3DRYKiqK\nXV3zOEwacZiYWTHV1WWC5e23M8Fy1FGlHSwOk0YcJmZWKl56KRMs77wD556bBMuRR5ZesJR0mEjq\nAdwD7AHUA7+IiJ9KGgtUppvtDKyMiH7pPlcBFwIbgBER8Xja3g+4C9gOmBwRIzdzTIeJmZWcF1/M\nBMvKlZlgOeKI0giWUg+TbkC3iJgjqQswCzgzIl7M2uYm4N2I+J6k3sD9wGFAD+AJYL+ICEkzgEsj\nYqakycDNEfFYE8d0mJhZSVu4MBMsq1Yly7mcdx58+tPFC5ZcwySvZUfE8oiYk75fCywEujfabDBJ\ngACcCYyNiA0R8QqwCDg8DaUdImJmut09wFn5rN3MLF9694Zrr02eGDl1Kuy4IwwbBj17wuWXw7PP\nQlv7/8QFy0BJPYG+wIystmOA5RGxOG3qDryWtduytK07sDSrfSkfDSUzszbnoINg1ChYsACmTIEu\nXeDLX06C5RvfgBkz2kawFGTFmfQW13iSPpC1WR99Hvhtax9v1KhRG99XV1dTXV3d2ocwM2t1Bx0E\n3/52Ei7z5ye3wS64AN5/P3Mr7LDDQC2+GZVRU1NDTU1N7l+UyvtoLkkdgUnAlIi4Oau9A8mVR7+I\neD1tuxKIiLgh/XkqcB2wBJgeEb3T9iHAoIi4uInjuc/EzNqNiOR22LhxyWv9+kywHHpo6wQLlHif\nSWoMsCA7SFInAgsbgiQ1ERgiqZOkfYB9geciYjmwStLhkgQMBSYUoHYzs6KS4OCD4bvfTUaETZiQ\nPNv+/POhVy/41rdg1qzi3wrL92iuAcDTwDwg0tfVETFV0p3AnyPiF432uQoYBnzIpkOD+7Pp0OAR\nmzmmr0zMrN2LgLlzM1cs9fXJ1crgwVBVtfVXLCU9NLgYHCZmVm4iYM6czHDjiCRUBg+Gvn2bFywO\nk0YcJmZWziKgtjYTLFImWA455KPBUl9fT21tLYceeqjDJJvDxMwsEQGzZ2eCpUOHTLD06QNz5szn\nwgvvoK6umnXrznGYZHOYmJl9VETSUd8QLB071rNmzUhWrBhNMhbLt7k24TAxM/t4EXDvvbMYNmwJ\nGzacnbaW/tBgMzMrIVIyQbJTp9b7ToeJmVkZqqqqorKyhmRB99w5TMzMylBFRQVjxlxE374j6dz5\nwZy/z30mZmZlzEODN8NhYma29drC2lxmZtbOOUzMzCxnDhMzM8uZw8TMzHLmMDEzs5w5TMzMLGcO\nEzMzy5lkbnUkAAAHzUlEQVTDxMzMcuYwMTOznDlMzMwsZw4TMzPLmcPEzMxy5jAxM7Oc5TVMJPWQ\nNE3SfEnzJH0t67P/krQwbf9h2ra3pHWSZqev27K27ydprqQ6SaPzWbeZmW2dfF+ZbAAui4iDgCOB\nSyQdIKkaOB04OCIOBm7K2uevEdEvfQ3Par8dGBYRlUClpJPzXHubV1NTU+wSSobPRYbPRYbPRevJ\na5hExPKImJO+XwssBLoDFwM/jIgN6WdvZ+32kfX0JXUDdoiImWnTPcBZ+ay9PfB/KBk+Fxk+Fxk+\nF62nYH0mknoCfYEZQCUwUNKzkqZLOjRr057pLa7pko5O27oDS7O2WZq2mZlZCehYiINI6gKMB0ZE\nxFpJHYGdI+IISYcB44BewBvAXhGxUlI/4BFJBxaiRjMza7m8P7Y3DY5JwJSIuDltmwzcEBFPpT//\nFfh0RLzTaN/pwOXA68D0iOidtg8BBkXExU0cz8/sNTNrgVwe21uIK5MxwIKGIEk9AhwHPCWpEtgm\nIt6RtBvw94iol9QL2BdYHBHvSlol6XBgJjAUuKWpg+VyMszMrGXyGiaSBgDnA/Mk1QIBXA3cCYyR\nNA/4gCQcAAYC35G0HqgHLoqId9PPLgHuArYDJkfE1HzWbmZmzZf321xmZtb+tZsZ8JI+I+nFdFLj\nt4pdTyFtbnKopJ0lPS7pJUmPSdqp2LUWiqSKdFTgxPTnsjwXknaS9EA6QXi+pE+X8bm4Kj0HcyXd\nJ6lTuZwLSb+WtELS3Ky2zf7u6blalP69Oak5x2gXYSKpArgVOBk4CPi8pAOKW1VBNTk5FLgSeCIi\n9gemAVcVscZCGwEsyPq5XM/FzSS3hXsDhwAvUobnQtLewFeBqojoQ3KL//OUz7m4k+Tfx2xN/u7p\nCNrBQG/gFOA2SVvsi24XYQIcDiyKiCUR8SEwFjizyDUVzGYmh/YgOQd3p5vdTZlM9JTUAzgV+FVW\nc9mdC0k7AsdExJ0AEbEhIlZRhucCWA2sB7ZPR5h+AlhGmZyLiPgjsLJR8+Z+9zOAsenfl1eARST/\nxn6s9hIm3YHXsn4u20mNWZNDnwX2iIgVkAQOsHvxKiuo/wW+STLgo0E5not9gLcl3Zne8vuFpM6U\n4bmIiJXAj4FXSUJkVUQ8QRmeiyy7b+Z3b/zv6TKa8e9pewkT46OTQ9n0H1Oa+LndkXQasCK9Uvu4\nS/N2fy5IbuX0A34WEf2A90hubZTj34tewNeBvYF/JblCOZ8yPBcfI6ffvb2EyTJgr6yfe6RtZSO9\ndB8P3BsRE9LmFZL2SD/vBrxZrPoKaABwhqTFwG+B4yTdCywvw3OxFHgtIp5Pf36QJFzK8e/FocAz\nEfH3iPgn8DBwFOV5Lhps7ndfBuyZtV2z/j1tL2EyE9g3XcK+EzAEmFjkmgqtqcmhE4Evpe8vACY0\n3qm9iYirI2KviOhF8vdgWkR8EXiU8jsXK4DX0onBAMcD8ynDvxfAS8ARkrZLO5OPJxmgUU7nQmx6\ntb65330iMCQd7bYPyeTx57b45e1lnomkz5CMXKkAfh0RPyxySQWTTg59GphHcqnaMDn0OZJ1z/YE\nlgCDsyaBtnuSBgGXR8QZknahDM+FpENIBiJsAywGvgx0oDzPxTdJ/vH8J1ALfAXYgTI4F5LuB6qB\nXYEVwHUkK5E8QBO/u6SrgGHAhyS3zR/f4jHaS5iYmVnxtJfbXGZmVkQOEzMzy5nDxMzMcuYwMTOz\nnDlMzMwsZw4TMzPLmcPErEAkrWnhfmeW2SrY1gY5TMwKp6WTus4iebRCs0nq0MJjmbWIw8TKRrrc\nzoJ09dy/SJoqadv0s+mS+qXvd5X0cvr+AkkPpw8RWizpUkmXp6vw/klS1yaOs7ukhyTNkVQr6YhG\nnw+S9GjWzz+VNDR9/8P0AU5zJN0o6UiSJcFvTI+5j6RekqZIminpqYblUtLVgW+X9Cxwg6SB6fFn\nS5olafs8nVqz/D4D3qwE7Qv8e0T8P0m/A84B7m9iu+yriINIlvXvDPwN+EZE9JP0E2AocEujfW8B\naiLi7HQdqC5b+H4A0iVfzoqIA9Kfd4yI1UqeFvloRDyUtj8BXBQRf5N0OHA7yVpTAN0j4oh0u4nA\n8Ij4c7r0/PtbODdmLeYwsXLzckTMS9/PAno2Y5/pEbEOWCdpJTApbZ8HHNzE9scBXwSIZL2i5vaV\nrAL+IelXwO+zjrNRenVxFPBA1tPvtsna5IGs988A/yvpPuChiCirlbStsHyby8rNB1nv/0nm/1Bt\nIPPfw3Yfs09k/VxP0/+HbEt9I9nH2ni8dGn0w0keJfBZYGoT+1YAKyOiX0RUpa9PZX3+3sYiIm4g\nWazvE8AzWasHm7U6h4mVm809MOsVkmdeAJyX4zGeBIYDSKqQtEOjYy8BDpS0Tdrncny6bWega0RM\nBS4D+qTbrwF2BIiINcDLks7d+AtJDdttQlKviJgfETeSPKbBI8IsbxwmVm42d9VwE3CxpFnALi3Y\nP9tI4FhJc4HngQOz942IpSTLnv8FGAvMTj/fEZgk6QWSRwp8PW0fC3wz7UTfBzgfGJZ20v+FpIO+\nqdpGSponaQ7J88+nNKN2sxbxEvRmZpYzX5mYmVnOHCZmZpYzh4mZmeXMYWJmZjlzmJiZWc4cJmZm\nljOHiZmZ5cxhYmZmOfv/XhQPqbb1GVsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10aa831d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "scores = []\n",
    "num_cluster_options = [5,10,20,50,100]\n",
    "\n",
    "for num_clusters in num_cluster_options:\n",
    "    kmeans = KMeans(num_clusters, n_init=10, max_iter=10)\n",
    "    kmeans.fit(X)\n",
    "    score = -1 * kmeans.score(X)\n",
    "    scores.append(score)\n",
    "    print('k=%d score=%g' % (num_clusters, score))\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(num_cluster_options, scores, 'bo-')\n",
    "plt.xlabel('num clusters')\n",
    "plt.ylabel('error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** How does error vary by initalization? **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "score=2866.97\n",
      "score=2860.65\n",
      "score=2868.09\n",
      "score=2866.2\n",
      "score=2882.51\n",
      "score=2863.72\n",
      "score=2859.34\n",
      "score=2867.88\n",
      "score=2883.09\n",
      "score=2869.8\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VfWd//HXO2ETqeLSVi0i2hpFQQhBrBvEIrbWPtSu\n4tLa4kxpdQS0nfm1jq2M7fxaH3Y62uk4dRBwr1bF2nHFLaj4cwHCIhTBOiqoqK0KIoLI/fz+OCdy\nwQAnITfnJvf9fDzyuDcn59z7CUve+a5HEYGZmdm2VOVdgJmZdQwODDMzy8SBYWZmmTgwzMwsEweG\nmZll4sAwM7NMShoYkvpIekjSQkkLJI1Ljw+T9JSkxvRxaHq8u6QbJc1Pr/lR0Ws9LGlxes0cSbuX\nsnYzM9tUlxK//gfA+RExV1IvYJak+4FLgAsjYrqk44FLgWOA0QARcYikHYBFkm6MiJfS1zs1IhpL\nXLOZmTWjpIERESuAFenz1ZIWA3sBrwK909N6Ay+nz1cAO0qqBnoC64BVRS/pLjQzs5yovVZ6S+oH\nNAADgF2BmUAAAo6IiGXpedcBnwd2AM6LiKvS4w8DuwPrgWkR8fN2KdzMzIB2+o097Y66FRgfEauB\nycC5EdEXOA+Ykp53BklQ7AHsB/wwDRqA0yJiIHA0cHR6rpmZtZOStzAkdQHuBO6JiMvTY6siYqei\nc96OiN6SrgBmRsQN6fHJ6XW3bvaaZwJ1ETGumffz5lhmZq0QEdra19ujhTEFWNQUFqmlkkYASBoJ\nLE2PLwZGpsd3BD4LLJZULWm39HhX4EvAM1t6w4goq4+LLroo9xpcU+eqyzW5prauKYuSDnpLOhI4\nHVggqZFkzOIC4LvAFZK6AWvTzwGuBCZLWkAytjE5Ip6R1BO4L22tVAMPAJNKWbuZmW2q1LOkZpL8\ngG/OYc2cvw74yNhERKwBhrZtdWZm1hKeptoO6uvr8y7hI1xTduVYl2vKxjVlk7WmdptW214kRWf7\nnszMSk0SUQaD3mZm1gk4MMzMLBMHhpmZZeLAMDOzTEq9W62ZmZWpQqFAY2P2DcAdGGZmFaixcSFj\nxlzJkiX1ma/xtFozswpTKBSoq5vA3LmXsXFkYtvTat3CMDOrEBs2wGuvwfTpjSxaVE9Lh7EdGGZm\n7aB4vKC2tpaqqradcxQBb70FL70Ey5YlH5s/f/VV2GUX2HXXJDxayl1SZmYltvl4QU1NA1OmjKW2\n9uDMr/Huuxt/+DcXBsuWQdeusPfeyUffvh993qcPdO/e+i4pB4aZWQk1/8O5wODBE5g9+zKqqqpY\nvx5eeaX5EGh6XLMm+YG/pTDYe2/YaaetVbKpjSE2AoA1a77mwDAzy9Ps2bMZPvxF1qz5yibHq6tv\no3//frz5Zh1vvAGf/OTWWwcf/zhoqz/OW664m2zo0KEe9DYzK0fV1TBhAowaBXvtBV1y+GlcVVVF\nXV1d9vNLWIuZWcWrra2lpqYBKBQdLXDQQTP4zndq6ds3n7BoDXdJmZmVWGPjQk499UqefXYEPXvC\n/vs3MHXq91o06F1qWbY3d2CYmbWDc84psH59I2PHlmZa7fZyYJiZlYE1a5KB68bGZBC7HPkGSmZm\nZeC22+Cww8o3LLJyYJiZldikSfB3f5d3FdvPXVJmZiW0eDEcc0yyAK9r17yr2TJ3SZmZ5WzyZDjz\nzPIOi6zcwjAzK5H3308Gux97DPbfP+9qts4tDDOzHN1xBxx0UPmHRVYODDOzErnqKvj7v8+7irbj\nLikzsxJ44QUYOhSWL4cePfKuZtvcJWVmlpPJk+H00ztGWGRV0sCQ1EfSQ5IWSlogaVx6fJikpyQ1\npo9D0+PdJd0oaX56zY+KXmtIenyJpMtKWbeZ2fb44AOYOrVzrL0oVuoWxgfA+RFxMHA4cLak/sAl\nwIURUQtcBFyanj8aICIOAYYCYyU1rY38L+CsiKgBaiR9vsS1m5m1yr33Jjc7Gjgw70raVkkDIyJW\nRMTc9PlqYDGwF/Aq0Ds9rTfwcvp8BbCjpGqgJ7AOWCVpD+BjEfF0et61wMmlrN3MrLUmTepcg91N\n2m0Xdkn9gMHAk8BSYKakXwECjgCIiPsknUESKDsA50XE25LqgOVFL7cc+FR71W5mltUrr8Cjj8IN\nN+RdSdtrl0FvSb2AW4HxaUtjMnBuRPQFzgOmpOedQRIUewD7AT9Mg8bMrEO4+mr4+tehV6+8K2l7\nJW9hSOpCEhbXRcQd6eHDImIUQETcKumq9PgRwO0RUQDekDSTZCzjMWDvopftw8ZurI+YOHHih8/r\n6+upr69vm2/GzGwrCoVk7cXNN+ddybY1NDTQ0NDQomtKvg5D0rXAXyPi/KJjs0kGw2dIGgn8MiIO\nTWdRDY6IMZJ2BJ4CvhERCyU9AYwDngbuAn4TEfc2835eh2FmuXjgAfjhD5P7XmirKxrKT+43UJJ0\nJPAIsACI9OMC4A3gCqAbsBY4OyIaJXUn6a4aRDK2MSUifp2+Vh1wNdADuDsixm/hPR0YZpaL0aPh\n6KPhnHPyrqTlcg+MPDgwzCwPf/0rfOYz8L//C7vsknc1LeeV3mZm7eTaa+HEEztmWGTlwDAz204R\nnW+jweY4MMzMttPjjyczpI46Ku9KSsuBYWa2nZru2d3RZka1lAe9zcy2w8qVsM8+sGQJfOITeVfT\neh70NjMrsRtvhFGjOnZYZOXAMDPbDp11o8HmODDMzFppzhx480049ti8K2kfDgwzs1aaNAnOOguq\nKuQnqQe9zcxa4d13Ye+9Yf785GZJHZ0Hvc3MSuSWW+DIIztHWGTlwDAza4WmtReVxIFhZtZCCxcm\nmwyecELelbQvB4aZWQtNngzf+Q50abebXJcHD3qbmbXAunXJuMUTT8CnP513NW3Hg95mZm3s9tth\n0KDOFRZZOTDMzFqgErYx3xJ3SZmZZfSXv8BnPwvLl0P37nlX07bcJWVm1oYmT4ZvfrPzhUVWbmGY\nmWWwfn2yjfmDD0L//nlX0/bcwjAzayN33w377dc5wyIrB4aZWQaVuLJ7c+6SMjPbhuXL4ZBDYNky\n2HHHvKspDXdJmZm1galTYfTozhsWWbmFYWa2FYVCMnYxbRoMGZJ3NaXjFoaZ2Xa6/37YbbfOHRZZ\nOTDMzLaikld2b85dUmZmW/D663DAAfDCC7DzznlXU1rukjIz2w7XXAMnn9z5wyKrkgaGpD6SHpK0\nUNICSePS48MkPSWpMX0cmh4/LT02J33cIOmQ9GsNkhYXfX33UtZuZpUtwt1Rmyv17T8+AM6PiLmS\negGzJN0PXAJcGBHTJR0PXAocExE3AjcCSBoA3B4R89PXCuDUiGgscc1mZjz6KFRXw+GH511J+Shp\nYETECmBF+ny1pMXAXsCrQO/0tN7Ay81cfipw02bH3IVmZu1i0qSkdaGt9upXlnYb9JbUD2gABgC7\nAjNJWg0CjoiIZZud/xxwYkQsSj9/GNgdWA9Mi4ifb+F9POhtZtvlrbdg333huedg9wrp/M4y6N0u\nd6RNu6NuBcanLY3bgXMj4o+SvgZMAUYVnT8MeLcpLFKnRcSrknYEpkk6IyKub+79Jk6c+OHz+vp6\n6uvr2/x7MrPO64Yb4Atf6Nxh0dDQQENDQ4uuKXkLQ1IX4E7gnoi4PD22KiJ2KjpnZUTsXPT5r4HX\nI+KXW3jNM4G6iBjXzNfcwjCzVouAwYPh17+GkSPzrqb9lMu02inAoqawSC2VNAJA0khgSdMXJAn4\nBkXjF5KqJe2WPu8KfAl4ph1qN7MKM2sWrF4NxxyTdyXlp6RdUpKOBE4HFkhqJBmzuAD4LnCFpG7A\n2vTzJsOBlyLihaJj3YH70tZKNfAAMKmUtZtZZZo0Cc46C6o8xeYjvNLbzCy1ejXsvTcsWgR77pl3\nNe2rXLqkzMw6hJtughEjKi8ssnJgmJmlrrrKd9XbGgeGmRmwYEFyZ70vfCHvSsqXA8PMjGSwe8wY\n6NIuq9M6Jg96m1nFe++9ZLB71izo1y/vavLhQW8zswymTYO6usoNi6wcGGZW8Zo2GrStc5eUmVW0\nJUvg6KNh2TLo1i3vavLjLikzs22YPBm+9a3KDous3MIws4r1/vvQty/MmJHcu7uSuYVhZrYVd96Z\nBEWlh0VWDgwzq1iTJnlld0u4S8rMKtKLL8KQIcnq7h12yLua/LlLysxsC6ZOhdNOc1i0hFsYZlZx\nNmxI7tn9P/8DgwblXU15cAvDzKwZ990He+zhsGgpB4aZVZyrrvLK7tZwl5SZVZQVK6B/f3jpJfjY\nx/Kupny4S8rMbDNXXw1f/arDojW887uZVYxCIemOuuGGvCvpmNzCMLOKMWNGMo122LC8K+mYHBhm\nVjGatjHXVnvqbUu2OegtqRoYFxH/3j4lbR8PeptZc/72N/j0p+H552HXXfOupvy0yaB3RGwATm2z\nqszMcnD99XDCCQ6L7ZFpWq2kfwe6AjcD7zYdj4g5pSutddzCMLPNRcDAgfDb30J9fd7VlKcsLYys\ns6QGp48XFx0L4HOtKczMrD098QSsWwcjRuRdSceWKTAi4phSF2JmVipXXZVsY+7B7u2TtUtqZ+Ai\nYHh6aAZwcUSsLGFtreIuKTMrtmoV7LMPLF4Mn/xk3tWUr7Zc6T0FeAf4RvqxCpiaoYA+kh6StFDS\nAknj0uPDJD0lqTF9HJoePy09Nid93CDpkPRrdZLmS1oi6bKMdZtZhfv97+Fzn3NYtIWsLYy5ETF4\nW8eauW4PYI+ImCupFzAL+DJwBfCLiJgu6Xjgnzbv9pI0ALg9IvZPP38S+IeIeFrS3cDlEXFfM+/p\nFoaZfejQQ+Hii+H44/OupLy1ZQvjPUlHFb3wkcB727ooIlZExNz0+WpgMbAX8CrQOz2tN/ByM5ef\nCtyUvt8ewMci4un0a9cCJ2es3cwqTKFQYPbs2dx442xee63AccflXVHnkHWW1PeAa9OxDIC3gDNb\n8kaS+pHMtnoSWArMlPQrQMARzVxyCnBi+vxTwPKiry1Pj5mZbaKxcSFjxlzJkiX1vP8+7L77Ncyf\nP5ba2oPzLq3D22ZgSKoCDoiIQZJ2AoiIVS15k7Q76lZgfESslnQ7cG5E/FHS10jGSEYVnT8MeDci\nFrXkfZpMnDjxw+f19fXUe+K1WUUoFAqMGXMlc+deRlMHyooVJzNmzARmz76MqirvhtSkoaGBhoaG\nFl2TdQxjVkQMbU1RkroAdwL3RMTl6bFVEbFT0TkrI2Lnos9/DbweEb9MP98DeDgi+qefjwZGRMT3\nm3k/j2GYVajZs2czfPiLrFnzlU2O9+x5G4880o+6urqcKit/bTmG8YCkH0raW9KuTR8Zr50CLGoK\ni9RSSSPSIkcCS4qKFslMrJuajkXECmBlOrtKwLeAOzK+v5l1ch98ANOnJ4Pba9bkXU3nlXUM45T0\n8ZyiYwHst7WL0sHx04EFkhrTay4AvgtcIakbsDb9vMlw4KWIeGGzlzsHuBroAdwdEfdmrN3MOqEN\nG+Cxx+Cmm+C222DffeHrX6/lueeuYdGik9n4+3CBmpoZ1NZ+Oc9yO4Usu9VWAYdHxMz2KWn7uEvK\nrPOKSLb5uPlm+MMf4BOfgNGj4RvfgP3SX183Dnon+4Dsv38DU6d+z4Pe25ClSyrrGEZjRNS2WWUl\n5MAw61wioLExaUn84Q/JDZBGj4ZTToEDD2z+mkKhQGNjIwC1tbUe7M6gLQPjV8D/A6aV+09jB4ZZ\n57BwYRISN92U3Fr1lFOSoBg40HtClUJbBsY7QE9gA8mYg4AonulULhwYZh3X0qVJd9NNN8HKlUlI\nnHIKDB3qkCi1ttzefGeSwet9I+JiSX2BPbe3QDOzF15Iuppuvhlefhm+/nW48ko4/HBwT1J5ydrC\n+C+gAHwuIvpL2gWYHhGHlrrAlnILw6z8vfIK3HJL0pJYuhS++tWkJTFiBFRX511dZWrLFsZhETEk\nnRpLRLyVTok1swqXdYD59deT6a833wzz5sFJJ8FPfwrHHgtdu7ZnxdZaWQNjvaRqknUUSPo4SYvD\nzCpY8b5NADU11zBlysZ9m956C6ZNS0LiqaeSHWPPOw8+/3no0SPHwq1VsnZJnU6yeG8IcA3wNeDC\niLiltOW1nLukrDMrp+mihUKBuroJm+zbBAUGDpzAD394GX/4QxWPPpq0IE45BU44AXbcMbdybRva\nbJZU+mIHAiNJZkg9GBF/3v4S254Dwzqrj/4237DJb/PtbUv7NsFtHHVUP8aOrePEE2GnsptLac1p\n08DoKBwY1hlt6bf5wYO3vgtroQDvvQfvvpvssfTuu233/J13ZrNu3YvApoGxww638eij3uivo2nL\nQW8zy9GsWY0sXlzPpvuFVrFgwQiOPrqRrl3rmv3hvnZtMlbQs2fSHbTjjtt+/vGPQ79+2z5vhx1q\nOfroa5g3b9N9mw44wPs2dVYODLMys3IlzJ+fzCRq+liwANat++i51dXJuoVBg5r/wd6zZynXMlQx\ndepYxoyZsMm+TVOmfM9bcXRS7pIyy0mhAM8/vzEUmkLijTdgwIAkBAYNgkMOgQEDChxzTMu7pNrn\n+yifgXhrPY9hmJWJ1auTVsLmrYZddtkYDE0fn/5084vXvAurlZIDw2w7tOY35wh46aVNg2HevGTL\ni4MO2rTVcMghsGvW25BtR01mWTgwzFopyxTW996DZ57ZdLxh/vxk++3NWw01NdDFI4ZWxhwYZq2w\npSms++47gbPOuowFC6qYNw9efDEJguJWw6BBySwjs47G02rNWqBQgGXL4M47G1m0qJ7Np7C++OII\nFi9u5IQT6rjgguTmPd28o5pVEAeGVZx33oElS+DZZ2Hx4uTx2WeTXVN32QX23DMJj8316AETJoDX\no1mlcmBY7koxkFsoJIPPm4fCs8/Cm28mXUkHHJB8nHRS8lhTAx/7GBQKtdTVXcPcuZsuSKup8YI0\nq2wew7Bcbe/+SO+8szEINm8t7LbbxlA48MCNz/fee9uL2TyF1SqNB72trGXdH2nDhuZbC4sXJ6ui\ni1sLTcFQUwO9em1/fZ7CapXCgWFlbUu7nXbrdhtnnNGPVavqePZZeO65pLVQ3Epoet6nj2/jadYW\nPEvKOqQNG5KVzl/5ShIMNTW+j4JZOXBgWG4KhVq6dLkG2HRweeDAGfzud192y8GszDgwrN09/TRc\nfDHMmVPF2LFjue++CTz3nHc7NSt3HsOwdvPkk/Av/5JsuvejH8FZZyVrGzy4bJY/D3pbWXj88aRF\nsWgR/PjHMGYMdO+ed1VmVixLYJT0VzlJfSQ9JGmhpAWSxqXHh0l6SlJj+ji06JpDJD0u6RlJ8yR1\nS48/LGlxes0cSbuXsnbbfo89BqNGwWmnJQPYzz0H3/++w8Ksoyr1GMYHwPkRMVdSL2CWpPuBS4AL\nI2K6pOOBS4FjJFUD1wGnR8QzknYB1he93qkR0Vjimm07zZiRdD298AL88z/DN7/pPZfMOoOSBkZE\nrABWpM9XS1oM7AW8CvROT+sNvJw+Pw6YFxHPpNe8tdlLunO7TEVAQ0MSFMuXJ0FxxhnQtWvelZlZ\nW2m3MQxJ/YAGYACwKzATCEDAERGxTNJ4oA74BLA7cHNEXJpe/3B6bD0wLSJ+voX38RhGO4qABx9M\nxihWrIALL0y6oHzvB7OOpWwW7qXdUbcC49OWxu3AuRHxR0lfA6YAo9J6jgSGAmuBByXNioiHgdMi\n4lVJOwLTJJ0REdc3934TJ0788Hl9fT319fUl/O4qUwRMn54Exd/+lgTF6NEOCrOOoqGhgYaGhhZd\nU/IWhqQuwJ3APRFxeXpsVUTsVHTO2xHRW9IpwBci4jvp8QuB9yLi3zZ7zTOBuogY18z7uYVRQhFw\n771J19OqVfCTn8A3vtH8PajNrOPIfZZUagqwqCksUksljQCQNBJYmh6/DxgoqUcaNCOARZKqJO2W\nnt8V+BLwTDvUbqkIuPNOOOww+Md/hPPOS9ZTnHqqw8KsUpS0hSHpSOARYAHJeEUAFwBvAFcA3Ui6\nns5umv0k6bT0nAJwd0T8SFLP9HW6ANXAAySzrz5SvFsYbSsC/vSnpOtp/Xr46U+TKbJeW2fWuXjh\nnrVaoQB33JEEBSRBcdJJDgqzzqpsBr2t4ygUYNo0+NnPkgHsiRPhxBNBW/1nZGaVwIFRYba0b9OG\nDXDbbUlQ9OgB//qvcMIJDgoz28iBUUE+ejvUa5g0aSxLlx7Mz36W3M/6kkvg+OMdFGb2UR7DqBBb\nuh1q9+4TqK29jIkTqzjuOAeFWaUql2m1VgYaGxvTlkXxX3kVMIL/+I9GPv95h4WZbZ0Do4IUCh89\nVl3toDCzbBwYFWDdOrj77lrWr28gWd7SpEBNzQxqa2tzqszMOhIHRif30EMwaBDMmlXFHXeMZfDg\nCfTseRs9e97GoEHjmTJlrO9wZ2aZeNC7k3rtNfjBD+DRR+E3v0kW3YFvh2pmzfOgdwXasAGuuAIG\nDIBPfSq5LWpTWABUVVVRV1dHXV2dw8LMWsTrMDqR2bOTW6D26AEPP5yEhplZW/GvmJ3AypUwblyy\nMvvss5NbpDoszKytOTA6sAi46SY46CB47z1YuBC+/W1PkzWz0nCXVAe1dGnSmnj9dbjlFjjiiLwr\nMrPOzi2MDmbtWrjoIjj88GTPp9mzHRZm1j7cwuhApk+Hc85J1lXMnQt9+uRdkZlVEgdGB/DKK8kt\nUZ9+Gn77W/jiF/OuyMwqkbukytgHH8Dllyctiv33h2eecViYWX7cwihTTz6ZrKno3TtZrX3ggXlX\nZGaVzi2MMvPWW0lQnHxysrXHgw86LMysPDgwykQEXHddsqZCSrb0OP10r6kws/LhLqky8Oc/J2sq\nVq2CO+6AYcPyrsjM7KPcwsjRmjVwwQUwfDh85Svw1FMOCzMrX25h5OSuu+Dcc5OAmDcP9tor74rM\nzLbOgVFCzd17YtkyGD8+mSJ75ZUwalTORZqZZeQuqRJpbFxIXd0Ehg9/keHDX2TIkAmcd95Camth\n8GCYP99hYWYdi++4VwKFQoG6ugnMnXsZGzO5QK9eE5g9+zJqapzTZlZefMe9nMyZ08jixfVs+sdb\nRaEwgnfeacypKjOz7VPSwJDUR9JDkhZKWiBpXHp8mKSnJDWmj0OLrjlE0uOSnpE0T1K39PgQSfMl\nLZF0WSnrbo1XXknWUZx5ZrJ9x9q1eVdkZta2St3C+AA4PyIOBg4HzpbUH7gEuDAiaoGLgEsBJFUD\n1wHfjYgBQD2wPn2t/wLOiogaoEbS50tc+1a9/Tb88Y/JTKf+/WHgwGQNxWc/C48+WsugQQ1AoeiK\nAjU1M6itrc2pYjOz7VPSWVIRsQJYkT5fLWkxsBfwKtA7Pa038HL6/DhgXkQ8k17zFoCkPYCPRcTT\n6XnXAicD95Wy/mJr18Ljj8MDDyTbdSxalNyT4thj4frrk4Hs6uqms6uYOnUsY8ZMYMmSEQDsv38D\nU6Z8j6oq9wKaWcfUboPekvoBDcAAYFdgJhCAgCMiYpmk8UAd8Algd+DmiLhUUh3wi4g4Ln2to4B/\niogTm3mfNhn03rAB5sxJwuGBB5LNAAcMgJEjk5A4/HDo3n3rr9HctFozs3KUZdC7XdZhSOoF3AqM\nT1satwPnRsQfJX0NmAKMSus5EhgKrAUelDQLWFXqGiNgyZKNLYiGBthzzyQcxo2DESNg551b9ppV\nVVXU1dWVpF4zs/ZW8sCQ1IUkLK6LiDvSw4dFxCiAiLhV0lXp8eXAI0VdUXcDQ4AbgL2LXrYPG7ux\nPmLixIkfPq+vr6e+vr7Z8155ZWML4sEHoaoqaUF89avwn/+ZBIaZWWfU0NBAQ0NDi64peZeUpGuB\nv0bE+UXHZpMMhs+QNBL4ZUQcKqk38ABwFMmA+T3Av0XEvZKeAMYBTwN3Ab+JiHubeb/YsGFDs90/\nb7+dtBwefDD5WLECPve5JCRGjkxuUuTdYc2sEmXpkippYEg6EngEWEAyXhHABcAbwBVAN5Kup7Mj\nojG95rT0nAJwV0T8OD1eB1wN9ADujojxW3jPGDz4XKZMGUv//gdvcaB65MjNB6rNzCpX7oGRB0kB\nG+jVawIRlzFwYNWHLYjDD4cePfKu0Mys/JTNoHf7q2L9+hHcd18jI0Z40NnMrC102nme1dXQq1fe\nVZiZdR6dNDC8qtrMrK11ysAYNGg8U6aM9UI5M7M21CkHvbc0rdbMzJpXsdubOyzMzNqef7KamVkm\nDgwzM8vEgWFmZpk4MMzMLBMHhpmZZeLAMDOzTBwYZmaWiQPDzMwycWCYmVkmDgwzM8vEgWFmZpk4\nMMzMLBMHhpmZZeLAMDOzTBwYZmaWiQPDzMwycWCYmVkmDgwzM8vEgWFmZpk4MMzMLBMHhpmZZeLA\nMDOzTEoaGJL6SHpI0kJJCySNS48Pk/SUpMb0cWh6fB9JayTNST+uKHqthyUtTq+ZI2n3UtZuZmab\nKnUL4wPg/Ig4GDgcOFtSf+AS4MKIqAUuAi4tuua5iBiSfpy92eudGhG16df+WuLa20xDQ0PeJXyE\na8quHOtyTdm4pmyy1lTSwIiIFRExN32+GlgM7AW8CvROT+sNvFx0mbbykh2yC60j/wNpT+VYE5Rn\nXa4pG9eUTdaaupS2jI0k9QMGA08CS4GZkn5FEhBHFJ3aT9IcYCXwk4h4rOhrV0taD0yLiJ+3S+Fm\nZga002/sknoBtwLj05bGZODciOgLnAdMSU99FegbEUOAHwA3ptcCnBYRA4GjgaMlndEetZuZWUIR\nUdo3kLoAdwL3RMTl6bFVEbFT0TkrI2LnZq59GPhBRMzZ7PiZQF1EjGvmmtJ+Q2ZmnVREbG1IoF26\npKYAi5rCIrVU0oiImCFpJLAEIJ359GZEFCTtB3wGeF5SNdA7Iv4mqSvwJeD+5t5sW9+wmZm1TkkD\nQ9KRwOnAAkmNQAAXAN8FrpDUDVibfg4wHLhY0vtAARgbEW9L6gncl7ZWqoEHgEmlrN3MzDZV8i4p\nMzPrHDrkNNXmSPpCurBviaT/k3c9AJImS3pN0vy8a2mypcWUOdfUXdKT6aLMhZL+b941NZFUlS4U\n/VPetQD8EHs+AAAFh0lEQVRIekHSvKZFr3nXAyBpZ0m3SPpz+vd3WBnUVFO0yLdR0soy+bf+4/TP\naL6kG9JelrxrGp/+LNjmz4NO0cKQVEUyDjISeAV4GhgdEYtzrusoYDVwbUQckmctTSTtAewREXPT\nGWizgZPK4M+qZ0SsScerZpJMdpiZZ01pXecBdcBOEXFiGdTzPMmEj7fyrqWJpKuBGRExNe027hkR\nq3Iu60Ppz4flwGERsSzHOvYBHgYOjIj3Jd0M3BUR1+ZY08HA74FDSRZa3wN8LyKeb+78ztLCGAYs\njYgXI2I9cBNwUs41ka4hKZv/2NDsYso/A5/KtyqIiDXp0+4k/y5z/3OT1Af4InBV3rUUEWX0/1bS\nTsDRETEVICI+KKewSB0L/CXPsEitAt4HdmwKVpJfcPPUH3gyItZFxAbgEeArWzq5bP7hbadPAcX/\nGJZTBj8Ey91miylzlXb9NAIrgIaIWJR3TcC/A/9IMlmjXARwv6SnJf193sUA+wJ/lTQ17f75b0k7\n5F3UZk4h+S06V2mr8N+Al0h2t3g7Ih7ItyqeIVnXtks6ueiLwN5bOrmzBIa1UDOLKXMVEYV0b7E+\nwHBJI/KsR9IJwGtpa0xsfcua9nRkurD1i8A5abdnnroAQ4D/TOtaA/wo35I2SqfhnwjcUga17Eey\nUHkfki2Sekk6Lc+a0q7oS0iWKdwNNAIbtnR+ZwmMl4G+RZ/3YdP9qaxI2hy+FbguIu7Iu55iaXfG\nXcDQnEs5EjgxHTP4PXCMpNz6mptExKvp4xvA7STdsXlaDiyLiFnp57eSBEi5OB6Ynf555W0oMDMi\n3ky7f6ax6bZIuYiIqRExNCLqgbdJ18U1p7MExtPAZ9Lt0bsBo4GymNVCef122qS5xZS5kbS7pJ3T\n5zsAo4C5edYUERdERN+I2I/k39NDEfGtPGuS1LNpqxxJOwLHkXQp5CYiXgOWSapJD40EyqE7scmp\nlEF3VOpZ4LOSekgSyZ/Vn3OuCUkfTx/7Al8GbtzSue22+WApRcQGSf8ATCcJwckRUQ5/ETcC9cBu\nkl4CLmoaHMyxpmYXU0bEvTmWtSdwTfqfqIqk5fNgjvWUq08Ct6fb33QBboiI6TnXBDAOuCHt/nke\n+E7O9QBJwJIMeH93W+e2h4iYl7ZSZ5N0+zQC/51vVQDcJmlXYD1w9tYmLXSKabVmZlZ6naVLyszM\nSsyBYWZmmTgwzMwsEweGmZll4sAwM7NMHBhmZpaJA8OszEh6WFI5rZY2AxwYZmaWkQPDLIN0W447\n05vxzJf0dUk/kfRU+vnvis59WNKv0x1lF0k6VNI0Sc9K+ll6zj7pDYeuT8/5g6QezbzvKEmPS5ol\n6eZ09bJZLhwYZtl8AXg5ImrTm2HdC/xHRAxLP++Z7nDbZF1EHAr8DrgD+B4wEPi2pF3Scw4AfhsR\nBwHvAGcXv6Gk3YALgZERMZRkS4kflO5bNNs6B4ZZNguAUZJ+IemoiHgHGCnpCSW34D0GOLjo/D8V\nXbcgIl6PiPeBv7DxfgMvRcQT6fPrgc23Kv8scBAwM93361tsuiuzWbvqFJsPmpVaRCxNB6K/CPxM\n0kPAOcCQiHhF0kVAcZfSuvSxUPQcks0et/T/bvON3QRMj4jTt/sbMGsDbmGYZSBpT+C9iLgR+BXJ\nPR8CeDPdcvxrrXjZvpIOS5+fBjy62defAI6U9Om0hp6S9m/VN2DWBtzCMMtmIHCppALJfZm/D5xM\ncj+KV4Gnis7d2hbQxV97luSueVOBhSTjHR+eExF/lfRt4PeSuqfHLwSWbvd3Y9YK3t7cLAeS9gHu\njIiBeddilpW7pMzy49/WrENxC8PMzDJxC8PMzDJxYJiZWSYODDMzy8SBYWZmmTgwzMwsEweGmZll\n8v8BpFruLvNlxhgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a9def28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scores = []\n",
    "for i in range(10):\n",
    "    kmeans = KMeans(20, n_init=1, max_iter=100, init='random')\n",
    "    kmeans.fit(X)\n",
    "    score = -1 * kmeans.score(X)\n",
    "    scores.append(score)\n",
    "    print('score=%g' % (score))\n",
    "     \n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(range(10), sorted(scores), 'bo-')\n",
    "plt.xlabel('sample')\n",
    "plt.ylabel('error')\n",
    "plt.show()"
   ]
  }
 ],
 "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.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}