pierre-haessig/logit_norm_imshow.ipynb

## logit_norm_imshow.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image plot with logit normalization - example of a transition matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pierre Haessig --  July 13, 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Objective**:\n",
    "\n",
    "> plot an 2D array (image) with `plt.imshow` using a [logit](https://en.wikipedia.org/wiki/Logit) normalization  \n",
    "> to **better distinguish** values that are *close to zero* AND values that are *close to one*.\n",
    "\n",
    "A closely related topic: logit scale\n",
    "\n",
    "* logit scale for plotting (y-scale), in particular cumulative probabilies: http://nbviewer.ipython.org/gist/pierre-haessig/7e3e6a818edeb6819708\n",
    "* Matplotlib Pull Request [#3753](https://github.com/matplotlib/matplotlib/pull/3753) which implements this scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.html.widgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example of the transition matrix of a dummy [Markov chain](https://en.wikipedia.org/wiki/Markov_chain#Example) with four states "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.9  ,  0.05 ,  0.04 ,  0.01 ],\n",
       "       [ 0.01 ,  0.979,  0.01 ,  0.001],\n",
       "       [ 0.   ,  0.001,  0.998,  0.001],\n",
       "       [ 0.   ,  0.01 ,  0.3  ,  0.69 ]])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = np.array([\n",
    "        [0.9, 0.05, 0.04, 0.01],\n",
    "        [0.01, 0.979, 0.01, 0.001],\n",
    "        [0., 0.001, 0.998, 0.001],\n",
    "        [0., 0.01, 0.3, 0.69]\n",
    "    ])\n",
    "T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  1.,  1.,  1.])"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Property of a transition/stochastic matrix:\n",
    "# each row sums to 1\n",
    "T.sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot with imshow, choice of colormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imopts = dict(interpolation='nearest',\n",
    "              cmap = 'BuPu_r')\n",
    "def plot_T(exponent=1):\n",
    "    plt.imshow(T**exponent, **imopts)\n",
    "    plt.xticks(range(4))\n",
    "    plt.yticks(range(4))\n",
    "    plt.colorbar()\n",
    "    if exponent == 1:\n",
    "        plt.title('linear norm')\n",
    "    else:\n",
    "        plt.title('power law with exponent $\\gamma$ = {:}'.format(exponent))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternate choice of colormaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Sequential\n",
    "#imopts['cmap'] = 'BuPu'\n",
    "\n",
    "### Diverging\n",
    "#imopts['cmap'] = 'RdBu_r'\n",
    "#imopts['cmap'] = 'RdYlBu_r'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear normalization (i.e. no normalization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observation:\n",
    "    \n",
    "* Good separation between three groups: values that are close to 1, close to zero, and in between\n",
    "* Bad separation at both ends:\n",
    "    * Values greater than 0.9 are very similar\n",
    "    * Values smaller than 0.1 are very similar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAEKCAYAAABUnUOUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFbhJREFUeJzt3XuQnFWdxvHvAyFoNihsBrlGXYxAkFuyggaQ5eYatUBK\nXSUuZeFld90q1NrSUkTL0rW84OIuoIgKiIgreEEhCMqighIIiOSKJOGmEkAwE0AJQUxmfvtHd4ZO\n093vOz19OT3n+aTequ7pt0+fhOThvKfPe36KCMzMUrFNvztgZlbLoWRmSXEomVlSHEpmlhSHkpkl\nxaFkZklxKPWRpN9JOqb6+HRJ5/e7T2b9NqXfHcjc2CKxiPhMPztilgqPlDInqWt/ByT5f3o2bg6l\nREj6hKRLqo9fLGlU0tsl/V7SOkmn15wrSadJukfSsKTvSNqp5vXvSfqDpMcl/ULSfjWvfUPSeZKu\nkbQBOKpBX26Q9J+SFkn6s6RrJc2oef0ESb+R9Jik6yXtW/Pa7yR9SNIK4AlJL6n+Xk6RdL+k9ZLe\nI+kQSSuqbXyx43+gNrAcSulodL/P4cDewLHAxyXtU/35+4ATgCOB3YDHgHNr3nc1MAvYGVgC/G9d\nuwuAT0XEdOCmJv1ZAJwCvACYCnwQQNLewLerfRgCrgGuqhsVnQS8FtgRGKn+7NBqn04CzgZOB44B\nXga8RdKRTfphmXEopUMNfvbJiHg6IlYAy4GDqj9/D/CxiHgoIjYBnwTevOVSLCK+ERFP1rx2kKQd\natq9IiIWV899usHnBnBRRNwTEX8BvgscXH3trcCPIuJnETECnAk8Fzis5r3nRMSDdW1/KiL+GhHX\nAU8A346I4Yh4CLgRmFP2D8omN1/zp+3hmscbgenVxy8CfihptOb1zcAukv4IfBp4M5WR0pZzhqiE\nQQAPjPOzn6r57N2B+7e8EBEhaS2wR835axu090hde/XPp2OGR0qD6n5gfkTsVHNMi4g/AG+jcml3\nbEQ8H/i76nsajcTa8SCVUKw0KgmYWf35Fu1sPeHtKgxwKA2qrwCfkfRCAEk7Szqh+tp04GngUUl/\nA9QvNSgbTs3O+x7weknHSNoO+ADwF+Dm8fwGxvF5lhmHUjqCrUcLrUYOZwMLgf+T9GdgMZWJZIBv\nAr+nMnK5o/pafbtlRiUN3xMRa4CTgS8C64DXA8dHxOaSbU3kHMuAvMmbmaXEIyUzS4pDycyS4lAy\ns6Q4lMwsKRNePCnJM+VmfRIRE1pKMd5/vxP9vDI6sqL7gcee7EQzpXzhc5/mA6d9tGef97EXvKdn\nnwWwdGQ5c7Y9qPjEDhmJkeKTOmj5yAoO2vbAnn3eE5s39OyzAFaxmtnsW3xiB1zBwo60s3FTub8D\n07bbtiOfV8S3mZhlL62LHYeSWeZSW6o4cKE074hX9bsLXbWrdul3F7pql0n++xtiqN9dGLdRh9LE\nHHbE5N52Z7dtdu13F7pq120mdyjtPICh5Ms3M0uKL9/MLCmJZZJDySx7iQ2VHEpmmUsrkhxKZtlL\nbKDkUDLLXWp7qjmUzDKXViQ5lMyyl9hAyaFklrtIbKzkUDLLXVqZ5FAyy11imeSdJ81yF1HuaETS\nfEmrJd0t6cMNXt9J0g8lLZd0q6SXFfXHoWSWuYgoddSTtC3wJWA+sB+wQNLsutNOB5ZExEHA26nU\nLGzJoWSWuSh5NHAocE9E/C4iNgGXAW+oO2c2cD2MFTJ9saSdW/XHoWSWuQlcvu0BrK15/kD1Z7WW\nA28EkHQo8CJgz1b9cSiZZa/tsVKZOfLPATtKWgqcCiwFWm4K7m/fzDLXbBL75kW/ZPGiX7Z664PA\nzJrnM6mMlmrajieAd255Lum3wH2tGlXRfS+S5gNnAdsCF0TEGXWvRy+rmfRar6uZ9Fqvq5n0Wq+r\nmfTSFSzsSImlBx7bWOrcPXeattXnSZoCrAGOBR4CfgUsiIhVNec8H3gqIv4q6V+AwyPilFaf03Kk\nVDO7fhyVVLxN0sLaDzWzAdfmfSYRsVnSqcC1VAYtF0bEKkn/Vn39q1S+lftGtb7cHcC7itotunwb\nm10HkLRldt2hZDZJTKRwQET8GPhx3c++WvN4MbDPeNosmuguM7tuZgNsAksCuqJopJTaCnQz67TE\ntgkoCqXC2XWolNLeYt4Rr5r0ZZDM+mEdwwwz3PF204qk4lD6NfBSSS+mMrv+VmBB/UkfOO2jHe+Y\nmW1tZ4a2qiu3hrs60m5iA6XWodRsdr0nPTOznhi47XAbza6b2eQxOmihZGaTW2KZ5FAyy93AXb6Z\n2eSWWCY5lMxy58IBZpYUj5TMLCkOJTNLyshE7sjtAu88aZa5KPmrkRLVTIYk/UTSMkl3SDqlqD8O\nJbPMtbtHd8lqJqcCSyPiYOAo4AvVzeGaciiZZa7dEkuUq2byB+B51cfPA9ZHxOZW/fGcklnmJjDR\n3Wi/tVfUnXM+8HNJDwE7AG8patQjJbPMTWBOqUycnQ4si4jdgYOBcyXt0OoNHimZZa7ZSOnXtyzi\n9ltuavXWMvutHQZ8uvI5cW+1msk+VLZFasihZJa50SZLAuYeejhzDz187PnXzvmv+lPK7Le2mkrh\nkZsk7UIlkFqWWHIomWWu3TmlktVMPgNcJGk5lemiD0XEo63adSiZZW4i976VqGYyDBw/njYdSmaZ\n820mZpYU76dkZklJLJMcSma5SyyTHEpmuWu2JKBfHEpmmfOckpklJbFM6kwonbrTs4rmThrfevry\nfnehq07e/k397kJXjTDS7y4kz3t0m1lSJuVIycwGl0PJzJLiiW4zS8qoQ8nMUpJYJjmUzHKXWih5\nO1yzzI0SpY5GSpRY+qCkpdVjpaTNknZs1R+HklnmulliKSLOjIg5ETEH+AhwQ0Q83qo/DiWzzLUb\nSpQrsVTrbcClRf1xKJllbjSi1NFAoxJLezQ6UdI04DVA4S0Snug2y1yzdUorb1/MyttvafnWcXzM\n8cCioks3cCiZZa/Zt2/7z53H/nPnjT2/7IKz6k8pU2Jpi5MocekGvnwzy16UPBoYK7EkaSqVEksL\n60+S9HzgSODKMv3xSMksc+3eZlKyxBLAicC1EfFUmXYdSmaZm8jiyaISS9XnFwMXl23ToWSWucQW\ndDuUzHLnPbrNLCneusTMkpJWJDmUzLKX2EDJoWSWO1++mVlS0ookh5JZ9lIbKRXeZiLp65IekbSy\nFx0ys94aHS139EqZe98uorKJk5lNQlHyV68UhlJE3Ag81oO+mFkfTGCTt67wnJJZ5hKbUupMKK1i\n9djjIYbYmaFONGtmNYZZz3rWd7zd1Ca6OxJKs9m3E82YWQtDzGCIGWPP7+LujrQ7kUiSNB84i8rW\nJRdExBkNzjkK+B9gO2A4Io5q1aYv38wyF23ekFtTzeQ4KrtQ3iZpYUSsqjlnR+Bc4DUR8YCkwsuo\nMksCLgVuBvaWtFbSO9r6HZhZkmI0Sh0NlKlm8jbg8oh4ACAihov6UzhSiogFhb8rMxtc7c8pNapm\n8oq6c14KbCfpemAH4OyIuKRVo758M8vcBOa5y7xzO2AucCwwDVgs6ZaIaDoh5lAyy1yzOaU1K37F\nmpW/avXWMtVM1lKZ3H4KeErSL4GDoPksvUPJLHdNhkr7HHAI+xxwyNjzH1365fpTxqqZAA9RqWZS\nP91zJfCl6qT49lQu7/67VXccSmaZa/fyrUw1k4hYLeknwApgFDg/Iu5s1a5DySx3E5lUKlfN5Ezg\nzLJtOpTMMhcjk3BFt5kNrkl5m4mZDa7EMsmhZJa9xFLJoWSWOV++mVlaerjVbRkOJbPMRS834C7B\noWSWucSu3hxKZrnznJKZpSWtTHIomeWu3Z0nu8WhZJY7X76ZWUoSyySHklnuUlsSUKZst5lNZqMl\njwYkzZe0WtLdkj7c4PWjJP1J0tLq8bGi7nikZJa5dpcElCmxVPWLiDihbLseKZllLqLc0UCZEksA\nGk9/HEpmuWs/lRqVWNqjvnXgMEnLJV0jab+i7nTk8m2EkU40k6QF25/Y7y501Xc3L+x3F7rqLVNK\nXzVkq9nl2z13LeHeu5a0fGuJ5pcAMyNio6TXAlcAe7d6g+eUzHLXZBJ71qy5zJo1d+z5dVd/vf6U\nwhJLEfFEzeMfS/qypL+NiEebdceXb2a5Gx0tdzzbWIklSVOplFjaaugtaRdJqj4+FFCrQAKPlMys\nzWVKZUosAW8G/l3SZmAjcFJRuw4ls9x1scRSRJwLnDueNh1KZrlL7D4Th5JZ7tK6y8ShZJY9j5TM\nLCkOJTNList2m1lSvPOkmSUlrUxyKJllz3NKZpYUX76ZWVLSyiSHkln2fPlmZknxkgAzS0piIyXv\np2SWu9EodzRQVM2k5rxDJG2W9Mai7jiUzHIXJY86NdVM5gP7AQskzW5y3hnATyhRRMChZJa79gsH\nlK1m8l7g+8C6Mt1xKJnlrv1ilIXVTCTtQSWozqv+qHACyxPdZrlrUrb7nvtXcO/ala3eWWaG/Czg\ntIiI6l7dhZdvDiWz3DWrZrLngcza88Cx59ct/nb9KYXVTIC/By6r1g4YAl4raVNENK3tVRhKkmYC\n3wReQCUZvxYR5xS9z8wGRPtLAsaqmQAPUalmsmDrpmOvLY8lXQRc1SqQoNxIaRPwHxGxTNJ04HZJ\n1zWoF25mg6jNUCpZzWTcCkMpIh4GHq4+3iBpFbA74FAymwwmsEd3UTWTup+/o0yb45pTqg7T5gC3\njud9ZpawxFZ0lw6l6qXb94H3R8SG2tfWcNfY4xnMYIgZHeugmVUMs571rO98w4MYSpK2Ay4HvhUR\nV9S/vg97d7pfZlZnqO5/+Hdxd2caHrQbcqtrCy4E7oyIs7rfJTPrqcRGSmVWdB8OnAwcLWlp9Zjf\n5X6ZWa+0v6K7K8p8+7YI345iNmlFYiMlr+g2y51DycyS0sNLszIcSma580jJzJIyaEsCzGyS80jJ\nzJKSViY5lMxyF4lVyPX6I7Pctb9Hd2E1E0lvkLS8uuj6dknHFHXHIyWz3LU5UKqpZnIclV0ob5O0\nsG6vtZ9GxJXV8w8AfgjMatWuR0pmuWu/7lthNZOIeLLm6XRguKg7HimZ5a79OaVG1UxeUX+SpBOB\nzwK7Af9Y1KhHSma5a3+kVCrNIuKKiJgNHA9cUnS+R0pmuWsSLfeuu5N7h1vuel2mmskzHxNxo6Qp\nkmZERNPd6hxKZplrtkvAXkOz2WvomSrcP139g/pTCquZSHoJcF+17tvc6ue13D7ToWSWu+5WM3kT\n8HZJm4ANwElF7TqUzHLXxWomEfF54PPjadOhZJY73/tmZknxLgFmlhSPlMwsKWllkkPJLHep7RLg\nUDLLnS/fzCwpaWWSQ6nItmzb7y501Tn/dGm/u9BV71x8br+70DVXzbu6Mw2NpFXOxKFkljvPKZlZ\nSlwh18zS4pGSmSXFIyUzS4pHSmaWlMRCydvhmmUuRkZLHY2UKLH0z9USSysk3STpwKL+eKRklrs2\n55RKlli6DzgyIv4kaT7wNeCVrdp1KJnlrv3Lt7ESSwCStpRYGguliFhcc/6twJ5FjfryzSxzEVHq\naKBRiaU9WnzUu4BrivrjkZJZ7pqMlO7bcDe/ffKeVu8sPcSSdDTwTuDwonMdSma5axJKe02bxV7T\nnqmw/fM/Xlt/SqkSS9XJ7fOB+RHxWFF3fPlmlruIcsezjZVYkjSVSomlhbUnSHoh8APg5IhoOeza\nwiMls8w1+7q/8H3lSix9HNgJOE8SwKaIOLRVuw4ls9xNYPFkiRJL7wbePZ42HUpmufO9b2aWlMRu\nM3EomWXOhQPMLC3h7XDNLCUeKZlZUlw4wMxSMlB7dEt6DvALYHtgKnBlRHykFx0zsx4ZpMu3iPiL\npKMjYqOkKcAiSUdExKIe9c/Mum2QQgkgIjZWH06lspT80a72yMx6K7HLt8IbciVtI2kZ8AhwfUTc\n2f1umVmvxGiUOnqlMJQiYjQiDqayY9yRko7qeq/MrHfa3yWgK0p/+1bdY/dq4OXADbWvreGusccz\nmMEQMzrVPzOrumPJYu5YckvH243RAVoSIGkI2BwRj0t6LvBq4JP15+3D3l3qnpltsf/ceew/d97Y\n8+9ceHZH2h2dwIruajGAs6jMN18QEWfUvb4vcBEwB/hoRHyhqM2ikdJuwMWStqFyqXdJRPysnc6b\nWZpGy+9qu5WS1UzWA+8FTizbbtGSgJXA3PF318wGRbQ/UipTzWQdsE7S68s26u1wzTI3SpQ6Ghhv\nNZNSfJuJWeaazSndv3kt94+sbfhaVVe+knMomWUuaBxKM6fswcwpzwx8btq0uP6UUtVMxsuhZJa5\nkRhp961j1UyAh6hUM1nQ5FyVbdShZJa50TYXRpapZiJpV+A24HnAqKT3A/tFxIZm7TqUzDI32uTy\nrYwS1UweZutLvEIOJbPMDdR+SmY2+U1kpNQNDiWzzLU7p9QtDiWzzI3Q9rdvXeFQMsvcRG7I7QaH\nklnmojsLs9vmUDLLnEdKZpaUZreZ9ItDySxz/vbNzJLidUpmlpQJ3JDbFQ4ls8z58s3MkpLaRLe3\nwzXL3GhEqaMRSfMlrZZ0t6QPNznnnOrryyXNKerPwIXSMOv73YWuWsdwv7vQVfeum9wFlu9Y8qzd\nGZNXbofuZ4+maqqZzAf2AxZIml13zuuAWRHxUuBfgfOK+jNwobR+kofS8CQPpfuGVxWfNMC6USyy\n2yYwUhqrZhIRm4At1UxqnQBcDBARtwI7StqlVX8GLpTMrLNGS/5qoEw1k0bn7NmqP57oNsvcBJYE\nlP3arn5/7pbv00R3nZOU1veJZhmJiNIb8jcy3n+/tZ8n6ZXAJyJifvX5R4DR2tLdkr4C3BARl1Wf\nrwb+ISIeafYZEx4pTfQPxcz6Z4L/fstUM1kInApcVg2xx1sFEvjyzczaVKaaSURcI+l1ku4BngTe\nUdTuhC/fzMw6aaC+fSuzUGtQSfq6pEckrex3XzpN0kxJ10v6jaQ7JL2v333qJEnPkXSrpGWS7pT0\n2X73aZANzEipulBrDXAclXLBtwELImJSLHyR9CpgA/DNiDig3/3ppGpBwl0jYpmk6cDtwImT5b8d\ngKRpEbFR0hRgEfDBiFjU734NokEaKZVZqDWwIuJG4LF+96MbIuLhiFhWfbwBWAXs3t9edVZEbKw+\nnEplfuXRPnZnoA1SKJVZqGWJq35TMwe4tb896SxJ20haBjwCXB8Rk/t+mi4apFAajOtMa6p66fZ9\n4P2taskPoogYjYiDqaxWPlLSUX3u0sAapFB6kK1rks+kMlqyASBpO+By4FsRcUW/+9MtEfEn4Grg\n5f3uy6AapFAaW6glaSqVhVoL+9wnK0GSgAuBOyPirH73p9MkDUnasfr4ucCrgaX97dXgGphQiojN\nVFaGXgvcCXxnkn17cylwM7C3pLWSCheZDZDDgZOBoyUtrR7z+92pDtoN+Hl1TulW4KqI+Fmf+zSw\nBmZJgJnlYWBGSmaWB4eSmSXFoWRmSXEomVlSHEpmlhSHkpklxaFkZklxKJlZUv4fEu8CXmCaTToA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa333d790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_T()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \"Square root\" type of linearization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Better sepation of values close to zero: plot of $T^\\gamma$ with an exponent $\\gamma < 1$\n",
    "\n",
    "(typically one gets $\\sqrt{T}$ for $\\gamma = 0.5$)\n",
    "\n",
    "Drawback: compression of the values close to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAEMCAYAAACCxKCJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGR5JREFUeJzt3Xm0HGWZx/HvDxIQRCQQRQlRQsIqCAEFwRiCG0FFOY7j\nEEVGRHEZ0KPjUXCBcRkZcEOOiCDCIBxhHFGIgiLKGpaAkEAghDWRKLIEghIiDOE+80dVbjqd7q66\nfbtvv33r98mpc7pvV731dPe9T9737er3UURgZpaK9XodgJlZLSclM0uKk5KZJcVJycyS4qRkZklx\nUjKzpDgpmVlSnJTMLClOSk1IWiLpTb2OA0DSHZKmt3g8mVjNhstJqbnIt56LiF0i4hoYTEBvrN+F\nRGIdaU1ej5E8/+aSfilpRR7LrIL9z5P0V0l/l/SApC+OVKz9wkkJkDSm1zEMQQDqdRAJ6fXrcSrw\nDPBS4P3AaZJ2brH/CcCkiNgUOBA4WtLM7ofZP3qWlPL/VY6RdKekJySdJWnD/LGdJF0laXk+dDko\n//nhkmbXtHGvpJ/V3F8q6dX57a0kXSjp0fx/pKMbnP9zkm4HnpLU8rXIY70v/x/uTkkHDzWmuvaG\n8lyWSHqTpHOBVwC/kvSUpM/WNDlV0m2SnpR0werXssF5G74ukiZLelzS1Jr9Hls9bGz1fuWPN3zP\nao7992bxtXqvWh1b8Hp0naQXAu8GvhwRKyPiOuBi4APNjomIOyPimZofrQIe7W6kfSYierIBS4Db\ngQnAOGAO8DVgDHAfcEx+e3/g78D2wCRgeX78VnkbD+b3twWeyG+vB9wCfClvYxJwP/DWuvPfmp9/\nwwbxLQbeWHP/PcDL8tvvBVYAW+bnLYypQfulnkt9LPVx1TyXG4GX5a/lQuCjDc7Z8nUBPgzcCWwE\nXAacVPR+5Y+NbfKebVcUX4mYWj63Rq/HMH4nfw0sb7LNbrD/VODpup99ptG+dfv8AHiaLCF9rFd/\ng6luvTtx9st0ZM39A/Nf7GnAX+v2/SlwfH77wfyX4RDg9PwXdgfgcOCifJ+9gT/VtXEscFbd+T9Y\nEF/TX3ZgHnBQ2ZiatFHqOIqT0mLgfTX3TwROa3C+Mq/LxcACYD4wtuj9ym+/oeA9axpfUUxFz63E\n+7Qp8LP8d+ub+c8mAtM78Dvc6Hl/BLiyxLECZgDLgL26+bfWb1uv51KW1tx+kKzHsFXdzwH+lP8c\n4GqyN3NKfvtJYD9gn/w+wCuBrSQtr2ljfeCaFudvSdJhwKeBbfIfbQKMH0JMjbR7XCMP19z+B2te\nr1plXpczyRLTRyLiubrjG71fUPyetYqvTExlnlszhwFHR8Qjkg6WtC+wdUT8rOjAElaQJb1aLwae\nKjowssx0laT/BWYBN3UgnlGh1xPdr6i7/VC+TZRUO3n5SuAv+e2ryYYHbwCuYs0f9n6s+UN+EFgc\nEeNqtk0j4h115y/1iZWkVwJnAP8GbB4R44A7WDPBWiamRto5rkzMzfZp+bpI2gQ4mSwxfUXSuLrj\nG71fUPyetbK0VUwlnlvR63FaRDwCEBEXAdOBPzfaUdJv8rmpRtslDQ65BxgjaUrNz3Yj+90oayzZ\nUM5yvUxKAj4haYKkzYEvAhcAc4GVwOckjZU0A3hH/his+UN+QUQ8RDa3MRPYnGxIBdn/Ok/lE9kb\nSVpf0i6SXtNmrC8k++VfBqwn6XBgl5rHy8TUSDvHPQJMLoi32adRRa/L94CbIuJI4BLgh3VtNnq/\noPg9a2Wo71X9c2v5ekTE83U/mhQR1zfZ98CIeFGT7e0N9n8a+AXwVUkbS5oGHASc2zBw6SWSDpH0\nwvx5HgD8M1nP1HK9TEpBNu/wO7KJzXuBr+dDhoPI5iweA74PfCAi7gGIiHvJusfX5vf/nh9/Xd4l\nJiIGyP4odgceyNs5g3W72uUCjVgIfBu4gWwosQtZAln9eGFMTdpt57gTgC/ln3J9plnTNOhBtHpd\nJL0LeCvw8Xz3zwB7aM11Nw3fr7zd/6PFe9YqvjxpDOW9qn9uZV6PWjeX2GcoPkH2wcCjwHlkE9d3\nrX5Q0qWSjsnvBvAxsp7a42Qf7HwgIjodU19Ti7+Z7p5YWgwcERFX9CQAG5LR8H7lva8tI6LRUMyG\nSNJZwNuBRyNi1yb7nEL2n9VKsg+WWo0cgN7PKZmNpFeRzd1ZZ5xNNt3QkKS3AVMiYjvgSOC0Mo06\nKVllRMQ5+TyQdUBEXEt2DVcz7wTOyfedC2wmacuidnt2SUBETOrVuW3o/H5ZGyaw9qUifwa2Jvtw\noin3lMysm+o/LS2cxB52T0lSJb+dbpaCiBjWl5GH+vc7xPP9hezq+dW2psS1ax0Zvv3tmfoLf7vn\nhK99lWO/fNyIne9Hx10+YucCuGzOeRww7dCRO+HKkXvvAC6bez4H7N1ydY+OWnXnshE7F8Dliy/m\nLZPeNSLnOubKIzrSzsrn6i/lamzjsesPtenZwFHABZJeBzy5+kLWVnr9NRMz67n2BjuSzif79sF4\nSUuB48muUCciTo+ISyW9TdJ9ZFetH16mXScls4pr91LFiCjs8kbEUUNtt++S0rTp+/U6hK6a/Ip1\nll4aVSZP2KV4pz627WY79DqEIUttUrjvktIb9hvdSWnKKE9KU7ZueOHvqDF53I69DmHIevWtjmb6\nLimZWWcllpOclMyqLrGc5KRkVnmJdZWclMwqLq2U5KRkVnmJdZSclMyqLrGc5KRkVnW+JMDMkpJY\nTnJSMqu6SGwA56RkVnVp5SQnJbOqSywnOSmZVV1qc0peDtes8qLkti5JMyUtknSvpM83eHycpF9K\nuk3SXEmvKorGScms4gai3FZP0vpkhUdnAjsDsyTtVLfbF4BbI2I34DCyKswtOSmZVVxEua2BvYD7\nImJJXtn6AqB+LeCdgCuz88TdwDaSXtIqHicls8pre/jWqITShLp9bgPeDSBpL+CVZAUEmvJEt1nF\nNZvovn7ONdww55qWh5Zo/r+A70maBywA5gEtKxUUJiVJM4GTgfWBMyPixBKBmFmfaJZZ9pk2nX2m\nTR+8/50Tv1G/S30JpYlkvaU1bUc8BXxo9X1Ji4EHWsXTcvhWciLLzPpZ+5NKfwS2k7SNpA2AfyEr\nqzRI0ovzx5D0EeDqiFjRKpyintLgRFbe6OqJrLtKPFUz6wPtXqYUEaskHQVcRjaS+nFE3CXpo/nj\np5N1Zv47L3p5B1BYrK4oKTWayNq7jfjNLFGNPu4vKyJ+A/ym7men19y+ARhSiZeipJTYtZ5m1nGJ\nXdJdlJQKJ7IgK6W92rTp+436MkhmvXD/8kU88OTdHW83rZRUnJQGJ7KAh8gmstapinnsl4/reGBm\ntrbJ43Zcq67cH5bMbrF3eYl1lFonpWYTWSMSmZmNiL5bebLRRJaZjR5ppSRf0W1WeQPD+fitC5yU\nzCqu74ZvZja6JZaTnJTMqs6FA8wsKe4pmVlSnJTMLCkDiWUlJyWzikstKXk5XLOKG8Ya3WWqmYyX\n9FtJ8yXdIemDRfE4KZlVXESU2uqVXATyKGBeROwOzAC+LanlCM1JyaziulzN5K/ApvntTYHHI2JV\nq3g8p2RWccO4TqnMIpA/Aq6Q9BDwIuC9RY06KZlVXLP5oj/eOIdbbryu5aElmv8CMD8iZkiaDFwu\nabe8oEBDTkpmFdcsKe259zT23Hva4P0zTvlm/S5lFoHcF/jP7Dxxf17NZAeytdoa8pySWcUNDESp\nrYHCaibAIuDNAJK2JEtILUssuadkVnHtzimVrGbyDeBsSbeRdYI+FxFPtGrXScms4oZz7WSJaibL\ngIOG0qaTklnFeT0lM0tKYjnJScms6hLLSU5KZlXn4ZuZJcWFA8wsKYl1lDqTlE478CedaCZJn/zd\nv/Y6hK465a3n9DqErnr21sW9DiF5XqPbzJIyKntKZta/nJTMLCn+9M3MkpLaGt1OSmYVl1hOclIy\nq7rUkpLXUzKruAGi1NZIiWomn5U0L98WSFolabNW8TgpmVVcu4UDylQziYhvRcTUiJgKHAtcFRFP\ntorHScms4rpczaTW+4Dzi+JxUjKruCj5r4FG1UwmNNpR0sbAAcCFRfF4otus4obxhdyhHHgQMKdo\n6AZOSmaV1+zTtwW33MAdt97Y6tAy1UxWO4QSQzdwUjKrvGbdnV323Idd9txn8P4FZ55cv8tgNRPg\nIbJqJrPqd5L0YmA62ZxSIScls4pr92smJauZABwMXBYR/yjTrpOSWcV1s5pJfv8coPQaOU5KZhWX\n2AXdTkpmVedVAswsKV6j28ySklZKclIyq7zERm9OSmZV5zklM0tKWinJScms8lLrKRWuEiDpLEmP\nSFowEgGZ2cgaxtIlXVFm6ZKzyRZxMrNRaCCi1DZSCpNSRFwLLB+BWMysB1LrKXlOyaziEptS6kxS\nunzxxYO3t91sByaP27ETzZpZjcWrlrBk1Z863m5qE90dSUpvmdRqWV4z64RJY7Zh0phtBu9f9ey1\nHWl3OClJ0kzgZLKlS86MiBMb7DMD+C4wFlgWETNatenhm1nFRZvffaupZvJmslUob5Y0OyLuqtln\nM+BU4ICI+LOk8UXtlrkk4HzgemB7SUslHd7WMzCzJMVAlNoaKFPN5H3AhRHxZ4CIWFYUT2FPKSLW\nWd7SzEaR9ueUGlUz2btun+2AsZKuBF4EfC8izm3VqIdvZhU3jHnuMkeOBfYA3gRsDNwg6caIuLfZ\nAU5KZhXXbE7p7ttv4u4FN7U6tEw1k6Vkk9v/AP4h6RpgN8BJycyaaNJV2mHX17LDrq8dvP/r839Q\nv0uZaiYXA9/PJ8U3JBvefadVOE5KZhXX7vCtTDWTiFgk6bfA7cAA8KOIWNiqXScls6obzqRSuWom\n3wK+VbZNJyWziovnR+EV3WbWv0bl10zMrH8llpOclMwqL7Gs5KRkVnEevplZWgZ6HcDanJTMKs49\nJTNLSrtLl3SLk5JZxbmnZGZpSSsnOSmZVZ2Hb2aWFg/fzCwlieUkJyWzykssK5Up221mo1g8H6W2\nRiTNlLRI0r2SPt/g8RmS/iZpXr59qSge95TMKq7dSwLKlFjKXR0R7yzbrntKZhUXUW5roEyJJQAN\nJR4nJbOqaz8rNSqxNKG+dWBfSbdJulTSzkXhdGT4dt2Vv+hEM0l6dvziXofQVZ97/Cu9DqGrTtri\n+F6HkLxmw7f77rmV+++5teWhJZq/FZgYESslHQhcBGzf6gDPKZlVXZNVAqZM2YMpU/YYvH/5JWfV\n71JYYikinqq5/RtJP5C0eUQ80SwcD9/Mqq794dtgiSVJG5CVWJpdu4OkLSUpv70XoFYJCdxTMrM2\nCweUKbEEvAf4uKRVwErgkKJ2nZTMqq6LJZYi4lTg1KG06aRkVnWJXdHtpGRWdV4O18yS4p6SmSXF\nScnMkuKy3WaWFK88aWZJSSsnOSmZVZ7nlMwsKR6+mVlS0spJTkpmlefhm5klxZcEmFlSEuspeT0l\ns6obiHJbA0XVTGr2e62kVZLeXRSOk5JZ1UXJrU5NNZOZwM7ALEk7NdnvROC3lCgi4KRkVnXtrzxZ\ntprJ0cDPgcfKhOOkZFZ1AyW3dRVWM5E0gSxRnZb/qHACyxPdZlXX/kR3mQNPBo6JiMjX6i4cvjkp\nmVVdk0sC7lt6O/cvXdDqyMJqJsCewAV57YDxwIGSnouI2TRRmJQkTQR+AryULDOeERGnFB1nZn2i\nSU9pyta7MmXrXQfvX37jT+t3GaxmAjxEVs1k1tpNx7arb0s6G/hVq4QE5XpKzwGfjoj5kjYBbpF0\neYN64WbWj9ocvpWsZjJkhUkpIh4GHs5vr5B0F7AV4KRkNhoMY43uomomdT8/vEybQ5pTyrtpU4G5\nQznOzBKW2BXdpZNSPnT7OfCpiFhR+9jd3DN4ewu2YDxbdCxAM8ssXrWEJav+1PmG+zEpSRoLXAic\nFxEX1T++A9t3Oi4zqzNpzDZMGrPN4P2rnr22Mw33W4ml/NqCHwMLI+Lk7odkZiNqIK2sVOaK7tcD\nhwL7S5qXbzO7HJeZjZT2r+juijKfvs3BX0cxG7WiH+eUzGwUc1Iys6SkNaXkpGRWee4pmVlSvEa3\nmSXFPSUzS0paOclJyazqIrEKub7+yKzq2l+ju7CaiaR3Sbotv+j6FklvLArHPSWzqmuzo1RTzeTN\nZKtQ3ixpdt1aa7+PiIvz/XcFfglMadWue0pmVdd+3bfCaiYR8XTN3U2AZUXhuKdkVnXtzyk1qmay\nd/1Okg4GTgBeDry1qFH3lMyqrv2eUqlsFhEXRcROwEHAuUX7u6dkVnVNUsv9jy3k/mUtV70uU81k\nzWkirpU0RtIWEfF4s/2clMwqrtkqAduO34ltx6+pwv37Rb+o36WwmomkycADed23PfLzNU1I4KRk\nZt2tZvJPwGGSngNWAIcUteukZFZ1XaxmEhEnAScNpU0nJbOq83ffzCwpXiXAzJLinpKZJSWtnOSk\nZFZ1qa0S4KRkVnUevplZUtLKSU5KRTbcY1KvQ+iq393S9FsBo8IRi4/tdQhdc/zmX+9MQx6+mVlS\nnk+rxpKTklnFuUKumaXFwzczS4p7SmaWFPeUzCwpiSUlL4drVnERUWprpESJpffnJZZul3SdpFcX\nxeOeklnVtXlJQMkSSw8A0yPib5JmAmcAr2vVrpOSWdW1P3wbLLEEIGl1iaXBpBQRN9TsPxfYuqhR\nD9/MKm4Yw7dGJZYmtDjVEcClRfG4p2RWdU16Sg+suJfFT9/X6sjSXSxJ+wMfAl5ftK+TklnVNUlK\n2248hW03XlNh+4pHL6vfpVSJpXxy+0fAzIhYXhSOh29mVRdRblvXYIklSRuQlViaXbuDpFcAvwAO\njYiW3a7V3FMyq7ho89O3kiWWjgPGAadJAnguIvZq1a6TklnVDePiyRIllj4MfHgobTopmVWdv/tm\nZklJ7GsmTkpmFefCAWaWlvDKk2aWEveUzCwpXqPbzFLSV2t0S3oBcDWwIbABcHFEjN6aNWZV1E/D\nt4h4RtL+EbFS0hhgjqRpETFnhOIzs27rp6QEEBEr85sbkF1K/kRXIzKzkZXY8K3wC7mS1pM0H3gE\nuDIiFnY/LDMbKTEQpbaRUpiUImIgInYnWzFuuqQZXY/KzEZO+6sEdEXpT9/yNXYvAV4DXFX72N3c\nM3h7C7ZgPFt0Kj4zy10/5xpumHNtx9uNgT66JEDSeGBVRDwpaSPgLcBX6vfbge27FJ6ZrbbvtOns\nO2364P3vnvSNjrQ7MIwruvNiACeTzTefGREn1j2+I3A2MBX4YkR8u6jNop7Sy4FzJK1HNtQ7NyL+\n0E7wZpamgfKr2q6lZDWTx4GjgYPLtlt0ScACYI+hh2tm/SLa7ymVqWbyGPCYpLeXbdTL4ZpV3ABR\namtgqNVMSvHXTMwqrtmc0oOrlvLg80sbPpbrykdyTkpmFRc0TkoTx0xg4pg1HZ/rnruhfpdS1UyG\nyknJrOKej+fbPXSwmgnwEFk1k1lN9lXZRp2UzCpuoM0LI8tUM5H0MuBmYFNgQNKngJ0jYkWzdp2U\nzCpuoMnwrYwS1UweZu0hXiEnJbOK66v1lMxs9BtOT6kbnJTMKq7dOaVucVIyq7jnafvTt65wUjKr\nuOF8IbcbnJTMKi66c2F225yUzCrOPSUzS0qzr5n0ipOSWcX50zczS4qvUzKzpAzjC7ld4aRkVnEe\nvplZUlKb6PZyuGYVNxBRamtE0kxJiyTdK+nzTfY5JX/8NklTi+Lpu6S0jMd7HUJX3b98Ua9D6KoF\nt6yzeuGocv2ca3odwpCVW6F73d5UTTWTmcDOwCxJO9Xt8zZgSkRsBxwJnFYUT98lpcdHeVJ64Mm7\nex1CV91x6429DqGrulEsstuG0VMarGYSEc8Bq6uZ1HoncA5ARMwFNpO0Zat4+i4pmVlnDZT810CZ\naiaN9tm6VTye6DaruGFcElD2Y7v69blbHqfhrjonKa3PE80qJCJKL8jfyFD/fmvPJ+l1wH9ExMz8\n/rHAQG3pbkk/BK6KiAvy+4uA/SLikWbnGHZPabgvipn1zjD/fstUM5kNHAVckCexJ1slJPDwzcza\nVKaaSURcKultku4DngYOL2p32MM3M7NO6qtP38pcqNWvJJ0l6RFJC3odS6dJmijpSkl3SrpD0id7\nHVMnSXqBpLmS5ktaKOmEXsfUz/qmp5RfqHU38GaycsE3A7Mi4q6eBtYhkt4ArAB+EhG79jqeTsoL\nEr4sIuZL2gS4BTh4tLx3AJI2joiVksYAc4DPRsScXsfVj/qpp1TmQq2+FRHXAst7HUc3RMTDETE/\nv70CuAvYqrdRdVZErMxvbkA2v/JED8Ppa/2UlMpcqGWJyz+pmQrM7W0knSVpPUnzgUeAKyNiYa9j\n6lf9lJT6Y5xpTeVDt58Dn2pVS74fRcRAROxOdrXydEkzehxS3+qnpPQX1q5JPpGst2R9QNJY4ELg\nvIi4qNfxdEtE/A24BHhNr2PpV/2UlAYv1JK0AdmFWrN7HJOVIEnAj4GFEXFyr+PpNEnjJW2W394I\neAswr7dR9a++SUoRsYrsytDLgIXA/4yyT2/OB64Htpe0VFLhRWZ95PXAocD+kubl28xeB9VBLweu\nyOeU5gK/iog/9DimvtU3lwSYWTX0TU/JzKrBScnMkuKkZGZJcVIys6Q4KZlZUpyUzCwpTkpmlhQn\nJTNLyv8Dy8IReQ8Ue+QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa320c650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_T(0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \"Squared\" type of linearization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $\\gamma > 1$, one gets a better separation of values close to one\n",
    "\n",
    "Drawback: compression of the values close to 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAEMCAYAAACCxKCJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF9NJREFUeJzt3XmwXGWZx/HvDwIIYhSSKRCIC0sYIipBBURBEJegooxD\nOcRRa3AtpwRLx0JcakpLHcd1wJFxYbEYLUALFMIg4sIaiAEhQCBkYRsiDEjCGsEZwn3mj3PuTdN0\n9zm3l9tv3/f3SZ2q7tvnvOfp7vTT73vO6fdRRGBmlorNhh2AmVkjJyUzS4qTkpklxUnJzJLipGRm\nSXFSMrOkOCmZWVKclMwsKU5KkyTpLkmHDTsOAEk3Szq4w+PJxGpWl5PS5EW5DF1E7B0RV8BEAnp9\n8yokEutUa/N6DCOOPST9RdKPhx3LqHBS6kDSjGHHMAkBaNhBJCSV1+Nk4Boy/XLoRnJJqfyGO0HS\nLZIelHS6pK3Kx/aSdJmkh8qhyxHl34+RtKihjTWSftZwf62kl5W3d5J0rqQ/SbpD0rEt9n+8pJuA\nxyR1fI3KWG+T9GgZ85GTjampvck8l7skHVZ+C78AuEDSY5I+1dDkfEk3SnpY0tnjr2WL/bZ8XSTt\nJmm9pPkN6z0wPmzs9H6Vj7d8zxq2/ad28XV6rzptW/F6TBlJRwMPAb8jjQQ5GiIiqQW4C7gJ2BnY\nDlgMfAmYAdwGnFDePhR4FJgLvBh4qNx+p7KNu8v7uwIPlrc3A64DPl+28WLgduBNTfu/vtz/Vi3i\nuxN4fcP9o4Ady9vvAjYAO5T7rYypRfu1nktzLM1xNTyX3wM7lq/lCuAjLfbZ8XUBPgjcAmwNXAx8\nver9Kh/bos17tkdVfDVi6vjcWr0ePfyf/C+K5NJqWdRmm5nAqvI9/ALw42F/tkZlSa6nRNHN/W5E\n3BMRDwFfARYCBwDPjoh/jYiNEXEpxX+WhRFxJ0WvZj5wMMUH515JewKvA64o234VMDsivly2cSdw\nKnB00/6/U+7/fyuDjTgnIu4rb/8MWAPsFxF31Iypub26z6WO8edyX/laXgDs02K9jq9LRJxKkVyu\noUi4n2vaR6v3C9q/Z++uEV/Ve1X3ubUkaaakn5W93G+Uf5vT6sRBRLwtIrZrs7y9zS6+BJwaEffi\nodukpHrMZG3D7bspvm12avo7wH+Xfwe4HDgE2L28/TDFh/jV5X2AFwI7SXqooY3NeeYHvXk/bUl6\nH/AJ4EXln7YFZk8ipla63a6V+xpuP8Gm16tRndflVOB84EMR8WTT9q3eL6h+zzrFVyemOs+tnfcB\nx0bE/ZKOlHQgsEv5xdITSfsAhwHzx//Ua5s5STUpvaDp9r3lMkeSouwfU/zHXVnevhx4O0Vy+ArF\nB/k9FN/W/16uczdwZ0TMrdh/rW82SS8Efgi8HlgSESFpGZv+E9aJqZVutqsTc7t1Or4ukrYFTqRI\nTF+U9POydzKu1fsF1e9ZJ2s7xdRC83Orej2+FxFPAUTEeZJOoE0vVNJFwGvbtHNFRLy16W+vo3jv\n7pYExRfV5pL2iohXVsSVvRSHbwL+UdLOkranGCqcDSwFHgeOl7SFpEOAt5WPQfFBPhR4VtllXgws\nALYHlpXrXEMxNDpe0taSNpe0t6Ru/6M8m+I//zpgM0nHAHs3PF4npla62e5+YLeKeNt9Y1e9LicB\n10TEh4ELge83tdnq/YLq96yTyb5Xzc+t4+sxnpAavDgirm6z7uER8Zw2S3NCguKLalfg5RRDyu9T\nvG5vbhePbZJiUgrgTODXFAc21wBfLocMRwCHAw8A3wXeGxGrASJiDfAYcGV5/9Fy+6vGv6UjYozi\nQ7EPcEfZzg8pDkpOPtCIFcC3gCUUQ4m9KRLI+OOVMbVpt5vtvgp8vjzL9cl2TdOiB9HpdZH0DuBN\nwEfL1T8J7Ctp/LhRy/erbPf/6PCedYqvTBqTea+an1ud16PRtTXWqSUinoiIP5XL/RQnP56IiPX9\n2sd0pg6fjaGQdCfwgYi4ZNixWLXp8H6Vva8dIuLCYcdiafaUzKbaS4DLhh2EFVI90G02ZSLijGHH\nYJskN3wzs7x5+GZmSel5+CbJXS2zIYmIni7MnOznt9f91dGXY0pnXXF7P5qp5ZzTT+Ko9398yvZ3\n5sHHTdm+AFaxmj2pe73g6PHz658L6M/JwsefbL5kq7Vttti8L/ur4gPdZtlLa7DjpGSWudTOdY1c\nUpo3f/9hhzBQs5g17BAGys8vPWNOSr2ZN/+AYYcwULNH8D/1ZPj5pSitrDRyScnM+svDNzNLSmI5\nyUnJLHuJdZWclMwyl1ZKclIyy15iHSUnJbPcpfajfCcls8yllZKclMyyl1hHyUnJLHeRWF/JScks\nd2nlJE/yZpa7qLm0ImmBpJWS1kj6dIvHt5P0C0k3Sloq6SVV8TgpmWUuot7STNLmFGWzFgDzgIWS\n9mpa7bPA9RHxcoqqxCdVxeOkZJa5iKi1tLAfcFtE3FXWZTwbeEfTOnsBl5b7WQW8SNJfdYrHScks\ncz0M33amKK8+7o/l3xrdCLwTQNJ+FGXbd+kUjw90m2Wu3SUBV195BVcvvqLjpjWa/1fgJEnLgOUU\nZec7zr/rpGSWvda55cCDDuLAgw6auP/tr32leZV7gDkN9+dQ9JY2tRzxGPD+8ftlReU7OkXj4ZtZ\n5ro90A38AdhD0oskbQn8HbCocQVJzy0fQ9KHgMsjYkOneCqTUtUpPzMbbd0eU4qIjcDHgIuBFcBP\nI+JWSR+R9JFytXnAckkrgTcDlaWIOg7fGk75vYGiq3atpEURcWtVw2Y2Inr4nUlEXARc1PS3HzTc\nXgLsOZk2q3pKdU75mdkIG4t6y1SpSkp1TvmZ2Qjr5YruQag6+5bYr2LMrO8SmyagKilVnvKDopT2\nuHnz95/2ZZDMhmEd61nP+r63m1ZKqk5KE6f8gHspTvktbF7pqPdXHlA3sx7NZtbT6sqtZk1f2k2s\no9Q5KUXERknjp/w2B07zmTez6WXkpsNtdcrPzKaPsVFLSmY2vSWWk5yUzHI3csM3M5veEstJTkpm\nuXPhADNLintKZpaU1JKS51Myy9xTY1FraaVGNZPZkn4l6QZJN0v6h6p4nJTMMhc1/zWrWc3kY8Cy\niNgHOAT4lqSOIzQnJbPM9TDzZJ2pjf4HmFnengmsLyeHa8vHlMwy18N1Sq2mNtq/aZ1TgEsk3Qs8\nB3hXVaPuKZllroeeUp1s9lnghojYCdgHOFnSczpt4J6SWebaXaf0h98v5g+/v6rTpnWmNjoQ+ApA\nRNxeVjPZk2IGkpaclMwy12709or9X8sr9n/txP0fnvSN5lXqTG20kmKO/6sk7UCRkDqWWHJSMsvc\nWJcTcLeb2mi8kklZQOBfgB9JupHicNHxEfFgp3adlMwy18vFkzWqmawDjphMm05KZpnzb9/MLCmp\n/czEScksc55PycySklhOclIyy11iOclJySx33V4SMChOSmaZ8zElM0tKYjmpP0npzIOP60czSTrt\n0XOGHcJAfWDmUcMOwYbM1ymZWVKmZU/JzEaXk5KZJcUHus0sKWNOSmaWksRykqfDNctdD9Ph1imx\n9ClJy8pluaSNkp7XKR4nJbPMjRG1lmZ1SixFxDcjYn5EzAc+A1wWEQ93isdJySxzAy6x1OjdwFlV\n8TgpmWWuh6TUqsTSzq1WlLQN8Gbg3Kp4fKDbLHPtzr4tv24Jy69b0mnTyRwiPwJYXDV0Aycls+y1\nu05p730PYO99D5i4f9YpJzavUqfE0rijqTF0Aw/fzLLXw/BtosSSpC0pSiwtal5J0nOBg4Hz68Tj\nnpJZ5rq9TKlmiSWAI4GLI+KJOu06KZllrpefmVSVWCrvnwGcUbdNJyWzzKV2RbeTklnmEstJTkpm\nufMc3WaWFE9dYmZJSSslOSmZZS+xjpKTklnuPHwzs6SklZKclMyyl1pPqfK3b5JOl3S/pOVTEZCZ\nTa2xsXrLVKnzg9wfUcwsZ2bTUNT8N1Uqk1JEXAk8NAWxmNkQ9DJH9yD4mJJZ5hI7pNSfpLSK1RO3\nZzGL2czqR7Nm1mAd61nP+r6328uBbkkLgBMppi45NSK+1mKdQ4B/A7YA1kXEIZ3a7EtS2pO5/WjG\nzDqY3fSFv5o1fWm325TUUM3kDRSzUF4raVFE3NqwzvOAk4E3R8QfJc2uatczT5plLsai1tJCnWom\n7wbOjYg/AkTEuqp46lwScBZwNTBX0lpJx1RtY2ajo4ekVKeayR7A9pIulfQHSe+tiqdy+BYRC6vW\nMbMR1v0xpTobbgHsCxwGbAMskfT7iGg79vTZN7PMtctJq266hlU3XdNp0zrVTNZSHNx+AnhC0hXA\ny6H9ATEnJbPMtRmaMXfvVzF371dN3L/gzJObV5moZgLcS1HNpHlkdT7w3fKg+FbA/sC3O8XjpGSW\nuy6Hb3WqmUTESkm/Am4CxoBTImJFp3adlMwy18vFkzWrmXwT+GbdNp2UzHKX2CXdTkpmmYunnJTM\nLCGpzafkpGSWucRykpOSWfYSy0pOSmaZ8/DNzNIyhVPd1uGkZJa5mMoJuGtwUjLLXGKjNycls9z5\nmJKZpSWtnOSkZJa7drMEDIuTklnuPHwzs5QklpNcOMAsdzE2VmtpRdICSSslrZH06RaPHyLpEUnL\nyuXzVfG4p2SWuy4vU6pTYql0eUS8vW677imZZS4iai0t1CmxBKDJxOOkZJa5iHpLC3VKLAVwoKQb\nJf1S0ryqeDx8M8tdm4xz26rruX319R23rNH69cCciHhc0uHAedC5pLaTUoUPzDxq2CEM1M82Lhp2\nCAP1rhm1D2Vkq90V3bvNnc9uc+dP3P/1hac1r1JZYikiHmu4fZGk/5C0fUQ82C4eD9/McjdWc3mm\niRJLkrakKLH0tG85STtIUnl7P0CdEhK4p2RmXc4SUKfEEnAU8FFJG4HHgaOr2nVSMstdDzOXVJVY\nioiTgWdUsezEScksd4ld0u2kZJY7JyUzS0paE086KZllzz0lM0uKk5KZJcVlu80sKZ550sySklZO\nclIyy56PKZlZUjx8M7OkpJWTnJTMsufhm5klJbFLAjyfklnuepgPt6qaScN6r5K0UdI7q8JxT8ks\nd10e6K5bzaRc72vAr6hRRMA9JbPcRc3lmepWMzkWOAd4oE44Tkpmuet++FZZzUTSzhSJ6nvje6sK\nx8M3s9x1P3VJnXHficAJERHlXN2VwzcnJbPctZmj+7a7b+L2tcs7bVlZzQR4BXB2WTtgNnC4pCcj\nom0ZHScls9y16SntvsvL2H2Xl03c/82SM5tXmahmAtxLUc1kYeMKEbHr+G1JPwIu6JSQoMYxJUlz\nJF0q6RZJN0s6rmobMxshXR5TioiNwHg1kxXAT8ermYxXNOlGnZ7Sk8AnIuIGSdsC10n6TfNpPzMb\nUT1c0V1VzaTp78fUabMyKUXEfcB95e0Nkm4FdgKclMymg1Geo7scO84Hlg4iGDMbglH97Vs5dDsH\n+HhEbGh8bBWrJ27PYhazmdW3AM2ssI71rGd9/xsexaQkaQvgXOAnEXFe8+N7MrffcZlZk9lNX/ir\nWdOfhhP7QW5lUioveDoNWBERJw4+JDObUon1lOr8zOQ1wHuAQyUtK5cFA47LzKbKWM1litQ5+7YY\n/0bObNqKxHpKvqLbLHdOSmaWlFG+TsnMpiH3lMwsKaN2SYCZTXPuKZlZUtLKSU5KZrmLxCrk+voj\ns9wNsMSSpHdIurG86Po6Sa+vCsc9JbPcddlRqlli6bcRcX65/kuBXwC7d2rXPSWz3I1FveWZKkss\nRcSfG+5uC6yrCsc9JbPcdX9MqVWJpf2bV5J0JPBV4PnAm6oadVIyy12bpHT7Ayu4/YGOE8zWymbl\ndEfnSToI+DGwZ6f1nZTMctcmtew2ex67zZ43cf+3K3/evEqdEkubdhNxpaQZkmZFRNvZ6nxMySxz\nEVFraWGixJKkLSlKLD2tfJKk3co52ZC0b7m/jtNnuqdklrsur+iOiI2SxkssbQ6cNl5iqXz8B8Df\nAu+T9CSwATi6ql0nJbPc9TBLQFWJpYj4OvD1ybTppGSWO//2zcyS4lkCzCwp7imZWVLSyklOSma5\nS22WACcls9x5+GZmSUkrJzkp5W7hjL8ZdggD9er93zbsEAbmgqUX9qehp9IqZ+KkZJY7H1Mys5S4\nQq6ZpcU9JTNLintKZpaUxHpKnk/JLHfdz9Fdp5rJ35fVTG6SdJWkl1WF456SWeaiy0sCalYzuQM4\nOCIekbQA+CFwQKd23VMyy133dd/qVDNZEhGPlHeXArtUheOkZJa77odvraqZ7NxhTx8AflkVjodv\nZpnr4Tql2htKOhR4P/CaqnWdlMxy1+Yg9h0b1nDnn2/rtGWtaiblwe1TgAUR8VBVOE5KZrlrk5R2\n3WZ3dt1mU4XtS/50cfMqE9VMgHspqpksbFxB0guAnwPviYiOGW6ck5JZ7gZbzeSfge2A75WVlp6M\niP06teukZJa5bi8JgFrVTD4IfHAybTopmeUusSu6nZTMcuffvplZUtxTMrOUuHCAmaUlPB2umaXE\nPSUzS4oLB5hZSkZqjm5JzwIuB7YCtgTOj4jPTEVgZjZFRmn4FhF/kXRoRDwuaQawWNJrI2LxFMVn\nZoM2SkkJICIeL29uSfH7lgcHGpGZTa3Ehm+Vk7xJ2kzSDcD9wKURsWLwYZnZVImxqLVMlcqkFBFj\nEbEPxTSWB0s6ZOBRmdnU6X463IGoffatnPj7QuCVwGWNj61i9cTtWcxiNrP6FZ+Zle54dDV3PLq6\nesVJirERuiRA0mxgY0Q8LGlr4I3AF5vX25O5AwrPzMbtOnMuu87c9Fm75J4L+9LuWA9XdJcVSk6k\nON58akR8renxvwZ+BMwHPhcR36pqs6qn9HzgDEmbUQz1fhwRv+smeDNL01j9qbafpmaJpfXAscCR\nddutuiRgObDv5MM1s1ER3feUJkosAUgaL7E0kZQi4gHgAUlvrduoSyyZZW6MqLW0MNkSS7X4ZyZm\nmWt3TOnujWu5+6m1LR8rDeSUnJOSWeaC1klpzoydmTNjU8fnqieXNK9Sq8TSZDkpmWXuqXiq200r\nSyw1UN1GnZTMMjc2wBJLknYErgVmAmOSPg7Mi4gN7dp1UjLL3Fib4VsdNUos3cfTh3iVnJTMMjdS\n8ymZ2fTXS09pEJyUzDLX7TGlQXFSMsvcU3R99m0gnJTMMtfLD3IHwUnJLHMxmAuzu+akZJY595TM\nLCntfmYyLE5KZpnz2TczS4qvUzKzpPTwg9yBcFIyy5yHb2aWlNQOdHs6XLPMjUXUWlqRtEDSSklr\nJH26zTrfKR+/UdL8qnhGLimtY/2wQxio6f78HmDdsEMYqEHUZRu0ejN0P7M31VDNZAEwD1goaa+m\ndd4C7B4RewAfBr5XFc/IJaX10/xDO92f3zonpeT00FOaqGYSEU8C49VMGr0dOAMgIpYCz5O0Q6d4\nRi4pmVl/jdX810Kdaiat1tmlUzw+0G2WuR4uCah72q55fu6O26nXWeckpXU+0SwjEVF7Qv5WJvv5\nbdyfpAOAL0TEgvL+Z4CxxtLdkr4PXBYRZ5f3VwKvi4j72+2j555Sry+KmQ1Pj5/fOtVMFgEfA84u\nk9jDnRISePhmZl2qU80kIn4p6S2SbgP+DBxT1W7Pwzczs34aqbNvdS7UGlWSTpd0v6Tlw46l3yTN\nkXSppFsk3SzpuGHH1E+SniVpqaQbJK2Q9NVhxzTKRqanVF6otQp4A0W54GuBhRFx61AD6xNJBwEb\ngP+MiJcOO55+KgsS7hgRN0jaFrgOOHK6vHcAkraJiMclzQAWA5+KiMXDjmsUjVJPqc6FWiMrIq4E\nHhp2HIMQEfdFxA3l7Q3ArcBOw42qvyLi8fLmlhTHVx4cYjgjbZSSUp0LtSxx5Zma+cDS4UbSX5I2\nk3QDcD9waUSsGHZMo2qUktJojDOtrXLodg7w8U615EdRRIxFxD4UVysfLOmQIYc0skYpKd3D02uS\nz6HoLdkIkLQFcC7wk4g4b9jxDEpEPAJcCLxy2LGMqlFKShMXaknakuJCrUVDjslqkCTgNGBFRJw4\n7Hj6TdJsSc8rb28NvBFYNtyoRtfIJKWI2EhxZejFwArgp9Ps7M1ZwNXAXElrJVVeZDZCXgO8BzhU\n0rJyWTDsoPro+cAl5TGlpcAFEfG7Icc0skbmkgAzy8PI9JTMLA9OSmaWFCclM0uKk5KZJcVJycyS\n4qRkZklxUjKzpDgpmVlS/h+JuwDoY7bWmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa2d7a090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_T(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAEMCAYAAACCxKCJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGBFJREFUeJzt3XuwHGWZx/HvD5LITQQSFwQiCiEYjECiglwMQVGDilBq\nuURdFFEsVlhLVxHU2l1rd3XxtmDposhlEUtYShTCcluUa7gEhIQEEkjCxQQxQAIoIbhLcp79Yzon\nwzDT3WfOzJl3Tv8+qa6ayXS//czMOc9533d63kcRgZlZKjbrdQBmZvWclMwsKU5KZpYUJyUzS4qT\nkpklxUnJzJLipGRmSXFSMrOkOCm1IOlRSe/sdRwAku6TNCPn8WRiNRsuJ6XWItt6LiKmRsTNMJiA\n3tG4C4nEOtJavB4jde5xks7NYvizpPmSZhUcc6OkFyQ9l21LRirefuGkBEga0+sYhiAA9TqIhPTy\n9RgDrABmRMS2wNeBSyTtlnNMAJ+LiFdm25SRCLSf9CwpZX9dTpV0v6SnJZ0n6RXZY1OyvyjPZEOX\nI7P/P07SnLo2lkm6pO7+Skn7ZLd3lnSppCclPSzp5CbnP0XSQuA5SbmvRRbr8uwv4v2Sjh5qTA3t\nDeW5PCrpnZIuBF4LXJH9lf1SXZPTJN0r6VlJF298LZuct+nrImkPSWskTavb76mNw8a89yt7vOl7\nVnfs37eKL++9yju24PXouohYFxHfiIgV2f0rgUeA6QWH+o9KnojoyQY8CiwEdgG2B+YC/0ztr89y\n4NTs9mHAn4HJwOuBZ7Ljd87aWJHd3x14Oru9GXA3tb9cY7LjHgLe3XD+e7Lzv6JJfI8A76i7/2Fg\np+z2R4C1wI7ZeQtjatJ+qefSGEtjXHXP5Q5gp+y1XAx8tsk5c18X4NPA/cCWwLXAt4ver+yxsS3e\nsz2L4isRU+5za/Z6DONn8r+BZ1psc0ocvyPwAjA5Z58bgCeBp7LX8NBe/Q6muvXuxLUfphPq7h+R\n/WAfAvyxYd9fAP+Y3V4BTAOOAX6S/cDuBRwHXJbtcwDw+4Y2TgPOazj/Jwvia/nDDswHjiwbU4s2\nSh1HcVJ6BPho3f3TgbOanK/M63I5sAhYAIwter+y228veM9axlcUU9FzK/E+bQtckv1sfSf7v4nU\nhlyd/HkeC/ym2evesN/+wNbZ/sdSS967j8TvXL9svZ5LWVl3ewW1HsPODf8P8Pvs/wFuAmYCk7Lb\nzwKHAgdm9wF2A3aW9ExdG5sDN+ecP5ekY4EvAK/L/msbYMIQYmqm3eOaWVV3+wU2vV71yrwu51BL\nTJ+JiBcbjm/2fkHxe5YXX5mYyjy3Vo4FTo6IJyQdLekgYNeIuKTowLKyof+FwF+Ak/L2jYg76+7+\nTNJs4L3ADzsVT7/r9UT3axtuP55tEyXVj7t3A/6Q3b6J2vDg7cCNbPrFPpRNv8grgEciYvu6bduI\neH/D+Ut9YpVNXJ4NfA7YISK2B+5j09xAmZiaaee4MjG32if3dZG0DXAGtcT0DUnbNxzf7P2C4vcs\nz8q8mEo8t6LX46yIeAIgIi4DZgCPNdtR0tV1n4o1ble2OEbAucCrgQ9FxIaCeKxAL5OSgL+VtIuk\nHYCvARcD84B1wCmSxkqaCbw/eww2/SJvERGPUxuXzwJ2oDakAriT2uT1KZK2lLS5pKmS3tJmrFtT\n++FfDWwm6Thgat3jZWJqpp3jngD2KIi31URq0etyJnBnRJwAXAn8uKHNZu8XFL9neYb6XjU+t9zX\no0mSeH1E3NZi3yNi06dijdv7WpziLOANwAci4n9bxQEg6VWS3iNpC0ljJH2M2h+ka/KOq5peJqWg\nNu/wP9QmNpcB/5INGY6kNmfxFLVu7d9ExFKAiFgGPAfckt3/c3b8rZEN2iNigNovxX7Aw1k7Z1Ob\nXxh6oBGLge8Bt1MbSkyllkA2Pl4YU4t22znuW8DXs0+5vtiqaZr0IPJeF0lHAe8GTsx2/yIwPRte\nbGzzZe9X1u7/kfOe5cWXJY2hvFeNz63M61HvrhL7lJL1oE8A9gVW1fWqZtftc5WkU7O7Y6l9mLNx\novtzwFERsbxTMY0Gyvmd6e6JpUeA4yPi+p4EYEMyGt6vrPe1Y9Q+urdE9XpOyWwkvZHa3J0lrNef\nvpmNmIi4oNcxWLGeDd/MzJrx8M3MkjLs4Zskd7XMeiQihvU9uqH+/g73fGV0ZE5p7pJVxTt1yLk/\n/A7Hn/TlETvf6VOOH7FzATzIUvZi8oiecyT5+XXOFXTmQ8R1L5a73nOrsZt35HxFPNFtVnlpDXac\nlMwqLrXPuvouKU3b/6Beh9BV4xnf6xC6ys8vPQNOSsMzff+Dex1CV03owx/qofDzS1FaWanvkpKZ\ndZaHb2aWlMRykpOSWeUl1lVyUjKruLRSkpOSWeUl1lFyUjKrutS+lO+kZFZxaaUkJyWzykuso+Sk\nZFZ1kVhfyUnJrOrSyklOSmZVl1hO8sqTZlUXUW5rRtIsSQ9IWibpK00e317SryXdK2mepDcWxeOk\nZFZxEVFqayRpc2o1/mYBewOzJU1p2O2rwD0RsS+1EupnFsXjpGRWcVFya2J/YHlEPJoVkb0YOKph\nnynADQAR8SDwOkmvzovHScms4oYxfNsFWFl3/7Hs/+rdC3wQQNL+wG7ArnnxeKLbrPKaZ5zbbrmZ\n2+bePPQDX+rfgDMlzQcWAfOB3EXBnZTMKq7VJPaBh8zgwENmDN7//unfbNzlD8DEuvsTqfWW6tqO\n54BPbbyflX9/OC+ewuFb0ey6mfW3Ycwp/Q7YU9LrJI0D/hqYU7+DpFdljyHpM8BNEbE2L57cnlLd\n7Prh1LLiXZLmRMSS3GdpZv2jze+ZRMR6SScB1wKbA+dGxBJJn80e/wm1T+X+M6svdx9QWLOsaPg2\nOLsOIGnj7LqTktkoMZzCARFxNXB1w//9pO727cBeQ2mzaPhWZnbdzPrYMIZvXVHUU0rtCnQz67TE\nlgkoSkqFs+tQK6W90bT9Dxr1ZZDMemE1a1jDmo63m1ZKKk5Kg7PrwOPUZtdnN+50/Elf7nhgZvZS\nExj/krpyS1nWkXYT6yjlJ6VWs+sjEpmZjYi+Ww632ey6mY0eA/2WlMxsdEssJzkpmVVd3w3fzGx0\nSywnOSmZVZ0LB5hZUtxTMrOkOCmZWVI2DOcbuV3g5XDNKi5K/mumRDWTCZKukbRA0n2SPlkUj5OS\nWcW1u0Z3yWomJwHzI2I/YCbwPUm5IzQnJbOKa7fEEuWqmfwR2Da7vS2wJiLW58XjOSWzihvGRHez\n9dYOaNjnp8D1kh4HXgl8pKhRJyWzims1X/S7O+byuztuzT+02FeBBRExU9IewHWS9s0KCjTlpGRW\nca16Sm8+4BDefMAhg/fPPvM7jbuUWW/tIOBfa+eJh7JqJntRWxapKc8pmVXcwECU2poorGYCPECt\n8AiSdqSWkHJLLLmnZFZx7c4plaxm8k3gfEn3UusEnRIRT+e166RkVnHD+e5biWomq4Ejh9Kmk5JZ\nxflrJmaWFK+nZGZJSSwnOSmZVV1iOclJyazqWnzc3zNOSmYV5zklM0tKYjmpM0np9CnHd6KZJF3w\n/K96HUJXfWLrD/Y6BOsxr9FtZkkZlT0lM+tfTkpmlhRPdJtZUgaclMwsJYnlJK+nZFZ17RYOgFLV\nTL4kaX62LZK0XtJ2efE4KZlV3ABRamtUpppJRHw3IqZFxDTgNODGiHg2Lx4nJbOKG0ZPqUw1k3of\nBS4qisdJyazihpGUmlUz2aXZjpK2At4DXFoUjye6zSpuGJ++DeXAI4G5RUM3cFIyq7xW1yktuvt2\nFt19R96hZaqZbHQMJYZu4KRkVnmtOkpTpx/I1OkHDt6/+JwzGncZrGYCPE6tmsnsxp0kvQqYQW1O\nqZCTklnFtTt4K1nNBOBo4NqIeKFMu05KZhU3nK+ZFFUzye5fAFxQtk0nJbOKS+2Kbicls4pLLCc5\nKZlVndfoNrOkeOkSM0tKWinJScms8hLrKDkpmVWdh29mlpS0UpKTklnlpdZTKly6RNJ5kp6QtGgk\nAjKzkTUwUG4bKWXWUzqf2spyZjYKRcl/I6UwKUXELcAzIxCLmfXAcNbo7gbPKZlVXGJTSp1JSg+y\ndPD2eMYzgfGdaNbM6qxmDWtY0/F2U5vo7khS2ovJnWjGzHJMaPiDv5RlHWl3OClJ0izgDGrrKZ0T\nEac32Wcm8O/AWGB1RMzMa9PDN7OKiza/kFtXYulwakvj3iVpTkQsqdtnO+BHwHsi4jFJE4raLXNJ\nwEXAbcBkSSslHdfWMzCzJMVAlNqaKFNi6aPApRHxGEBErC6Kp7CnFBEvW3PXzEaR9ueUmpVYOqBh\nnz2BsZJuAF4JnBkRF+Y16uGbWcW1ykkPLryTBxfemXtoiebHAtOBdwJbAbdLuiMiWk6IOSmZVVyr\nOaXJU9/K5KlvHbx/xS9+1LhLmRJLK6lNbr8AvCDpZmBfaD1L7wq5ZlXX/tWTgyWWJI2jVmJpTsM+\nlwOHSNo8q5J7ALA4Lxz3lMwqrt0ppTIlliLiAUnXAAuBAeCnEeGkZGY5ul9i6bvAd8u26aRkVnGx\nYRRe0W1m/WtUfs3EzPpXYjnJScms8hLLSk5KZhXn4ZuZpWUEl7otw0nJrOJiJBfgLsFJyaziEhu9\nOSmZVZ3nlMwsLWnlJCcls6prd+XJbnFSMqs6D9/MLCWJ5SSvp2RWdTEwUGprRtIsSQ9IWibpK00e\nnynpT5LmZ9vXi+JxT8ms6tq8TKlMNZPMTRHxgbLtuqdkVnERUWprokw1EwANJR4nJbOKa3813KbV\nTHZpbB44SNK9kq6StHdRPB6+mVVd+zPdZQ68B5gYEeskHQFcBvkltZ2UCnxi6w/2OoSuumR94zrv\no8tHxpSeyqisVld0L196Dw8tvSfv0MJqJhHxXN3tqyX9h6QdIuLpVo06KZlVXYuJ7kmTpjNp0vTB\n+9ddeV7jLoPVTIDHqVUzeUnxWkk7Ak9GREjaH1BeQgInJTNrc5WAMtVMgA8DJ0paD6wDjilq10nJ\nrOqGsXJJUTWTiPgR8LIqlnmclMyqLrFLup2UzKrOScnMkpLWwpNOSmaV556SmSXFScnMkuKy3WaW\nFK88aWZJSSsnOSmZVZ7nlMwsKR6+mVlS0spJTkpmlefhm5klxZcEmFlSEuspeY1us6obiHJbE0Ul\nlur2e6uk9ZIKl3J1UjKruii5NagrsTQL2BuYLWlKi/1OB66hRGUTJyWzqmu/nEnZEksnA78EnioT\njpOSWdUNlNxerrDEkqRdqCWqs7L/KpzA8kS3WdW1WKN7+YqFPLRyUd6RZWbIzwBOzQoHiBLDNycl\ns6prVc1k132YtOs+g/evu/0XjbsUllgC3gxcXMtHTACOkPRiRLSs7VWYlCRNBH4G/BW1zHh2RPyg\n6Dgz6xPtXxJQWGIpInbfeFvS+cAVeQkJyvWUXgS+EBELJG0D3C3puohYMrT4zSxJbSalkiWWhqww\nKUXEKmBVdnutpCXAzoCTktlo0MUSSw3/f1yZNoc0p5R106YB84ZynJklLLEruksnpWzo9kvg8xGx\ntv6xB1k6eHs845nA+I4FaGY1q1nDGtZ0vuF+TEqSxgKXAj+PiMsaH9+LyZ2Oy8waTGj4g7+UZZ1p\nuN++kJtdW3AusDgizuh+SGY2ohLrKZW5ovtg4OPAYZLmZ9usLsdlZiOl/Su6u6LMp29z8ddRzEat\nSKyn5Cu6zarOScnMkjKCQ7MynJTMqs49JTNLSr9dEmBmo5x7SmaWlLRykpOSWdVFYhVyff2RWdW1\nv0Z3YTUTSUdJuje76PpuSe8oCsc9JbOqa7OjVFfN5HBqq1DeJWlOw1prv4mIy7P93wT8GpiU1657\nSmZV137dt8JqJhHxfN3dbYDVReG4p2RWde3PKTWrZnJA406Sjga+BbwGeHdRo+4pmVVd+z2lUtks\nIi6LiCnAkcCFRfu7p2RWdS1Sy0NPLeah1bmrXpepZrLpNBG3SBojaXxEtFytzknJrOJarRKw+4Qp\n7D5hUxXu3zzwq8ZdCquZSNoDeDir+zY9O1/u8plOSmZV191qJh8CjpX0IrAWOKaoXScls6rrYjWT\niPg28O2htOmkZFZ1/u6bmSXFqwSYWVLcUzKzpKSVk5yUzKoutVUCnJTMqs7DNzNLSlo5yUmp6k7c\n8pO9DqGrDv3iib0OoWuu+P6VnWloQ1rlTJyUzKrOc0pmlhJXyDWztLinZGZJcU/JzJLinpKZJSWx\npOTlcM0qLjYMlNqaKVFi6WNZiaWFkm6VtE9RPO4pmVVdm3NKJUssPQzMiIg/SZoFnA28La9dJyWz\nqmt/+DZYYglA0sYSS4NJKSJur9t/HrBrUaMevplVXESU2ppoVmJpl5xTHQ9cVRSPe0pmVdeip/Tw\n2mU88vzyvCNLd7EkHQZ8Cji4aF8nJbOqa5GUdt9qErtvtanC9vVPXtu4S6kSS9nk9k+BWRHxTFE4\nHr6ZVV1Eue3lBkssSRpHrcTSnPodJL0W+BXw8YjI7XZt5J6SWcW1+ri/8LhyJZb+AdgeOEsSwIsR\nsX9eu05KZlU3jIsnS5RY+jTw6aG06aRkVnX+7puZJSWxr5k4KZlVnAsHmFlawsvhmllK3FMys6S4\ncICZpaSv1uiWtAVwE/AKYBxweUScNhKBmdkI6afhW0T8RdJhEbFO0hhgrqRDImLuCMVnZt3WT0kJ\nICLWZTfHUbuU/OmuRmRmIyux4VvhF3IlbSZpAfAEcENELO5+WGY2UmIgSm0jpTApRcRAROxHbcW4\nGZJmdj0qMxs57a8S0BWlP33L1ti9EngLcGP9Yw+ydPD2eMYzgfGdis/MMstXLuShlYs63m4M9NEl\nAZImAOsj4llJWwLvAr7RuN9eTO5SeGa20aSJ+zBp4qZiINfdcVFH2h0YxhXdWTGAM6jNN58TEac3\nPP4G4HxgGvC1iPheUZtFPaXXABdI2ozaUO/CiPhtO8GbWZoGyq9q+xIlq5msAU4Gji7bbtElAYuA\n6UMP18z6RbTfUypTzeQp4ClJ7yvbqJfDNau4AaLU1sRQq5mU4q+ZmFVcqzmlFetXsmLDyqaPZbry\nkZyTklnFBc2T0sQxuzBxzKaOz60v3t64S6lqJkPlpGRWcRtiQ7uHDlYzAR6nVs1kdot9VbZRJyWz\nihto88LIMtVMJO0E3AVsCwxI+jywd0SsbdWuk5JZxQ20GL6VUaKaySpeOsQr5KRkVnF9tZ6SmY1+\nw+kpdYOTklnFtTun1C1OSmYVt4G2P33rCicls4obzhdyu8FJyaziojsXZrfNScms4txTMrOktPqa\nSa84KZlVnD99M7Ok+DolM0vKML6Q2xVOSmYV5+GbmSUltYluL4drVnEDEaW2ZiTNkvSApGWSvtJi\nnx9kj98raVpRPH2XlFazptchdNVof35/HFjV6xC6avnKhb0OYcjKrdD98t5UXTWTWcDewGxJUxr2\neS8wKSL2BE4AziqKp++S0ppR/ks72p/fqnii1yF0VTeKRXbbMHpKg9VMIuJFYGM1k3ofAC4AiIh5\nwHaSdsyLp++Skpl11kDJf02UqWbSbJ9d8+LxRLdZxQ3jkoCyH9s1rs+de5yGu+qcpLQ+TzSrkIgo\nvSB/M0P9/a0/n6S3Af8UEbOy+6cBA/WluyX9GLgxIi7O7j8AHBrRehw/7J7ScF8UM+udYf7+lqlm\nMgc4Cbg4S2LP5iUk8PDNzNpUpppJRFwl6b2SlgPPA8cVtTvs4ZuZWSf11advZS7U6leSzpP0hKT+\n+0y5gKSJkm6QdL+k+yT9Xa9j6iRJW0iaJ2mBpMWSvtXrmPpZ3/SUsgu1HgQOp1Yu+C5gdkQs6Wlg\nHSLp7cBa4GcR8aZex9NJWUHCnSJigaRtgLuBo0fLewcgaauIWCdpDDAX+FJEzO11XP2on3pKZS7U\n6lsRcQvwTK/j6IaIWBURC7Lba4ElwM69jaqzImJddnMctfmVp3sYTl/rp6RU5kItS1z2Sc00YF5v\nI+ksSZtJWgA8AdwQEYt7HVO/6qek1B/jTGspG7r9Evh8Xi35fhQRAxGxH7WrlWdImtnjkPpWPyWl\nP/DSmuQTqfWWrA9IGgtcCvw8Ii7rdTzdEhF/Aq4E3tLrWPpVPyWlwQu1JI2jdqHWnB7HZCVIEnAu\nsDgizuh1PJ0maYKk7bLbWwLvAub3Nqr+1TdJKSLWU7sy9FpgMfBfo+zTm4uA24DJklZKKrzIrI8c\nDHwcOEzS/Gyb1eugOug1wPXZnNI84IqI+G2PY+pbfXNJgJlVQ9/0lMysGpyUzCwpTkpmlhQnJTNL\nipOSmSXFScnMkuKkZGZJcVIys6T8P9Z0ybinLqtDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa2cca710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plot_T, exponent=[0.2, 5, 0.1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logit normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the logit transform\n",
    "\n",
    "$$ a \\mapsto \\log (\\frac{a}{1-a})$$\n",
    "\n",
    "expands values that are close to zero and close to one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def logit(a):\n",
    "    '''logit, with base 10 logarithm'''\n",
    "    return np.log10(a/(1-a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_T_logit(vmax=5):\n",
    "    # cut values that are zero\n",
    "    T_sat = np.where(T>0, T, 1e-16)\n",
    "    lT = logit(T_sat)\n",
    "    plt.imshow(lT, vmin=-vmax, vmax=vmax, **imopts)\n",
    "    plt.xticks(range(4))\n",
    "    plt.yticks(range(4))\n",
    "    plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The logit normalization works best with a diverging colormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Diverging\n",
    "imopts['cmap'] = 'RdBu_r'\n",
    "#imopts['cmap'] = 'RdYlBu_r'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAEACAYAAAD4GBC1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEZVJREFUeJzt3X2QXXV9x/H35y4JDw0WZlASIBaE2FGgZCNNMhBMGMWJ\nVKXOtGOZgQ7IdGwtNDotWmSmFh9HHRWEsaMYqWiVTlWgDEakCCQrEnnYyEOC3TjBCQ8BkUgJEUj2\nfvvHXtYl5O4999xz9vzuPZ/XzBn2uOfhO8h+5vs7v/OgiMDMrGqNqgswMwOHkZklwmFkZklwGJlZ\nEhxGZpYEh5GZJWGfqgsws/4l6WHg/4BxYFdELM57LIeRmfUigBUR8XSvB/Iwzcx6pSIO4jAys14E\n8D+S7pb0N70cyMM0M+vFyRHxuKRXAzdLeigi1uU5UM9hJMkPt5lVJCJ6GiJ1+/e75/ki4vHWP38t\n6VpgMVBNGAHs3nxnEYfJ5JLLruSjq3rqBruyjqNn7FwA/37ZZzln1Ydm7Hwv7B6fsXMBfPPyz3H2\nBRfO2Pkeeuq5GTsXwJrVl/H281bNyLk+sKyY/zbfrz/KtN2X41cvW5d0ADAUEc9K+gPgbcAleevw\nMM2s5oaUsbl6ZQ91KHCtJvbfB/iPiPhR3jocRmY1N5RzoBcRW4CFRdXRd2G0fMmiqkso1cIlJ1dd\nQqn+ZPFJVZdQqmOGl1RdQtdmNzKmUckj+r4LoxVL31R1CaVauHSww+iEAQ/bBYuWVl1C1zIP00rW\nd2FkZsXKO0wrmsPIrObcGZlZElJ5DMNhZFZz7ozMLAm+ZmRmScg8tV8yh5FZzXmYZmZJ8DDNzJLg\nzsjMkuDOyMyS4DAysyR4mGZmSeh1al/SEHA38EhEvDPvcRxGZjVXwDBtFbAROLCXg6TyWIqZVWRI\nyrTsjaQjgNOBr9HjJ4vcGZnVXI+d0ReBC4FX9VqHw8is5tp1Pfe/8BwPvNj+gwaS3gE8GRGjklb0\nWkfHMJK0ErgUGAK+FhGf6fWkZpaORpswOmG/OZyw35zJ9Wue+/Wem5wEvEvS6cB+wKskXR0Rf52r\njul+2bpKfgWwEngjcKakN+Q5kZmlSUPKtOwpIj4SEfMj4ijgr4Af5w0i6NwZLQY2R8TDAJKuAc4A\nNuU9oZmlZWj2UFGH6umDrp3C6HBg65T1R4D++/yBmbW1t66nWxFxO3B7L8foFEb+dLXZgGsk8jxI\npzB6FJg/ZX0+E93Ry1xy2ZWTPy9fsmjgPydkVoWxe+9k8+j6wo+rRhq3GyqiffMjaR/gF8BbgMeA\nnwFnRsSmKdvE7s13ll1nZdZRzPfMU/XC7pK/zFexh55qPzXd7z6w7Ggioqe2RlLccvyfZtr2Lfff\n1fP5pjNtZxQRuyWdD9zExNT+6qlBZGb9r4hrRkXoeJ9RRKwB1sxALWZWgQJn03riO7DNak5+Ib+Z\npaAxlMYFbIeRWc31zTUjMxtsDiMzS4KHaWaWBHdGZpaEoVme2jezBPTLs2lmNuDyDtMk7cfEk/r7\nArOB6yPiorx1OIzMak45L2BHxPOSTo2Ina3nWEckLYuIkTzHcxiZ1Vwvw7SI2Nn6cTYTz68+nbuO\n3FWY2UBQQ5mWve4rNSRtAJ4Abo2IjXnrcGdkVnPt7jP62banuOuJ30y7b0Q0gYWS/hC4SdKKiLgt\nTx0OI7Oaa7R5an/paw9l6WsPnVz/8n1jbY8REc9IuhE4EbgtVx15djKzwaFGI9Pyiv2kQyQd1Pp5\nf+A0YDRvHe6MzGquh8dB5gHfkNRgorH5ZkTckvdgDiOzmuthav9+YFFRdTiMzGoubxgVzWFkVnOp\nfB3EYWRWcxryg7JmloChWWnEQBpVmFllfM3IzJIwUGF0xbZXF3GYJL3v2e9XXUKpvnLgW6suoVQ3\njD5WdQnJ8wVsM0vCQHVGZta/HEZmlgR/HcTMktDw1L6ZpcDDNDNLgmfTzCwJjZyPg0iaD1wNvAYI\n4KsR8aW8dTiMzGquh2HaLuCDEbFB0hzgHkk3R8SmPAdzGJnVXA/vM9oGbGv9vEPSJuAwwGFkZt0r\nYjZN0pHAMLA+7zEcRmY11+tsWmuI9l1gVUTsyHsch5FZzbWbTVu7cQvrNm6Zfl9pFvA94FsRcV0v\ndTiMzGpOjb3Ppi0/7hiWH3fM5Pqnv3/ry/eTBKwGNkbEpb3WkcYNBmZWncZQtuWVTgbOAk6VNNpa\nVuYtw52RWd3lvOkxIkYosKFxGJnVnN+BbWZp2Gd21RUADiOz2vOzaWaWhjazaTPNYWRWdw4jM0uB\nh2lmlgZ3RmaWhH4JI0lfB/4MeDIiji+/JDObSZo1q+oSgGx3T14F5L7F28wSl/9xkEJ17IwiYl3r\nXSVmNoDaPSg703zNyKzuBmk2bc3qyyZ/PmZ4CQsWLS3isGY2xfaxUbaPjRZ+3IHqjN5+3qoiDmNm\n0zh4wTAHLxieXN+y5qpiDpwzjIqe3EqjPzOz6jQa2ZZXKnRyq2MYSfoOcAfweklbJZ1b1MnNrHqa\nNTvTsqeIWAdsL6qOLLNpZxZ1MjNL0CBdMzKz/uVn08wsDW06o9vW38vt6++dsTIcRmZ1p713RiuW\nnsiKpSdOrn/s8tWllpFGf2Zm1VEj27LnbgVPbrkzMqu5aNMZddyv4Mkth5FZ3Xk2zcyS4Nk0M0tB\n3mFa0RxGZnXnMDKzJDiMzCwJDiMzS0E00oiBNKows+pIVVcAOIzMzMM0M0uBp/bNLA2+6dHMkuDO\nyMyS4DAysxSkMrWfRiSaWXVyvs8IQNJKSQ9JGpP04V7KcBiZ1Z2UbXnFbhoCrmDic0VvBM6U9Ia8\nZTiMzOouf2e0GNgcEQ9HxC7gGuCMvGWkMVg0s8r0cJ/R4cDWKeuPAEvyHqyQMLrw/E8UcZgk3XDe\ne6suoVRXbvmXqkso1Q1H/W3VJaSvTRitXbuWtWvXTrdnFFmGOyOzmos2z6adsnw5pyxfPrn+yU99\nas9NHgXmT1mfz0R3lIvDyKzmxpu5G5y7gQWSjgQeA94D5H5Jv8PIrObyRlFE7JZ0PnATMASsjohN\neetwGJnVXP7GCCJiDbCmiDocRmY1F1HodejcHEZmNddLZ1Qkh5FZzSWSRQ4js7pzZ2RmSRj3NSMz\nS0EiWeQwMqs7D9PMLAme2jezJDSrLqDFYWRWc4k0Rg4js7prJpJGDiOzmhtPI4scRmZ1l0hj5Hdg\nm9Vdk8i0dEPSX0p6UNK4pEVZ9nEYmdVcRLalS/cD7wamfW/tVB6mmdVcGTc9RsRDAGrzStu9cRiZ\n1Vwq14wcRmY1l/dBWUk3A3P38quPRMQN3R7PYWRWc+3uM7r7jhHu/ulI2/0i4rQi6+gYRpLmA1cD\nr2HiPUxfjYgvFVmEmVVnvM3zIMNLlzG8dNnk+le++Jm8p8h04SjLbNou4IMRcSywFPj7Xr6nbWZp\naUZkWroh6d2StjKRGTdK6vjS/o6dUURsA7a1ft4haRNwGJD7kyRmlo4yXq4WEdcC13azT1fXjFof\naxsG1nezn5mlq++eTZM0B/gusCoidkz93fjjo7/fbs5cGgfOK6xAM5uwfWyU7WOjnTfsUrtrRjMt\nUxhJmgV8D/hWRFy35++H5g0XXZeZ7eHgBcMcvOD3f2tb1lxVyHF3NdNIoyyzaQJWAxsj4tLySzKz\nmdRPr509GTgLuE/SSz3iRRHxw/LKMrOZMp5IGmWZTRvBD9SaDay+u4BtZoPJL1czsyS4MzKzJPTN\nNSMzG2y7HEZmlgIP08wsCU13RmaWAs+mmVkSyhimSfoc8A7gReCXwLkR8cx0+/hmRrOaG4/ItHTp\nR8CxEXEC8L/ARZ12cBiZ1VyzGZmWbkTEzRHx0hO464EjOu3jYZpZzc3A1P57ge902shhZFZzZX4d\nRNLFwIsR8e1Ox3MYmdVcuzuwx+69k82j7V/q2unrIJLOAU4H3pKlDoeRWc21C6PXLVzC6xYumVy/\n6arsHwWStBK4EFgeEc9n2cdhZFZzJT2bdjkwG7i59Ynrn0bE+6fbwWFkVnNlhFFELOh2H4eRWc35\nqX0zS8KLu/vkhfxmNtjcGZlZEhxGZpYEh5GZJWG3w8jMUuDOqE+8c/iwqkso1QFnf6XqEkr1g/Fn\nqy6hNPteUMxxXhz3bJqZJcCdkZklwWFkZklwGJlZEsabvmZkZglwZ2RmSXAYmVkSXijhQVlJHwfe\nBQTwG+CciNg63T7+OohZzY03I9PSpc9GxAkRsRC4Dvhopx3cGZnVXEkvV5t6t+kc4KlO+ziMzGqu\nrGtGkj4JnA3sBJZ22t5hZFZz7cJo+9go28dG2+7X6VNFEXExcLGkfwa+CJw7XR0OI7OaizZhdNDR\nCzno6IWT61vWXPXy/Tp8qmiKbwM/6LSRw8is5rr9dHUWkhZExFhr9QygfYvV4jAyq7lmOU/tf1rS\nHwPjwC+Bv+u0g8PIrObK6Iwi4i+63cdhZFZzkcajaQ4js7qL8OMgZpaAMoZpeTiMzGqu3dT+THMY\nmdWcw8jMkjDuF/KbWQr6ojOStB9wO7AvMBu4PiIumonCzGxm9MUF7Ih4XtKpEbFT0j7AiKRlETEy\nQ/WZWcn6Zmo/Ina2fpwNDAFPl1qRmc2oVG567PimR0kNSRuAJ4BbI2Jj+WWZ2UxpNiPTUraOYRQR\nzdarI48A3ixpRelVmdmMiWZkWsqWeTYtIp6RdCNwInDb1N+NP/77twNozlwaB84rqj4za7l95A7W\n/uSOwo/bF1P7kg4BdkfEbyXtD5wGXLLndkPzhksqz8xesnzZSSxfdtLk+ic++/lCjltm1yPpH4HP\nAYdExLTXmzt1RvOAb0hqMDGk+2ZE3FJMmWaWgrLCSNJ8JhqYX2XZvtPU/v3AogLqMrNElXhx+gvA\nh4Drs2zsO7DNaq6M+4wknQE8EhH3Scq0j8PIrObaDdN+99gD/O6xB9vuN83XQS4GLgLeNnXzTnU4\njMxqrt0wbd+5x7Lv3GMn1397z3+97Pftvg4i6TjgKODnra7oCOAeSYsj4sl2dTiMzGquufvFQo8X\nEQ8Ah760LmkL8KZeZ9PMbMBFc7z0U2TZyGFkVnMxXm4YRcTrsmznMDKruRnojDJxGJnVnMPIzJLg\nMDKzJBQ9m5aXw8is5prujMwsBR6mmVkSHEZmloSy7zPKymFkVnPujMwsCQ4jM0tCc/euqksAHEZm\ntefOyMySkEoYdfxumpkNtmZzPNPSDUn/KukRSaOtZWWnffoujJrPPl51CaUau/fOqkso1U/Wra26\nhFLdPlL8d83KFuPjmZZuDwt8ISKGW8sPO+3Qd2EUO7ZVXUKpNo+ur7qEUt0xsq7qEkpVxkcWyxbN\n8UxLDtnexN/Sd2FkZsUqMYwukPRzSaslHdRpY1/ANqu5vE/td/g6yL8BH2utfxz4PHDetMfr9ZtJ\nksr7Nq6ZTSsiuhoK7anbv98855N0JHBDRBw/3XY9d0a9/ssws+qU9fcraV5EvDTb9G7g/k77eJhm\nZmX4jKSFTMyqbQHe12mHnodpZmZF6KvZNEkrJT0kaUzSh6uup0iSvi7pCUkd29l+I2m+pFslPSjp\nAUn/UHVNRZK0n6T1kjZI2ijp01XX1I/6pjOSNAT8Angr8ChwF3BmRGyqtLCCSDoF2AFc3elCX7+R\nNBeYGxEbJM0B7gH+fFD+vwOQdEBE7JS0DzAC/FNEjFRdVz/pp85oMbA5Ih6OiF3ANcAZFddUmIhY\nB2yvuo4yRMS2iNjQ+nkHsAk4rNqqihURO1s/zgaGgGk/5Wyv1E9hdDiwdcr6I63/zfpIa5p3GBio\nW80lNSRtAJ4Abo2IjVXX1G/6KYz6YzxpbbWGaN8FVrU6pIEREc2IWAgcAbxZ0oqKS+o7/RRGjwLz\np6zPZ6I7sj4gaRbwPeBbEXFd1fWUJSKeAW4ETqy6ln7TT2F0N7BA0pGSZgPvAf674posA0kCVgMb\nI+LSquspmqRDXnr2StL+wGnAaLVV9Z++CaOI2A2cD9wEbAT+c8BmY74D3AG8XtJWSedWXVOBTgbO\nAk7t5v02fWQe8OPWNaP1TDz6cEvFNfWdvpnaN7PB1jedkZkNNoeRmSXBYWRmSXAYmVkSHEZmlgSH\nkZklwWFkZklwGJlZEv4fhZ5DadgtaIUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa33e7290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_T_logit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS0AAAEACAYAAADm0SAGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFVZJREFUeJzt3XuwnVV5x/Hv7xwSLhNo2mEkIQlyC95JQpFkuCUotCko\nYB21VEpFbR0dFO20VVBrrSKFakWkOIiARhRaQTEMNyNEQkROuSQhkFADQpsEEkSu4SI5Zz/9Y78c\nTnb2Pvvdt7PP2u/vM/NO9rv3etf7DMx5Zq31rnctRQRmZqno63YAZmaNcNIys6Q4aZlZUpy0zCwp\nTlpmlhQnLTNLipOWmTVF0k6SBiStlLRG0tk1yp0vaZ2kVZLmtHrfHVqtwMyKKSJeknRURLwgaQdg\nuaTDI2L5K2UkHQvsHxEzJc0FvgXMa+W+bmmZWdMi4oXs40SgH3iyosjxwPeysgPAZEl7tHJPJy0z\na5qkPkkrgc3A0ohYU1FkGrB+xPkGYHor93TSMrOmRUQpImZTTkRHSlpQpZgqL2vlni2PaUnyy4tm\nXRIRlQmhIY3+/da6X0Q8I+k64GDgFyN+2gjMGHE+PfuuaW0ZiH9xyaXtqCaXLy+6hs+dcuKY3e+8\n/sPH7F4AP//u+Rz9gU+M2f2ee2lwzO4FsOzyCzjy5NPG7H43371hzO4FsOGWRUx/2yljcq+Bz/9J\nW+r5mF6bq9yF8b/bnEvaHRiMiKcl7QwcA3yx4rLFwGnAlZLmAU9HxOZW4vXTQ7OC61fOxtr2bbKp\nwPck9VEeavp+RNws6SMAEXFRRFwv6VhJDwLPA6e2Gq+TllnB9TfZwYyI1cBBVb6/qOK8rU3r5JLW\nkbNe3+0QOmrf2XO7HUJHvfbAQ7odQkftts+sbofQsIl9ObPWUGfjyMtJa5xx0kpbikkrd/dwnEgu\naZlZezXbPewWJy2zgnNLy8ySktoMcycts4JzS8vMkuIxLTNLSu4pD+OEk5ZZwbl7aGZJcffQzJLi\nlpaZJcUtLTNLipOWmSUlte5hapNhzazNJvYp11FJ0gxJSyXdL+k+STVXr5T0VkmDkv681Xjd0jIr\nuBa6h1uBT0XESkmTgLslLYmItSMLSeoHzgFuZPv14hvmlpZZwfVLuY5KEbEpIlZmn7cAa4E9q9zi\n48BVwG/bEa9bWmYF146BeEl7A3OAgYrvpwEnAG8D3kqLO/GAk5ZZ4dUaiF/9++e57+Xn616fdQ2v\nAk7PWlwjnQd8JiJCkmhD97Bu0pK0MLtxP/CdiDin1Zua2fjRVyNpzdppErN2mjR8fuXz2/fuJE0A\nrgYuj4hrqlTzx5R34gHYHfgzSVsjYnGz8Y6atLIBtAuAoynvVXanpMWVA21mli412T/MWk6XAGsi\n4rxqZSJi3xHlLwOubSVhQf2W1iHAgxHxSHbTKyn3T520zHpE/8T+Zi89DDgZuFfSiuy7M4G9YPtd\nedqlXtKaBqwfcb4B6O2dF8wKptmWVkQsp4EZCBHR8p6HUD9pect7sx7Xl9h7PPWS1kZgxojzGZRb\nW9v48qJXx9+OnPX6nt/my6wbnn14Fc8+vKrt9aovrema9ZLWXcDMbA7Go8D7gJMqC33ulBPbHpiZ\nbWu3fWZts6/ixqWXt6XenmppRcSgpNOAmyhPebjETw7NekuzY1rdUneeVkTcANwwBrGYWRe08PSw\nKzwj3qzg5I0tzCwlff29NRBvZj2u58a0zKy3OWmZWVLcPTSzpLilZWZJ6Z/gKQ9mlpDUZsSn1Zk1\ns7ZTv3Id210nXSpps6TVVeuVdpd0o6SV2W49H2hHvE5aZgWn/r5cRxWXAQtHqfo0YEVEzAYWAF+T\n1HLvzt1Ds4JrtnsYEbdliynU8hhwYPZ5N+B3ETHY1M1GcNIyK7gOvsZzMXCLpEeBXYH3tqNSJy2z\ngqs1T+u/Nz3BnZt/10rVZwIrI2KBpP2AJZJmRcRzrVTqpGVWcH01VnmYt9cezNtrj+HzC+9d12jV\nhwJnAUTEQ5IeBl5HeZ2+pnkg3qzg1NeX62jCA5R38kLSHpQT1m9ajdctLbOCa/Y1HklXAPOB3SWt\nB74ATIDhnXi+AlwmaRXlBtI/RsSTrcbrpGVWcDWmM9QVEdstvV7x+xPAO5uqfBROWmYF12zS6hYn\nLbOC67XdeMysx6nfL0ybWUL6J6SVBtKK1szazmNaZpaUQiatBXdMb0c149Llj3++2yF01Mmv+VC3\nQ+ioe66+otshjHseiDezpBSypWVm6XLSMrOkeDceM0tKn6c8mFlK3D00s6Sk9vQwrWjNrO36+vtz\nHZVy7MbzfkmrJN0r6ZeSDqxWruF421GJmaWrg7vx/AY4MiIOBL4EfLsd8bp7aFZwLaynNepuPBHx\nqxGnA0BbZqE7aZkV3Bg9PfwQcH07KnLSMiu4Tj89lHQU8EHgsHbU56RlVnC1nh4uW/Mwt615uLW6\ny4PvFwMLI+KplirLOGmZFZz6qi8COP/N+zP/zfsPn5/946WN1SvtBfwYODkiHmwhxG04aZkVXY2k\nVU+O3Xj+CfhD4FuSALZGxCGthuukZVZ0TU4uzbEbz4eBDzdV+SictMwKzmvEm1ladpjY7Qga4qRl\nVnCpvXvopGVWdE0OxHeLk5ZZ0TlpmVlK3D00s7S4pWVmSem1pCXpUuA44PGIeEvnQzKzsaQJE7od\nQkPydGbrLfRlZinr6893jBN1W1r1Fvoys7TVemF6vPKYllnRFfHp4YZbFg1/3m2fWey2z6x2VGtm\nI5See4zYsqnt9RaypTX9bae0oxozG0XfrlNh16nD56XNK9tUcfNJS9JC4DygH/hORJxTpcwC4OuU\nl615IiIWNH1D3D00sya7h5L6gQuAo4GNwJ2SFkfE2hFlJgP/AfxpRGyQtHvL4eYI7ArgduAASesl\nndrqTc1s/NCEibmOKg4BHoyIRyJiK3AlcEJFmb8Ero6IDQAR8USr8eZ5ejjqQl9mlrjmu4fTgPUj\nzjcAcyvKzAQmSFoK7Ap8IyK+3+wNwd1Ds8Jr4d3DyFFmAnAQ8HZgF+BXku6IiHXN3tRJy6zoarS0\nfjFwD7cO3DPalRuBGSPOZ1BubY20nvLg+4vAi5KWAbMAJy0za5Kqt7QWzDuYBfMOHj7/l29eUlnk\nLmBmNvn8UeB9QOVw0k+BC7JB+x0pdx//vZVwnbTMiq5G0qonIgYlnQbcRHnKwyURsVbSR7LfL4qI\nByTdCNwLlICLI2JNK+E6aZkVXDSZtAAi4gbghorvLqo4/yrw1aZvUsFJy6zoijgj3swSVsR3D80s\nXa10D7vBScus6Jy0zCwpTlpmlhQnLTNLSfSllQbSitbM2k/qdgQNcdIyKzp3D80sJZ7yYGZp8eRS\nM0uKW1pmlhQnLTNLSWpTHtJKsWbWfurLd1S7VFoo6QFJ6yR9ukaZ87PfV0ma02q4TlpmRSflO7a7\nbHgLsYXAG4GTJL2hosyxwP4RMRP4W+BbrYbrpGVWdM23tPJsIXY88D2AiBgAJkvao5VwnbTMCi7U\nl+uootoWYtNylJneSrxtGYG75+or2lHNuPRGdul2CB111rqPdTuEjrpn5hHdDmH8qzFetWzZMpYt\nWzbalXm2EAOo7Fvmva6qtB4bmFnbRY13D4+YP58j5s8fPj/rK1+pLJJnC7HKMtOz75rm7qFZwQ2V\nItdRxfAWYpImUt5CbHFFmcXAKQCS5gFPR8TmVuJ1S8us4Jrtq+XcQux6ScdKehB4Hji11XidtMwK\nrnojKp+cW4id1vwdtuekZVZwES2Ni485Jy2zgmulpdUNTlpmBZdYznLSMis6t7TMLClDHtMys5Qk\nlrOctMyKzt1DM0uKpzyYWVJK3Q6gQU5aZgWXWEPLScus6EqJZS0nLbOCG0orZzlpmRVdYg0tr6dl\nVnQlItfRCEl/JGmJpF9L+pmkyTXKTZZ0laS1ktZka26NyknLrOAi8h0N+gywJCIOAG7Ozqv5BnB9\nRLwBOBBYW69iJy2zgitFvqNBw7vwZP+eWFlA0h8AR0TEpVBeVDAinqlXsZOWWcF1qKW1x4hllTcD\n1bYN2wf4raTLJN0j6WJJdXeScdIyK7ihiFxHpWzManWV4/iR5aI85b5a2tsBOAi4MCIOorwcc61u\n5DYXmVmB1Zqnddfty7nrV8trXhcRx9T6TdJmSVMiYpOkqcDjVYptADZExJ3Z+VW0I2lJmgEsAl5D\nOVt+OyLOr3edmaVhqMZ7PHPmHc6ceYcPn1/09XMaqXYx8NfAOdm/11QWyBLaekkHRMSvgaOB++tV\nnKeltRX4VESslDQJuFvSkoioO8pvZuNfh2bE/yvwX5I+BDwCvBdA0p7AxRFxXFbu48APsi3IHiLH\nbj11k1ZEbAI2ZZ+3SFoL7EmOR5NmNv51YhHAiHiScsup8vtHgeNGnK8C3tpI3Q2NaUnaG5gDDDRy\nnZmNXz377mHWNbwKOD0itoz8beixFa+WmzSFvl2nti1AMysrPfcYsWVT2+utNaY1XuVKWpImAFcD\nl0fEdgNq/VPntDsuM6vQt+tUGNEgKG1e2ZZ6t5bSylp5nh4KuARYExHndT4kMxtLvbjc8mHAycC9\nkl7pB54RETd2LiwzGytDiWWtPE8Pl+OZ82Y9q2cH4s2sN3kRQDNLiltaZpaUnhvTMrPettVJy8xS\n4u6hmSWl5JaWmaUktaeHnn9lVnCliFxHIyS9R9L9koYkHVSjzAxJS7Ny90n6RJ663dIyK7hOLE0D\nrAbeBVw0Spmm1upz0jIruE6MaUXEAwDlV5drlmlqrT4nLbOCGw9THhpZq89Jy6zgmu0eSloCTKny\n05kRcW0D9dRcq68aJy2zgqs1I37dPXfw4IraDZ/RduPJq95afdU4aZkVXK2kte/suew7e+7w+U2X\nNb0JV9WBrWbX6vOUB7OCGypFrqMRkt4laT0wD7hO0g3Z93tKui4r9spafUdJWpEdC+vV7ZaWWcF1\n4oXpiPgJ8JMq3w/vxtPsWn1OWmYF51UezCwpLw/22MYWZtbb3NIys6Q4aZlZUpy0zCwpg05aZpYS\nt7R6zEHvPqnbIXTUez76jW6H0FGf3Lq52yF0zMSpl7WlnpeH/PTQzBLilpaZJcVJy8yS4qRlZkkZ\nKnlMy8wS4paWmSUltaTl9bTMCu73g6VcRyPybCGWlTsjK7da0g8l7Vivbicts4LrxCKAvLqF2LJa\nBbLNLP4GOCgi3gL0A39Rr2J3D80KrkOLANbdQgx4lvLeh7tIGgJ2ATbWq9tJy6zgujWmFRFPSvoa\n8H/Ai8BNEfHzetc5aZkVXK2k9dS6FTy1bkXN61rdQkzSfsAngb2BZ4AfSXp/RPxgtOuctMwKLmok\nrcn7zWbyfrOHzx++Ydt3HduwhdjBwO0R8TsAST8GDgVGTVoeiDcruFIpch0tqDWw9QAwT9LO2XZi\nRwNr6lXmpGVWcKWhUq6jEXm2EIuIVcAi4C7g3uzSb9er291Ds4JrsRVVVZ4txLLzc4FzG6nbScus\n4CKtVw+dtMyKLiKt13ictMwKrhPdw05y0jIruFpTHsYrJy2zgnPSMrOkDHljCzNLSU+1tCTtBNwK\n7AhMBH4aEWeMRWBmNjZ6aiA+Il6SdFREvCBpB2C5pMMjYvkYxWdmHdZzUx4i4oXs40TKi3Q92dGI\nzGxMpTa5tO67h5L6JK0ENgNLI6LuC41mlo4xeGG6reomrYgoRcRsYDpwpKQFHY/KzMZMlCLXMV7k\nfnoYEc9kb2cfDPxi5G9Dj726UJgmTaFv16ntis/MMrfePsCttw+0vd6emvIgaXdgMCKelrQzcAzw\nxcpy/VPndCg8M3vF/EPnMv/QucPnX/7aN9tSbydaUZL+DXgH8DLwEHBqRDxTo2w/5eVpNkTEO+vV\nXa97OBW4JRvTGgCujYibGwnezMa3DnUPfwa8KSJmAb8GRpsqdTrlxf9y3aTelIfVQM09y8wsfR1a\nT2vJiNMB4N3VykmaDhwLnAX8XZ66vXKpWcFFRK6jBR8Erq/x29eBfwByD6z5NR6zgqvV9Xvx0ft4\n8dH7a16XZzceSZ8FXo6IH1a5/h3A4xGxopFZCU5aZgVXq3u445Q3seOUNw2fP333j7b5vd5uPJI+\nQLnr9/YaRQ4Fjpd0LLATsJukRRFxymj1OmmZFVxp8OW21ylpIeVu3/yIeKlamYg4EzgzKz8f+Pt6\nCQs8pmVWeFEaynU06JvAJGCJpBWSLoRtd+OpFkqeit3SMiu4GGo4IdWvM2Jmje+32Y1nxPe3Ul5R\npi4nLbOCa6IV1VVOWmYF56RlZklx0jKzpHTi6WEnOWmZFVzJLS0zS4m7h2aWFCctM0tKJ+ZpdZKT\nllnBuaVlZklx0jKzpJQGt3Y7hIY4aZkVnFtaZpaU1JKWl6YxK7hSaSjX0QhJX5K0StJKSTdLmlGj\n3EJJD0haJ+nTeepOLmmVnnus2yF01LMPr+p2CB11x/Lbuh1CR3ViX8JOi6GhXEeDzo2IWdlGz9cA\nX6gskG0ddgGwEHgjcJKkN9SrOLmkFVs2dTuEjur5pPVLJ63xphOLAEbEcyNOJwFPVCl2CPBgRDwS\nEVuBK4ET6tXtMS2zguvUmJaks4C/Al4A5lUpMg1YP+J8AzC3SrltJNfSMrP2Kg2+nOuoJGmJpNVV\njncCRMRnI2Iv4LuUtwqr1NS+ZGpxPzMktX+nRzPLJSLUyvWN/v02cz9JewHXR8SbK76fB/xzRCzM\nzs8AShFxzmj1tdw9bPU/mpl1T6f+fiXNjIh12ekJwIoqxe4CZkraG3gUeB9wUr26PaZlZp1wtqTX\nAUPAQ8BHobwbD3BxRBwXEYOSTgNuAvqBSyJibb2KW+4empmNpaQG4puZiJYKSZdK2ixpdbdjaTdJ\nMyQtlXS/pPskfaLbMbWTpJ0kDWQTKddIOrvbMfWyZFpa2US0/wGOBjYCdwIn5WlOpkDSEcAWYFFE\nvKXb8bSTpCnAlIhYKWkScDdwYq/8vwOQtEtEvCBpB2A55d2Sl3c7rl6UUkurqYloqYiI24Cnuh1H\nJ0TEpohYmX3eAqwF9uxuVO0VES9kHydSHp95sovh9LSUkla1iWjTuhSLNSl7UjQHSG/q+Cgk9Ula\nCWwGlkbEmm7H1KtSSlpp9GOtpqxreBVwetbi6hkRUcres5sOHClpQZdD6lkpJa2NwMg3xWdQbm1Z\nAiRNAK4GLo+Ia7odT6dExDPAdcDB3Y6lV6WUtIYnokmaSHki2uIux2Q5SBJwCbAmIs7rdjztJml3\nSZOzzzsDx1B9MqW1QTJJKyIGgVcmoq0B/rPHnj5dAdwOHCBpvaRTux1TGx0GnAwcJWlFdizsdlBt\nNBW4JRvTGgCujYibuxxTz0pmyoOZGSTU0jIzAyctM0uMk5aZJcVJy8yS4qRlZklx0jKzpDhpmVlS\nnLTMLCn/D4Yy7RbmL3ysAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa314c6d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_T_logit(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interactive choice of vmax: 4 (which corresponds to minimum probability of 0.0001) seems good for this example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAEACAYAAAD4GBC1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEZVJREFUeJzt3X2QXXV9x/H35y4JDw0WZlASIBaE2FGgZCNNMhBMGMWJ\nVKXOtGOZgQ7IdGwtNDotWmSmFh9HHRWEsaMYqWiVTlWgDEakCCQrEnnYyEOC3TjBCQ8BkUgJEUj2\nfvvHXtYl5O4999xz9vzuPZ/XzBn2uOfhO8h+5vs7v/OgiMDMrGqNqgswMwOHkZklwmFkZklwGJlZ\nEhxGZpYEh5GZJWGfqgsws/4l6WHg/4BxYFdELM57LIeRmfUigBUR8XSvB/Iwzcx6pSIO4jAys14E\n8D+S7pb0N70cyMM0M+vFyRHxuKRXAzdLeigi1uU5UM9hJMkPt5lVJCJ6GiJ1+/e75/ki4vHWP38t\n6VpgMVBNGAHs3nxnEYfJ5JLLruSjq3rqBruyjqNn7FwA/37ZZzln1Ydm7Hwv7B6fsXMBfPPyz3H2\nBRfO2Pkeeuq5GTsXwJrVl/H281bNyLk+sKyY/zbfrz/KtN2X41cvW5d0ADAUEc9K+gPgbcAleevw\nMM2s5oaUsbl6ZQ91KHCtJvbfB/iPiPhR3jocRmY1N5RzoBcRW4CFRdXRd2G0fMmiqkso1cIlJ1dd\nQqn+ZPFJVZdQqmOGl1RdQtdmNzKmUckj+r4LoxVL31R1CaVauHSww+iEAQ/bBYuWVl1C1zIP00rW\nd2FkZsXKO0wrmsPIrObcGZlZElJ5DMNhZFZz7ozMLAm+ZmRmScg8tV8yh5FZzXmYZmZJ8DDNzJLg\nzsjMkuDOyMyS4DAysyR4mGZmSeh1al/SEHA38EhEvDPvcRxGZjVXwDBtFbAROLCXg6TyWIqZVWRI\nyrTsjaQjgNOBr9HjJ4vcGZnVXI+d0ReBC4FX9VqHw8is5tp1Pfe/8BwPvNj+gwaS3gE8GRGjklb0\nWkfHMJK0ErgUGAK+FhGf6fWkZpaORpswOmG/OZyw35zJ9Wue+/Wem5wEvEvS6cB+wKskXR0Rf52r\njul+2bpKfgWwEngjcKakN+Q5kZmlSUPKtOwpIj4SEfMj4ijgr4Af5w0i6NwZLQY2R8TDAJKuAc4A\nNuU9oZmlZWj2UFGH6umDrp3C6HBg65T1R4D++/yBmbW1t66nWxFxO3B7L8foFEb+dLXZgGsk8jxI\npzB6FJg/ZX0+E93Ry1xy2ZWTPy9fsmjgPydkVoWxe+9k8+j6wo+rRhq3GyqiffMjaR/gF8BbgMeA\nnwFnRsSmKdvE7s13ll1nZdZRzPfMU/XC7pK/zFexh55qPzXd7z6w7Ggioqe2RlLccvyfZtr2Lfff\n1fP5pjNtZxQRuyWdD9zExNT+6qlBZGb9r4hrRkXoeJ9RRKwB1sxALWZWgQJn03riO7DNak5+Ib+Z\npaAxlMYFbIeRWc31zTUjMxtsDiMzS4KHaWaWBHdGZpaEoVme2jezBPTLs2lmNuDyDtMk7cfEk/r7\nArOB6yPiorx1OIzMak45L2BHxPOSTo2Ina3nWEckLYuIkTzHcxiZ1Vwvw7SI2Nn6cTYTz68+nbuO\n3FWY2UBQQ5mWve4rNSRtAJ4Abo2IjXnrcGdkVnPt7jP62banuOuJ30y7b0Q0gYWS/hC4SdKKiLgt\nTx0OI7Oaa7R5an/paw9l6WsPnVz/8n1jbY8REc9IuhE4EbgtVx15djKzwaFGI9Pyiv2kQyQd1Pp5\nf+A0YDRvHe6MzGquh8dB5gHfkNRgorH5ZkTckvdgDiOzmuthav9+YFFRdTiMzGoubxgVzWFkVnOp\nfB3EYWRWcxryg7JmloChWWnEQBpVmFllfM3IzJIwUGF0xbZXF3GYJL3v2e9XXUKpvnLgW6suoVQ3\njD5WdQnJ8wVsM0vCQHVGZta/HEZmlgR/HcTMktDw1L6ZpcDDNDNLgmfTzCwJjZyPg0iaD1wNvAYI\n4KsR8aW8dTiMzGquh2HaLuCDEbFB0hzgHkk3R8SmPAdzGJnVXA/vM9oGbGv9vEPSJuAwwGFkZt0r\nYjZN0pHAMLA+7zEcRmY11+tsWmuI9l1gVUTsyHsch5FZzbWbTVu7cQvrNm6Zfl9pFvA94FsRcV0v\ndTiMzGpOjb3Ppi0/7hiWH3fM5Pqnv3/ry/eTBKwGNkbEpb3WkcYNBmZWncZQtuWVTgbOAk6VNNpa\nVuYtw52RWd3lvOkxIkYosKFxGJnVnN+BbWZp2Gd21RUADiOz2vOzaWaWhjazaTPNYWRWdw4jM0uB\nh2lmlgZ3RmaWhH4JI0lfB/4MeDIiji+/JDObSZo1q+oSgGx3T14F5L7F28wSl/9xkEJ17IwiYl3r\nXSVmNoDaPSg703zNyKzuBmk2bc3qyyZ/PmZ4CQsWLS3isGY2xfaxUbaPjRZ+3IHqjN5+3qoiDmNm\n0zh4wTAHLxieXN+y5qpiDpwzjIqe3EqjPzOz6jQa2ZZXKnRyq2MYSfoOcAfweklbJZ1b1MnNrHqa\nNTvTsqeIWAdsL6qOLLNpZxZ1MjNL0CBdMzKz/uVn08wsDW06o9vW38vt6++dsTIcRmZ1p713RiuW\nnsiKpSdOrn/s8tWllpFGf2Zm1VEj27LnbgVPbrkzMqu5aNMZddyv4Mkth5FZ3Xk2zcyS4Nk0M0tB\n3mFa0RxGZnXnMDKzJDiMzCwJDiMzS0E00oiBNKows+pIVVcAOIzMzMM0M0uBp/bNLA2+6dHMkuDO\nyMyS4DAysxSkMrWfRiSaWXVyvs8IQNJKSQ9JGpP04V7KcBiZ1Z2UbXnFbhoCrmDic0VvBM6U9Ia8\nZTiMzOouf2e0GNgcEQ9HxC7gGuCMvGWkMVg0s8r0cJ/R4cDWKeuPAEvyHqyQMLrw/E8UcZgk3XDe\ne6suoVRXbvmXqkso1Q1H/W3VJaSvTRitXbuWtWvXTrdnFFmGOyOzmos2z6adsnw5pyxfPrn+yU99\nas9NHgXmT1mfz0R3lIvDyKzmxpu5G5y7gQWSjgQeA94D5H5Jv8PIrObyRlFE7JZ0PnATMASsjohN\neetwGJnVXP7GCCJiDbCmiDocRmY1F1HodejcHEZmNddLZ1Qkh5FZzSWSRQ4js7pzZ2RmSRj3NSMz\nS0EiWeQwMqs7D9PMLAme2jezJDSrLqDFYWRWc4k0Rg4js7prJpJGDiOzmhtPI4scRmZ1l0hj5Hdg\nm9Vdk8i0dEPSX0p6UNK4pEVZ9nEYmdVcRLalS/cD7wamfW/tVB6mmdVcGTc9RsRDAGrzStu9cRiZ\n1Vwq14wcRmY1l/dBWUk3A3P38quPRMQN3R7PYWRWc+3uM7r7jhHu/ulI2/0i4rQi6+gYRpLmA1cD\nr2HiPUxfjYgvFVmEmVVnvM3zIMNLlzG8dNnk+le++Jm8p8h04SjLbNou4IMRcSywFPj7Xr6nbWZp\naUZkWroh6d2StjKRGTdK6vjS/o6dUURsA7a1ft4haRNwGJD7kyRmlo4yXq4WEdcC13azT1fXjFof\naxsG1nezn5mlq++eTZM0B/gusCoidkz93fjjo7/fbs5cGgfOK6xAM5uwfWyU7WOjnTfsUrtrRjMt\nUxhJmgV8D/hWRFy35++H5g0XXZeZ7eHgBcMcvOD3f2tb1lxVyHF3NdNIoyyzaQJWAxsj4tLySzKz\nmdRPr509GTgLuE/SSz3iRRHxw/LKMrOZMp5IGmWZTRvBD9SaDay+u4BtZoPJL1czsyS4MzKzJPTN\nNSMzG2y7HEZmlgIP08wsCU13RmaWAs+mmVkSyhimSfoc8A7gReCXwLkR8cx0+/hmRrOaG4/ItHTp\nR8CxEXEC8L/ARZ12cBiZ1VyzGZmWbkTEzRHx0hO464EjOu3jYZpZzc3A1P57ge902shhZFZzZX4d\nRNLFwIsR8e1Ox3MYmdVcuzuwx+69k82j7V/q2unrIJLOAU4H3pKlDoeRWc21C6PXLVzC6xYumVy/\n6arsHwWStBK4EFgeEc9n2cdhZFZzJT2bdjkwG7i59Ynrn0bE+6fbwWFkVnNlhFFELOh2H4eRWc35\nqX0zS8KLu/vkhfxmNtjcGZlZEhxGZpYEh5GZJWG3w8jMUuDOqE+8c/iwqkso1QFnf6XqEkr1g/Fn\nqy6hNPteUMxxXhz3bJqZJcCdkZklwWFkZklwGJlZEsabvmZkZglwZ2RmSXAYmVkSXijhQVlJHwfe\nBQTwG+CciNg63T7+OohZzY03I9PSpc9GxAkRsRC4Dvhopx3cGZnVXEkvV5t6t+kc4KlO+ziMzGqu\nrGtGkj4JnA3sBJZ22t5hZFZz7cJo+9go28dG2+7X6VNFEXExcLGkfwa+CJw7XR0OI7OaizZhdNDR\nCzno6IWT61vWXPXy/Tp8qmiKbwM/6LSRw8is5rr9dHUWkhZExFhr9QygfYvV4jAyq7lmOU/tf1rS\nHwPjwC+Bv+u0g8PIrObK6Iwi4i+63cdhZFZzkcajaQ4js7qL8OMgZpaAMoZpeTiMzGqu3dT+THMY\nmdWcw8jMkjDuF/KbWQr6ojOStB9wO7AvMBu4PiIumonCzGxm9MUF7Ih4XtKpEbFT0j7AiKRlETEy\nQ/WZWcn6Zmo/Ina2fpwNDAFPl1qRmc2oVG567PimR0kNSRuAJ4BbI2Jj+WWZ2UxpNiPTUraOYRQR\nzdarI48A3ixpRelVmdmMiWZkWsqWeTYtIp6RdCNwInDb1N+NP/77twNozlwaB84rqj4za7l95A7W\n/uSOwo/bF1P7kg4BdkfEbyXtD5wGXLLndkPzhksqz8xesnzZSSxfdtLk+ic++/lCjltm1yPpH4HP\nAYdExLTXmzt1RvOAb0hqMDGk+2ZE3FJMmWaWgrLCSNJ8JhqYX2XZvtPU/v3AogLqMrNElXhx+gvA\nh4Drs2zsO7DNaq6M+4wknQE8EhH3Scq0j8PIrObaDdN+99gD/O6xB9vuN83XQS4GLgLeNnXzTnU4\njMxqrt0wbd+5x7Lv3GMn1397z3+97Pftvg4i6TjgKODnra7oCOAeSYsj4sl2dTiMzGquufvFQo8X\nEQ8Ah760LmkL8KZeZ9PMbMBFc7z0U2TZyGFkVnMxXm4YRcTrsmznMDKruRnojDJxGJnVnMPIzJLg\nMDKzJBQ9m5aXw8is5prujMwsBR6mmVkSHEZmloSy7zPKymFkVnPujMwsCQ4jM0tCc/euqksAHEZm\ntefOyMySkEoYdfxumpkNtmZzPNPSDUn/KukRSaOtZWWnffoujJrPPl51CaUau/fOqkso1U/Wra26\nhFLdPlL8d83KFuPjmZZuDwt8ISKGW8sPO+3Qd2EUO7ZVXUKpNo+ur7qEUt0xsq7qEkpVxkcWyxbN\n8UxLDtnexN/Sd2FkZsUqMYwukPRzSaslHdRpY1/ANqu5vE/td/g6yL8BH2utfxz4PHDetMfr9ZtJ\nksr7Nq6ZTSsiuhoK7anbv98855N0JHBDRBw/3XY9d0a9/ssws+qU9fcraV5EvDTb9G7g/k77eJhm\nZmX4jKSFTMyqbQHe12mHnodpZmZF6KvZNEkrJT0kaUzSh6uup0iSvi7pCUkd29l+I2m+pFslPSjp\nAUn/UHVNRZK0n6T1kjZI2ijp01XX1I/6pjOSNAT8Angr8ChwF3BmRGyqtLCCSDoF2AFc3elCX7+R\nNBeYGxEbJM0B7gH+fFD+vwOQdEBE7JS0DzAC/FNEjFRdVz/pp85oMbA5Ih6OiF3ANcAZFddUmIhY\nB2yvuo4yRMS2iNjQ+nkHsAk4rNqqihURO1s/zgaGgGk/5Wyv1E9hdDiwdcr6I63/zfpIa5p3GBio\nW80lNSRtAJ4Abo2IjVXX1G/6KYz6YzxpbbWGaN8FVrU6pIEREc2IWAgcAbxZ0oqKS+o7/RRGjwLz\np6zPZ6I7sj4gaRbwPeBbEXFd1fWUJSKeAW4ETqy6ln7TT2F0N7BA0pGSZgPvAf674posA0kCVgMb\nI+LSquspmqRDXnr2StL+wGnAaLVV9Z++CaOI2A2cD9wEbAT+c8BmY74D3AG8XtJWSedWXVOBTgbO\nAk7t5v02fWQe8OPWNaP1TDz6cEvFNfWdvpnaN7PB1jedkZkNNoeRmSXBYWRmSXAYmVkSHEZmlgSH\nkZklwWFkZklwGJlZEv4fhZ5DadgtaIUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3aa2d820d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(plot_T_logit, vmax=[1, 7]);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}