{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GLM,  fixed versus random effects in fMRI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Remind what is the issue with random and fixed effects\n",
    "2. Construct some data set from one bold series\n",
    "    1. Extract a time series   \n",
    "    2. Add some signal (3 conditions)\n",
    "3. Solve the GLM for one subject\n",
    "4. Solve the GLM for 7 subjects\n",
    "    1. Fixed\n",
    "    2. Random\n",
    "5. Testing. What is the distribution of our contrast of beta values ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - \n",
    "from __future__ import print_function, division\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# - set gray colormap and nearest neighbor interpolation by default\n",
    "plt.rcParams['image.cmap'] = 'gray'\n",
    "plt.rcParams['image.interpolation'] = 'nearest'\n",
    "\n",
    "# - import numpy.linalg with a shorter name\n",
    "import numpy.linalg as npl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - \n",
    "np.set_printoptions(precision=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. The issue with random and fixed effect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See http://www.fil.ion.ucl.ac.uk/spm/doc/books/hbf2/pdfs/Ch12.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Construct some data from a BOLD time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - nibabel\n",
    "import nibabel as nib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - change this to the path of this image on your system\n",
    "my_path_to_file = \"/home/jb/code/pna2015/pna2015/day10/ds114_sub009_t2r1.nii\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(64, 64, 30, 173)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7efeb592f950>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmUZUd1p/ttKbMya55UKg2lklCpBEhCDJowNlBgBrUH\nhkUvAw8EBrrdBjd+pllgwO8ZuR/LGIOx8WuDXxuLbjCooW2sBpvZptxgWUjIAgok2aUBpNJQg2pS\njVlS7fdHROSJe26c4Z47593fWrky77nnnht5h/jF/u0dEaKqGIZhGJPJKcNugGEYhjE8TAQMwzAm\nGBMBwzCMCcZEwDAMY4IxETAMw5hgTAQMwzAmmFIREJHrRGSniGyLjv0PEbnN/9wrIrdF971bRLaL\nyJ0i8qLo+GUiss3f95H+/CuGYRhGp1RFAp8Aro4PqOqrVPXpqvp04K/8DyJyEfBK4CL/mI+KiPiH\nfQx4k6puBjaLSMs1DcMwjOFQKgKq+i1gX+o+38H/EnC9P/RS4HpVPaGqPwbuAq4SkTOB5ap6sz/v\nk8DLetB2wzAMo0u6yQk8G9ipqnf722cBO6L7dwBnJ44/4I8bhmEYQ6YbEXg18JleNcQwDMMYPFNN\nHiQiU8DLgWdEhx8Azolub8BFAA/4v+PjDxRc1xYyMgzDaICqSvVZ6QeW/gDnAdtyx64Gvpk7dhHw\nPWAR8ATgbkD8fd8BrgIE+BJwdcFzaVV7Bv0DXDvsNlibFla7rE3Wpj60S5s+tqpE9HrgRuBCEblf\nRN7g73olWUI4iMntwOeA24EvA29R3zrgLcDHge3AXar6lWp5MgzDMPpNqR2kqq8uOP6GguO/C/xu\n4vitwFOaNNAwDMPoHzZjuJqtw25Agq3DbkCCrcNuQAFbh92ABFuH3YAEW4fdgARbh92ABFuH3YBe\nI5ljM3xERLVpcsMwDGNC6abvtEjAMAxjgjERMAzDmGBMBAzDMCYYEwHDMIwJxkTAMAxjgjERMAzD\nmGBMBAzDMCYYEwHDMIwJxkTAMAzDI8JaEVYPux2DxETAMAwj423Abwy7EYOk0X4ChmEYC5SlwLJh\nN2KQWCRgGIaRMQOcMexGDBITAcMwjIxZYP2wGzFITAQMwzAyLBIwDMOYYGaB9SJMzJL2JgKGYRgZ\ns7h90lcNuyGDwkTAMAwjY8b/nhhLyETAMAwjYxZ4DBMBwzCMiWQGuJ8JqhAyETAMw8iYBX6CRQKG\nYRgTiYmAYRjGBDODEwGzgwzDMCaQWeDHWCRgGIYxkZgdZBiGMcFMnB0kqjrsNswjIqqqEzNd2zCM\n0UGEU4DHcUJwGJhV5fHhtqoe3fSdFgkYhmE4FgHHVZkDDgBrh9yegWAiYBiG4ZgFjvu/H2ZCLCET\nAcMwDMcscMz/vZMJSQ6bCBiGYThmyCKBPcBpQ2zLwDARMAzDcMSRwBwwPcS2DAwTAcMwGiPCb4nw\nsmG3o0fMkInACSZEBKaG3QDDMMaaC4FDw25Ej4gTw48xISJgkYBhGN0w638WArEddIIJGSSXioCI\nXCciO0VkW+74W0XkDhH5oYh8wB+bFZHrReQHInK7iLwrOv8yEdkmIttF5CP9+VcMwxgCs2S7cY07\ncWJ4YuygqkjgE8DV8QEReR7wEuBSVb0E+JC/61UAqnopcBnwH0Rko7/vY8CbVHUzsFlEWq5pGMbY\nspAjARMBVf0WsC93+M3A+1X1hD9ntz/+ELBURE4FluKy6wdF5Exguare7M/7JCyYRJJhTDomAmNO\nk5zAZuA5InKTiGwVkcsBVPWrwEGcGPwY+KCq7gfOBnZEj3/AHzMMY/xZzITZQSKcIsLG1H3jSJPE\nxxSwWlWfKSJXAJ8DzheR1+I+EGcCa4BvicjfdXpxEbk2urlVVbc2aKNhGINhIUcCiwvOuxz4CPBT\ng2hUChHZAmzpxbWaiMAO4PMAqnqLiJwUkdOAZwF/raqPA7tF5B9xuYFvAxuix2/ARQNJVPXaBm0y\nDGM4LLTEcBCBshLRpf5naPjB8dZwW0Te2/RaTeygG4Dn+ye+EJhW1T3AndHxpcAzgTtV9WFcbuAq\nERHgGn8NwzDGn4UWCcR2UNEgeYaF8z9XloheD9wIXCgi94vIG4DrcPbPNuB64PX+9P8PWOSP3wxc\np6o/9Pe9Bfg4sB24S1W/0vt/xTCMIbDQRKBOYnhBiUCpHaSqry6465rEuceB1xZc51bgKR23zjCM\nUWfiEsP+vIXyP9uMYcMwusIigTHHRMAwjEaIMIVzExbKqLiuCCykZLiJgGEYjQkd4UIZFXdkB4mw\nIPZDNxEwRhoRniQyGTM3x5DZ3O9xJ44EykpEZ3K/xxoTAWPU+TPcHBRj9AiTqRZEZ0j7fgJlJaLx\n77HGRMAYdZawQL5sC5BZ4DALKxKoaweF88ceEwFj1JllQhbyGkNmgf0sHJHupDoonD/2mAgYo84s\nsGjYjTCSzAIHWCCdIZ0lhuPfY42JgDHqWCTQY0TYIMKf9uBSQQRaKmVEmBLhCz24fhsi/D8iXFbz\n3H8vwis6uLxFAoYxgpgI9J6zgZ/pwXUW43ICc7RGa+uAXxTpS/9yBfDEmue+EHhSB9euu9F86Pwt\nEjCMAbAYs4N6zTS96cDCyPk4raPiM/zvfrxv08DymudeELdBhLeIsLrk/PxG81XVQRYJGEY/8RaD\nRQK9p9cicCx3vaGLgP/sbKL1s/NrwM+WPMzsIMMYMaYBwUSg1/RDBFKRQMtz+FzEhV0+5zSwosZ5\np/nz4s/OIuDKksd0lRgW4WIR1ka3nz0Os4pNBIxRJnQsZgf1ln7bQev97/z79kHgrV0+5yLq2UGb\n/O+4I5/G5RSK6CQSSM2P+G3g5dHtv8HttDjSmAgYo0z4klkk0Fum6I2wLgaOUmwHzR8T4QzgFT14\n3ro5gQv87/j5poHLShLWnYjAAdqFdAWwEkCEU/3toi0qRwYTAWOUCV8gE4HeMozE8L8DjtAbEahj\nB20CdtNuBynF1UWd2EGp+RHzIkAmVEtqtHWomAgYo4zZQf1hGjjFLwXdDVWJ4RmYX3L6P+DWger2\nvaxrB10A3EG7HXQTibyAb6PgqoKgexEIv00EDKMLzA6qiQgvFOF9NU8Pr2dpNCDCN0RYV3JKUWJ4\nPa2j/mcDO3Hbzg7KDtpEWgT+kXReYAY4por62y0loiL8vch8h15kB60EVkV/g4mAYXSFiUB9LqT+\nxKhKERBhEfB84Kkl15nF5QRSdtB90fXX+NsnGJwdlIoEFlEsAvEcAWiPBH4K5ucYzAIHKY8EghhY\nTsCYHERYL8I7enhJs4Pqs5L6HU5SBER4jQiX+5sbcPbIxSXXWUzODhJhFlgKPEz2voWZuPmZxU2o\ntINEWOHbcH94Pl+qOYWLRp6S2KMiTgpDJAL+sTP+mpCwg/w5ZgcZE88FwOt6eD1LDNdnJfUnLxVF\nAtcAV/u/N/rfZSKQSgyvx1k/x2mdVHUcJwLdvpd17KBNwN2555sCHlflUVx556rcY+KkMLRGAiFf\nUGYHLfXnmAgYE80U9af012Ei7CARThHhF7u8TCciELzuvAhcQFZaeS6wA7io5DqpxPB6XBQQj/pD\nB9uLSKCOHbQJuIvWjnyRvw3OyslfoyUSUOVxmC/1DP9bHAnk7aC8DRRumx1kTBT9EoGFbgedBfz3\nLq+xii4iAW+PnEsmAhuBLwMXl8x6TSWGz6A9Eui1HXSqSGlS+2ycgMUiMO2fH2qIgCc8PiUC+Uhg\nBbAHiwSMCWcaWNHDqfKzlO/1ulBYRNbBNKWJHRR3yBtxnWSYaXsucCuuIyya9ZpKDJ9BeyQQ20Ed\niYAIS0RY4/8+Bddn7ad8sHEarkOO7aBpWiOBlbnH5O0gaBeBJb4N08CjtL7eK3A5iPD5X4X77JoI\nGBPFlP/p1RK7i3Ff2IUuAjPAIl+R05RuE8MXAN8BVoqwDCcKPwFup9gSaksMk84JdGMH/TLwO1G7\n50iP5GPWAo/QWo2UF4E6kUAoE40jgUW+DUdpF4HwnEtw78dOzA4yJozQufTKEprFjbgWuh0U/r9u\nooFuE8ObgO3APcD5uEjgPuBHFCeHi+ygVE6gqR20hux1Cb7+o5R/xmIRiKOeYAcdoF0E4r0EAik7\nKAhaLHLgXv+D/tqr/O2HsUjAmDBCwrFOHXcdQj32JEQCAMu6uEZXOQFcJHAXrqpmMy4SqCsCcYcY\nRKCoOqhTEVhO+2i+SgSCHZTPCVRFAkV2UHhdl5CJRX6C3Ap/3f04AVgJPIiJgDFhBBHoZSQwCSIQ\nOrluRKAXkcDdOCF4JnBElUM4ESiyg8oSw72KBJZH7WxiB9XNCdRNDAexyC+VEUTgANnM4YcwETCa\nIMI7RWpvoTdKmB3UjK4iAb/uzWJgtmZSPlUiGkcCz8flA6A8J5BKDMc5gV6UiK6IHlPXDkolhmM7\nqNNIYD4xTKsdlI8EDpCJQLCDLCdgNOLFwJOH3YgG9NoOWoz7wk9KJNA0J7AC9zrVXZahJRLwFS9P\nIIsEno6zgsB1pssLktapxPDpZJFAr+ygOBI4gevE6+YEUonhVE5gDbA3d6wsJ1BkB8U5AYsEhokI\nK0R4xbDb0ZBFjMGHJ8FA7CBfNvjKHj3HKNBtTmAlzovOd0yI8G9F2mbH5u2gs4ADqhzGCYHgIwG/\noNo+SO7N22IH+SUjZnEdYT4S6MYOijvyOZzgJQcaYZE3VY7QWU5gHbArd6xKBFKJYcsJjBCXAO8Z\ndiMaMsMYfHgS9MMOOkh7x/FEsrLBhUC3OYFVuI63TQSA/xt4Tu7YNG7phPC8YYYtuM7/MbJIANwI\neU3iefOJ4XXAbi8ccSQQOs4TuFLYTuaRrIiuU8cOClYQFFcHpUTgdNz+AzH5EtEqOyhEAmv9fXsw\nO2io9LJefdAsYgw+PAn6UR2UsoOmGFGRFEFEStfbSdGLSKBIBNaQTQALTAOHoue9ABcBoMpjOCH4\nSXT+XnKRQDRpai563nVkHWkcCcwCx/1SDAqc2sH/lqoOKrODghUEnSWG47YH4kjgKNV2UMgJbPTP\ncZgR/ZzGLGQROJX61RKjxrhGAoOqDpqm+xm2/WIzcGOHo91uI4FYBPKDhzVkS0EE8iJwFm6ZhcBH\ngH+Kbu+jPRKI198PHWI8ms5HAsei453keFLVQYV2EC4SCCJQNGM4lROosoP2UT5PII4ENvrfRxiD\n7/FCFoFxjwRG/sOTYBp4nN6JQJgxnLeDphnC6+NH+VWj2A24DmFDB5fuhQjsJzeL1a+vs4S0CBym\n1es+FO5U5f9V5YHo/JQdFJLC0GoHhY40tWxE/ngdOq0OWkurHRQ/tlM7KBaBvbTaQflIIM4JbCR7\nP0Y+oi8VARG5TkR2isi23PG3isgdIvJDEflAdPxSEfknf/wHIuLX8pbLRGSbiGwXkY/0519pY5xF\nYJwjgX30vkQ0FQnMlmwY3i9eD3yw4pyz/e+y1TfzdGsHFeUEgoWTt4OmaI0EluJGrUWkRCCuqy+y\ng4oigVoi4N/fMPqGDu2g3EqgdRLDRSIwi3sN6s4TOJcFFAl8gmx9cQBE5HnAS4BLVfUS4EP++BTw\nKeBX/PHnku3X+THgTaq6GdgsIi3X7BNTjK8dNM45gb30f8ZwuD3o1+g03Be8jCACneQFFuE6lqYW\nV1FOYA1wL3BObhOVfCSwxN8uokoEQpL0dNKRQLw4WyeRQBDF2tVBtNpBkHXkhTmBXFUTicfm7aBj\nuL7t1Gif5jgnsIyFIgKq+i3cPx/zZuD9qnrCnxPU80XAD1R1mz++T1VPisiZwHJVvdmf90ngZb36\nB0qwSGDwTOM6jH5PFgsd2qBfozDaLeMs4E46E4EZXMfV68TwGlyt+kNkm8RAe05gKdUikC8RDRPF\nIBsVF0UCTe2gFbRbOp3YQZB15PN2kKoTr2g56riqicRj2+ygKBcy4/M/y327gpCEMtnpGhbiUGkS\nTm8GniMiN4nIVhG5PDquIvIVEblVRMI2g2Ft78ADZKOlfjJF5+Voo8K45gQGaQfBcETg9Ipzzga+\nTmd20CJcJ9OLeQJxdLTaX/cuWvMCKREos4NSieEqOyi1bET+eBXLceLYiR2UjwRCcjiOBKDVEkpZ\nQdBaIhrsoDiqCeK3DDjq7af9/r79XihGPi/QRASmgNWq+kzgHcDn/PFp4GeA/8P/frmIPB/a1HVQ\nFO2eNNJE+5mO9AengBY7SISfFuFnurhe0VLS4b0ddIVQ3UjgG8BFHQxAuo0EQk4gv7zxGlwHXiUC\nTeygVGI4toNSS0lD5yKwh87soFQksIhqEchXBoXHxpFAXgSCDRbyAdAaCcAYWEJT1ae0sQP4PICq\n3iIiJ0XkNNyGCv9bVfcCiMiXgGcAf0FrpcQGaKk8aEFEro1ublXVrQ3aCK0ikF8YapQ5ldb9TMeJ\nvB30KtwX6dsNrxcigSkRJArXhxkJrBFhytfTpzgb+AGuQ96A+15UESKBKoEposwO2ov7vsXJ4SAC\na/3tOnZQN5FAN3bQHuASL6h17aBUTiCuDoJ6kUAsAo/ivpfLaI8EQu4K2sWgLyIgIluALb24VhMR\nuAG3wNQ/iMiFwCJV3SMiXwPeKSKLcS/ec4EPq+rDInJQRK4CbsZtZv3HRRdX1WsbtClF8OFSCZ9R\nJh6djRtTuC9T+IKeS3tOqROC7xx2Fwtf4mGJQIjOTsMtDtaCr2Y5A+fBh9U364hAiASe0LBdVSJw\nF/Ds6PgU7SWiVdVBqwFEWOnbmc8J9Lw6CPc5OoArO54iG80fxXnt06oto3uolxiG1uRwHMHExCJw\nHPearSazfOJI4AC4iiQRDpH1OX2xg/zgeGu4LSLvbXqtqhLR64EbgQtF5H4ReQNwHXC+Lxu9Hnid\nb9Q+4MPALcBtwK2q+mV/qbcAH8dtWnGXqn6laYM7YCztILIvyDjaQSESCCOsjWSjzSaE0WZ+gtEw\nIwEoHrGvw63Bcxy3+mbd5HCIBPpRHbQXNxu4KidQFgnEOYGXA//LPyZ07I+RTc4Mnd8cLic3BUgU\nOXUqAgfJBGUamPMRYVE0UJQYzotAPGGsTiQQRGANrZFA3g4CJxJBKMbbDlLVVxfcdU3B+Z8GPp04\nfivwlI5b1x3jLgIj/cEpYIrMvgkblzey4kLpnSqPicz7uqGjGlUROJvM6vwRcEXN6wbPuZucQFli\n+G7gCSKcospJOs8J7ANW+UjnApy4v8I/H6qoCMeAfZFlFzru2EOHzu2gR6PHBDsIMhGYX/nTV/sE\n6yZ+vjp20D2J58+LwBHaRWCGbKJYIJSKwhiIwEKfMQzjN1dg3O2g4NmejeucmkYCi8nshngNGBiu\nCBymuEIoFoHvAs+pOaGtF9VBhYlhvzroYTLx6qg6yI/iD+M6zU3A14B/S6vAH6fVUgkddz4nV3e5\na8jKLuNIIIjAftrLVtcCe3OlnnUSw1V2UMhp5COBVGIYXHXYdv/3gqwOGhfGORI4yHiKwDTOGjiI\ns0J20VwE4sTjqNhBi3Eef1EkcBZu+WCAW3Ed68/WuO589UmnJc1+nf9p/1xFdhD+/nDf/GQx/3xV\nkQBkyeELgPf723HnfoxWS6UXkUCwg8Jj4rzQg7jXOyZvBUG9nECdEtEqO2g+76jK21S519+0SGCI\njKsIzOBGOSM9eiggjgQuwVXJrIhmVXZCLAL50eMwS0TLRGA+EvCj0T/B5cOqWITrLOao8b6L8HWR\n+f99JS4PES/kFohFIE7UxpHADHAiLLFQQkgOb8JZXX9KawI2LwKh485v29jJAnKxHTRDqx30AGkR\nKNoYJm8HNckJHMG9BnFV1AzudS4qPjERGCLjagctwn2gFo/hRLcp3OgpiMC9pMP2OuRFYBQigSAC\ndewggM/gLKGNBecHQidziApLyH8mngec6Q+FfACk1w7am7gvXjaiqjIosA8XBZyKG22/l9Y9Hcrs\noG4igWAH5S2dB2mfdNqyEJ6nLBKoawcFS6vIDnoimf2T5wgjPqAbWREQ4X1dTrce50jgKO7LEq8I\nuUyE/2torapHbAddgluT/hFc2V6nxCWIKTvoBMMRgfuoZwfhN2r/C+CNFdcNo9RKEcB1KKeSCVHI\nB0CUGPa5iDCTONwXvgvxAnJVlUGBvbhE912qqConfZI5UGQH5ffu7TQxHOyg+eogf19q5YEltAta\nbCW1iUCUTD5IO6nqoCnaE8MX4aKjFEexSKBzvM/5W3S3/EAQkHETgdAh5EcQ5wBvG0qL6hPbQU/G\ndZh7aJYXiGek5u2gaVznNiwRqBsJgCuXPq/iumEBucNUW1zhO7Eu+h063zgxvBI4FNk88U5YwQ4K\ny5N0IgJ3F9x/PGpHWMFTaS0lhe4igSo7KC4mCFRNFlsH7EmsGxQ/NraDoDUSWI6LkO4s+B/MDmpI\nGDl204GPqx0Uf+DiD88i/GzVobSqHrEdNEMWCTQRgSo76AA9+HKJcGkNuyZQJzGcF4HURi95Zqgf\nCQQRCEJ0BtnEtdjyyW+cnreD4kigjh20F7iMbBvKPMdot1SO4zrabu2gOBIos4NSkUDZPIGVuNcv\nZQXFj40jAWiNBC4GHlBtE5+AiUBDwpesFyIwrpFAPowM/0c3k6/6TbCDQp12iASa2kFl1UE9EQHg\n13HljnXblLSDvHW5mtZkKbSXbaboxA6Kk5kA64Gd/u8qEZiJLNajdG4HLaU4EvgQ7cuDzNGdCAQ7\nKM4JdGoHxSKQjwROA/4A+ELB89cRgWdQbAWBlYg2ppciMM6RQPzhCV+cqlUsh0mwgw4CJ3Ff1MJI\nQKT0f6mqDjqAt05EWOqXM2jCGuqLySxuBLoytz4/+EqWnE8O6X1/89RODNNuB+UjgfCZiZPCkNlB\nYUQcPPu6dlBY/iMZCajyRdU2AQx2Sbd2UGqy2C7cBLb4WospFoH4sZCVMe+nNcEdky8RTdlBT6Vc\nBCwSaEjoHCY5Esh/eOK1z0eV2A560K/rkkwM+077ARF+oeBacWK4yg76NeA9DdtcSwSi1V2P4DrX\nvLDFCdqYOiLQSSTQiR0Ur9sUkphBBB7Dff9XUN8OguJIIEW3kUDhZDEvtjvJqqTAvY95W6YoMfwg\n8OfAaxPCHUhNFoPWSGAWE4G+EDq6bkbxIYs/biIQRh1FdtAoi0BsB/3EHytKDK/EfUGvE+GSxP1x\nYjhlB8UT6tbTfLZt3UhgEa6e/iQuAZp/H7oRgU4igRW4zj22g+rkBFoiAZ8IncNFDHXtoGNE1U81\naJwT8KK7jNZIIG/p5C2h2naQKkdU+XeqLUtM5KljB4FbJ6oIKxFtSK/soMOMgR0kwmtF+C1/s6g6\naCTtIBFuiyYuBTvoJ7iqGCi2g5b7896Fm3iUp8wOylcHnUbzL1pdEYhFaRft70OZCFS1LY4E5quD\nRHiXCP8xd+5y3Do3sR0UcgJx/qEoMRyiNXCd2RrqicB9wE0lo+YUIRJoYgctBY75KqMwmMtbOvkK\noTIRyD+2DnWqg05SXBkEViLamF7ZQfFyuaPMOWRlhKFSZOTtIG/pPI1s9DoFPKbK36ryVn+sKDG8\nHNfpfQ04P3F/J9VBa2ku9qup9yWN25OKBOJJWzF1EsNFkcALcJ5zzHKcJRO+I0WRQD4nkLeD8M+5\nmhp2kCr3qvK8qvNydFMdFKyg+DEpSyeOBMpKRPOPrcMJf83HvRilIoF7SiqDwOygxvQiEjiV1kWy\nOkaEq0V4eRdtqMsSsg9KqBkfBzsolFbGyxHkN1spigRCqL8LWJeYGFg1WSwWgUaRgN9gPH7ty4hF\n4GHgfSJ8S2R+ddxGdlC0WcocrpNZ5o+fgqvLz5evrsBHAr79y8i8/zjqWEtrTiCfGA7H6kYCTUhF\nAnUXkIsXZWtZSjo6p7EdVJMTtG4iE16n8P8cArZVXMPsoIasw41ihm0HPQu62h6xLovJOqKiSCAs\nLDdKdtC5/nf4Ugc7KKbMDjqkyhyu88yfU2UH5SOBJl+0sJxFnTWIYlH6HeCXce/R0/yxpjmBKdxI\n8yStkcBm3P90bu785TjhPObP2RlZNPFzxQnjcF9RJNAvEQjVQb2KBJrYQUVRRB3yIpC3g/4a+JWK\na5gd1JDTcZNyhm0HLaYPKi7CrAivzD1PPhLIjyBmcFt7jlsksBdYnVhSOUQC4Dzt9bn78zOG40ig\npUSU5nZQ2Cilo0hAlUdU+TbOCw5WV9OcQDyTNRaBK3BLEm/MrSEVOsfduD06dkb3HQNm/fn52ctF\nkUAtO6ghRdVBdRaQi0UgtZQ0pO2gXuYEHiMdCcwBqHJMtW3V0jxmBzVkHa7D61YEurKD6JMIAJcD\nv5t7nlgEinICoyYCqUigRQR8mehhaKvjj7/kD+NGrjFVk8UO4hbZm8Z1ZE0jgeN0nhgOxFFOUU6g\nKhKIF1iLReBK4Ju4kWScUwkrawYRmB/t+3X/T+Leh7wIhHYM2w6qGwmcRfa/dVIdVJYTaGIHLaVV\nBOYKlpgowkSgU/zkj2W4D8Cw7aBZ+iMCF9D6wVhCNqotKhFdhBOBUbKD8pFAyg6CdHJ4GdmKj1Ui\nkLKDjuO+1GdG53fKGlxH0mlOIBCXvxZFAieAU0qW+yiLBG7BVVDFeYGwxv4uciLgOYYTt6W0L/Uc\nRtNNqoOa0E1i+AKyiWllkcBZUaTUybIRdThBq511hNb/pQ6WE2hA2Cg6XvWwCb2yg/pRYrqJVh+6\nKBLI20EP4Warjsr6QefiOpTwpU7ZQZDOC+QjgbwdFI/qUtVBIcLYiNuIvMkXbQ3OdmwqAvFEuKQI\nROv8F30O40jgX4CnivAs4FLgn3GlmXFeoMwOwj/X+bjJevGINdhBsVAPyw6qIwKbyEQgmRNQ5SCu\nACS8f/0oEY3toF3AH3V4jUO4iG5kGUURCKsixqseNqEXkUC/7KALgCXRCCafGE4tIBeWmN7H6Kwf\ntBG3Z0BZYhhcZ3muyPy6N9AaCeykPRKIcwZFS0kfwZXXPkTzSGAH9UUgbzXUiQSg3BKajwRUuQ/4\n98CXgXv9RKZ8JBDbQRtJRwIXkF7ILj+ansO9zv2OBJrYQReQzU5OLSUdiDeHSeUEiqykOoQI9DiA\nKnOq/HYZdpBhAAAgAElEQVQnF1DlqCq/1OHzDpRRFYFQ/TAKkUC/ROBUso4tFQmk7KCwccfQLSFv\n252OE4HYDkpFAj/ATQjbK8Iz/LGqnEB8f8oOCiKwEdeRN80J7KB+dVBZJFCUE4Dy5PB8JwOgyg24\n6qMb/KFUJBDsIEiLwPm0z+wtyglAf3MCp9B9JJDaVCYQi0BVTqBJJADt7/uCYhRF4HSySKDbeQKj\nKgKbcAm8OIxNRQJ5O2iO9ESlYXA2bgR+BFjko5qkHaTKu1RZDXyVrKSvKicQ319kB4VIYAfNI4EH\ncRU1Vd+FqsRw00ggvK/zqPJh1fkNhFI5gRAJQLsIHMV9vvKRQFF1EPTPDjqe+w01RECE5WR5wfCY\n1IxhaN0hrB92EHSeBxgrRlEEYjuoF9VBI2EHifAJEU4XYQ1OoHaRdfyLcZUup1BeHRQ27qglAiKs\nFOF/1ujgmnAuroMKX9BToG23qTxhs3KozgnkI4F8iWgvIoE1OEtn3noU4c9Ekq9vUSSw1gtgmQiU\nzRrOb3aSJx8JxHYQpHMCKREomicA/Y0EwnPHx9pEQIQpEW7wO31tAu6OchqppaQDYYewKdz3Kn9/\nN9VBcQJ9wTKqIrCL3ojAKEUCvwj8HFmYe5hWEQDXUZRVB3VqB63HrZX/gm4aXsBGXAcVvqBFVlBM\nLAJVOYHl0f35jiNEHEEEdoLrSDr8H8L6OrHgvoT2pRogIQKqHMFFdEvpLhIo62TmIwHfQapqy36+\nKTtoE+12UBwJ5Du3fuYE4t9QHAlcDLwU+Gla8wHxY1KWTogEFgNHE+WbJ0oeW4VFAkOi15HAqIjA\nCuDFZB/wI2RedGwLlVUHtdlBImwR4ekFzxkvtdxr8pFAUWVQTFEksAe3Nnw82o8Tw1V20B6abd4R\n1teJBXklrhPNk0oMh7afgXsNijrTNhEQ4U0+UV4VCezGjXSX0Pqa7fbtya+Cecy3p5NIoJ/VQeG5\n42MpEbgCtx3li2nNB0DxAnLgRGAlaSsoPF9Y/6eT+n4wERgap9PbSGDo8wT8CG4KeCFwIe4DHo8+\nF5PNgC2rDkrZQa/BRRkplgLfB35GpG35gW5JRQJVI629ZEs1zI/0/eJce2j9v8rsoLhEdA1ZSXGn\n73VLJODX4pnBCXWeVE4A/9znAwdKOplUYvjdwCVURALeXrsfJ3bBCgJngb0o8ZyhjUWJ4XyJKIxG\nJHAF8FmygdJdiceUVQelKoPA/a9LEo+rg4nAkHgyvYsERsUOWo4r7dwJ/BJZJBCLwCO0RgJ17aDT\nKN6+cQnutfwU8Ks9+D9i8pFAHTtoH612UDyKnU8Oh92i/LpCkLaDQiQA7rUrjAREuKqgPWHjlRCV\nhVnNKRFI5QTAidf5FFtBkBaoNf6nKhIA9zqfS1YZhCrql67IE6KVuonh49Fm9L0m/F9xJ1q0gNyV\nwH/B/Z8/TasdVDRZDDI7qCgSCLN+O7WCwmPD8y9YRlEEvg38kNERgZnc2i1NCCO4rwIXEUUC3gIR\nsgXRypaSTlUHraV43kD4Yvwver8Q3uk4UYuTdk3tIGjNC+TvK7ODwHXEyUjAV5rclJujEPYEXoEr\n6wyvdRCBIjuoKBLYRLkItCSGfaJ+Fe61qMoJgFs19ALaX5cUx4D9Pl+RP56aJ9AvKyhcPzx3fKxF\nBERYDDwRuBX4e9xAMB8JzOISv3nBikUgZdeFSMBEoICREwG/288+hjxZzHf8i3EfhG5nDYcR3Ff9\n7bvJfOiwBnroiMoWkEvZQXVE4CD1auE7ISR240igjh20xneC+b1t4wqhOGkM5dVBUB4JBGG5KHd8\nJW5P4MdpFYHtwKaE8FeJQNEcAWgXqJU44V9NvUjgR779sR1U9lyp3b+K5gn0ywoK149/Q3oBuacB\nt6tyDPiKb9/9uessJb1uT1VOIEQCZgcVMHIiENGLeQJHgVMblkguwo1sD9O9JRRGcP8b+BLuSxo6\nnuBlhtthxN8rO+gw7VFFL1jqr92kOmgpcCRXThrPFciPeMuqg8J1q0Tg4tzxeNOVWATu9/9Xai2j\nosRwp3bQmuh3nUjgdlz75+2giufKW0GQtlT6LQIpOyiVE7gSt04SwN8C/zNnUYWZzalBRlVOYI7m\ndlD4PNtksSHRixnDJ2guJmGE3qTqJM8K3KjzqCo/7zu/4EOH5wmRQejsTwASVcwEcdiLXz/Ij1bX\nUB0JHKb3kcBSWiOBTuyg/Egf2u2gZCTg/+dTyUTggF+ptCgxXBQJxNsvxiJwAGdF5PMCZYnhKjso\nnxiORaBuJBBEoE4kkBKBOBKIS0T7aQfl9+KFtAiExfJQ5UFVXpO4TpEI9M0O8t/Tk1gkMDQKO28R\nflWE11Y8PoxMm9pKvRSB1AgujgRiO2iGLOyNLaFFZEm8sH7QCpytMC8CInxWZH5lzaXUFAEvKv/g\nK2RK8R1xKhKo+qLtx3W0KVsjbwcV5QSmcVtYqn/+sFJm0fu0HlfFlI8EQlIYMgFehevM76Y9L1CW\nGF5GfyOBh3Cv8SaqRWAPrUnVwGM48ZxluJHA47hVVeOd5C4DvltxnWWkxbJuYriJHRQebyIwJMpG\n8KHUsoxYBIYdCaRGcKHTDx/efE4AWi2huLMIltBp+NmyoaIGeBGwwf+dv3YZVwHPIVuauYxZ4IQX\npLh8rzQS8OvdH8YtOZGPBOLlpvORQDx6jO2MI/5xUB4J/B1pEYgjgVAdtJ90JFCWE4AOEsM4K+oQ\nNXMCXvB+hHuPquygPwDeX3CNY7RaIwPPCfh2zOcF/AS/JwD/WnGdGcojgbIS0VMKHlsHE4EhUtZ5\nX1xyXyCIQJP6ccg84J7ZQbljqcTwUlo7hbjzjteYCcnhtf7vsHxBGM0u9+cFETgGTCf28Y15sf+d\n98JThCgAOisRBdfxnkv76xEvxlYVCYQv9CEyESjLCXwHl5CON7ZZTRYJ5O2gTiKBIAKdJIbX4ISm\nbiQATgSuoCISUOVkScnnMVptlaO0i3EvmcNFbfn2xKJ+DrDLJ4XLrgPFIlCVGC56bB0WvAiMyrr0\nKZIi4CdeXZC6L8eoRQKd2EHxTM5g48QRQhCBGbL5BWvJNiRpEQFVVGT+WkUjyRf7a3UqAp3YQeBE\nYCPpSCDYWoU5AVpF4MvAbf7vIrFfj0vE34ErPbzJH19Fuwgsw9XkF0UCRYlhqLaD4oquIAJPxb12\ndTri22kXx045jnttw3pDfwPc3MX16jxfqgONRSA/O7joOuFxeUJiuCgnMJf73SkLXgTGMRK4EOdt\n1hGBx0uuU0VIBPbTDlpK9uENt/Mj3fyOY9BqB+0hG0WHrfaCGMRlmPNRhQhPEuG80BAR1gJPAr5I\n+0JuKUJSGDpLDIPreDeSjgSCCOQ7u7jTmBcbVY6oco8/XhYJ7CRLrgZWkAliKhLoJDEMnecEQiTQ\nspR0CT/yv6vsoDKO4T6L4fU7rMq9XVyvijnSr1n8fubXCSq6DqQHGY9SnRMoemwdJlsEROQ6Edkp\nIttyx98qIneIyA9F5AO5+zaKyCEReXt07DIR2SYi20XkIzXbVtR5hy9yv+2gMEJv+viYlB2ULxE9\njBudnohKJ+OEbpEd9AjZKDos05y3g/LX+lVoSay/AFe+ej/1IoF4I5JeRQJHcNVQ+TVyIJEYTly3\nLCfwMFmZZSAlAmFPgD3AVM4+KrKDQoVUpyJwj3++WeqNUm/3v7uNBJZRT6x7Qa8igcKO3FeGzeE+\n//0QgeAmLFiqIoFPAFfHB0TkebiVFi9V1UuAD+Ue82FcrW/Mx4A3qepmYLOIXE01RVU9F+G+QFUi\nEMoIx8kOWkVrhxB33Ck7KGzFWRQJhOqg/LWW05oovho3kS21rn+KnucEfMIwRANlk8WKVoNse598\nFVOY2VwWCYT8zEqyNYDySzgnRSBqd1lOIJUY3oV7b06jngg8iBOabkSgJRIYAHNUi0BlJBAlk4te\np4O4z21RiWh4ziZMdiSgqt8i800Dbwber6o+pNSwrjki8jJcB317dOxMYLmqBu/xk8DLarStLBL4\n54L7/HMi9MYOGkR1UCwCq2n9wB0GlvkKCo0SbMEOWktmB63FicCjpCOBOMm8gtaS0cuBG2kmAp0s\nGwFZJJDqzML/UTZZrEgEUpHAauCwX3p5O63J3pW0RgL55aDzm7kURQIA18O8LZUiFQns9T9nUqOT\n8R3hf/X/R1PyieF+swP468TxvAhURQLhMUXtPoD73PbLDrLJYjk2A88RkZtEZKuIXA4gIsuAdwLX\n5s4/G/dhCDxANmIt4wTp2b4X45KBZR17vMFJt3ZQv6uD4jLOfCQQcgL53afK7KB/oTUnUBQJxCKw\nHDeS3UnznEAndtAM6WRoKBOtmxiOSb1PIR8QnndVdN8Ksg4/nxOAdCSQGmmiyttV5xPEKVKTxYII\nnEHNUaoq71RNLglRl7JJVz1Hlb2qvD1x1wmy3ejOpzonANkEyhQhEiiaMRyeswmTHQkUMAWsVtVn\nAu8APuePXwv8oaoega4XXAsjn5ZRvK8MOg+3wFzZFnWxNdE0EigsEfWTqtYkH5Wmrh2UigRC2Wh8\nPLaD8onhf6E6J5C3g8LIvm4kUJQTqJsYhvJIoG6JaExK7EM+APxEtWhNoLLEMLRu5hLWkWraGaQi\ngX10EAn0iEHbQUWESOBM4JBqLYurjh2UEoEQPTe1g3aQ7eK2IGlSIroD+DyAqt4iIidF5DTc+h+v\nEJHfx424TorIUX/uhujxG0hPawdARK7Nbn3jcfjZGbIR2IW4jc0fpbxj74UIlEUCl+I+dHvzDyqg\nrDpoMW7UG0QgFovQcedryYMdNEdriejZuAW4nuzPS1YH0W4HBRE4DJwhglRswJHKCXRiB0FnkUAd\nOyj1Pq3Hi4AqJ0TmR8KhoiRvB4UZw+AigbDD2DRuU5KmCdV5EYiW+ggisJzmHVSTdgzSDioiTBbb\nRL0oAKojgWSJqC+PPlHy2FJU+fkmj+s3IrIF2NKLazURgRuA5wP/ICIXAotUdQ9utmlo4HuBR1X1\no/72QRG5CleTfA3wx0UXV9Vrs+vwZlo78Itw+Yaqjj0WgV7YQfklF+Jlh+tQVR20G9eprqZ11HEY\n96XN20F7/fOfwHWcS3Cd55m4SODKqJ2ldpCfQDaLX9BNhMeoXqis23kC0DwSKIo4iiKBeA/efbiO\n/lHacwJrcRObgtjGOYGyfEAd4sRw2OXqqMh8VDSoSOA42cq4wySI+kbq5QPCY8pEAIrXQWosAqOK\nqm4Ftobbvs9tRFWJ6PW4hOGFInK/iLwBuA4435eNXg+8rsbzvAX4OC6pdZeqfqVm+/Kd/Qbcl7OO\nCIQwsKfzBHynuYhWf7mKMjsozgnklxBI2kHR+kHryaqDnogbPe8ms4Py1UEhEojtoCW4vVlDWWpq\nv988RTmBbiOBODHcq5xAvAfvfmC1H43nI4EzaS3zvI/eiUAsUPFyFeH3ICMBGFyJaBFBBOrMEQgc\np/h1inM7KUIZqZGgNBJQ1VcX3HVNxeN+J3f7VuApnTUNaO/Aw8iuk0ig9gJyIiwFXqjKDbgO5SDt\nnUv4OxkJiPBkYLmqm4npVwGdpj1UTS0bAcWJ4fxocRew1I8o9+BEaZt/zDLf0cUh8hFgadQBhqgg\nHtVDlhcoW8tlKdkkqU6rg8pyAntwSyOkIoGwamqnInBHdHs/7nWawRUOhNc0bD4Ui8CDwOl+Tabw\nHjUlTgzHC9cFERhkJADDHxXHIvCFDh5TFQkUvUdlj514RnnGMLR34MHjrRKBMEcAOluS+mnA7/u/\nYzsobkO8IXmKlwF/GiUgl+OSX3mPPSxstYJWEUglhvN2ELgRf6hI2Y/bpPtBshLRaVxHFz788bWm\nqBaBMuLEcPhCd2oH1Y4E/Gv3GMwLat3E8HrSdlAcBUD22s/X+nv//yFcnmWhRQLD7hDncNbnCyC5\nRWaKqpwATJAd1EvGQQTykUCbCIiwKVdK2jQxvJJsjZeixHCVCMwAT4f5fW2T/nq0VPRaiiOBkBNI\nLS0QFo6L7aEH8JEA7dPoQ+QRVw5Bq7UD9cpEu5ksdtT/L0U5gdNxHWZ+dctgCXVtB9GaDwiPhfZZ\nv2Fv336LwCCrg2D4HeIc8HrgH1X5SQePKasOArODGjGOIpCygz4LPDO63cgOwo0SV+UsgHyNdx0R\neBj4NX+7bEvAIAIhJxDaGyiqDgJnB8W16XtwIhAigbwIhAqYFTjBCJFAPKqHepFA48liXvz+kNYR\nevw/nIub4JWPnOIlq+tGAmfhRvOBZCTgR/2ppR9CXqCXieF4R7NgCw2qgxolO2gL8CcdPKYsEqiT\nExj2/zyyjLoI5K2cIjtoGa3r4OergzqJBMBV2hTNEwh/FyWGZ3DJ818UYR3lu0GFZQOO+tF8PvlV\nNFkMokjA8wjODgrCkbd5wvHluNcw5AcK7aDEPruBOHo4QXnn3IYq7y5YOvgR3Eg5ZRWFSKAo4sgn\n8JfgOvtd0TkhEognigWOJI71KhKYwy3lfQrDzQmMUiRwF/D1Dh9TZQcVvUcmAiWMugjMd/Z+6YQ1\nuM4vv0XdLK2j127sIHCWUDd20IO4hNcvUV5ueQT3P8XJ2zqTxcBNmPtBdPtWYJsf1R7HiUuRHbQX\nt23eNGkRWC/CW3HlwCnmowdfVXQC9xp1W3Vy0F8jJZrhfawbCZwD3J/bxzgkhvM5AUiLQIgEnker\nrdQR0cTHWUYjJzDs6qB7gA/m3psqyqqDDgLHSq5XZiVNPKO8nwC0duDrgEdUedxP+pmJJjUtplwE\n6tpBvRKB47hVOZ+DG4kWRQKHcUnseH2fVE6gzQ5S5S9zt98a3XwUFzUV2UGPRrfzIrATeBZuLkhR\nhVD+Md1s5j2Pn9jzCOlIIOQ66uYEzoU2v3kfbgJYPicA7v/JLwL3E+B3cJ+lq+iOIFJronZNZHWQ\nKh+oPquNqkigbK9kiwRKGJtIgGjij7dOTpKJ2GJak5nxPIEjZMnQKoLFczrFIrCYbK/cFKHDvhk3\naavKDoLqSCBlB5VxiHYRiO2gg7TaRnGnuwP3P/4a2RpEefIicJweiIDnEdKvV1z1VCcS2IgbycfE\ndlCdSOBuf+7LVYtnudcktG8tWed/FCdMg1qgbFTsoCaU5QT2UL6Kq4lACeMkAvNLACTuK4sEbgF+\nusTfjlmJ+1KGSCC1qcwSXLKxKCcw69t2O25y2wbK7SBoFYFUTqDuxiOBVCQQ20GPRrdbEsN+gbL1\nuOUnykQgv6zDUnpjM+whHQmEvWS7jQRq20GqbAfOUOWWuo0vISSHL8SvBOqj2E2qpaPYXjISkUBD\nCi0dL9CXlzzWqoNKGCcRyJf7BUtoCtfpxyIQzxPYjvsAXAwgwktECjepX+nPr7KDHqIiEvDe/G24\nKoi6kcBhWjv7o/56S+hMBA7hopkmdhCqHCSzX1LkK4pCJNALEWgUCYR1ffznAaojgTqJ4fBa9IKw\neNtm3NIe4fr5pdr7yUKNBKpeR4sEShgHEQghfn4dmCAQoYPO20GPwfxo66vAi33p538Ffq7g+Vbi\nqhbydlBsMyymhgj4v2/B+eud2EHzI5ZoLsFqOhvJhEigqDoo2EFLSIhA1JYl+aW8/e2wG1qgJzkB\nT1EkUGUHQatgn0u7CJRFAr8J/F2D9tblGG4g8sAAR/6pNsB4doj/BbdMTRNsxnAJ4yACVXZQWIXz\njMjyyZcRfhW3kfrL/XXideJjVlEvEngYWJ7Y6wDaRWCG8sTw8aiqIZ8TCOfkl5iuIpUTCGsVxXZQ\nKicAzOddjtK65DR4myxXiTGISCDYQWWT0uK8wEba7aBQHdSWGFbl1ppLGjflGHAZ2V7BwyB8hoZd\nHdQxqmzvIi9jdlAJ4yQCSTsI1yntw3VYwafPdxR/D/wU8J9w+x/EO0bFhEhgHdk8gVDjfao/ZwmZ\nnZKyS2IRCLupleUEjuZu5z+sh3EVJb3ICaTsoLy1k3/u/AqqqchhDve69GK09UX/k6d2JODfq/xm\nRuDEbgnu9eyVzVOXIAK3V53Y5zbA5I2KP4MroTYSjLoIxBO98iIQlisII/Z4uYMWEVDlAPA9XJL2\njymOBEJOYD1+gpC3ZOIRZni+A6STw7EI3IOrBCmzg/IikO/sD+E6rU7toHxO4BjZ6qd17KDw3MvA\nLYTnF8NLRQ49SwyrcqMq30zcVUcEwvt0BrA3PyHNRy8HcHMIyjaG7wfHcMuJDDMSmEgRUOWzqrVX\nK504xmGeQBiJpnICi3C7mB3DhfpnAHeStgz+3J9/N4lIwI8el+I67rOBucjyCJZQ6DhDEjGVF5jx\n7Ql175+meM30eLkIcHsn/zh3TpNI4BCus5+/tm/LEZzA5e2gShEA3oMry70hcX4v7aAiDuJmhR8m\nnTOA7H06jfZ8QGA/bhAw6EjgKC4KGwU7aKJEwChnHESgKicwhfuCxWvetImAKv8N5hObK0RYkkvQ\nhZUrH8GJRTx6j/MCQQSK5gqEEtHwvL9e8v+1RAKqybVUDtP5FoSh7fkE5GHca5SvDirqVGMRWO9/\nvk7aDupVYriIEAmUJfnC+5TKBwT24fa1HYYdpESVQUNgIiMBo5xRt4NCGegMrjPal7+PYjvocRL4\n0f0OnCUQswo44O/fQ+sEnrwIBDtopQjnivDJ6NzUYm9F5O2gFCES6HSyWLh+/lpn0GoHleUEYhFY\nhatuSUUOx/21+hkJdGIHpSqDAmFS0TBE4J4hVgaBe59Odrhcg7HAGQsRwHXuuxIVKbEIxJFAPE8g\nRbx1YGAlWQexm9bOORaBUB4Z7KArcUnnQCcicJjy6e6Q5QSaRAL5zjrMni6dJ5B77iACK3GbgKQW\neAuJ4X7bQUEEip6nTiQwTBEYphUU2mBRgNHCuIhAPh8Q35cSgaq17e+jPTm8kixZuItiEYhzAmF0\nHJdR9loEDtP5jOGySAA6ywmEnMxKnLg+I3H+ILzmsDl82eY1YdnviymOBPbhynIHtV5P4AjDrQyC\nrNrNMOYZh5zALG5ryvwXKF5VsrQ6KEFRJBBEYDetfn88YSzYQSEncDGt8wg6EYGv4xLZZcSbt9Sl\nLCcA7dVBRTmBsIAdZJVTVwL35s4LbRu2HXQUeBtOLL5RcM5+Bh8FAPwRw1+4bb8Izx5mG4zRY1wi\ngSugbf2WXkcCq2gVgbp20LwI+Mlq+c3iC1HlqGrpXr7QunlLXYpE4AguV3KMzu2gVcCNuDVaiiKB\nQdlBZZHAhcBLVAv/p30MQQRUeUC1ZW+DoaDK94fdBmO0GBcRuJJs4lUgP0+gExGoigTq2kHrcJUm\nM77qaBFwoseJtyYiUGYHPernPhzG/c/TFK9imc8J3IgbZadyAjAYO6hMBD4BXF2xZeGwIgHDGElG\n3Q46hut8noSb7BWTjwR2Aaf5zrhpTqBOYjiuDrqcTFBm/fP2elng0OH2yg46GP29DjiS2Moxfu6z\nozkUN0WPjRlEJBAEaREFIqDKt2pcZx+DnyhmGCPLOEQCTwXuTGxHGCaLzeK2ZzyB65xWUS0C95N1\nboE4EvgusDW67xDZngTBDtpPlqsIItEyR6BHdBMJpKqDYoFYR3E+IFxnGdkcijtxdlJqngD0UQSi\nHdNW0l3EsQ34Wk8aZRgLgFGPBEJHn1rPPb92EGSdVuE8AQBVjomwD2cfhUWpVuJmC6PKd3FCENiD\nm4UKrXaQ4Mr+rvTtEEZDBMoigbh8dD1uK8wiQnXQStwcijkRtjOc6iBwUcwauhAbVX6I25rTMAzG\nIxKA9nxAuC+2gyAre6yaJwBu68Sro9txYjjPbmCdT/zGaweBE4Ej/ngnlUF1aVIdNIdbHiPfWeft\noMWJc2KCqK4is8puwi2lnX8+6P/qlI/iVlS1WnfD6BHjEAlAcSSwklYRCCPXKjsI4D8CfyfCD1X5\nDq12UJ5duD1mZ/BrConMnxvsoDBZqtciEOya2tf1Hv8FibvydhCUi0AoEZ1/bVR5Q+K8QYrA2ZgI\nGEbPGPVI4AiuI0pNsimKBIIdVNohqfID4I3A50Xm15gv2qd0N84/D1YQOAvqBG4tmJAT6Gck0Ivr\nHiT7H8N16+QEygQSBmsHLR7A8xjGxDDSkYAqD4vwlLB1YI4yO2iKGvaJKl8U4TXA6yjv6HbjlmYO\nlUGo8qgIT/T5hSAC/cwJ9GKm56dxq4DG161jB1WJwCAjATARMIyeMdIiAKDaNjs1kJ8nAK0iUHeh\nrj8B/oxqOygfCcRtCzkBZYQjAT+B6rD/+4QIj1E/JzAKkYCJgGH0mJEXgRJSkUAnOYHAt3GCcgbF\nHV2oDlpKWlxCTuBxRlgECq5dNxIosspgcJFASGqP3faIhjGqjHpOoIyucgIBn0T9KG7DlKQ/rsqc\nv/ZZtE4iC8TzBEZhslhd6ojAfIloyXmD2rvWIgHD6DELIRIIewFDqx3USYf0F8CTSmbOgrOENlIc\nCSz2ber1iP048BH60/EdoTwxfBz3Wq7Fz6EoYBDLRoCJgGH0nHGPBBZRPE+gcLJYHlUOqfIbFaft\nxi01kRKBIzg7qOfVQaqoKr9RIVBNKY0EojWGzmY0IoFgB5kIGEaPGHcRCHZQsGCaRgJ12A2cR7kd\n1I8S0X5SZQeBixQ2YNVBhrEgWSgiECeGO8oJdMAuiiOBcRWBMA+jjEO4SKAsMWzVQYYxpiw0Eeh3\nJFBmB42jCOyjdd/mFIeoLhENkUBtC64hVh1kGD2mVARE5DoR2Ski23LH3yoid4jID0XkA/7YC0Xk\nuyLyA//7edH5l4nINhHZLiIf6VHbq+YJ9EMEzqS8RHTcROBXgL+pOCckjqtyAoPYwNwiAcPoMVWR\nwCdoXWQN37m/BLhUVS8BPuTv2g38gqpeCrwe+FT0sI8Bb1LVzcBmEWm5ZkOO46yfk9GM4n6LACyg\nnIAq+wtmY8fUEYE5BtMxmwgYRo8pFQFV/RbtdsGbgfer6gl/zm7/+3uq+rA/53ZgsYhMi8iZwHJV\nDUQl6TAAAApJSURBVCuBfhJ4WQ/afhxnU8Sdcr9zAlCeE+jHPIFhc4iSORSeOQZj0QQ7qN+2k2FM\nDE1yApuB54jITSKyVUQuT5zzCuBWLxRnAzui+x7wx7rlOG67wVgEBhEJLKScQB0O4/YSKCtRPc5g\nROBR3Pad/SiXNYyJpMlksSlgtao+U0SuAD6H22cXABG5GPg94IVNGiQi10Y3t6rq1oJTj+MWbEuJ\nQEfzBGpSZQeNY06gDoeo3o7xEeADA2jLfuA/D+B5DGOkEZEtwJZeXKuJCOwAPg+gqreIyEkRWauq\nj4jIBn/fNaoaFld7AFdnHthAtptXG6p6bc12hM52FCKBscwJ1KRSBPzWnu/vd0NUeRx4X7+fxzBG\nHT843hpui8h7m16riR10A/B8/8QXAou8AKwC/hb4TVX9p6ixDwEHReQqERHgGrLljLuhSAT6khPw\n6wcdwETAMIwFRFWJ6PXAjcCFInK/iLwBuA4435eNXo9bix/cTl2bgPeKyG3+J+zL+xbg48B24C5V\n/Uq3DfdVLSdpF4ElwDT98ah3k7aDFnJO4BDlE8UMwxhjSu0gVX11wV3XJM59HwWhuqreCjyl49ZV\nM0dUjaPK4yIcB5bTHxG4D+d/51nIOYHdwM5hN8IwjP4wzquIgutw8yPzw7ilj/shAj9PupOPS0QX\nmgh8DvirYTfCMIz+sFBFYBV9EAHVwjkACzYn4JOxVpdvGAuUcV47CNIicAg3f2CQ68vES0kvtMli\nhmEsYBaiCBymP/MEyjiBey2XssAiAcMwFjYLVQRggJGAn8F6FJeLMBEwDGNsMBHoHUeA1ZgIGIYx\nRixEEQgLnQ1aBI7iJqqZCBiGMTaMuwjMMTqRQGiHiYBhGGPDuIvAKNlBJgKGYYwdC0EE8iWZw8wJ\ngImAYRhjxLiLwBHaNzuxSMAwDKMm4z5j+O3A3tyxYSaG52zDE8MwxomxFgHV5L4EIRIY9FIHR7Eo\nwDCMMWPc7aAUw8wJmAgYhjFWLGQRsEjAMAyjgoUoAoeAk6qcHPDzmggYhjF2LEQROMzgrSAwO8gw\njDHERKB3WCRgGMbYYSLQO0wEDMMYOxaiCBwgm707SI5iG8oYhjFmLDgRUOUR4GlDeGrLCRiGMXYs\nOBEAUGXnEJ7W7CDDMMaOsZ4xPGJ8A7h32I0wDMPoBFEdnaVuRERVVYbdDsMwjHGim75zQdpBhmEY\nRj1MBAzDMCYYEwHDMIwJxkTAMAxjgjERMAzDmGBMBAzDMCYYEwHDMIwJxkTAMAxjgjERMAzDmGBM\nBAzDMCYYEwHDMIwJplQEROQ6EdkpIttyx98qIneIyA9F5APR8XeLyHYRuVNEXhQdv0xEtvn7PtL7\nf8MwDMNoQlUk8Ang6viAiDwPeAlwqapeAnzIH78IeCVwkX/MR0UkLGj0MeBNqroZ2CwiLdccZURk\ny7DbkMfaVJ9RbJe1qR7WpsFQKgKq+i1gX+7wm4H3q+oJf85uf/ylwPWqekJVfwzcBVwlImcCy1X1\nZn/eJ4GX9aj9g2DLsBuQYMuwG5Bgy7AbUMCWYTcgwZZhNyDBlmE3IMGWYTcgwZZhN6DXNMkJbAae\nIyI3ichWEbncHz8L2BGdtwM4O3H8AX/cMAzDGDJNNpWZAlar6jNF5Argc8D5vW2WYRiGMRBUtfQH\nOA/YFt3+MvDc6PZdwGnAu4B3Rce/AlwFnAHcER1/NfCnBc+l9mM/9mM/9tP5T1VfXvTTJBK4AXg+\n8A8iciGwSFX3iMgXgM+IyIdxds9m4GZVVRE5KCJXATcD1wB/nLqw7SpmGIYxWEpFQESuB54LrBWR\n+4HfBq4DrvNlo3PA6wBU9XYR+RxwO/AY8BbN9q58C/DfgMXAl1T1K334XwzDMIwOGak9hg3DMIzB\nMhIzhkXkaj/BbLuI/OaQ2nCOiHxTRH7kJ8H9uj++RkS+LiL/KiJfE5FVQ2jbqSJym4h8cYTatEpE\n/tJPGrxdRK4adrv8ZMUf+YmJnxGRmUG3KTXBsqwNRRMsB9CmD/r37vsi8nkRWTnINhW1K7rv7SJy\nUkTWDLJdvZog2+82iciVInKz7xdu8UU6zdrUNJnQqx/gVFxy+TxgGvge8OQhtOMM4Gn+72XAvwBP\nBn4feKc//pvA7w2hbf8J+DTwBX97FNr034E3+r+ngJXDbJf//NwDzPjbnwVeP+g2Ac8Gnk5rMUWy\nDbiJld/zn/vz/PfglAG16YXhuYDfG3Sbitrlj5+DKyy5F1gzAq/V84CvA9P+9roRaNNW4MX+738D\nfLNpm0YhErgSuEtVf6xuAtr/wE08Gyiq+rCqfs//fQi4A5fgfgmuw8P/HuhENxHZAPwc8HEgJM6H\n3aaVwLNV9ToAVX1MVQ8MuV0HgRPAEhGZApYADw66TZqeYFnUhtQEyysH0SZV/bqqnvQ3vwNsGGSb\nitrl+TDwztyxob1WdDZBdlBtegg38AJYhZt/1ahNoyACZwP3R7fDJLOhISLn4ZT3O8B6Vd3p79oJ\nrB9wc/4QeAdwMjo27DY9AdgtIp8QkX8WkT8TkaXDbJeq7gX+ALgP1/nvV9WvD7NNEUVtKJpgOWje\nCHzJ/z3UNonIS4EdqvqD3F3DbFenE2QHwbuAPxCR+4APAu9u2qZREIGRykyLyDLgr4D/U1Ufje9T\nF28NrL0i8gvALlW9jSwKaGHQbfJMAc8APqqqzwAO4z6UQ2uXiGwCfgMXAp8FLBOR1w6zTSlqtGGg\n7ROR3wLmVPUzJacNpE0isgR4D/De+HDJQwb1Ws1PkMUNyD5Xcu6g2vTnwK+r6kbgbbiqzSJK2zQK\nIvAAzgMMnEOrkg0MEZnGCcCnVPUGf3iniJzh7z8T2DXAJj0LeImI3AtcDzxfRD415DaBe392qOot\n/vZf4kTh4SG263LgRlV9RFUfAz4P/NSQ2xQoer/yn/0NZGF93xGRX8ZZja+JDg+zTZtwIv59/5nf\nANwqIuuH3K4duM8T/jN/UkROG3KbrlTVv/Z//yWZ5dNxm0ZBBL6LW1n0PBFZhFuJ9AuDboSICE5d\nb1fVP4ru+gIuwYj/fUP+sf1CVd+jqueo6hOAVwF/r6rXDLNNvl0PA/eLmywI8ALgR8AXh9iuO4Fn\nishi/16+ADdnZZhtChS9X18AXiUii0TkCfgJloNokLiVfN8BvFRVj+XaOpQ2qeo2VV2vqk/wn/kd\nwDO8lTa0dpFNkEWiCbJDbtNdIvJc//fzgX/1f3fepl5nshtmv/8NrhrnLuDdQ2rDz+B89+8Bt/mf\nq4E1wDf8i/w1YNWQ2vdcsuqgobcJeCpwC/B93Chp5bDbhUsm/gjYhkvATg+6TbiI7UHcRMr7gTeU\ntQFnf9yFE7EXD6hNbwS2Az+JPusfHWSbcu06Hl6r3P334KuDhvBazbfJf44+5T9XtwJbhvz+vQEX\n+X7H91f/BDy9aZtssphhGMYEMwp2kGEYhjEkTAQMwzAmGBMBwzCMCcZEwDAMY4IxETAMw5hgTAQM\nwzAmGBMBwzCMCcZEwDAMY4L5/wEW5mVFwhx7EgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efebaef23d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load this image, get the 32,32,15 voxel time series without the 5 first points.\n",
    "# plot the time course\n",
    "to_rm = 5\n",
    "img = nib.load(my_path_to_file)\n",
    "print(img.shape)\n",
    "(x,y,z) = (32,32,15)\n",
    "plt.plot(img.get_data()[x,y,z,to_rm:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.9044132678 1647.08928571\n"
     ]
    }
   ],
   "source": [
    "# - store the previous time series at the x,y,z voxel in 'bold'\n",
    "# compute bold mean and std\n",
    "bold = np.asarray(img.get_data()[x,y,z,to_rm:])\n",
    "bold_std = bold.std()\n",
    "bold_mean = bold.mean()\n",
    "print(bold_std, bold_mean)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.965607232455 100.0 1.93121446491\n"
     ]
    }
   ],
   "source": [
    "# - Normalized to have mean 100\n",
    "bold = bold * 100. / bold_mean\n",
    "bold_std = bold.std()\n",
    "bold_mean = bold.mean()\n",
    "sig_mag = 2*bold_std\n",
    "print(bold_std, bold_mean, sig_mag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(42, 3)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7efeb587c710>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGdJJREFUeJztnX+s5XdZ519vZlqURahNTV3bIWVxMOUPsPwYUBe5suwy\nbZSajbFWMQoaq0l1dzVSYBM7yWZ3dRXWmK7sLBZiiFINsDgkLbUGbvwF2GpbinR0RqjOtFKpoqvG\n1Znts3+c76Wnh7n3nPvc82vevF/JpOf7Pd8555nXeT7Pec73+5xTVRUhhBD8eMqqAwghhLAYUuBD\nCMGUFPgQQjAlBT6EEExJgQ8hBFNS4EMIwZSpBV7SOyQ9KumBHY75OUknJN0v6ar5hhhCCKHDLB38\nO4HD290p6Rrgq6vqIPADwNvmFFsIIYQ9MLXAV9VvAZ/b4ZDXAL84HPsx4CJJl84nvBBCCF3mcQ7+\nMuDU2PZp4PI5PG4IIYQ9MK+LrJrYzu8fhBDCitk/h8d4GDgwtn35sO9JSErRDyGEBlU12UTP/Ben\n/gGuAB7Y5r5rgNuH2y8DPrrNcTXLc3X+QL0W6pd2//c4sqiY+v+WJ2KC+jGoty7Q28uhftvB1cS/\n69ug3rtAb8+COnW+ezrHv+slUHcvNpYqKK2jpzV+/ar7d6d28JLeDbwCuETSKeBm4ILhWY9W1e2S\nrpF0Evh74HWtd5q98f/wnOl/CqN/26KItx7x1sfV3VoytcBX1fUzHHPjfMJp8zieSfMURv+2RRFv\nPeKtj6u7tcRFdDdpNuccxzzYHLu9ToVqc4FxdNncZv8qvW0u8Hm7bM543DIL/OaCn6fL5qoDmCdf\n1AW+qjbnH8remIhpbQr8eeBqnJV5O888TbK0Ar+OnmA9X7+98EVd4M8D1qbAn2fEW4+cojHDRfTj\nwL5VB7EA9rH4QhVvuyfe+ri6W0ucCrzLv2WcdKI94q1HOngzXES7Jk0KVY9Fj/s5e0uBN8JFtOts\nbea5eyy6UDl7yxy8ES6iXbuCdPA94q1HOngzXES7Jk0KVY9465ECb4aLaNekSaHqEW89UuDNcBHt\nmjQpVD3irUfGJM1wSVLXpMk8d49FeytA0hf8fxDOd9LBm+Ei2jVp0on2WKi3KopRkXdzlwJvhoto\n19GrjEn2yLhfj3gzw0W0a1eQDr5HOtEe8WaGi2jXpEmB75FC1SPezHAR7Zo0KfA9Uqh6xJsZLqJd\nkyYFvkcKVY+MSZrhkqCuSZMxyR4pVD3yxmiGi2jXpEkH3yOFqke8meEi2jVpUuB7LGPcz9FdCrwZ\nLqJdZ2szB99jGYXK0V3m4M1wEe3aFaSD75FOtEe8meEi2jVpUuB7pFD1iDczXES7Jk0KfI8Uqh7x\nZoaLaNekSYHvkULVI+OlZrgkqGvSLHTBDb+KiGSTB1ukUPXIG6MZLqJdkyYLrkfGJHsk38xwEe06\nepWxtR4Zk+yRfDPDRbRrV5COqke89Yg3M1xEuyZNFlyPeOsRb2a4iHZNmiy4HvHWI97McBHtmjRZ\ncD3irUemj8xwSVDXpMmC6xFvPfLGaIaLaNekyYLrkTHJHsk3M1xEuyZNFlyPeOsRb2ZMFS3psKTj\nkk5Iuukc918i6YOS7pP0CUnfu5BId8Z1tjZzyT0yB98j+WbGjqIl7QNuAQ4DzwOul3TlxGE3AvdW\n1dcCG8BbJO1fQKw74doVpKPqEW894s2MaaIPASer6qGqOgPcBlw7ccyfA88Ybj8D+MuqOjvfMKfi\nmjRZcD3irUe8mTGt074MODW2fRp46cQxbwc+JOkR4MuAb59feDPjmjRZcD3irUemj8yYVuBrhsd4\nM3BfVW1Ieg5wl6QXVNXfTh4o6cjY5mZVbc4c6c7U6OHR1i8kmpBC1SPeesTbGiBpg9Hp7j0zrcA/\nDBwY2z7AqIsf5+uB/wxQVX8i6dPA1wD3TD5YVR1pR7oDVZREsZyLRMskHVWPeOuRAr8GDI3v5ta2\npJu7jzVN9D3AQUlXSLoQuA44NnHMceBVQyCXMirun+oGtAccEycLrkfm4Hsk38zYsYOvqrOSbgTu\nZNSt3FpVD0q6Ybj/KPBfgHdKup/RC/eGqvqrBcd9LhzHrzK21iNjkj2Sb2ZMHWesqjuAOyb2HR27\n/RjwLfMPbdc4dgbpqHrEW494M8NJtGPiZMH1iLce8WaGk2jHxMmC6xFvPeLNDCfRjomTBdcj3npk\n+sgMpwR1TJwsuB7x1iNvjGY4iXZMnCy4HhmT7JF8M8NJtGPiZMH1iLce8WaGk2jH+drMJffIHHyP\n5JsZTqIdO4N0VD3irUe8meEk2jFxsuB6xFuPeDPDSbRj4mTB9Yi3Hpk+MsMpQR0XXApVj3jrEW9m\nOIl27AzSUfXYx3LGJN28pcCb4STaMXGy4HrEW494M8NJtOP4VcbWemRMskfyzQwn0Y6dQTqqHvHW\nI97McBLtmDhZcD3irUe8meEk2jFxsuB6xFuPXNQ3wylBHRdcClWPeOsRb2Y4iXZMnCy4HvHWI97M\ncBJt9dFPGr02VdSCn8rK20Dm4HukwJvhJNotcZax2MDPG6RQdYk3M5xEu83XLmMmGfy8Qebgu2QO\n3gwn0W6dQTr4PulEe8SbGU6i3RInBb5PClWPeDPDSbRb4qTA90mh6pE5eDOcEtRtwaXA90mB7xFv\nZjiJdusMltFNgZ83yJhklxR4M5xEuyVOOvg+KVQ94s0MJ9Fu41cZk+yTMckeGZM0w0m0W2eQDr5P\nOtEe8WaGk2i3xEmB75NC1SPezHAS7ZY4KfB9Uqh6ZEzSDKcEdVtwKfB9UuB7xJsZTqLdEicFvoGE\nAFWlUDVIgTfDSbTbR7/MwfdY5hujkzdIgbdjqmhJhyUdl3RC0k3bHLMh6V5Jn5C0OfcoZ8MtcdLB\n94i3PinwZuzf6U5J+4BbgFcBDwN3SzpWVQ+OHXMR8D+AV1fVaUmXLDLgHXCbr80cfI9lFXg3b5A5\neDumiT4EnKyqh6rqDHAbcO3EMd8JvLeqTgNU1WPzD3Mm3DqDdKI94q1POngzpom+DDg1tn162DfO\nQeBiSR+WdI+k755ngLvALXFSqHrEW5+MSZqx4ykamOn/B3oB8ELgXwFPAz4i6aNVdWKvwe0StwWX\nQtUj3vqkgzdjWoF/GDgwtn2AURc/zingsar6B+AfJP0m8ALgCwq8pCNjm5tVtbnbgHfALXFSqHos\n69qFmzdIgV8LJG0AG/N4rGkF/h7goKQrgEeA64DrJ475NeCW4YLsU4GXAm8914NV1ZE9xDoNt49+\nGZPsEW99UuDXgKHx3dzalnRz97F2LPBVdVbSjcCdjJL51qp6UNINw/1Hq+q4pA8CH2f04r29qj7Z\nDWgPuCXOMjv4C5bwPMsin3z6pMCbMa2Dp6ruAO6Y2Hd0YvtngJ+Zb2i7xm38KmOSPTIm2SdjkmY4\niXbrDNKJ9oi3PungzXAS7ZY4KVQ94q1PCrwZTqLdEieFqke89ckcvBlOCeq24FKoesRbg61f4WS2\n777sBStv646TaLfESaHqkTn4HgKqKgXeCSfRbh/9Ms/dI956pKEwxEm0W+JkwfWItx7xZoiTaLf5\n2szB98gcfI/kmyFOot06g3RUPeKtR7wZ4iTaLXGy4HrEW49cuzDEKUHdFlwKVY946xFvhjiJdkuc\nLLgeGZPskXwzxEm020e/fGTuEW89UuANcRLtljhZcD3irUe8GeIk2m38KmNrPTIm2SP5ZoiTaLfO\nIB1Vj3jrEW+GOIl2S5wsuB7x1iPXLgxxSlC3BZdC1SPeesSbIU6i3RInC65HxiR7JN8McRLtljhZ\ncD3irUe8GeIk2u3cXs6J9oi3HinwhjiJdkucLLge8dYj3gxxEu02X5u55B6Zg++RfDPESbRbZ5CO\nqke89Yg3Q5xEuyVOFlyPeOuRaxeGOCWo24JLoeoRbz3izRAn0W6JkwXXI3PwPZJvhjiJdvvol4/M\nPeKtRwq8IU6i3RInC65HvPVYircqCkBCi36u4JWgWXA94q1HvPVxc7e2OEl2m6/NXHKPzMH3WFa+\ngZ+7tcVJsltXkE60R7z1WNa1C/C7frG2OCWo24JLoeoRbz1yisYQJ8luSZNC1SNjkj1S4A1xkuyW\nNCnwPeKtRwq8IU6S3c7rZZ67R7z1SIE3ZKpkSYclHZd0QtJNOxz3EklnJf3b+YY4M25Jk060R7z1\nSIE3ZEfJkvYBtwCHgecB10u6cpvjfgr4IKzsCwxuo1cZk+yRMckeGZM0ZJrkQ8DJqnqoqs4AtwHX\nnuO4HwbeA3x2zvHtBreuIJ1oj3jrkTFJQ6Yl6GXAqbHt08O+zyPpMkZF/23DrppbdLvDbcGlUPWI\ntx45RWPINMmzFOufBd5YVcXo9MyqTtG4JU0KVY+MSfZIgTdk/5T7HwYOjG0fYNTFj/Mi4DZJAJcA\nV0s6U1XHJh9M0pGxzc2q2txtwDvgljQp8D3irUcK/JogaQPYmMdjTSvw9wAHJV0BPAJcB1w/fkBV\n/YuxwN4JfOBcxX049sgeYp2GW9KkUPVYlrcCJKGtX0g8z0mBXxOGxndza1vSzd3H2rHAV9VZSTcC\ndzK6KHJrVT0o6Ybh/qPdJ14AbhduMs/dYyneqiiJYrnTJ4skBd6QaR08VXUHcMfEvnMW9qp63Zzi\n6uCWNOnge6yiUKXA7w63nFtbnCS7zdZmDr7HMguVk7t9LHcO3ulT49rikpzg1xWkg++RTrRHvBni\nJNktaVLgeyzzlImTuxR4Q5wkuyVNCnyPFKoe8WaIk2S3pEmB75FC1SPeDHGSnHG/HvHWx8ldCrwh\nTpLdkiYdfI8Uqh7xZoiTZKeRNciYZJeMSfbIzwUb4iTZrStIB98jnWiPnNoyxCU5wWuxQQp8lxT4\nHvFmiJNkt6RJge+ROfgeKfCGOEl2S5oU+B4pVD3izRAnyW5Js5QFt/VTt9LK/kct8yaFqke8GeIk\n2e3CTS569Yi3HinwhjhJdkuaLLge8dYj3gxxkuw2W5u55B6Zg++Rnws2xCU5wa8rSEfVI956xJsh\nTpLdkiYLrkfGJHsk3wxxkuyWNFlwPeKtR7wZ4iTZLWmy4HrEW494M8RJstPIGmTcr0u89UiBN8RJ\nslvSZMH1iLce8WaIk2SnkTXImGSXjEn2yJikIS7JCX5dQTqqHvHWI94McZLsljRZcD0yJtkj+WaI\nk2S3pMmC6xFvPeLNECfJBSi/itjCacHFW494M8RG8vCztwUp8A2cFly89Yg3Q9wkO80lZ567R7z1\nSIE3xE2yU+JkwfWItx7xZoibZKe55MzB98gcfI/MwRvikpxbOHUG6ah6xFuPeDPETbJT4mTB9cgc\nfI/kmyFukp0SJwuuR7z1iDdD3CQ7JU4WXI946xFvhrhJdhpby7hfj3jrkQJvyEySJR2WdFzSCUk3\nneP+75J0v6SPS/odSc+ff6gz4ZQ4WXA94q1HvBkyVbKkfcAtwGHgecD1kq6cOOxTwDdW1fOB/wT8\nr3kHOiNOY2sZk+yRMckeGZM0ZJbkPAScrKqHquoMcBtw7fgBVfWRqvqbYfNjwOXzDXNmnDqDdFQ9\n4q1HvBkyi+TLgFNj26eHfdvxfcDtewlqDzglThZcj4xJ9ki+GbJ/hmNq1geT9E3A64Fv2Ob+I2Ob\nm1W1Oetjz4hT4mTB9Yi3HvG2JkjaADbm8VizFPiHgQNj2wcYdfGTQT0feDtwuKo+d64HqqojjRh3\ng1PiZMH1iLce8bYmDI3v5ta2pJu7jzWL5HuAg5KukHQhcB1wbPwASc8C3ge8tqpOdoOZA06JkwXX\nI956xJshUzv4qjor6UbgTkZXvm+tqgcl3TDcfxT4CeDLgbdJAjhTVYcWF/a2OM0lZ567R7z1SIE3\nZJZTNFTVHcAdE/uOjt3+fuD75xtaC6fEyYLrEW89lv3GeMGSnuuLGpfk3MJpLjlz8D0yB98j+WaI\nm2SnjiqdaI+MSfZIvhniJtkpcbLgesRbj3gzxE2yU+JkwfWItx7xZoibZKfEyYLrEW894s0QN8kW\niSMhgKrZv0W8Ryy8DaRQ9Yg3Q9wku8wlL3NkDXy8Qebgu6TAG+Im2WX8apmTIODjDTIm2SU/F2yI\nS3Ju4dIZLLNIgY83yJhkl3TwhrhJdkmcFPg+KVQ94s0QN8kuiZMC3yeFqke8GeIm2SVxUuD7pFD1\niDdD3CS7JE4KfJ8Uqh7xZoibZJextYxJ9smYZI8UeEPcJLskTjr4PilUPfLGaIhLcm7hMpecOfg+\nmYPvkZ8LNsRNsktHlQ6+T+bge+STjyFukl0SJwW+TwpVj3gzxE2yS+KkwPdJoeoRb4a4SXZJnBT4\nBsOvcCq/wtkiBd4QN8kuiZMC30PEW5cUeEPcJLuMX2UOvke89cmYpCFuBd5l/Cpjkj2W/cnHxRtk\nTNISN8kuH/1yiqbHst8YXbxBTtFY4ibZJXFS4HvEW58UeEPcJLskTgpVj3jrkwJviJtkl8RJoeoR\nb31S4A1xk+ySOClUPeKtTwq8IW6SXcavMu7XI976pMAb4ibZJXHSifaItz6ZgzfEJTm3cJmvzRx8\nj3jrkzl4Q9wku3RU6UR7xFufnKIxxE2yS+KkUPWItz4p8Ia4SXZJnBSqHvHWYOtXOCG/wumGm2SX\nxEmh6hFvPQRUfmbZj6mSJR2WdFzSCUk3bXPMzw333y/pqvmHOTMuiZNC1SPeesSbKTtKlrQPuAU4\nDDwPuF7SlRPHXAN8dVUdBH4AeNuCYp2FXY1fSdpYXCg9hpjWbp57jV2Ns3Jv54mnSVbg7fZnLPH5\nZmYdX7+9MO1d9BBwsqoeqqozwG3AtRPHvAb4RYCq+hhwkaRL5x7pbOx2/GpjQXHshQ3Wc9xvYwlx\n7JaNie118LaxxOeflY0p96/A24eevsTn2w0bqw5gnkxb1JcBp8a2Tw/7ph1z+d5Da+Hy0S8fmXvE\nW48VeJOW+HxftOyfcv+sF10mX6xlXayZ5CzwPRJfN9vhP/hciRctNKJd84PPBa4GHlvik54Frpb4\nwPaHrKeriZi+jOV2omeBFz3Z23nhaZL9LN3blzxz53xbFa3X7zur+NuFhLNHVLV9LZb0MuBIVR0e\ntt8EPF5VPzV2zP8ENqvqtmH7OPCKqnp04rFWVfRDCOG8pqpan3imdfD3AAclXQE8AlwHXD9xzDHg\nRuC24Q3hryeL+14CDCGE0GPHAl9VZyXdCNzJ6Er7rVX1oKQbhvuPVtXtkq6RdBL4e+B1C486hBDC\nVHY8RRNCCOH8ZeETALN8UWoZSDog6cOS/lDSJyT9yLD/Ykl3SfpjSb8u6aIVxLZP0r2SPrAOMUm6\nSNJ7JD0o6ZOSXroGMb1peO0ekPTLkp667JgkvUPSo5IeGNu3bQxDzCeG/P83S47rp4fX735J75P0\nzGXGda6Yxu77MUmPS7p4HWKS9MODq09IGr++uJKYJB2S9HtDTbhb0kvaMVXVwv4wOq1zErgCuAC4\nD7hykc+5QyxfCXztcPvpwB8BVwL/DXjDsP8m4CdXENuPAr8EHBu2VxoTo+81vH64vR945ipjGvLn\nU8BTh+1fAb5n2TEBLweuAh4Y23fOGBh9MfC+Ie+vGNbBU5YY17/eej7gJ5cd17liGvYfAD4IfBq4\neNUxAd8E3AVcMGx/xRrEtAm8erh9NfDhbkyL7uBn+aLUUqiqz1TVfcPtvwMeZDTD//kvag3//dZl\nxiXpcuAa4Bd4Ytx0ZTENnd7Lq+odMLoOU1V/s8qYgP8DnAGeJmk/8DRGF/2XGlNV/RbwuYnd28Vw\nLfDuqjpTVQ8xWoyHlhVXVd1VVVuz7R/jie+mLCWubVwBvBV4w8S+Vcb0Q8B/HeoTVfXZNYjpzxk1\nVQAXAQ93Y1p0gZ/li1JLZ5gKuopR4l9aT0z9PAos+1u4/x34cZ78RZNVxvRs4LOS3inpDyS9XdI/\nW2VMVfVXwFuAP2NU2P+6qu5aZUxjbBfDVzHK9y1WmfuvB24fbq8sLknXAqer6uMTd63S1UHgGyV9\nVNKmpBevQUxvBN4i6c+Anwbe1I1p0QV+7a7gSno68F7g31XVk76cUKPPQUuLWdI3A39RVffyhV8W\nW0lMjE7JvBD4+ap6IaPJqDeuMiZJzwH+PaOPpV8FPF3Sa1cZ07mYIYalxyfpPwL/VFW/vMNhC49L\n0tOANwM3j+/e4a8sy9V+4Mur6mWMGq1f3eHYZcV0K/AjVfUs4D8A79jh2B1jWnSBf5jRObctDvDk\nd6ClIukCRsX9XVX1/mH3o5K+crj/nwN/scSQvh54jaRPA+8GXinpXSuO6TSjLuvuYfs9jAr+Z1YY\n04uB362qv6yqs8D7gK9bcUxbbPdaTeb+5TzxUXspSPpeRqf/vmts96rieg6jN+j7h3y/HPh9jX63\napWuTjPKJ4acf1zSJSuO6VBV/e/h9nt44jTMrmNadIH//BelJF3I6ItSxxb8nOdEkhi9M36yqn52\n7K5jjC7YMfz3/ZN/d1FU1Zur6kBVPRv4DuBDVfXdK47pM8ApSc8ddr0K+EPgA6uKCTgOvEzSlw6v\n46uAT644pi22e62OAd8h6UJJz2Z0KuD3lhWUpMOMOtJrq+r/TsS79Liq6oGqurSqnj3k+2nghcPp\nrVW6ej/wSoAh5y+sqsdWHNNJSa8Ybr8S+OPh9u5jmvdV4XNcJb6a0cTKSeBNi36+HeL4l4zOc98H\n3Dv8OQxcDPzGIPHXgYtWFN8reGKKZqUxAS8A7gbuZ9TdPHMNYnoDozeaBxhdzLxg2TEx+pT1CPBP\njK4tvW6nGBidkjjJ6A3q1UuM6/XACeBPx3L955cZ11hM/7jlauL+TzFM0awypiGP3jXk1e8DGyuK\naTynXszo+uB9wEeAq7ox5YtOIYRgisNPnYYQQjgHKfAhhGBKCnwIIZiSAh9CCKakwIcQgikp8CGE\nYEoKfAghmJICH0IIpvx/mBRAoPEveA4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efeb825d850>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# - divide the total duration time in 12, into 12 intervals of \"cond_duration\" TR\n",
    "# we will fill the time with 3 conditions, (eg A,B,C)\n",
    "# A, B will be our experimental conditions, C is for control condition (baseline)\n",
    "\n",
    "cond_duration = bold.shape[0] // 12 # integer division\n",
    "\n",
    "# make a signal for two conditions: A, B. Each block of A or B lasts for cond_duration, \n",
    "# and these blocks occur regularely, eg if the time is divided in 12 intervals,\n",
    "# 'A' blocks will occur at intervals 0, 3, 6, 9, 'B' at intervals 1, 4, 7, 10, and the \n",
    "# remaining intervals will be baseline \n",
    "\n",
    "ones = np.ones((cond_duration,1))\n",
    "a = np.kron(np.eye(3), ones) \n",
    "\n",
    "print(a.shape)\n",
    "signal_A = np.kron(np.ones(4), a[:,0]).flatten()\n",
    "signal_B = np.kron(np.ones(4), a[:,1]).flatten()\n",
    "baseline = np.ones(shape=(bold.shape[0]))\n",
    "plt.plot(signal_A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(168, 3)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD/CAYAAAAUnaZMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADG5JREFUeJzt3V+MnNddxvHvQ5xc9I8UrFRu/rg4EolKpEjJRQwiirIX\nELVFSnKB0oIQFqqiiqqAuCEpIJw7okggxA1CkEiuaF0soJEjJTQOykq5IW1E2hoct7FUS3bkbEGk\n0JAbh/642JeyWby7MzuzO35//n6klc+8c/zOOT6eZ8+cmXlPqgpJUk8/tugGSJJ2jiEvSY0Z8pLU\nmCEvSY0Z8pLUmCEvSY3NPeSTfCzJ6SSvJ3lk3ueXJE0u8/ycfJKrgG8DPwe8AXwd+KWqem1uDyJJ\nmti8Z/IHgTNVdbaqLgJfBh6Y82NIkiY075C/ETi35vb54ZgkaQH2zPl8W679JPE6CpK0DVWVaf/O\nvGfybwD719zez+ps/j2qqu3P4cOHF96GnfzRuHX+/9m5b7M89+Yd8q8AtyQ5kOQa4JPA8Tk/hiRp\nQnNdrqmqd5N8DvgqcBXwZPnJGklamHmvyVNVzwHPzfu8Y7G0tLToJkgb6vz/s3PfZjHXz8lP9IBJ\nubY7XsnU7/voMuJzb7ySUJfBG6+SpMuIIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRny\nktSYIS9JjRnyktSYIS9Jjc39KpST8CJX4+UFrqRxcSYvSY0Z8pLUmCEvSY0Z8pLUmCEvSY0Z8pLU\nmCEvSY0Z8pLUmCEvSY0Z8pLUmCEvSY0Z8pLUmCEvSY0Z8pLUmCEvSY0Z8pLUmCEvSY0tZGcojZe7\neo2bO3tdeZzJS1JjhrwkNWbIS1JjhrwkNWbIS1JjhrwkNWbIS1JjhrwkNWbIS1JjhrwkNWbIS1Jj\nhrwkNWbIS1JjhrwkNWbIS1JjhrwkNWbIS1JjhrwkNbaQ7f/cgmy83P5PGhdn8pLUmCEvSY0Z8pLU\nmCEvSY0Z8pLUmCEvSY0Z8pLU2LY/J5/kLPCfwH8DF6vqYJK9wF8DPwGcBR6qqu/PoZ2SpG2YZSZf\nwFJV3VlVB4djjwInqupW4B+G25KkBZl1uWb91x/vB44M5SPAgzOeX5I0g1ln8i8keSXJw8OxfVW1\nMpRXgH0ztU6SNJNZrl1zd1VdSPIh4ESS02vvrKpKcsmL1Dz22GM/Ki8tLbG0tDRDMySpn+XlZZaX\nl2c+T+ZxsbAkh4G3gYdZXad/M8n1wItV9dF1dcsLlI2XFygbN59745WEqpr6Cbit5Zok70vywaH8\nfuA+4CRwHDg0VDsEPL2d80uS5mNbM/kkNwNfGW7uAb5YVX84fITyGPARNvgIpTP5cXMmP24+98Zr\nuzP5uSzXTPWAhvyoGfLj5nNvvHZ1uUaSNA6GvCQ1tpDt/3zJP16+3JfGxZm8JDVmyEtSY4a8JDVm\nyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtS\nY4a8JDW2kJ2hNF7u6jVu7ux15XEmL0mNGfKS1JghL0mNGfKS1JghL0mNGfKS1JghL0mNGfKS1Jgh\nL0mNGfKS1JghL0mNGfKS1JghL0mNGfKS1JghL0mNGfKS1JghL0mNGfKS1NhCtv9zC7Lxcvs/aVyc\nyUtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtS\nY4a8JDVmyEtSY4a8JDVmyEtSY5uGfJKnkqwkObnm2N4kJ5J8J8nzSa5dc9/nk7ye5HSS+3ay4ZKk\nrWWzrfiS3AO8DXyhqm4fjj0B/FtVPZHkEeDHq+rRJLcBXwLuAm4EXgBuraofrjune/+NmFs3SouR\nhKqaev/NTWfyVfUS8Na6w/cDR4byEeDBofwAcLSqLlbVWeAMcHDaBkmS5mc7a/L7qmplKK8A+4by\nDcD5NfXOszqjlyQtyJ5Z/nJV1RbLL762l6RtWF5eZnl5eebzbLomD5DkAPDMmjX508BSVb2Z5Hrg\nxar6aJJHAarq8aHe3wOHq+rldecz+EfMNXlpMXZkTX4Dx4FDQ/kQ8PSa459Kck2Sm4FbgK9t4/yS\npDnZdLkmyVHgXuC6JOeAPwAeB44l+TRwFngIoKpOJTkGnALeBT5bTvskaaG2XK6Z+wO6XDNq/t6W\nFmM3l2skSSNhyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSYzNdalhXnmTq\nb1XrMuJlKa48zuQlqTFDXpIaM+QlqTFDXpIaM+QlqTFDXpIaM+QlqTFDXpIaM+QlqTFDXpIaM+Ql\nqTFDXpIaM+QlqTFDXpIaM+QlqTFDXpIaM+QlqTFDXpIaW8j2f25BNl5u/yeNizN5SWrMkJekxgx5\nSWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5SWrM\nkJekxgx5SWrMkJekxhay/Z9byI2XWzdK4+JMXpIaM+QlqTFDXpIaM+QlqTFDXpIaM+QlqTFDXpIa\nM+QlqbFNQz7JU0lWkpxcc+yxJOeTvDr8fHzNfZ9P8nqS00nu28mGS5K2ls2+wZjkHuBt4AtVdftw\n7DDwg6r643V1bwO+BNwF3Ai8ANxaVT9cV8+vTI6Y33iVFiMJVTX15QI2nclX1UvAW5d6vEscewA4\nWlUXq+oscAY4OG2DJEnzs901+d9I8s0kTya5djh2A3B+TZ3zrM7oJUkLsp2Q/zPgZuAO4ALwR5vU\n9bW9JC3Q1FehrKrv/W85yV8Czww33wD2r6l603BMkjSl5eVllpeXZz7Ppm+8AiQ5ADyz5o3X66vq\nwlD+beCuqvrlNW+8HuT/3nj9yVr3AL7xOm6+8SotxnbfeN10Jp/kKHAvcF2Sc8BhYCnJHawuxXwX\n+AxAVZ1Kcgw4BbwLfHZ9wEuSdteWM/m5P6Az+VHz97a0GDsyk5fWc1evcfOX9JXHyxpIUmOGvCQ1\nZshLUmOGvCQ1ZshLUmOGvCQ1ZshLUmOGvCQ1ZshLUmOGvCQ1ZshLUmOGvCQ1ZshLUmOGvCQ1ZshL\nV5B5bCd3uerct1kY8tIVpHMQdu7bLAx5SWrMkJekxtzjVZJGYjt7vO56yEuSdo/LNZLUmCEvSY0Z\n8pLU2I6HfJK9SU4k+U6S55Ncu0G9s0m+leTVJF/b6XbNKsnHkpxO8nqSRzao86fD/d9Mcudut3EW\nW/UvyVKS/xjG69Ukv7+Idm5HkqeSrCQ5uUmdMY/dpv0b+djtT/Jikn9J8s9JfnODeqMcv0n6N/X4\nVdWO/gBPAL8zlB8BHt+g3neBvTvdnjn16SrgDHAAuBr4BvBT6+p8Anh2KP808I+Lbvec+7cEHF90\nW7fZv3uAO4GTG9w/2rGbsH9jHrsPA3cM5Q8A32723Jukf1ON324s19wPHBnKR4AHN6k79ceDFuQg\ncKaqzlbVReDLwAPr6vyo31X1MnBtkn2728xtm6R/MJ7xeo+qegl4a5MqYx67SfoH4x27N6vqG0P5\nbeA14IZ11UY7fhP2D6YYv90I+X1VtTKUV4CN/rELeCHJK0ke3oV2zeJG4Nya2+eHY1vVuWmH2zUv\nk/SvgJ8dXg4/m+S2XWvdzhvz2E2ixdglOcDqK5aX193VYvw26d9U47dnTo05werLjPV+7z0tq6pN\nvgx1d1VdSPIh4ESS08OM5HI06ZcL1v+2HcuXEiZp5z8B+6vqnSQfB54Gbt3ZZu2qsY7dJEY/dkk+\nAPwN8FvDjPf/VVl3e1Tjt0X/phq/uczkq+rnq+r2S/wcB1aSfHho+PXA9zY4x4Xhz38FvsLqksHl\n6g1g/5rb+1mdLWxW56bh2Bhs2b+q+kFVvTOUnwOuTrJ395q4o8Y8dlsa+9gluRr4W+CvqurpS1QZ\n9fht1b9px283lmuOA4eG8iFWf+u8R5L3JfngUH4/cB+w4ScfLgOvALckOZDkGuCTrPZzrePArwIk\n+Rng+2uWrS53W/Yvyb4kGcoHWf329L/vflN3xJjHbktjHruh3U8Cp6rqTzaoNtrxm6R/047fXJZr\ntvA4cCzJp4GzwEND424A/qKqfoHVpZ6/G9q9B/hiVT2/C23blqp6N8nngK+y+kmUJ6vqtSSfGe7/\n86p6NsknkpwB/gv4tQU2eSqT9A/4ReDXk7wLvAN8amENnlKSo8C9wHVJzgGHWf0U0ejHDrbuHyMe\nO+Bu4FeAbyV5dTj2u8BHoMX4bdk/phw/r10jSY35jVdJasyQl6TGDHlJasyQl6TGDHlJasyQl6TG\nDHlJaux/ACCWASBhosDlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efeb5951290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# - construct the design matrix with condition A, B and baseline using np.vstack\n",
    "X_ = np.vstack((signal_A, signal_B, baseline)).T\n",
    "plt.imshow(X_,aspect='auto',interpolation='nearest')\n",
    "print(X_.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7efeb57e9650>,\n",
       " <matplotlib.lines.Line2D at 0x7efeb5946490>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmYXFd1rv8uqeduqVttyZIseZRHDAYbbEYbAQYbQww4\ncYCYmYRcCIkhBHD43cRWbiBAQgYSSH4JMTgMDkMMGELAZhBwmQcbbPCAZ8losobW2K2WtO4fe++u\nXafOOXVq6Bq61/s8/XTVqWnXqVPf+erba+8tqophGIYxt1nQ7gYYhmEYs4+JvWEYxjzAxN4wDGMe\nYGJvGIYxDzCxNwzDmAeY2BuGYcwDcsVeRK4VkS0iclu0bVxEbhaRu0XkJhEZSzzmOBHZKyJvma1G\nG4ZhGLVRzdl/GLg4se0q4GZVPRX4mr8e87fAfzeneYZhGEYzyBV7Vf02sDOx+VLgOn/5OuCF4QYR\neSFwH/DLJrbRMAzDaJB6MvvlqrrFX94CLAcQkRHgbcA1zWmaYRiG0Swa6qBVN9dCmG/hGuDvVHU/\nIA22yzAMw2giPXU8ZouIrFDVzSKyEtjqt58H/KaIvBcYA46IyAFV/WDyCUTEJuQxDMOoEVWt20jX\nI/Y3Aq8E3uP/f8434oJwBxG5GtiTJvSBRho9G4jINap6TbvbEWNtKkYntgk6s13WpmJ0aJsaMsnV\nSi+vB74LnCYiG0Tk1cC7gWeLyN3AM/11wzAMo4PJdfaq+tKMmy6s8rh1dbfIMAzDaDo2grbE+nY3\nIIX17W5ACuvb3YAU1re7ARmsb3cDUljf7gaksL7dDUhhfbsb0GykHYuXiIh2WmZvGIbRyTSqm+bs\nDcMw5gEm9oZhGPMAE3vDMIx5gIm9YRjGPMDE3jAMYx5gYm8YhjEPMLE3DMOYB5jYG4ZhzANM7A3D\nMOYBJvaGYRjzABN7wzCMeYCJvWEYxjzAxN4wDGMeYGJvGIYxDzCxNwzDmAeY2BuGYcwDTOwNwzDm\nASb2hmEY8wATe8MwjHmAib1hGMY8wMTeMIyqiPBWEYbb3Q6jfnLFXkSuFZEtInJbtG1cRG4WkbtF\n5CYRGfPbzxORW/zfz0XkxbPdeMMwWsYfAce3uxFG/VRz9h8GLk5suwq4WVVPBb7mrwPcBjxeVc8G\nngN8QEQWNrOxhmG0jSFgsN2NMOonV+xV9dvAzsTmS4Hr/OXrgBf6+x5Q1SN++yAwoaqHm9hWwzDa\nxyAw0O5GGPVTT2a/XFW3+MtbgOXhBh/l/AL4BfDHTWifYRhtRgTBib05+y6mp5EHq6qKiEbXfwic\nKSKnA18WkfWqOpH2WBG5Jrq6XlXXN9IWwzBmjeDoTexbiIisBdY26/nqEfstIrJCVTeLyEpga/IO\nqnqniNwLnAz8JO1JVPWaOl7bMIzWM+T/m9i3EG+A14frInJ1I89XT4xzI/BKf/mVwOd8Q04QkR5/\n+XjgFOBXjTTOMIyOIIi8ZfZdTK6zF5HrgacDS0VkA/DnwLuBT4nIa4EHgN/2d38acJWITAPTwOtU\ndfdsNdwwjJYxmPhvdCGiqtXv1ewXFVFVlZa/sGEYNSPCY4FbgTer8vftbs98pVHdtBG0hmFUw5z9\nHMDE3jCMasy5DloRjvElpfMGE3vDMKoxFzto/wt4arsb0UpM7A3DqMacc/bAIuC4djeilZjYG4ZR\njbmY2Q8Cq9rdiFZiYm8YRjWGgP3MLbEfAFa3uxGtxMTeMIxqDAI7mFuZvTl7wzCMBEPAdhLOXoQe\nEVa2p0kNY2JvGIaRIDj7ZIzzXODfW9+cxvAllxbjGIZhJEh19rhqlm5cqnAAOAwsF2HeLLBkYm8Y\nRjWyMvvVKdu6gQFgL+49Hd3mtrQME3vDMKqRFeOsAvqTdxbhLBEWtaJhdTIIHAAeZh5FOSb2hmFU\nY4hssS9z9iIMAjcBz2xN0+oiFvt500nb0EpVhmHMC7Kc/WqgL7HtNbilSiscfwcxCEwCGzFnbxiG\nMUOFs/cVLWUxjgi9wNuAu+h8sZ93zt7E3jCMaqR10C7GnQTibb+DW53uO3S22A9gYm8YhlHBELAb\nZtw7OJHcSLmon4tbtnSKzhZ766A1DMNIYRA3N84BSlHOauBeyp19uF83iH3I7M3ZG4ZheIZwQh+L\n/SrgIUBFZgo9BnAiepDOF/uZGGe+LGJiYm8YRjViZx+cfIhxpqJtQUSnqKzS6SQGgAOqLprC9T/M\neUzsDcOoRnD2k5THOA9THtkEZ98NMc4Bf3kvc2vq5kxM7A3DyCSaNCwtxnkYJ+5pzr5bxH6azv4V\n0jRM7A3DyKMfOKjKESrFPsQ43ejsJ/3lg5jYO0TkWhHZIiK3RdvGReRmEblbRG4SkTG//dki8mMR\n+bn//4zZbLxhGLNOWKUKyjP7EOPEzj78AugGsY+dfW/OfecMRZz9h4GLE9uuAm5W1VOBr/nrANuA\n56vqWcArgY82q6GGYbSFWBgPAIMi9AOjwFbKhT045k6vxgknJTBnX0JVvw3sTGy+FLjOX74OeKG/\n762qutlv/yUwKCLz4qxpGHOU2NmHDtqVwGYf7WQ5+04W0PgEdhBz9rksV9Ut/vIW3MRHSX4T+Imq\nTtf5GoZhtJ8KZ48T+01+W5qz77YYJ/XEJILMpRr8hme9VFUVEY23iciZwLuBZ2c9TkSuia6uV9X1\n9by+CC8BBlT5SD2PNwwjl2RmnxT7bs3si3TQXol7/+9qRaOSiMhaYG2znq9esd8iIitUdbOIrMRl\ndwCIyGrgBuDlqnp/1hOo6jV1vnaS04GRJj2XYRjlJJ39AD7G8du60dnHmX1eB+0y2liD7w3w+nBd\nRK5u5PnqjXFuxHXA4v9/zjdmDPhv4O2q+r1GGlYDPcyTQRFG5yDCpSL8Vbvb0QLSMvsVJGIcv5Zr\nD84pd7rYJzP7LGc/wBzSliKll9cD3wVOE5ENIvJqfEQjInfjVqR5t7/7G4E1wNUicov/WzpLbQ/0\n0p3rYBrdzeOprFKbi2Rl9sHZhxhnAJhURen8apyiHbRzSuyrxjiq+tKMmy5Mue9fAn/ZaKNqxJy9\n0VRE+FfgU6p8NeduxwNnitCvylSLmtYOwlQJ+P9LSO+gjaORbqjGCZl93gjaOSX2c2EErTl7o9kc\nD5xW5T7H4Y69M2e/OW0lTIIGpcw+jnGCs48F1GKcDmQuiL05e6PZhFw6j+OBHwBnz35z2krS2Sdj\nnCxn38liX7SDdpA5ZCTngtibs58DiHSUQw6ClooIC3DTBXyeuS/2sbOfBIZxVSpbom3m7LuAuSD2\n5uy7HBEGgFva3Y6IIXLEHuf6dwHfY36IfezsjwN2qRIGS3ajs5+XHbRzRezN2Xc3Q0CvL9/rBKrF\nOMcDDwK3Amd1ULubhgij/iScHFR1IqW8HtKdfafPNzMvO2gbHkHbAfQyhz6QecqQ/99HyXG1kyHy\nDcTxwIOq7BJhG3AycFdLWtY63gkosJByZ38M8Ivofl3l7BPz80P1GGfOGElz9kYnEE7WneIGB4Gj\ncxz7cThnDy5+OqclraoREU4U4Yw6Hz4G/C7OyceZvVDqnA3b0px9r+/bqAsRnlt0XhoRnuV/hRSh\nDzikymF/vVoH7ZwxknNB7M3Zdz+xs+8EhnDL1WUNCDwet9g2wJ3AKa1oVB28BPi9Oh87jOuEvZhy\nZw/lMU6Fs48GVjXyeX4OV9NfhA8DZxW8b5zXg3XQdhXm7LufjnH2IvTivhcPkd1JGzJ7cB21nbpg\ndR/1xynDuAnADlGe2UO5sw8LjsfOPmyv67W9o+/DzZlf7b4jwLF4dy7CQhF+J+chaWJvHbRdgpVe\ndj/B2XdCzhtKDTeR3Ukbi/1uOlfsG/luDOPWpHgT8HO/Lc3ZT1KZ2UNjuX0Q36pij5sIEUpG4Wjg\nYyKZvwqSJ6VqHbT9jcRRnURb34QI/SKc1+DT9NBZlRxG7XSMs6fk/DaR7ezjzH4PnSv2fTQm9vtU\n+YDqzHvNi3GSItpIjBMeV0TsQ59EX/RfgKdm3D95UqoW40DipCXC6/wvQERYJDIzKWRH0+4z1vnA\nPzX4HMEFmLvvXjopsw+lhptJcfYijOG+N7v8pk529o2I/QiwL7EtiHlaB20znX04DsYK3DdN7AEu\nyLh/oRjHR0n9wASVUc5f437dATwGeEeBdraddov9ChoX6VA+OmeytXnIrDt7ES4KbqxAW/KcfSi7\nDAv27AYWNaeVTafRGCcp9rU4+2aIfVFnf4hKsT8/4/5Jsc+KcfpxJ4L9RNriE4TFlDqPl+D2VcfT\nbrFfSeNib86++2mFs/834FEF25In9vEkYDB3nf0wriIpZhp4VWJ7Jzj7OykX+7twg93SRDgtbkoz\nAQP+fmEO/0A4AY1H/03sC7CSxh25OfvupxWZ/VKKfSlDB21qjINzcjui6y3N7EW4QqRwXX9VZy/C\nuAjXJLYJKc5eFVXluuhXDbTR2YvQB5yAG+QVi/0E8DPgiSkPy3T2IiwT4e1+exD7MNNnIJyAzNnX\nSDNinF7cSD9z9t3LrDp77/AGKbZ8ZTVnP0652Lfa2V8GhYsaijj7c4CrEhHXAG7g0aECrzEbzj48\nrlqMczKuRHYPpWMnxC/fJj23z+ugPQV4RXS/IPaxkQxiHzv73oIRYVtpt9g3I8bpwf2sNGffIkS4\nRITLmviU4bObrdLLMDiqWc5+Z3R9N7Co6GjPJrCE4sd6qtiL8I8iM0K6Brff45G2aXl9FrNZjVMt\nxjkDuCPxWn2UxD6tIievg3aAkiEI7ycp9sHRjyeud7y7b7fYN8vZ72nC8xjFOQ94ehOfb7ZjnGX+\nf2Fnr8oecKV1idvLnL1fpUpp3RiBcUq/hKpREeP4aQXeQGm2zpP9/3j2zrS8PoswqGo2Mvtqzj5P\n7B8k/ZdZXgdtLPbh/WQ5+6ToFzm22kq7xX4l0CPS0IRsPTixN2ffOnop1nlWlNnuoA3OvsgXMp6/\nfRPwcRH+NRrHkXT2UDC3F+GpIvxxgTbk0aizPwX3vQ9Ofg3wYyrFvqizjwdVtboa5zRc52ya2O8i\n/RjN66AdoFRZldVBuwQ4jDn74ojMTDK0n8ZceRB7c/ato5diZXFF6Thn7y+/HLge+C3cyEyozOyh\nePnl6cDTCtwvj3GKi31aB+0Zif9rgM9Qv9jH0yU009nvo7qhOAbYSKXYT+FOyGmjaKs5+17f8ZvX\nQfsglc7exD6HFbhcNLkzayXEOObsW0ezxX6I8k62ZlNPZo8q31flepygLPe3Z4l9kU7agYL3S8V3\nAo5QPMZJmxvnDNxMnWf4foYg9o+LpgVIG1CVxUGc4RqmAWcvwl+J8KSo3duofowdDWwl3dnvw011\nkDym8jpog4aMkN9Bex/lzn4nJva5hHUsQ29+vZizbz3NjnEGceVys+nsd1M8xknOqb+VkrNPi3GK\nin1/wftlEdxkIzHOGcAN/v/RwKQq9+L2/0n+PoUze1+GOYUT5qSI1uLsz8NNpxzaXUTsl5Mh9r5d\naVFOtQ5acL/S8jpo76e8GmcDJva5hMEp5uy7j9lw9ruYXWd/P8VjnP2JbVvIj3GK1tr309ho2yAw\ntXTQLkz0iZ0BfAUnWo8D7vXbb6EU5dQS44ATxTEqnX0tn+cYpZNDP07EMw2F70NZAjxCudiH0ktw\nJ+VqYp+McaDk7A9QmdkHZ7/E/zKaG2IvIteKyBYRuS3aNi4iN4vI3SJyk4iMRdu/ISJ7ROQfC7z2\nSpzYm7PvPmYjs9/F7FW0LAMeoHiMk+bsQ4yT5eyLiHhDMQ71OfvwukEgT8FVsNwFPI/miP0UTgQb\nyexHKX2Hw8CoBTmLkhyFWwv3EOkxDmQ7+7wOWsiPcZbg9lkYORv6B7q+GufDuMULYq4CblbVU4Gv\n+evgdsz/Bv6k4GuHzD555qyVrnX2IvyqGwZjpNALjDaxtrwrnL0XnoUpt7cqxhnHOdFanD2UROwE\n4BFV9uIE/zcoif3PgUf7y/U4+0arcUaj+4dO1l1km4rluM8FssU+rZM2+d7ynH1WB+2vcdp5jH+N\nfXS7s1fVb1PpYi4FrvOXrwNe6O+7X1W/g/uQitCws/cdSgtwO7urnL3/aX0yXeAIUujFiV6zDvDU\nzF6EBU2aujo4+0Yz+yXAzsR0AVBbB+1IA/OjLwEeJsXYZJiGpIiFunT8/xMoif0mSr9eaumghdJ3\nvi5n701DHOMEwZ4gW+xD5yzU5uzjk0TysWG/LiK/g3YnLspb4//vpdvFPoPlqhp21hZKB0gg+UXI\nohkdtD24Ge8azf3bQTiwGz5IRLhYhL9t9HlqoJbFJYqQ5ex/D/g/TXj+pURiL8JRInwxpy1ZMU5a\nXg+1ZfZQ/wl+HCf2Zc7eTwexMeXEGIQvTezv9P+D2D9CedVS0UFVUBL7qs5ehCeIzJjFwCDlK841\nKvahPWnOPiQKgawYJ6+Ddpd/7jV0kbNvZDATqqoiUlTcyxl81bn0fGsvh89bw6hcJOv+84Gan+P3\nHzvAZz51mCX3jjA1ukLWPfUU3ICH+/Vqra9draMmsZd1chSlDrpyTv7SU9ix5o/kuR/4Ak96/8YG\n2nSfXq2Hq99t5ssxhhOfRgnOPrkvllDqGK0LeeMZSxj//BiPukG56/lLZN1jTuFp7zqDX/7W8+R/\nXf4YVv5ssuwBK356FMt/PibrXlVaV/bJf9PHvc85luO/fRYPPXW/rHtc+Zqz/Tv3MTV2VIHmBDFZ\njPs1UCtB7B+X2L4ct59OBO6Jtvf614nF/gf+chD9IPbbKI1HGGZo6x5Zt7zY2rpH/1w5NAB/cOax\nsm7abVvz5cUcHFku655W/hwnfP3l7F51hbz0T/6J07/wS71a91ES9KSzzxoYFd5zVowT9m3a40Oi\nEEjGOIcpj3HSOmh3Ue7s91Had01DRNYCa5v1fPWI/RYRWaGqm0VkJaWza2287lvjLN5wHvsnl9K3\n9wrgBTU/x/LbFnDFJb0MPfL7HBoYAJ4MrAKeDXynrna1jvAFLOYIdqy5j9EH97DwUNJ1wmUvW8Lk\nkgX07f0s7ktbD8uB/wV8osB9m+3sQwdt0oX10qhjGnvwE7z66cLQ9vfz+P9/FfAlLvjLQc6+FhZv\n+G+SseNLL13N4I6TiN39hX/aw3kfWM3QtvcxtXgx8KXoEUdx2Su+zvU37inQmiBm9eb2IcZ5SmJ7\ncORnUC72fTjnGY61E4GP+8u/wi3qHYRvAhgUoR8Y5iWXXQS8xW/P54rnreZwn7JwurRffusloxzu\n76to64tesYrDfbBo01dxx9rrKQlys5x9aowTDeSMo+mks99GKcbZTZQa+D6bBX5bcPbbcGJ/QkY7\n60ZV1wPrS+2Xqxt5vnrE/kbglcB7/P/PJW4v1mn3D/cpbrjzh4CvqvIftTZEhKNwB+1bgItVeZms\nkxspHfydTG0xTt/eYT79ydfpHZd9MnmTCH+JO0BfDjxdlV/X2hhZJ/9E8f3Wi3NADdfa+7w2q86+\nj0b7NHTBCm746Abuv/CJwG2qnCLCb+IGEl2lWn5yE+FHwB+o8sNoWxCIdwDPUC0tQyfr5K0Mbzmf\n4pk91F9+OQ7cSmVmH1zlGcAXou3JeaMW4cVblYPAi8IdVVGRmShnhMHtA8Cb9GqtevIX4avAOaqc\nEm17LfBUVV4TbevFOeErOe+fXsUlfxiOt6SzD+WT1cT++/5yXunlCdFjVgCbE30u00CfPw6D2GdV\n44zhKoBUhB3AE3BVTfvogr63aqWX1wPfBU4TkQ0i8mrg3cCzReRu4Jn+erj/A8D7gFeJyEMicnrK\n0wbu8gdcI5l9L+7Dip8jrwe/k4jzwer07VnApnOyqjBGcNUmnwauqLM9tey3XlzG24z9HL6cYY6V\nmMbFHsaYWryd8k60EIedkXL/isxelQO+jSdQWbCwi4GJXooJeLOcffI4WIo7+SbfT4g0YrHPy+K3\nEeb979vbT2npxWpMUZ7Xh23Jz/McXI36x9mx5lEc7gmfQziOks6+WoxTawdtMsJBlcO4fsaFOFF/\nhGyxj8tuu66DNtfZq+pLM266MOP+JxR9YdWZMq9mddCGD2SC5o7unC0KO3tZJ730LBR2H5slrqGc\n7AFc/XE9TJA+pW8avRQb4ViEIK5pg3DC9AD1I0cWM7X4Dv8ag9FgnA2ki308EVrMVtzcNr9MbJ+g\nf3cfxb7sAzS22Mk4ruxvQASJHOoy4CdE78c71aSzH/HXs3jEP9cwvfsHKBLhOIIoxqSJ/QXAt1TZ\nJScv+BmTY8f57eH7miy9bDTGSXbQJjtnAyHKGcDtg6wRtCGvD8/dTxd10LZ71ktorJJmLjj7IgfJ\nYqYWgy5MVj4FRnDuYi/1i2M9zr4ZJ9UgrmnznxcV0WwWTA+zf+lmVY741xnCieZ3KejsPUHsK529\nE8ai1ThbqT/GiUeMxt+Zpbj528N8N5BeqRaOkyxKzr5nKvSjFKGosz/ftxMWTH+NIz0hfhql/BdI\nox20oR+mqrP3hE7aIPZZHbSx2O+I/pvYFyTT2YtUHbgzL5w9U8NLmByDqDIlWnwCahB7EUQkdf/U\nst/KYhwRHpMy4VRRgrimiX1Dzl7WSS8LDvdyYDx8wcP+WYLLe09KmV47y9lvwfUxJUsvJ+iZHKZ4\nZr+t4H3TCKWfyXLAZbhSyilKc7j3ERkh/z2qVj9fcvYLDw7TmLMvmxvHjy14GkHsp4fupmcyvIdR\n3P4tVGfv30vRDtqks08T+9jZx5l9mM8+6FMc48T/TewLkjeCdj3krrfZw3xw9vtWHM3kKHixF+FY\n3IjHQCz21Vzjc4Evp2yv1dnHMc7HqH8xk2rOvpEYZzHTQwdhwSP+evhShjjk17jcNSbP2Q9QKfa7\nWDA9TPHMvi6xj+Zh2YnbX/F3ZilOqO+g9GulF8r6xAZxE4TlLTUYnP0ICw6N0Lizjz/PU3Cdm05s\nJ46/l979cQnvFopX4wwDqjpz4sqLcZLOPi3GiZ19XI2TF+OYs6+DVGfvKyAeQ/4Xo5f54OynB4/2\nzj7EOGsorz8Pmf0eqovj5aSXidXj7MeiaXKzIqZq5Dn7RmOcMQ4OH6JUjho7+52Ui2OoFllASSxi\nQmSQjHEmkCOLKTYyNoh9PTHOME6sp3D7K+6kXeafN34/Zc6e6nk9hIFVCw4Og/ZT2wIm1TL7R+EW\nBnfsXPMgILJOBsh29lkxTt4o2IoO2igdyIpxwuNr6aCN/zcSn7aMjhV7nNAvzLgtkObs8zp1Ooni\n1ThHepYyVXL2wHG4n+bhAC8U43gx+w1gaUrsUst+i2Oc5TghqnfwU+zskxlvL64srt6IaJSpxUco\nfTHD/glxSJnY+7YcSJkOAUqRQWWMIzPT+1Y7MTUS48Sjd4s4+yB6ocqpWl4PYWBV/+5hYHcNAxOL\nZPbx6F2ArUyNwu5VS6gU+2qll3GEAxmll77a7yClz6WWGKdaB605+zrI6qANM/DliX1w9nEUlNep\n00nUUmd/FNNDB3CTcQlwvN8eRCNV7EU4JeE2L8CVvm2icn3OWvZbHOOEGGS2nD3U/0UaY2r0CCXX\nmevsyc7roSQuSWc/CQg9B4rMj1N3jEO5qyzi7JMxTrWySwiZ/cCuMH1FUcKvjeS2TLFXZZqpxYd5\n4II1uOMuxGRQHuOkHZNHU8zZQ/lxXSTGSeugDe1Kiv0RXMeyiX1Bspx9EbEPzj4+YXSTsz9EoYNE\nx5ke3I87uEZIF/sQ48QRwaeBx0fXL8MtXPEwbqRxzAQwKuukyKC4uBpnjW9Xo84+q/QS6v8ijTI5\neoSS64wz+x24Ts3TEm1Jy+shI8bx7neCwR37qR7PDFB/NU6qs49Wr9qF+1xD+Ww9Mc424DgGdk0h\nhTtnoSSKMdWcPUwPHuTA+ImkxzhT+AnoUoo04hp7yBf7nbi55xfiTorxSSJ+fOzsk5l96OCOT7jb\ngZercsT/gqCBX6AtoVPEPq2D9hzc0O+57Ox3UETIFhwZ5+CiKUqLaASxDye1MHFVMsYZx58QvMN/\nEdliH76sRcpg4xhnDfAz6hf7Is6+3jx0jMkxpfTe9uKOjWGcI7uX8g7arM5ZcOKyJ6ODcxdD2w7Q\nRGcvwogIl0ebspz9UcAOX1oaRw5JZ18kxnkEOI6h7WF64aKkOfuZWM4fe6eTFPtDg/s50nM8GR20\nquzxz50c2Z0X48Sll1DSg6W4DuLplPanOfsQ6R3xt/dTmgQNVTQx+rrj3X2niH2ZwPhyuDNxkzYV\ncfZTuIWCF1CbQ20n/Th3UP0AkcNjHByZpDT74nF40fBuInQqJsV+LLq+Alioyt24NVVXxy8x41CL\n/SpKxjjfo/4Yp1o1DtQv9qNMLhHKxX41MOG/xFtx65SG95wX49wH/EXGbRMMPzJJjoh7Z9iHE5Mi\nMc7ZwDuj67Gzj0U95PVQbnqSzr5ojNPD4I4QoRTlG1ROmxL/UlsN7FZNPOfhvr0405HVQQtuce/j\nKCdMMxwoEuNk5fXh8SPAtHfph3ybwnET9vfJuIGLaXT8KNqOFHucC3gYd/BVdfbROpgDerVO4oY/\nd/qUxwOkiL0IZ4pwUdk9FxwaY2pkP+VifzvugBwG9vp9sB8nXgv9z9ZFlIRyZl4U0p09FP9V1Etp\nkeVTcQOUGolx8ursD9GIsz8wHov9PuBYvGj6fXYfJXef6exVmVLlbzJeZxfDWw6SL+JBfCcoFuOM\nUt6vEsouoTQ4DEp5PZSfBOIO2kIxjhe63Qxtn6YGZ6/KLar8T2JzHONURjgAh3t3seDwSkpin8zs\nwYn98YlHJscLHAJ6vNlLjXHIrsTB338x5aZgKLo+ifu+hFWq0uj4+XE6QezTOmjPxi2TVm0qheDs\nk8/TDbl9cPbJA+QZUJo8CoCF04uYWrwf94V4NO7A2oQ7QGdW3vHiFQ66ICjx//BlzxL7Wpz9lH++\nx+CcfVq2WoSwMlSWsw8nlXoY5cCShZR/iY+l3BXGUU6es89jgqFth6u0c4DSPivi7Edx5ZzxL7PQ\nuZjl7OPvQD0xDsA2BrcfojZnn0Z1sdcFO5AjR+OOzUdId/YPUSn2ZatN+eM+5O5Zzj6rcxachsRi\nH74n8aJhxeMRAAAgAElEQVQsTwZu9b8I07AYpwBpmX0s9nmr3QTnl3yebsjtwwCd5AEyQLJSZsH0\nIg4u3otz9ufi3E4Q5uSXOEQ54f3Hzj4W+7IYx1N1v0Wrgx32bTiCm4QtfGFqJXb2aaWXjazvOcaB\n8aTYH0d5+WQs9nmZfR67GN4G+Utj9vt2FF3VKnwOocN1FaW1A5LOPoh9iDMXUvolMUXxGAfgEYYe\nOUJtmX0aU7hfmUKW2KNbkSMrfbv2UzpRxTNXpjn7tGUTg1mIHwslZ38G+c5+EeXHyaGofyaI/S0Z\njwcT+0KkufdTcJUS88HZp4l9+YRkCw+OcGBsD07sz8N9AYJoZIl9eP9pYr+R+p19Ly7bVJwg3Ocv\nh87jWqnm7HfQSGZ/YEkPlTFO7OzvoSnOfruQvzZscPYHcIJcbe3h8DmEE38s9klnvw1mHG74zqQ5\n+yJz7m9jcEfov6kbP5vkEdx3NF3sF05vRg4v8691CGb66+JO1rTMPm0lrXD8pDn7y3Gzwf5bRnOD\nUYlLdOOT/gHgqZjYN0yaoK/CCVI1sU86+26aMiHL2feTdPYLDw4xOb6b0jKQD1Hu7GOXE8ov88T+\nYeCYlNilyH4Lk8/h2xAyzLBOaxkinCHCX+SMLi3L7BNtKotxRHiTSMX0BtkooxwYD5ETlDrRks7+\nZH+5fmc/sGMBBZy9F+RkiWwaeWKf5eyh9As3rfSymLMf2AmNO3soTdiWLvZ9ezfSOzmM6zCf6Xej\nemaf5+yT1TiP4L4zz1HNzNvTMvu4lPQArm8qT+ytg7YAaYK+GndgF3H28U+tbpoyYcbZJwRugPKs\nFnqmhti3bIJSuVns7JMuJxnjVGT2fn72/VSWtBXZb3lin1aRcy7wZ8A/ZGT6Q8B+n4Ueonza7WSM\n8zu4Kq2CyBhTo0e8y4TSfmp+Zj+4M8yHnkU/JREqEuWM4URohT9RHkMVZx/dFgtmHZn9rlDV1ihT\nuLEdG0nLyxdvfJCBXVA6sYTYtqbM3pPl7D8FPFo1mqqhkqTY76Fc7Cf9fZLTW8eYsy9AWQetXxYt\njKibEXtfYfLmxGNj4el4Zy/C20Vm5vEfwH35lPL4IryHFQCyToSFUwPsWbmLcrGvltmPUhqEBeXO\nHtKjnNT9JsJqEd7rr8b7fBtupR7IjnEW4yZKezJuCbok8UCmZJQTnH14D0eTH5WUozLG5Gj8pQ37\nKXb2G3Cdy/005Owneqq0LQzSASf2M85ehL8RqZivaBS3CttKnKDvUZ15fDxdQtLZhxNBcvrvopn9\nJgZ2Cc1x9lPAnwPvTJ2CYviRrQzsCH0/4f5Jsd+KM0SxkKbN3nkQ975DfxIAqkyqpg6kikl20KY5\n+9szavQDVo1TgGQH7Urc0mGHKRfwo4D3JbLObnP2z8FNCAWlDrukIwidlCG3H0IXHGF6ZB+l0X95\nmX2YDG0U1yGVJfZZo2jT9tu5uNkyoVzs3wx8xF/OcvajOEH9IK6/IUnI7KFS7MMydiP+V8FyahF7\n0TGmysQ+CMSM2PtOuA249Vnrd/ZuAZOizn6mIsefZP4At49jxnD9Visoj3CgfFBVmrMPMU49mf0/\nc9RdO2iesz+IG8iXxgRD2w9R7uzLYhx/kniI8tw+y9mHyeKKzukTPzbZQZsU+59WeQ5z9gWYxtXI\nLvTXQ4QD5WI/hFvfNu68jDto63b2Ivy1CHlLKDaLIUriGzrskgdJeA8hqx110/RyACdS05Rn9skD\nP0xzPIZz72mll+ArckR4scjMCvZZ++1kSiehGbFXZSIMFScjs8eJ2m6yyz2LOPthSvOVFBJ7N6hO\nFzM5Fjv1tBgHSlHOKuoT+1307gu/DLJIOvsQ4zzB33Zy4v6jOLFfSakPKxAbm6PJdvY1xziqTLLg\nyCKa4+wPAH+VU664y8c4ec4eKnP7LLEfSTyuKMkO2mSMsxNKaxJnYGJfjaiCIIhJfGDHAh52ZFwy\nmFV6WauzfxqV85rPBkG0INvZhyHbQezHmB6eppRrn6PKNopV42wk39n/HnA9Jdeetd/WUD7gJe3n\nbFaMM+qft2LUrieOTpLll3FmH3415LnnxPPKIQ73V4txwFXkXAVcBHy84PPHTPjVqmrJ7MNJ+Hyc\nsCaPv1jsYwMEvoPWj55eSnlJYeyO64lxwB0DzXD2v0X+/pygf3fcPxC+wwuhbFqKmdze/8Jrttin\nddDGJuFK4Noqz2EdtAWJRT3+yZp09uH2QOzs44moas3sh6glHiiICL8twlujTcOUDogsZ9+PG5Id\nfsGMMrUozNmPKrf77VmZfYhxxnD7MUvsg1v6ICWXmbXf1pDi7BNkxThFnH1FjBNNLxDEPpxIin5O\noxxZGOZWCQSBSDr7u3Hv8VmqbCj4/DG76Jms9qsjy9lfAHyCSrEfw1WwZMU4g7iTwKbEfD0NxTh+\nmpFwgm4IVW6POsfT2E3PZA8LpsPJagq3X5JRTOzs+4DDKfl5iGKaJfYzJsHn/nnvA8zZFybupM2L\ncaBcMGJn/32YiSNqdfazIva4L/KjoutxjJPn7O8ndvZTiw9T2XEYO/u0GGeUfLH/OG407vcpCU8R\nZ58n9nnOfgJYKFJRhRLHOPF8KgtxndcTlM+XX/RzGuNI7z7Kf45nOfsPAY9V5Z6Cz51kwi/jV1Nm\n76PLpwDXke7s78FNk3A85TFOKL08DieEMWkdtEXns8c/77RerfWIZk3o1XoYObKXty37D79pRuwT\nd41r7bOWVgzOfirltmqkjaBNzuJZDRP7gsQRTFaME77kcRQQO/vPAc/zP21rHVQ1W2J/BuU99MkY\nZ4rKycsGSDr7ybF4TvbAbqpX42Rm9n6uly2U13xX7DffIX4s1Z19XjXObu/U0iqAhknvoA2Zc6hy\nCL8aijv7Q337qexo20tymmJlv4/G6mU3Cw4NIoeLOvtfAxfiZnbdhOv8CxVBYZ/34T7j7f5+ac7+\neFzEQcptSWdfNMZpiqsvjDDBwEQ45iZJd+cbcMcgpEc40HiME3fQbqByv1bDxL4gsajHzn6Kysw+\n1dn7tS3vwM0tU+t0CXG80kzOCM8bZY1xjJNVjfMAsbN30/SWib3vGJ3Gld5l1dlvpFTJknT2gThS\nSNtvx+PEKQx9DyKcZCewSISjRMrGRiwmfwK2Rb4NUDl74XT0fsKCFcWd/aHBMrH3J5wT/DiDpqFX\n62F0wQEGduUdQ7Gzfz/uu3cD8K1ERRA4wQ0nyE2446gis8d9NklnnxxBO40zRUMUW2YwXqCjFcTH\nXJaz3477hQOzI/bTRH1HqnxFlT+s8Tk+i+uj6FhyxV5ErhWRLSJyW7RtXERuFpG7ReQmERmLbvtT\nEfmViNwpIs+poR1FM/ukM4ydPbgvz2/SAc7eT5u7kpJrD1MRJ519WozzALGzdzM3pgnUBG6wTVbp\n5TZcrX0fxcR+NzAi6yQ+Ltbg6r3DnN+pzt53Ht/q2/5wNIBqlJKYl3XS+sFCcbvSnH0Q++W4eKu4\ns58eqlhUQ5XtBR9fG7pwDwO78uqsQx8Nfh3Zy3AjMm/0t8cjeWN3vRl33KRV46SJfezsp6MCiP05\nVTExrXX25d/VLGcfnxDSpkqAxp19eP268IuY1Fry2VKqOfsPAxcntl0F3KyqpwJf89cRkUcBL8Zl\n1BcDHxSRor8cJnGrwSzACWSW2N9NdmYP7uz6Ah45dTcFnb3/ydxLk8TeL5QOzo2FAxBKgh6XXmZl\n9g8BR/l5QsZyxH43TuyzMvuJ6HrWz/iZGEev1sP+ueJc/WScEIXsNyvGQZVzVVmEE5uwP0MHLVQ6\n+2GcCIXOr7R5ycP+ORp3IilajRPWAKj7C1wTKhN+Ob8sQh+Nu7uLji5V5b/9pngkb+yuN+E++9ht\nF83sYxGrpRKn05z9TsrFfjZiHGjVsdImcsVYVb9NZeXCpbgOJfz/F/rLLwCuV9VpVX0A17mUNogm\njdBBu4zykYJJsf8VsCpyjWXOXpX7gK38649OodKhZhHEo2Gx96Mg7/Edb2fgnG4Q8uT/LGff77dt\nxwlcmKY3rf57Aud4s2KcCUpOv4izD88Z/ypagxOiEKllin3EI8Ay/znliX3SRcall8kYZzlO7Is7\n+6nFaQthzxYT9E8U7aBNIxb7pLN/OOEaCzt7v60Wse9EZ78XGIyWYJyNGCe8/pylp/pdKliuqmEk\nZ5iYC5zD/H50v6yZFdMIHbTJErNknf0jfts4TgyTzh7gAQ4uDgL4Glkn+fnsG08d45t/DqMPnSXr\n3nFFwfam87LnnMTPXnEMZ37yajaffTZHenaze9Upsu61V/DaJ63kh2+E/t3Hybo3XMFZ/zHE8/7g\nRfzwjSehskbWvcu99uOuHeNp776E7715ijM//WrgLA4ctZBsZy9UxjjjuM92P6Wl+PpIP2GUDd3H\nOa2XyjpxJYhPet/TWXHrt3hw7QKe/LcvZvNjT+DXTzha1v1x9r56wj8f5oz/ehkrf/prvvz3h7ns\nFb8t64AXvewkNp3zxJnHvuYpq7j1lci633fXn/S3yzjmJxfKuo+P83tPWM5Pfr+fS3/vcm74D2Hh\n1OmsvPUetp1xvKx7Y5HP6clMLgkdlLOP6A5O++KwrHtWedsmVg8zunEfz/jfj2fhwUlZ9970tl/2\nO6v49bnny7o3X8GLXvYENp89IuvecgUXX7mKHWumZd0flR73Z73C564dRPQELn7TU2TdjtI6wxe+\n7VQO9/ZxaGCEka2LZd0/XsHZH+phwXSPrHt9kf0W+rxaxS7gWbJOlAv+4niO9PYxOTZc1tZrgM9+\n+ADPuOY13P+MM9l65lGy7q3l7+X8vzyOQwOLONIzIOveVNv3+Pm/fw4Png/HfvcsWffBRjRgg16t\n32rg8bNKPWI/g6qqiOTlVJm3icg1pWufWgSXD1A5UjCeE3sIJ2ThJLKdysweXFndElwn2DOqvonh\nbSOc8iUY2XICcEnV++dx1K9WcMqXYHDnb7P6+yMs+vUGtp05AlzC2ANLOPWLR+iZXIrKJZz85YX0\n7Xk2a246kUMDgzOvffKXhxjdsJbTP9/L2IMXAg/w0FMfTXZmD5UxTlh6T0XYi4vG9mZkigdws032\n+trlD+F+kT0GgOO/tYZjv7+Zvr19jGx6NsdMjzG8bWnuvjrt8/0sue+ZLJzexWk3Hpm574lfH6d/\n9wkz18ceXMppX+ifuX7Kl5YyuuE8YAVjD45yyn+7/XLqFw9zZOEyjv3+iSzesCz3tUvs465LN5G/\nHkLz0AWfYHjzRSiXzPzuPDjcz8YnPY/Rz9zAid84lZ4D+8lq+5qbR+ndfxJwCStuPYnFD48Dl/Co\nzxzF3uWTZY9bOA2nfkmRQ0cY2vGssuc56atrODQ4xOHeXoa2jwB9nPrFfhZOk/nalXyi+l2axo3A\nK4BLOPkrq5keGuZQf+n7EDjtRujb8zyO+fEyxh6sPP5O/spqDg4vQvRIxW3VOP5bJ9C3B1bc+iiK\nx4Rp/BhomtiLyFpK5eSNo37l3Kw/4ATgtuj6ncAKf3klcKe/fBVwVXS/LwNPzHhOLb+uHwN9Gejr\nQf81cdsU6ADoP4BeCfpl0Ev8bdeCvjZx/78DfXO19xXd/zGgCvr5oo/Jea7LQe8AfQj0XtBHgx4G\nFdDz/faHQAdBp/xjXg/6L9Fz7AcdAv1weG+gk6BDKa/3Ed/2JdG2ZX7bvf76V0FfA7oxp907QcdT\ntotvzyLQ2/37eQnoJ6vsh4+Bvhz0dNC7ou3HgG6Jrl8E+pXo+sdBr/CXHwd6q7+8we/HR4HemfGa\nvaDHJLa9FfSvG/1ca/j8D4L2R9fP8p9FL+gHQP8g57GD/nNeCPom0H+o8lq7QH+csv0NoP8M+h+g\nr/TbbgH9n1bthwb237tAbwL9ZsptPwY9F/SNoB9Iuf19oN8BvaGO173cf06Xt3sf5LcTbeTx9ZRe\n3gi80l9+JaWFhm8EXiIifSJyIm4BkmrzSQRCXHMS6aVkYd6R/ZTnvmnOfielMq0iDCX+N8IYbj3W\nA7i64LsorUw/jBt4NEJ5Z91MZu9/wYSqjU3ASt9pXda5F5Hm7EMuvyu6fgz5oyeTUU5gHDcHexiJ\nWlNmT2X+uwVY4sdCkHJ71sLRe3GVRWFt0DQuBj6d2BbXtreCeM4aKE0hfRRVMnt15aBbceZqlOpR\nyn4qvytxG5IdtEUmQWs3k6R30EKpk7Yjq3G6gWqll9fjxOs0EdkgIq8G3g08W0TuBp7pr6Oqv8TN\nHf1L4H+AN6g/HRUgdNA+FbeeaUw4EYTBN7HYp2X2O6hN7IcpzZjXKEG8bgDuUReLhEFBQzixCxN6\nhS9+3EEbFlA/jBf7cF9NL5vbjRtaHh/gU5SWDARmYpy8L3vWuqhx52qYpKqI2G/DCV38ePz72kJp\nDEHZ7aTX2Yf3sJXyRTuSHAc8LppQD1ov9sn2LYv+F2nLLbglOYvMTXOAfLGvt4O2nUyRPeXBLlw8\nax20dZKb2avqSzNuujDj/u8C3lVHO0Kn62Mp7+QNt4V5R/bhMvtQ5ZOX2RdlCCdOzXD2wZFdDzNz\nrIRqkmGcOwkDnCqcPeUngc24vC6eTiDJBIkDX3Ump69F7LMW04greOLBOkWc/QmUD6gKhPVvH6S8\nBh/Kq3FiZ7oPt1/jedyTrMZ9hqdSWhUprAbWKrKc/VKqV+NASexHYWYOpCz2kz7KMxinBdRXetlO\n8px9KNEcpryII2DOvgqdNIJ2LfBz1YqKkVjsizj7emKcZon9GK5j9FeqfNBvi8V+n/87inJnn5wv\nB0rOPm+O9d1k186HGKCo2KfFOLHY1+Psk2IO5VVayZNBNWc/iSvBS1vxahXuF83Z0bZudPbnUKzW\nPcvZh8q25MI+3RDj5Dn7OMaZLWff1JHVnUYnif1TSe/JTor9ZkqjS7Ocfa1i/wjNdfYxwbmHXyZ7\ncWKf5ezD9vA+qzn7NLGPnX2RzD4rxokHYlUdVBURMvtkTANu6oVj/OU0Z58cQYtvw1YfZcVTaMSs\nwh0/50TbOiGzP0J9zr5ajPNB4Ds5bYj7PLrF2YdBVWn7KcQ4eWI/lPHYapizbyGTuIPz2xm3xZl9\n7ELjlaoCtcY4oX6/mZl9TNLZB7FPLoINlTHOSvKXyttN+oG/l3Jnv4LGY5xaOmi3kd5BC+WzYyad\nfTzrZdLZh7EdWbn9auALtN/Zx2K/DDdYKjj7akK0Abd/T6eKs1flI6ozy1TGNDqCtp1M4mLOep19\n/L8WTOxbyAFcTX6aU0lm9vEsjR0Z4yS2ZcU4k4nbIYpxVNmLe28ryBb724HPpGxPZvZ91BfjxHOg\n1+rsKzpoPbHY5zn72JleB3zeX64Qex/rrAK+CJwdxTztcPZx25bi+g+Cs89tiyqKmwHzGOofxRov\nOB4+p4/gSqE7nXAyrNZBm1WNk/XYasyLDtpOEftJ4DbViqkZoOQoQ4wTi31ajLMLGPUli0UYwp0g\nevxcNI2QFuME5x7an+bsw/tJitMm3EyIqWKvygZV3pNy0x7KY5z4fxp5MU49zn4Hbl+MUyla8VTI\neZn9jDNV5duqM52uac5+MXBElV/5doa5zzvB2d/p/xeJccBFOVC/2IfMPl7H9Zt+33Q64bOq1kFr\nzr4OOkXsfwi8N+O2ZGa/HzeqtocUZ69uuti9pItXGuEXQ16lR1HSYovQAZuMcWJnP+RPTsmf+ptw\nYw9q7Tj6FKUS1uCCmlGNU8jZ+xLLXbgTVSPOPu110j6neFrskHtDZzn7om1pVOyTi5d0E3nOfjZj\nHOugbRXqli/LWquyLLP3P3VD9JHm7KG2KCf0Beyn8dy+SIxT5ux9p+N+KgdbgcvtaxZ7VT4WOeGi\nYl+kGqeoswcXja0hPbMP8ymlOfu00suYpKBC+ZxKP6V9Yj9zIvJRUhD7Wp39pGpdHY2Q3kHbLVRz\n9tU6aLMeW43wmHr3eVfQEWJfhUncB9xH+bJhi0jP7KG2ipz4F0PduX20QEie2O/HHahLKT+wwvvJ\ninGySi+LUETsi8Q4tWT24HL71VQ6+zjGqdfZp4l9mFPpPtoX48QnorDM3cPU5uzvws0s20gbajkp\ndxJTif8xRZ19vdU4h1RTtWTO0C1iv4SSq4eSOGY5+1rFPnSczoiICCMivLGGdo7gHFnygIkz+7QY\nB8rFPj5Y63L2CYJY51VjzJazh/RYS0QYpnqdfZpLSxP7OMbZjtu/UKBTtMnEpZfLcPsgdFYXqcbB\nT2NycwNtCCflfuaWsw/rI89GB+1B5niEA90j9uOUu9tqzn4nxcsvs5z9KcDbamhn1hD3rMw+zdkn\nxWmTf95GDsRmOftaBlWBEzlIOHt/wt6Cc+P9lH+ucellVoyT5eyD2O+gJPbtHFS1DHhE3doMB3Gf\n46y3xe/fsKZqt4l9Zmbv9+MR8qdLSH1sNXxhyLOq3rHL6RaxP4p0sW+Gs8/K7IdreA7InrwqGeNU\nc/ZJsYfZF/umdtB6spw9uNz+ZErrrAaaEePEzr6Qm24isbNfSmkfhBNfq9pyAHfMdVuMk+fswZm4\nKS2tbBbTiLNHlR/V87huolvEfpz0pfeamdmXxTh4gY5maKxG1qjHIPZxXJQUobwYB1oj9s2OcYLA\npf3k3or75ZTM87NG0MbUEuO0s/RyGaV9sI3SBHetIBwvc8bZe3aRvWh6Q2I/H+gmsa/F2efGOCJc\nJDKzEnxWjBNcfurziPB2ER4XbcqazyRk9nGMA+nOPi3GgQbE3vchVJsbZTY6aLfhlphMm61zC9XF\nPsvZJ6ckgMoYZ0lUytquDtqks2/1LwyYm87exL5Oul3s63X2FwDP9ZfjOvs0sc96nouB/xNdz3L2\ncWYfYhwoVo2zHff+Gu08eh4ll5lGs6dLwL9eVq14cPbJ24uUXpZ9TiL04060WwGiaaVH6Sxnb2Jf\nHXP2s0i3iH2zM/sllCZTy8vsw33TGAAuEeGx/npejLOI8mocKJDZe1e8hcZKL1Hl64lsPHn7FKBe\nOGMacfYPk32CmcnsE9vrqcY5Btic+AWxHV/L3+Jyujxn38qTziQwnfeZdyjVyid3kV1VNi9q5Ruh\nW8Q+mdk3Wo0zTmkBjazMvpqzHwD+E3iHv54X44zjvnyHo/eRdPZhUFXyYN1Ea8rCyqIcP24gnhun\nJmevym24xW3S2IKb7z7N2dfaQXsSrrY+Zjsu2mn18PdOcvbd5upDJdEUFuPMCt0i9iM039knxb6m\nzB73pf574FkiHEt+jHM0pYO0sLP33I8Tr9kmGeUM4k5Q4WRaq7MPJW1pbMUde0lnH09fXLSDdg1u\nZsmY7TjH32qx7xRnf4DuFb08sbcYpwEanfirFQRHlCX29WT248DRPrYQVQ6KsB/nzgMj0X3TGMB9\nmb8OPMM/Nm3loDDr5P7oOqRn9qEzNeZltMalJSty4ggHyoW40faEqXmTJ8e4DUWdfZbYm7PvTsK4\nhDTM2TdAtzh7SI9xssS+2tw447j3fjwlEU7L7LeTn9lP4ubgP5/sOvt9Gf+zBlWViYIqB1uUvSYr\nctLEviZnn0OYmz7p7Cdw+xHy58aJq3E6SewPUJrUbjmlk9r9lMpoW9WObhW9PGf/EG7xmzRM7KvQ\nTWKfdPZLgMMZQrgfWCCSOYvlElwH4snR86Zl9hvId/aTuNWRLiAjxvE5/QHyY5zQidvq6pGYZIyT\nFPta1qCtxnbc+gVp8wgN+BlNizr7k+kcsQ/O/kRgp6p7f6r8XHWm+qsVTDIHnb0qn1blTRmPM7Gv\nQreLfeoB7U8AYU6SMrzrWoybZ/wkyp19Uuw3kp/ZT+IWEDkaOI3s1YVCaWe4DMVKL1vJTkrLBUL5\nkoTQRGfvT4CPkD6Vwm7cibNqNY7vRO40Zz+Im3Xzpy1+7WQ7ulX07qLU11EYf0yFZSuNFLpZ7MdJ\nj3ACD+DEPMkoTsQ24oQiiG+a2D9EirP3J4w+SkO3v+NfK6uufG94HVUO4r6IWYOq2nWwfgy4UoSF\n/nqas29WjAMuyknbXyHKKdJBuww4mNIR3NYYB7cO7i1V7jvb7ehKsVflN1RTF1IvQvhuGSl0k9gn\nM/sl5Iv9vTgxTzKO68DdRHmMk8zsR3AxTpqz78fNcBkipLB2blWx9+yj85z9V3Cu+nJ/PS67hNoH\nVVXj/1JZMgml2Q2LxDhprh6c2K+gvc6+3WLfrTFOI3yT7lhrty3ULfYicqWI3CYit4vIlX7bY0Xk\neyLycxG5UUTS5luplSxnv5j8AzpL7JfgIoswfXBeZr+R9Mw+KcpB7LNinL2Ut38v6c6+bWLvT1xX\nA1d7dz+rzl6V1/ta/CS7KO7s88R+IS3el370rgLn0V6xz6tombOocrGfHdNIoS6xF5FHA78LnAs8\nFni+iKwBPgS8TVXPAj4LvLUJbcwSe2jc2VfL7LM6aJOi/GPcYtdZriJMgBa4ifKqgk6IcQBuxrnC\n88kuvZztRTFCjNOIs9/h/7fjix+Opw1teO3AfHX2Rg71OvvTgR+o6qSqHsb9fPpN4BRVDS73q35b\no6SJ/RRO6PMO6HvIFvvg7AfIz+w34ibVEhEGRWaqVULnLDBTHvkbOSWSZTGOKr+rWtYJtR8nbiO0\nL8YJ7v7/Ao8nIfZ+cJXiB1vNYjPiGKda6WWes4f27MsDwC1tnqqgazN7Y/aoV+xvB84XkXERGQIu\nwU0ze7uIvMDf53Lg2Ca0sSKz91+kPdTn7JdQcvZQHuMk6+x34r40I8BbcDEH1B63JGOcMqJ1dZfW\n+LyzQViwO+nswZ1kZ3ue9GoxzgFgMKcSB0rHRrucfTsrccDE3kihrhG0qnqniLwHF0fsA24FDgOv\nBd4vIn8G3EjOASci10RX16vq+oy7po2gBfeFzhOdR4AeEcZVZ37WQ8nZJ8U+zdnvw0+ZCzyB0ojI\nAWqbrybZQZvGHlwFSbtLx24B/gTX5mQHapinqG0xjipHRGYipbQae1RREbbTRmffhtdNtsFinC5H\nROvg544AAA1FSURBVNYCa5v1fHVPl6Cq1wLXAojIu4CHVPUu4CK/7VTc1LpZj7+m2OtwSIRDpIt9\n3uNUZMbdx2K/BNiiyl6Rsvr3uH67DxDcySqMxj0b+K6/b63O/iaq1w7Hs0u2k1/iBgXdTbqzh9l3\n9qvIdvbgPqtrcCffTRn3aZfYfwL4RhteN+ZW0hejMboIb4DXh+sicnXmnQvQSDXO0f7/ccCLgE+I\nyDK/bQHwv4F/bqRxEZNUiv1eqotOWpQTnD04oQiO+yDul0AP3tX7eGUHzkEeR8n5l2X21VDlBtWZ\nip0sOkLs/TiAO4Gn0B6xr9ZBC+5Y+C3gORmLo0CbxF6Vd6q2dGqEtDbcocq/trMNRufRSJ39Z0Tk\nF7i45g2quht4qYjcBdwBbFTVjzShjQBfoHLmx2qZPaSLfcjswYn9fpjJzUP5ZYhw8Pd9Ji6mCpn+\nbJRIxlMJt5uf4kYFJ8U+vOfZXF4vFvssZ/8+4FmqM6tTpdEuZ28YHUkjMc4FKdveD7y/oRalvha/\nk7J5D+UdqmncAzw5sS129psp/8UQopxY7MPK87dQLvbNnmN+D264dysX28giZM5pzn62F8XYhavG\nySzxVOXvCzyPib1hRHTDCNos6nX2oc4e3LiAr0e3TeCcf9LZn4orSQwxzmw5+8k2l+wFssS+FRNs\nFXH2RfgYbtyAYRh0x3z2WVSrxgEn9qeKsCDKdsMIWlS5KXH/Tbhh9ocod/bgxP43/OVZE/smP2e9\n/BxmyltjpmiN2Oc6+yKoljq2DMOY+85+A66z8V98XTaUO/skm3ErWI1QGg27A+cwf0wpxqmpg7Yg\ne+iMvB5V9gKvpHIO9lY4+124NYcP53S+GoZRI90u9rnC48XiUuAs4D0iDOB+zWTVvAdnn8zsb8c5\nzlmPcZr8nHWjykf9jJ4xrXL2S7BBQYbRVLo9xqnamanKHhEuwc2T/QXcohJZufgmnLPfTknsv4qb\nx2YfMOx/IcxGB21ycrROZNadvSqTIhyc7dcxjPnGnHb2AT+C9oPAu8mOcKAU48w4e1W2q/JtP6Ph\nEVzH4Ww5+46IcXJohbMHF+WYszeMJtLNYv9N4Loa7v9+3AydyYUuYtJinJhQhz/nY5wMWiX2E5jY\nG0ZT6doYR5X7cQs5F73/dhH+GTgj524Vzj5BWOBkNjpoJ2h+NNRsWrW26S7cLyjDMJpE14p9nazD\nzc6ZRcjsR3DL5iUJM2POhrNfj1sGsZNppbMfa8HrGMa8oZtjnJpRZa8qd+bcZQdO6MfJjnGC2DfV\nhasypcodzXzOWaBVzn6iRa9jGPOGeSX21fClmltwK1hlxTizldl3A9ZBaxhdiol9JWEh8jxnPxuZ\nfTdgzt4wuhQT+0o24+ZTrxbjzEext2ocw+hSTOwrCYthWIxTySStEWGLcQyjyZjYVxLmg9mbctus\nddB2Ca109hbjGEYTMbGvJM/Zz/cYZzv5g9KaxcMteh3DmDfMtzr7IgRnnxfjzNcO2i8CX2rB63yN\n9q/jahhzChP7Sqo5+0XMU2fvJ5CbzSUJW/o6hjGfsBinEotxDMOYc5jYV7IZ+KFqajVIXI0zHzto\nDcPoUizGSeBF/okZN+/DTacwQOdPR2wYhjGDOfva2IebN+egLZlnGEY3YWJfG0HsLa83DKOrqFvs\nReRKEblNRG4XkSv9tvNE5IcicouI/EhEzm1eUzuC/bjFsC2vNwyjq6hL7EXk0cDvAufiVn96vois\nAd4L/Jmqng38ub8+l9iHE3tz9oZhdBX1OvvTgR+o6qSqHsYtEXgZbmHuUX+fMdxIyLlEqLM3sTcM\no6uotxrnduCdIhLy6+cBPwSuAr4rIn+DO5E8uSmt7Bz2+/8m9oZhdBV1ib2q3iki7wFuwrndW4Aj\nwL8Df6iqnxWRy4FrgWenPYeIXBNdXa+q6+tpS4sJA61M7A3DmFVEZC2wtmnPp6qNP4nIO4GNwHtU\ndbHfJsAuVR1Nub+qqjT8wi1GhF7c1LvrVXlGu9tjGMb8oVHdbKQa52j//zhcXv8J4B4Rebq/yzOB\nu+t9/k5ElWnc1Lvm7A3D6CoaGUH7GRE5Cid+b1DVCRF5HfABEenHlSe+rhmN7DD2YWJvGEaXUbfY\nq+oFKdt+TPZUA3MFE3vDMLoOG0FbO/sxsTcMo8swsa+dfdgIWsMwugwT+9qxGMcwjK7DxL52LMYx\nDKPrMLGvHXP2hmF0HSb2tWNibxhG12FiXzt7KM2RYxiG0RU0ZbqEml+0S6dLABBhFbBPlV3tboth\nGPOHRnXTxN4wDKMLaNvcOIZhGEb3YGJvGIYxDzCxNwzDmAeY2BuGYcwDTOwNwzDmASb2hmEY8wAT\ne8MwjHmAib1hGMY8wMTeMAxjHmBibxiGMQ8wsTcMw5gHmNgbhmHMA0zsDcMw5gF1i72IXCkit4nI\n7SJypd/2SRG5xf/dLyK3NK+phmEYRr3UJfYi8mjgd4FzgccCzxeRNar6YlU9W1XPBv7L/3UFIrK2\n3W1IYm0qRie2CTqzXdamYnRimxqlXmd/OvADVZ1U1cPAN4HLwo0iIsBvA9c33sSWsbbdDUhhbbsb\nkMLadjcghbXtbkAGa9vdgBTWtrsBKaxtdwNSWNvuBjSbesX+duB8ERkXkSHgecDq6PbzgS2qem+j\nDTQMwzAap6eeB6nqnSLyHuAm3ALctwBHoru8FPhE480zDMMwmkFTliUUkXcBD6nqv4hID7AROEdV\nf51x/9avhWgYhtHlNLIsYV3OHkBEjlbVrSJyHPAi4In+pguBO7KEHhprsGEYhlE7dYs98BkROQqY\nBt6gqrv99hfTXR2zhmEYc56mxDiGYRhGZ9PSEbQicrGI3CkivxKRt7fytaM2HCsi3xCRX/gBYX/k\nt4+LyM0icreI3CQiY21o20I/IO0LHdSmMRH5jIjcISK/FJEntrtdIvKn/vO7TUQ+ISL9rW6TiFwr\nIltE5LZoW2YbfJt/5Y//57SwTX/tP7uficgNIjLa7jZFt71FRI6IyHgr25TXLhH5Q7+/bvdFKC1r\nV8bnd56I/NDrwo9E5Ny626SqLfkDFgL3ACcAvcCtwBmtev2oHSuAx/nLI8BdwBnAe4G3+e1vB97d\nhrb9MfBx4EZ/vRPadB3wGn+5BxhtZ7v88XMf0O+vfxJ4ZavbhCsvPhu4LdqW2gbgUf547/XtvwdY\n0KI2PTu8FvDuTmiT334s8GXgfmC8lW3K2VfPAG4Gev31Ze3eV8B64CJ/+bnAN+ptUyud/XnAPar6\ngKpOA/8JvKCFrw+Aqm5W1Vv95b3AHcAq4FKcsOH/v7CV7RKR1cAlwIeA0IHd7jaNAuer6rUAqnpI\nVSfa3K7duH6iIV/5NQT8utVtUtVvAzsTm7Pa8ALgelWdVtUHcF/M81rRJlW9WVVDWfQPKI2HaVub\nPH8LvC2xrSVtymnX64G/8vqEqm5rZbsy2rQJZ7AAxoCH621TK8V+FbAhur7Rb2sbInIC7kz6A2C5\nqm7xN20Blre4OX8HvJXy8QrtbtOJwDYR+bCI/FRE/k1EhtvZLlXdAbwPeAgn8rtU9eZ2tikiqw3H\n4I73QLuO/dcAX/KX29YmEXkBsFFVf564qd376RTgAhH5voisF5EndEC7rgLeJyIPAX8N/Gm9bWql\n2HdUT7CIjODm7rlSVffEt6n7ndSy9orI84GtqnoLJVdfRqvb5OkBzgE+qKrn4AbQXdXOdonIGuBN\nuJ+uxwAjIvKydrYpjQJtaGn7ROT/Aw6qat5gx1lvkx9x/w7g6nhzzkNauZ96gCWq+iSc8fpUzn1b\n1a5/B/5IVY8D3gxcm3Pf3Da1UuwfxuV0gWMpPzO1DBHpxQn9R1X1c37zFhFZ4W9fCWxtYZOeAlwq\nIvfjylafKSIfbXObwH0+G1X1R/76Z3Div7mN7XoC8F1V3a6qh4AbgCe3uU2BrM8reeyvpvRzfNYR\nkVfhIsIros3tatMa3In6Z/54Xw38RESWt7FNgY244wl/zB8RkaVtbtd5qvpZf/kzlKKamtvUSrH/\nMXCKiJwgIn24evwbW/j6wMwkbf8O/FJV/z666UZcRx/+/+eSj50tVPUdqnqsqp4IvAT4uqq+vJ1t\n8u3aDGwQkVP9pguBXwBfaGO77gSeJCKD/rO8EPhlm9sUyPq8bgReIiJ9InIiLi74YSsaJCIX41zq\nC1R1MtHWlrdJVW9T1eWqeqI/3sNo+y3talPE54BnAvhjvk9VH2lzu+4Rkaf7y88E7vaXa2/TbPR0\n5/Q2PxdX/XIP8KetfO2oDU/D5eK34ub0uQW4GBgHvup35k3AWJva93RK1ThtbxNuCusfAT/DuZ7R\ndrcL17H3C+A2XEdob6vbhPsF9mvgIK4v6tV5bcBFF/fgTlYXtahNrwF+BTwYHesfbFObpsJ+Stx+\nH74ap1VtymqXP44+6o+rnwBr2/z5vRr3S/YHXq++B5xdb5tsUJVhGMY8wJYlNAzDmAeY2BuGYcwD\nTOwNwzDmASb2hmEY8wATe8MwjHmAib1hGMY8wMTeMAxjHmBibxiGMQ/4f8+gFEPGmgZYAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efeb58888d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# - add some signal with a signal to noise ratio of 2, \n",
    "# sig_mag = 1.0*bold_std\n",
    "\n",
    "sig_mag = 1.0\n",
    "c = [sig_mag*bold_std, 0, 0]\n",
    "Y = X_.dot(c) + bold\n",
    "time = np.arange(X_.shape[0])\n",
    "plt.plot(time, Y, time, X_.dot(c)+100.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efeb839c990>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD/CAYAAAAUnaZMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD3NJREFUeJzt3V+oXNd1x/HfT2PHrZNQIxxdS7JaFWqjCgIytGqpW3wf\nHJMmYOeh2HUpEY1bQkOaUgq1k5bKT6kbqCnNQwmpHZQ2UWra2ijYxpJCZATCTkXsxq0syQILZFW+\nsp04/eMH68/qwxyZ0VT3zty5Z886s+/3A5d7ZmZrZrH31dpr1pyZcUQIAFCnNdkBAADKIckDQMVI\n8gBQMZI8AFSMJA8AFSPJA0DFWk/ytj9q+6jtV2zf3/b9AwDG5zbPk7fdk3RM0u2STkv6V0n3RsTL\nrT0IAGBsbVfy2yWdiIiTEXFO0rck3dXyYwAAxtR2kt8o6dTA5dea6wAACa5q+f5G9n5s8zkKADCB\niPBy/03blfxpSZsGLm9Sv5q/TER0/mfnzp3pMYzzg3bNyrrPSpySZLvTP7MQ46U4J9F2JX9Y0k22\nN0v6T0n3SLq35ccAMCNsa82abp+pffHixc7HKEkXLlyY6N+1muQj4rztz0p6RlJP0iPBmTXAqraS\nKnRaZiHGSbVdySsinpb0dNv3O23z8/PZISDBrKz7rMTZ6/VmokqehRgn1ep58mM9oB30kttTcwWS\ngb/Ndl1zzTXZIVTj3XffVUzwwmvrlTwAXEIRko8kD6CYmtsgs4IkD6AYKvl8bLMAUDEqeQDF0K7J\nR5IHUAztmnwkeQDFUMnnYwUAoGJU8gCKoV2TjyQPoBjaNflI8gCKoZLPxzYLABVLqeTZ3dvDB2qh\ny2jX5KNdA6AYCrp8JHkAxVDJ52MFAKBiVPIAiqFdk48kD6AY2jX5SPIAiqGSz8c2CwAVo5IHUAzt\nmnwkeQDFkOTzsQIAUDGSPABUjHYNgGJ6vV52CKseSR5AMZxCmY8kD6AYXnjNxwoAQMWo5AEUQ7sm\nH0keQDEk+Xwk+RnHf6J28U1b7aInn48VAICKUckDKIZKPh8rAAAVo5IHUAzveM1HJQ8AFaOSB1AM\nZ3/lI8kDKIYXXvOR5AEUQyWfj20WACpGJQ+gGCr5fCR5AMXQk89HkgdQDEk+HysAABUjyQNAxWjX\nACiGjzXIR5IHUAxn1+SjXQMAFaOSB1AMZ9fkS0nyfMVae3g6jC4jyeejkgdQDEVIPrZZAKgYlTyA\nYmjX5CPJAyiGdk0+kjyAYqjk87ECAFAxKnkAxdCuyTdxkrd9UtJ/Sbog6VxEbLe9VtI/SvoZSScl\n3R0Rb7cQJ4AZRLsm30pWICTNR8QtEbG9ue4BSfsi4mZJ32kuAwCSrHSbHX4udqekXc3xLkmfWOH9\nAwBWYCU9+ZC03/YFSV+JiK9KmouIheb2BUlzKw0QwOzio4bzrSTJ3xoRZ2x/SNI+20cHb4yIsH3F\nD6l58MEH3zuen5/X/Pz8CsIA0FW88Dq5N998U2+99daK72fiJB8RZ5rfb9h+XNJ2SQu2b4iI122v\nl3T2Sv92MMkDqBcvvE5u3bp1Wrdu3XuXjx8/PtH9TLQCtq+1/cHm+P2S7pD0kqQ9knY0w3ZIemKi\nqAAArZi0kp+T9HjzVOwqSd+IiL22D0t6zPZ9ak6hbCVKADOJdk2+iZJ8RLwqadsVrv+hpNtXGhSA\nOpDk8/GOVwDF0JPPxwoAQMVSKnmewrWHr1JEl1HJ56NdA6AYCrp8JHkAxVDJ52MFAKBiVPIAiqFd\nk48kD6AY2jX5SPIAiqGSz8c2CwAVo5IHUAztmnwkeQDFkOTzkeQBFENPPh/bLABUjEoeQDG0a/KR\n5AEUQ7smH0keQDFU8vlYAQCoGJU8gGJo1+QjyQMohnZNPpL8jKNSahfftIXasM0CQMWo5AEU0+v1\nskNY9UjyAIqhnZiPJA+gGF54zccKAEDFqOQBFEO7Jh9JHkAxJPl8JHkAxdCTz8cKAEDFqOQBFEMl\nn48kD6AYevL5SPIAiqGSz8cKAEDFqOQBFEO7Jh9JHkAxtGvykeQBFEMln49tFgAqRiUPoBjaNflS\nkjxfsdYeng6jy0jy+ajkARRDEZKPbRYAKkYlD6AY2jX5WAEAqBiVPIBier1edgirHpU8AFSMSh5A\nMZxdk48kD6AYXnjNR5IHUAyVfD62WQCoGJU8gGJo1+QjyQMohnZNPpI8gGKo5POxAgBQMSp5AMXQ\nrslHkgdQDO2afCR5AMVQyedbcpu1/ajtBdsvDVy31vY+28dt77V93cBtn7f9iu2jtu8oGTgAYLRR\nlfzXJH1Z0tcHrntA0r6I+JLt+5vLD9jeKukeSVslbZS03/bNEXFx+E7Z3dvDVymiy2jX5FsyyUfE\nQdubh66+U9JtzfEuSQfUT/R3SdodEecknbR9QtJ2Sc+1GC+AGUJBl2+SnvxcRCw0xwuS5prjDbo8\nob+mfkUPYJWiks+3ohdeIyJsL9UvoJcAABM4duyYjh8/vuL7mSTJL9i+ISJet71e0tnm+tOSNg2M\nu7G5DsAqRbtmclu2bNGWLVveu/zkk09OdD+TJPk9knZI+svm9xMD13/T9sPqt2lukvS9iaICUAXa\nNfmWTPK2d6v/Iuv1tk9J+nNJD0l6zPZ9kk5KuluSIuKI7cckHZF0XtJnglM/gFWNSj7fqLNr7l3k\nptsXGf9FSV9caVAAgHbwjlcAxdCuyUeSB1AMST4fSR5AMfTk87HNAkDFqOQBFEO7Jh8rAAAVo5IH\nUEyv18sOYdWjkgeAilHJAyiGs2vykeQBFMMLr/lI8jOOSqldfNxSu/j7zMc2CwAVo5IHUAztmnwk\neQDF0K7JR5IHUAyVfD5WAAAqRiUPoBjaNflI8gCKoV2TjyQPoBgq+XxsswBQMSp5AMXQrslHkgdQ\nDO2afCR5AMVQyedjBQCgYlTyAIqhXZOPJA+gGNo1+UjyAIqhks/HNgsAFaOSB1AM7Zp8KUmer1hr\nD0+H0WUk+XxU8gCKoQjJxzYLABWjkgdQDO2afKwAAFSMSh5AMb1eLzuEVY9KHgAqRiUPoBjOrslH\nkgdQDC+85iPJAyiGSj4f2ywAVIxKHkAxtGvykeQBFEO7Jh9JHkAxVPL5WAEAqBiVPIBiaNfkI8kD\nKIZ2TT6SPIBiqOTzsc0CQMVSKnl29/bwVYroMto1+WjXACiGgi4fSR5AMVTy+VgBAKgYlTyAYmjX\n5CPJAyiGdk0+kjyAYqjk87HNAkDFqOQBFEO7Jt+SSd72o5I+LulsRHy4ue5BSb8r6Y1m2Bci4unm\nts9L+pSkC5I+FxF7C8UNYAaQ5PONquS/JunLkr4+cF1IejgiHh4caHurpHskbZW0UdJ+2zdHxMUW\n4wUwQ+jJ51tym42Ig5J+dIWbrrRyd0naHRHnIuKkpBOStq84QgDAxCbtyf+B7U9KOizpjyPibUkb\nJD03MOY19St6AKsU7Zp8k6zA30r6WUnbJJ2R9FdLjOXTswAg0bIr+Yg4e+nY9t9J+nZz8bSkTQND\nb2yuA7BK9Xq97BBm1qFDh3To0KEV349HfVSt7c2Svj1wds36iDjTHP+RpF+MiN9qXnj9pvp9+I2S\n9kv6uRh6ANtU9y3io4bRZWfOnMkOoRobNmxQRCz7lexRp1DulnSbpOttn5K0U9K87W3qt2JelfRp\nSYqII7Yfk3RE0nlJnxlO8ABWF86uyTeykm/9AankW8U+ii47e/bs6EEYy9zcXPuVPLqPSqldbJrt\n4u8zH+c3AUDFqOQBFMN58vlI8gCKoV2TjyQPoBgq+XysAABUjEoeQDG0a/KR5AEUQ7smH0keQDFU\n8vnYZgGgYlTyAIqhXZOPJA+gGNo1+UjyAIqhks/HCgBAxajkARRDuyYflTww4MCBA9khjGVW4jx4\n8KDWrFnT6Z9ZiHElbS+SPDBgVpLnrMT57LPPynanf2YhxpU8IyLJA0DF6MkDKMZ258+wmYUYV4Lv\neAWAGTHJd7xOPckDAKan3ucoAACSPADUjCQPABUrnuRtr7W9z/Zx23ttX7fIuJO2f2D7BdvfKx3X\nwON+1PZR26/Yvn+RMX/T3P5vtm+ZVmxDMSwZp+152z9u5u8F23+WEOOjthdsv7TEmC7M5ZJxdmQu\nN9n+ru3/sP3vtj+3yLjU+Rwnzo7M50/Yft72i7aP2P6LRcZlz+fIOJc9nxFR9EfSlyT9SXN8v6SH\nFhn3qqS1peMZesyepBOSNku6WtKLkn5+aMzHJD3VHP+SpOemGeMy4pyXtGfasQ3F8GuSbpH00iK3\np8/lmHF2YS5vkLStOf6ApGMd/dscJ870+WziuLb5fZWk5yT9atfmc8w4lzWf02jX3ClpV3O8S9In\nlhg77Q+62C7pREScjIhzkr4l6a6hMe/FHxHPS7rO9tx0wxwrTmn683eZiDgo6UdLDOnCXI4Tp5Q/\nl69HxIvN8f9IelnShqFh6fM5ZpxS8nxKUkS80xy+T/3C6YdDQ9Lns3nsUXFKy5jPaST5uYhYaI4X\nJC02aSFpv+3Dtn9vCnFJ0kZJpwYuv9ZcN2rMjYXjGjZOnCHpV5qnmU/Z3jq16MbXhbkcR6fm0vZm\n9Z95PD90U6fmc4k4OzGfttfYflH9PPTdiDgyNKQT8zlGnMuaz1be8Wp7n/pP24b96WWRRcQSb4a6\nNSLO2P6QpH22jzYVV0njvklgeNec9psLxnm870vaFBHv2P51SU9IurlsWBPJnstxdGYubX9A0j9J\n+sOmUv5/Q4Yup8zniDg7MZ8RcVHSNts/JekZ2/MRcWBoWPp8jhHnsuazlUo+Ij4SER++ws8eSQu2\nb5Ak2+slnV3kPs40v9+Q9Lj6LYrSTkvaNHB5k/q791Jjbmyum6aRcUbEf196mhcRT0u62vba6YU4\nli7M5UhdmUvbV0v6Z0n/EBFPXGFIJ+ZzVJxdmc+BeH4s6UlJvzB0Uyfm85LF4lzufE6jXbNH0o7m\neIf6u85lbF9r+4PN8fsl3SFp0TM0WnRY0k22N9t+n6R7mngH7ZH0ySa2X5b09kD7aVpGxml7zu5/\nVJ3t7eq/m/lKvbxMXZjLkbowl83jPyLpSET89SLD0udznDg7Mp/Xuzmzz/ZPSvqIpBeGhnVhPkfG\nudz5nMYHlD0k6THb90k6KenuJrgNkr4aER9Xv9XzL03cV0n6RkTsLR1YRJy3/VlJz6j/AscjEfGy\n7U83t38lIp6y/THbJyT9r6TfKR3XJHFK+g1Jv2/7vKR3JP3mtOO0vVvSbZKut31K0k71zwbqzFyO\nE6c6MJeSbpX025J+YPvSf/IvSPrpS3F2ZD5HxqluzOd6Sbtsr1G/uP37iPhO1/6vjxOnljmffHYN\nAFSMd7wCQMVI8gBQMZI8AFSMJA8AFSPJA0DFSPIAUDGSPABU7P8AblSkRer58B4AAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efeb57f39d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Add a trend to X_, call this new matrix X\n",
    "trend = np.arange(X_.shape[0])[:,np.newaxis]/X_.shape[0]\n",
    "X = np.hstack((X_, trend))\n",
    "plt.imshow(X,aspect='auto',interpolation='nearest')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, recall the GLM and how to test a contrast "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking again at our time series, we should add a trend..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0.51    -0.4    100.285] [  0.825  -0.243  99.189   1.89 ]\n"
     ]
    }
   ],
   "source": [
    "# - estimate the parameters of the model X_ and X using the OLS (see below)\n",
    "\n",
    "pinvX_ = npl.inv(X_.T.dot(X_)).dot(X_.T) # this is to match the formula we have seen, \n",
    "pinvX = npl.inv(X.T.dot(X)).dot(X.T) # in practice do not use this, use npl.pinv(X2)\n",
    "\n",
    "beta_ = pinvX_.dot(Y)\n",
    "beta = pinvX.dot(Y)\n",
    "print(beta_, beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note our estimates 'beta' $\\hat{\\beta}_1$, $\\hat{\\beta}_2$, $\\hat{\\beta}_3$, etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, we need to _test_ if our signal of interest $\\hat{\\beta}_1$ is significant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need 2 things: signal, and standard error of the signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First, what is our observed signal of interest ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.825176582124\n"
     ]
    }
   ],
   "source": [
    "# express our signal of interest as c.dot(beta) - print the signal magnitude. \n",
    "c = np.asarray([1,0,0,0]).reshape(4,1)\n",
    "signal = c.T.dot(beta[:,np.newaxis])[0,0]\n",
    "print(signal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Second, what is our observed standard error ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall the linear model for one observation i :\n",
    "\n",
    "$$ Y_i = X \\beta + \\epsilon_i $$\n",
    "\n",
    "and in matrix form:\n",
    "\n",
    "$$ Y = X \\beta + \\epsilon $$\n",
    "\n",
    "with $ \\epsilon_i $ independant gaussian variable with variance $\\sigma^2$, such that the vector $\\epsilon$ has variance $\\sigma^2 I$,  $I$ the identity matrix.\n",
    "\n",
    "Our OLS estimates are:\n",
    "\n",
    "$$ \\hat\\beta = (X^tX)^{-1}X^t Y $$\n",
    "\n",
    "Our fitted data : \n",
    "\n",
    "$$ \\hat{Y} =  X\\hat\\beta = X(X^tX)^{-1}X^t Y = M_X Y$$\n",
    "\n",
    "And our residuals:\n",
    "\n",
    "$$ r = \\hat{\\epsilon} = Y - \\hat{Y} = Y - M_X Y  = (I - M_X) Y = R_X Y $$\n",
    "\n",
    "$R_X$ or $R$  when there is no ambiguity, is the __residual forming__ matrix. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Let's now  compute $ \\mathrm{var}(c^T \\hat\\beta) $. We have\n",
    "\n",
    "$$ \\mathrm{var}(c^T \\hat\\beta) = \\mathrm{var}(c^T (X^tX)^{-1}X^t Y) \\qquad  (A)  $$ \n",
    "\n",
    "Let's call $ X^+ = (X^tX)^{-1}X^t $.\n",
    "\n",
    "We have assumed that $\\mathrm{var}(Y) = \\sigma^2 I$\n",
    "with $I$ the identity matrix.\n",
    "\n",
    "$$ (A) \\;\\; = c^T X^+ \\mathrm{var}(Y) X^{+T} c  = \\sigma^2 c^T X^+ X^{+T} c = \\sigma^2 c^T (X^T X)^- c $$ \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How do we know $\\sigma$ ? We don't, we have to estimate it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will need the residual sum of square:\n",
    "\n",
    "Our residuals, the estimated noise, are $ r = \\hat{\\epsilon} = Y - \\hat{Y} $, with $ \\hat{Y} = X \\hat{\\beta}$\n",
    "\n",
    "The residual sum of squre is :\n",
    "\n",
    "$$ \\textrm{RSS} = \\sum_{i=1}^{n} (Y_i - \\widehat{Y_i})^2 $$\n",
    "\n",
    "But the mean sum of square is not the simple division by the number of time points, we have to correct for the number of estimation we have done when computing $\\widehat{Y}$, so that will be \n",
    "\n",
    "$$\\hat\\sigma^2 = \\textrm{RSS} / (n - p)$$\n",
    "\n",
    "Where p here is 4, and n = 196.\n",
    "\n",
    "This estimation is also called the **Mean Residual Sum of Square** (MRSS), and $\\hat\\sigma = \\textrm{MRSS} = \\textrm{RSS} / \\it{df} $  where $\\it{df}$ is the degrees of freedom $(n-p) = 196 - 4$, the number of observations (here, 196) minus the rank of the design matrix (ie the number of independant columns, here 4). The MRSS is the estimation of our noise level.\n",
    "\n",
    "Coming back to our t-statistics, we can now calculate the **Standard Error** (SE) of $c^T \\hat\\beta$ with $\\textrm{SE} = \\sqrt{\\hat{\\sigma}^2 c^T (X^t X)^+ c} $\n",
    "\n",
    "The t test is then : $t = {c^T \\hat\\beta}/{\\textrm{SE}} $\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102.523017949 0.625140353349\n"
     ]
    }
   ],
   "source": [
    "# - Compute the estimate of sigma: compute the residuals, take the sum of square, and \n",
    "# divide by n-p, \n",
    "# n is the number of observations, p is the number of (free) parameters that we estimate\n",
    "\n",
    "n,p = X.shape # ok because X has no linear dependencies for its columns !\n",
    "r = Y - X.dot(beta)\n",
    "RSS = (r**2).sum()\n",
    "sigma2_hat = RSS/(n - p)\n",
    "\n",
    "print(RSS, sigma2_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.038  0.019 -0.025  0.013]\n",
      " [ 0.019  0.036 -0.022  0.006]\n",
      " [-0.025 -0.022  0.043 -0.044]\n",
      " [ 0.013  0.006 -0.044  0.076]] \n",
      " 0.153752432606\n"
     ]
    }
   ],
   "source": [
    "# - what can you say of this compared to bold_std ? \n",
    "# Compute the standard error of the signal\n",
    "\n",
    "V_1 = npl.inv(X.T.dot(X))\n",
    "\n",
    "SE = np.sqrt( sigma2_hat * c.T.dot(V_1).dot(c))[0,0]\n",
    "print(V_1, \"\\n\", SE)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.36691724573\n"
     ]
    }
   ],
   "source": [
    "# - compute the t statistics, that follows a student distribution with n-p df\n",
    "\n",
    "t_observed = signal/SE\n",
    "print(t_observed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.82517658212413614, 0.62514035334885676, 0.15375243260563373)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# - \n",
    "signal, sigma2_hat, SE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How far from the t distribution centered on 0 ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - \n",
    "import scipy.stats as sst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.34762696233e-07\n"
     ]
    }
   ],
   "source": [
    "# compute the p of the t statistics, use sst.t(df) to define a t_variable\n",
    "# then find how to compute the p-value (see cdf or sf functions of t)\n",
    "df = n-p\n",
    "t_variable = sst.t(df)\n",
    "p_ttest = t_variable.sf(t_observed)\n",
    "print(p_ttest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  F statistics and model selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Review the F test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The F-test is best understood in terms of model selection.\n",
    "The question we are asking is : should all the regressor be in the model ? \n",
    "\n",
    "To answer this, we construct a reduced model X0 from X. From this model we look at the RSS of model X (full model) and model X0 (reduced). The residuals are always going to be larger for the reduced model (how do we prove this ? Let's see in class). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{eqnarray} \n",
    "F_{\\nu_1, \\nu_2} & = & \\frac{(\\hat\\epsilon_0^t \\hat\\epsilon_0 - \\hat\\epsilon^T\\hat\\epsilon)/ \\nu_{1} }{\\hat\\epsilon^T\\hat\\epsilon/\\nu_{2}} \\\\ \n",
    " & = & \\frac{(Y^TR_{X_0}Y - Y^TR_{X}Y)/ \\nu_{1} }{Y^TR_{X}Y/\\nu_{2}} \\\\ \n",
    " & = & \\frac{(Y^T(I-P_{X_0})Y - Y^T(I-P_{X})Y)/ \\nu_{1} }{Y^T(I- P_{X})Y/\\nu_{2}} \\\\ \n",
    "& = & \\frac{(\\textrm{SSR}(X_0) - \\textrm{SSR}(X))/\\nu_1}{\\textrm{SSR}(X)/\\nu_2}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "What are $\\nu_{1}$ and $\\nu_{2}$ ?\n",
    "\n",
    "They are the same things that we had for the t-test. The difference $ \\hat\\epsilon_0^t \\hat\\epsilon_0 - \\hat\\epsilon^T\\hat\\epsilon $ will have the difference of each term df. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's apply it in our simple case :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# compute the F test for regressor 0 in our previous model\n",
    "# what is the difference with the t-test\n",
    "\n",
    "X0 = X[:,1:]\n",
    "# plt.imshow(X0, aspect='auto')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# -  define a glm function that returns beta, Yfitted, residuals\n",
    "\n",
    "def glm(X, Y):\n",
    "    \"\"\" a simple GLM function returning the estimated parameters and residuals \n",
    "    \n",
    "    parameters:\n",
    "    -----------\n",
    "    X: np array\n",
    "        The desing matrix\n",
    "    Y: np array\n",
    "        The data\n",
    "    \n",
    "    returns:\n",
    "    --------\n",
    "    betah : np array\n",
    "            estimated beta\n",
    "    Yfitted: np array\n",
    "            fitted data \n",
    "    resid:  np array\n",
    "            residuals\n",
    "    \"\"\"\n",
    "    \n",
    "    betah   =  npl.pinv(X).dot(Y)\n",
    "    Yfitted =  X.dot(betah)\n",
    "    resid   =  Y - Yfitted\n",
    "    return betah, Yfitted, resid\n",
    "\n",
    "\n",
    "# define an F test function:\n",
    "\n",
    "def F_test(X, X0, Y):\n",
    "    \"\"\"\n",
    "    parameters:\n",
    "    -----------\n",
    "    X: np array\n",
    "        The design matrix\n",
    "    X0: np array\n",
    "        The reduced design matrix\n",
    "    Y: np array\n",
    "        The data\n",
    "        \n",
    "    returns\n",
    "    --------\n",
    "    F: float\n",
    "        the F statistics\n",
    "    nu1, nu2: floats\n",
    "        the degrees of freedom\n",
    "    \"\"\"\n",
    "    betah, Yfitted, resid = glm(X,Y)\n",
    "    betah0, Yfitted0, resid0 = glm(X0,Y)\n",
    "    \n",
    "    nu1 = X.shape[1] - X0.shape[1] # assume X full rank (no linear dependencies)\n",
    "    nu2 = X.shape[0] - X.shape[1]\n",
    "\n",
    "    numerator = ((resid0**2).sum() - (resid**2).sum())/nu1 \n",
    "    denominator = (resid**2).sum()/nu2\n",
    "    F =  numerator / denominator\n",
    "    return F, nu1, nu2\n",
    "\n",
    "# define a t test function:\n",
    "\n",
    "def t_test(X, c, Y):\n",
    "    \"\"\"\n",
    "    \"\"\"\n",
    "    assert c.ndim == 2\n",
    "    assert c.shape[1] == 1\n",
    "    betah   =  npl.pinv(X).dot(Y)\n",
    "    signal = c.T.dot(betah)[0,0]\n",
    "    RSS   =  ((Y - X.dot(betah))**2).sum()\n",
    "    p = npl.matrix_rank(X)\n",
    "    n = X.shape[0]\n",
    "    df = n-p\n",
    "    V_1 = npl.pinv(X.T.dot(X))\n",
    "    SE = np.sqrt( (RSS/df) * c.T.dot(V_1).dot(c))[0,0]\n",
    "    return signal/SE, df \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28.8038007225 1 164\n",
      "2.69525392467e-07 2.69525392466e-07\n"
     ]
    }
   ],
   "source": [
    "# compute the F test for regressor 0 in our previous model \n",
    "# use the F_test function\n",
    "# use sst.f to find the p value, compare it to the p_ttest\n",
    "# why is p_ttest half of the p_Ftest ?\n",
    "\n",
    "f_obs, nu1, nu2 = F_test(X, X0, Y)\n",
    "print(f_obs, nu1, nu2)\n",
    "\n",
    "F_variable = sst.f(nu1, nu2)\n",
    "p_Ftest = F_variable.sf(f_obs)\n",
    "\n",
    "print(p_Ftest, p_ttest*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - F test can be computed by stacking several one dimensional contrasts \n",
    "# into a \"contrasts matrix\" - no demonstration provided here but see reference\n",
    "# Poline and Brett, Neuroimage 2012."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.19458741707e-11 2.69525392467e-07\n"
     ]
    }
   ],
   "source": [
    "# - now, test A and B :\n",
    "X0 = X[:,2:]\n",
    "f_obs, nu1, nu2 = F_test(X, X0, Y)\n",
    "\n",
    "F_variable = sst.f(nu1, nu2)\n",
    "p_FtestAB = F_variable.sf(f_obs)\n",
    "\n",
    "print(p_FtestAB, p_Ftest)\n",
    "# is this expected ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's create  5 subjects with the same design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efeb31d8310>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFNJJREFUeJzt3W2MXOd53vH/Ra5JxaoVSnC7pEiCUlPSCIOidttQaN1A\nDsAGStKQDhBEFuCCsISggOoo+VKQbJFU+UJIAqKmRdGgaRyBcU1abFqrBFrUpFzISFCoshzGMUKr\nMosSFhnv0olDR0qcii93P8yRtGJJ7i5nzsvO/n/AYs+cOTP7cDh77T33Oc85qSokSSvbmr4HIEka\nn2EuSVPAMJekKWCYS9IUMMwlaQoY5pI0BVoJ8yQPJHklydeT7G/jZ0iS3pFJH2eeZC3wv4DdwHng\nS8BDVfW1if4gSdLb2qjMdwFnqupsVV0CPgvsbeHnSJIabYT5ZuC1BbfPNeskSS1pI8w9P4AkdWym\nhec8D2xdcHsro+r8bUkMfEm6BVWV661vI8xfBrYnuQf4Q+BB4KHrDOimT/L444/z+OOP3/D+5Lr/\nnmUZd+fvHXfcMfYYXn/99bGfQ5ImHuZVdTnJJ4HPA2uBT03rkSzr1q3rewiSBLRwaOKSfugE2ixD\nOHXvpk2bxn6Oubm5CYxE0mrRZZtl1bAylzQUhvkYDHNJQ7Fiw3wIO0ANc0lDsWLDfAje85739D0E\nSQIM87EY5pKGwjAfw/r16/segiQBPYb5uP3qSfTMx2XPXNJQWJmPwTaLpKEwzMdgZS5pKAzzMViZ\nSxqK3sJ83J73EKbzuwNU0lBYmY9hZsaXT9IwmEZjsDKXNBSG+RjcASppKFZsmHtuFkl6x4oN8yGw\nZy5pKEyjMViZSxoKp/OPwR2gkoZi0TBP8hvAjwMXquqvN+vuAp4FtgFngZ+uqovNfQeBh4ErwGNV\ndaKdoffPylzSUKxZwjbPAA9cs+4AcLKqdgBfaG6TZCfwILCzecy/SbKUn7EizczMjP0lSZOwaJpU\n1W8nueea1XuA+5vlw8ALjAJ9L3C0qi4BZ5OcAXYBL05ovINiZS5pKG61NJytqvlmeR6YbZbv5t3B\nfQ7YfL0nmIbp/Ia5pKEY+3N+VVWSmyVr/6nbEsNc0lDcapjPJ9lYVXNJNgEXmvXnga0LttvSrJtK\nnjVR0lDcapgfB/YBTzbfn1uw/kiSpxm1V7YDL407yKGyMpc0FEs5NPEoo52d70/yGvCLwBPAsSSP\n0ByaCFBVp5McA04Dl4FHq6XmttP5JekdSzma5aEb3LX7BtsfAg6NM6iVwjCXNBQe6DwGe+aShsIw\nH4OVuaSh8NwsY7AylzQUVuZj8ERbkobCMB+DlbmkoegtzJ3OL0mTY2U+BsNc0lAY5mMwzCUNhWE+\nBnvmkoZixYa50/kl6R0rNsyHwDCXNBSG+RgMc0lDYZiPwWt4ShqKVZ1GQzglwCSOlx/C/gMY/4/b\nlStXxh6DtFqt6jDXZI3bdvrud787oZFIq49hrokxzKX+GOaaGPchSP1ZymXjtgK/CfwVoIBfq6p/\nleQu4FlgG82l46rqYvOYg8DDwBXgsao60c7wV76h9LsnwaN7pP5ksSBIshHYWFW/l+QvAV8GPgp8\nAvijqnoqyX7gzqo6kGQncAT4QUYXdX4e2FFVVxc85zDSZ0oMJczvvffesR5/9uzZyQxEmmJVdd0K\ncCnXAJ0D5prlN5J8jVFI72F0oWeAw8ALwAFgL3C0qi4BZ5OcAXYBL475b9DAWZlL/VlWkzPJPcCH\ngP8JzFbVfHPXPDDbLN/Nu4P7HKPwV0uG0qqxZy71Z8m/fU2L5T8CP1dVry8MkKqqRVonw+gDqFVW\n5lJ/lhTmSd7DKMg/XVXPNavnk2ysqrkkm4ALzfrzwNYFD9/SrNOUM8yl/izlaJYAnwJOV9WvLLjr\nOLAPeLL5/tyC9UeSPM2ovbIdeGmSg9YwGeZSf5ZSmX8Y+Djw+0lONesOAk8Ax5I8QnNoIkBVnU5y\nDDgNXAYeraEcbjFAQ5nOPwme313qz1KOZvkdYM0N7t59g8ccAg6NMS6tQFbmUn88/EATY5hL/THM\nNTGGudQfw7xnQzlGfBLsmUv9Mcw1MVbmUn8Mc02MYS71xzCfAkNp1dhmkfpjmGtirMyl/hjmmhjD\nXOqPYa6JMcyl/hjmPZum6fyGudQfw1wT4w5QqT+GuSbGylzqj2GuiTHMpf4Y5j0byjHik2CYS/0x\nzDUxXgNU6o+/fZqY9evX9z0EadUyzDUxtlmk/hjmU2AofXfDXOrPTcM8yW3AF4H1wDrgP1fVwSR3\nAc8C22iu/1lVF5vHHAQeBq4Aj1XVifaGryGxZy7156a/fVX1F0l+uKr+PMkM8DtJ/h6wBzhZVU8l\n2Q8cAA4k2Qk8COwENgPPJ9lRVVdb/ndoAKzMpf4s5YLOf94srgPWAn/CKMzvb9YfBl5gFOh7gaNV\ndQk4m+QMsAt4cbLDnh7TNJ3fHaBSfxYN8yRrgN8Fvg/41ar6gySzVTXfbDIPzDbLd/Pu4D7HqELX\nKmBlLvVnKZX5VeCDSb4X+HySH77m/kpys/JyGDNa1Dp75lJ/lvzbV1XfSfJfgL8FzCfZWFVzSTYB\nF5rNzgNbFzxsS7NOq4CVudSfxY5meT9wuaouJvke4O8DvwQcB/YBTzbfn2sechw4kuRpRu2V7cBL\nLY19KgzlsMJJMMyl/ixWmW8CDjd98zXAp6vqC0lOAceSPEJzaCJAVZ1Ocgw4DVwGHq2hJI1aZ5hL\n/UkfWbtIj13LNJS/l9/4xjfGevy2bdsmNBJpelXVdT/Ou8dKE2NlLvXHMJ8CQ+m7G+ZSfwxzTYxh\nLvXHMNfEeA1QqT+GuSbGylzqj2EuYBjndxnKeWomMY5JfEq5fPny2M+h1cMwl1owiU8phrmWwzCX\nWuD+A3XNMNdgDKVFMgnuP1DXDHOpBYa5umaYSy2wzaKuGeZSC6zM1TXDXFNlKH13w1xdM8ylFhjm\n6pphLrXAMFfXDHOpBV4PVV3zHafBGMp0/klYv35930PQKrOkME+yFngZOFdVP5HkLuBZYBvNZeOq\n6mKz7UHgYeAK8FhVnWhj4NKQ2WZR15Zamf8co+t6vq+5fQA4WVVPJdnf3D6QZCfwILCT0QWdn0+y\no6quTnjc0qAZ5uramsU2SLIF+DHg14G3PsPuAQ43y4eBjzbLe4GjVXWpqs4CZ4BdkxywtBLMzMyM\n/SUtx1LeMf8C+CfAHQvWzVbVfLM8D8w2y3cDLy7Y7hyjCl1a1FCOEZ8EK3N17aZhnuQfABeq6lSS\nj1xvm6qqJDf7DRrGb5fUIXeAqmuLVeZ/F9iT5MeA24A7knwamE+ysarmkmwCLjTbnwe2Lnj8lmad\ntKpYmatrN+2ZV9U/raqtVXUv8DHgv1fVPwSOA/uazfYBzzXLx4GPJVmX5F5gO/BSO0OX/n9Jxv6a\nhLVr1479JS3HcveyvNUyeQI4luQRmkMTAarqdJJjjI58uQw8WkNpYkodss2iri05zKvqi8AXm+Vv\nA7tvsN0h4NBERietULZZ1DWPf5JaYJira4a5BmOapvN7cQp1zTCXWmBlrq4Z5lILDHN1zTCXWmCY\nq2uGuQZjmqbz2zNX1wxzqQVW5uqaYS61wMpcXTPMNVWG0qpxBqi6ZphLLbAyV9cMc6kF9szVNcNc\naoFhrq4Z5hqMaZrOb5ira4a51AJ75uqaYS61wMpcXTPMpRYY5uqaYa7BGMox4pNgmKtrhrnUgpkZ\nf7XUrSW945KcBf4UuAJcqqpdSe4CngW20VwHtKouNtsfBB5utn+sqk5MfujScDkDVF1bavlQwEea\na3++5QBwsqqeSrK/uX0gyU7gQWAnsBl4PsmOqro6yYFLQ2abRV1bzmfBaxuae4D7m+XDwAuMAn0v\ncLSqLgFnk5wBdgEvjjdUaXFD6bsb5uracirz55NcAf5tVf07YLaq5pv754HZZvlu3h3c5xhV6NKq\nYc9cXVvqO+7DVfXNJH8ZOJnklYV3VlUluVk5M4xDDKSOWJmra0sK86r6ZvP9W0k+x6htMp9kY1XN\nJdkEXGg2Pw9sXfDwLc06aUUYyikBxm33DKXlNJTXc9plsf+sJO8F1lbV60luB04AvwTsBv64qp5M\ncgDYUFVv7QA9wijwNwPPA3+tFvygRap4SUxPmE/iyJ4333xz7OeYFlV13f/YpVTms8DnmjfGDPCZ\nqjqR5GXgWJJHaA5NbH7Q6STHgNPAZeDRGspMDkmdm8T+A8N8cYtW5q38UCtzaVHTUpnfeeedYz/H\nxYsXx36OaTFOZS6pB+OG8VA+ELszuBuGuaRWGebdMMwltcpj7rvhqyypVZ6nphuGuTSlhrID1Ksu\ndcMwl9QqK/NuGOaSWmVl3g3DXFKrPJqlG4a5NFBDmDQ0CYZ5NwxzSa2yzdINw1xSq6zMu2GYSwPl\ndH4th2EuqVWGeTcMc0mtMsy7YZhLapU7QLthmEtTaijT+a3Mu2GYS2qVlXk3DHNJrfLcLN1Ys5SN\nkmxI8ltJvpbkdJL7ktyV5GSSV5OcSLJhwfYHk3w9yStJfqS94UsaupmZmbG/tLilvkr/EvivVfVT\nSWaA24F/BpysqqeS7AcOAAeS7AQeBHYCm4Hnk+yoqqstjF+aWtMynd/KvBuLhnmS7wV+qKr2AVTV\nZeA7SfYA9zebHQZeYBToe4GjVXUJOJvkDLALeHHyw5c0dO4A7cZSKvN7gW8leQb4G8CXgZ8HZqtq\nvtlmHphtlu/m3cF9jlGFLmkVMsy7sZQwnwH+JvDJqvpSkl9hVIG/raoqyc0+Ew5jXrGkztnz7sZS\nXuVzwLmq+lJz+7eAg8Bcko1VNZdkE3Chuf88sHXB47c06yQtg+dm0XIsGuZNWL/W7MR8FdgN/EHz\ntQ94svn+XPOQ48CRJE8zaq9sB15qY/CShs8doN1Y6uefnwU+k2Qd8L+BTwBrgWNJHgHOAj8NUFWn\nkxwDTgOXgUdrKCWCpM5ZmXdjSWFeVV8BfvA6d+2+wfaHgENjjEvSmIYynd+eeTd8lSW1ysq8G4a5\npFYZ5t0wzCW1yjDvhmEuDdS0TOf3rIndMMwltcrKvBuGuaRWGebdMMwltcow74ZhLg3UtEznt2fe\nDcNcUquszLthmEtqlZV5NwxzaUoNZTq/J9rqhmEuqVVW5t0wzCW1yp55NwxzSa0yzLthmEu6oaGc\nEmASvfshHOq5Zs2a1sZgmEtSR8bdf/Dmm2/e8D7DXJI6Mm7LaawwT/IB4LMLVv1V4BeAfw88C2yj\nuWxcVV1sHnMQeBi4AjxWVSduceySNDXa3H+Q5fSBkqwBzgO7GF0X9I+q6qkk+4E7q+pAkp3AEUaX\nmdsMPA/sqKqrC55nGPOMJa0aQzi9waZNm8Z6/NzcHFV13eb/ctssu4EzVfVakj3A/c36w8ALwAFg\nL3C0qi4BZ5OcYRT+L97K4CVpWrR5PdTlPvPHgKPN8mxVzTfL88Bss3w37w7uc4wqdEla1dqcDbvk\nME+yDvgJYP+191VVLdI66f/zjST1rM3ZsMupzH8U+HJVfau5PZ9kY1XNJdkEXGjWnwe2Lnjclmad\nJPVmCMeZL7cyf+ONN3jjjTeWtO1ywvwh3mmxABwH9gFPNt+fW7D+SJKnGbVXtgMvLePnSNJUWm7P\nfMOGDWzYsOHt2/Pz8zfcdknPnOR2Rjs/f2bB6ieAY0keoTk0EaCqTic5BpwGLgOP1hB2I0tSz3pv\ns1TVnwHvv2bdtxkF/PW2PwQcGnt0ksQwpvNPwm233dbaczsDVJI60uakIcNckjqydu3a1p7bMJek\njgziOHNJ6stQLoE3LtsskjQFDHNJmgK9H5ooSRqflbkkjWkI0/kNc0maAoa5JE0Be+aSNAWszCWt\natNybhYrc0maAs4AlaQpYGUuaVVzOv/iDHNJ6ohhLklTwDCXpCnQa888yUHg48BV4KvAJ4DbgWeB\nbTTX/6yqiwu2fxi4AjxWVSdaGbkkLcOqns6f5B5GF3H+/qr6v0meBT4G/ABwsqqeSrIfOAAcSLIT\neBDYCWwGnk+yo6qutvYvkKQVos82y58Cl4D3JrkCvBf4Q+AgcH+zzWHgBUaBvhc4WlWXgLNJzgC7\ngBcnP3RJWll6C/Oq+naSXwa+AXwX+HxVnUwyW1XzzWbzwGyzfDfvDu5zjCp0SVr1Zmba2025WJvl\n+4CfB+4BvgP8hyQfX7hNVVWSmzWT+j+4U9KKNi3T+Zc7A/TUqVOcOnVqSdsu9mfibwP/o6r+GCDJ\nfwL+DjCXZGNVzSXZBFxotj8PbF3w+C3NOkla9ZbbZrnvvvu477773r79zDPP3HDbxcL8FeAXknwP\n8BfAbuAl4M+AfcCTzffnmu2PA0eSPM2ovbK92V6SVr0+e+ZfSfKbwMuMDk38XeDXgPcBx5I8QnNo\nYrP96STHgNPAZeDRGsIcWkkagDZ75ukjaxfpsUvSxA2hrnz11VfHevwHPvABquq6zX9ngEpSRzwF\nriRNAc/NIkljGsJ0/t6OM5ckTY6VuSRNAcNckqaAYS5JPRvC6QBuppfjzCVJk7Wm7wFIksZnmEvS\nFBhcmCd5IMkrSb7eXMVIY0hyNsnvJzmVxJOeLUOS30gyn+SrC9bdleRkkleTnEiyoc8xriQ3eD0f\nT3KueX+eSvJAn2NcyQYV5knWAv8aeIDRpeceSvL9/Y5qxSvgI1X1oara1fdgVphnGL0XFzrA6JKJ\nO4AvNLe1NNd7PQt4unl/fqiq/lsP45oKgwpzRpeYO1NVZ5tLz32W0aXoNJ5h74YfqKr6beBPrlm9\nh9GlEmm+f7TTQa1gN3g9wffnRAwtzDcDry247WXnxleMLqz9cpKf6XswU+BGl0zUrfvZJF9J8inb\nVrduaGHucZKT9+Gq+hDwo8A/TvJDfQ9oWjTn6vc9O55fBe4FPgh8E/jlfoezcg0tzK+97NxWRtW5\nblFVfbP5/i3gc4xaWbp180k2AlxzyUTdgqq6UA3g1/H9ecuGFuYvA9uT3JNkHfAgo0vR6RYkeW+S\n9zXLtwM/Anz15o/SIo4zulQivPuSiboFzR/Et/wkvj9v2aCm81fV5SSfBD4PrAU+VVVf63lYK9ks\n8LlmGvIM8JmqOtHvkFaOJEeB+4H3J3kN+EXgCa5zyUQt7jqv5z8HPpLkg4zaVf8H+Ec9DnFFczq/\nJE2BobVZJEm3wDCXpClgmEvSFDDMJWkKGOaSNAUMc0maAoa5JE0Bw1ySpsD/A/jvvDd2PkdyAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efeb3741090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# - \n",
    "N = 5\n",
    "XN = np.kron(np.eye(N), X)\n",
    "plt.imshow(XN, aspect='auto',interpolation='nearest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# - create N subject data, add noise to the previous data\n",
    "YN = np.kron(np.ones(N), Y)[:,np.newaxis]\n",
    "YN += np.random.normal(0, 0.7, size=YN.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 1)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# - create contrast\n",
    "c = np.asarray([1,0,0,0]).reshape(4,1)\n",
    "cN = np.kron(np.ones(shape=(1,N)), c.T).T\n",
    "cN.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.20941620009 820 1.27693693883e-12\n"
     ]
    }
   ],
   "source": [
    "# - compute t stat for fixed effect analysis\n",
    "t_Nfix, df = t_test(XN, cN, YN)\n",
    "print(t_Nfix, df, sst.t.sf(t_Nfix, df)*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a random effect model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.424  0.739  0.771  0.691  0.631]\n"
     ]
    }
   ],
   "source": [
    "# compute signal (c.T.dot(beta)) for each subject\n",
    "signals = np.zeros((N,))\n",
    "nTR = Y.shape[0]\n",
    "pinvX = npl.pinv(X)\n",
    "\n",
    "for sub in range(N):\n",
    "    Y_sub = YN[sub*nTR:(sub+1)*nTR]\n",
    "    signals[sub] = c.T.dot(pinvX.dot(Y_sub))\n",
    "\n",
    "# print the signals\n",
    "print(signals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.5857974999 4 0.000450662649435\n"
     ]
    }
   ],
   "source": [
    "# What is the random effect t-test ? Create a design matrix and a contrast, use t_test()\n",
    "X_random = np.ones((N,1))\n",
    "c_random = np.ones((1,1))\n",
    "t_random, df_random = t_test(X_random, c_random, signals.reshape(N,1))\n",
    "print(t_random, df_random, sst.t.sf(t_random, df_random)*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10.585797499881789, 0.00045066264943452163)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# - check that we get same between t_test and scipy stat ttest_1samp\n",
    "sst.ttest_1samp(signals, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compare the two statistics: which one is \"better\"  ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}