Last active
August 16, 2024 13:45
-
-
Save fladd/9d3707f7a457cf8272042e742952c950 to your computer and use it in GitHub Desktop.
Comparing GLM modelling approaches of SPM and FSL
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Comparing GLM modelling approaches of SPM and FSL" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 112, | |
"metadata": { | |
"jupyter": { | |
"source_hidden": true | |
} | |
}, | |
"outputs": [], | |
"source": [ | |
"import numpy as np\n", | |
"from nipy.modalities.fmri.glm import GeneralLinearModel\n", | |
"from nipy.modalities.fmri.hrf import spm_hrf_compat as hrf\n", | |
"from matplotlib import pyplot as plt" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Let's simulate some data\n", | |
"\n", | |
"Imagine an experimental design with 20 volumes of a task in condition A, and 20 volumes of a task in condition B, each surrounded by 20 volumes of rest. \n", | |
"\n", | |
"Here is a simulated timecourse of a single voxel that gets differently positively activated by the two task conditions (note how it is delayed by 6s):\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 113, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"[<matplotlib.lines.Line2D at 0x1f5246c9438>]" | |
] | |
}, | |
"execution_count": 113, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABSl0lEQVR4nO29eZgcV3n/+zm9b9Ozb5qRPNps2bJlG8s2u8EQjLGBhAABB3AIiUOAB0Juwg8nv9/lXgK54SaXLYsDBgIEggkEgsP6M4uxA3iRkTd50S6NRtLsW+/buX9Unerq7uqZnrVnus/nefRourqmu2qq+633fN/veY+QUqLRaDSa5sBV7wPQaDQazfqhg75Go9E0ETroazQaTROhg75Go9E0ETroazQaTRPhqfcBLERXV5ccGhqq92FoNBrNpuKRRx6ZkFJ2Oz1Xc9AXQriBA8CIlPJmIYQAPgK8AcgDd0gpP21u/xTwKiAB/J6U8tfma9wK/E/zJT8ipfzSQu85NDTEgQMHaj1EjUaj0QBCiFPVnltKpv8+4Gkgaj7+PWArsEdKWRBC9JjbbwR2m/+uBe4ArhVCdAAfAvYDEnhECHG3lHJ6Cceg0Wg0mhVQk6YvhBgEbgI+Z9v8x8CHpZQFACnlmLn9tcCXpcEDQJsQoh+4AbhHSjllBvp7gFeu0nloNBqNpgZqLeR+EvgAULBt2wn8jhDigBDiB0KI3eb2AWDYtt8Zc1u17SUIIW4zX/PA+Ph4jYen0Wg0mlpYNOgLIW4GxqSUj5Q95QdSUsr9wJ3AF1bjgKSUn5VS7pdS7u/udqxDaDQajWaZ1JLpvwB4jRDiJHAXcL0Q4isYmfq3zH2+Dewzfx7B0PoVg+a2ats1Go1Gs04sGvSllLdLKQellEPAm4CfSinfAvwn8FJzt+uAw+bPdwNvEwbPBWallOeAHwGvEEK0CyHagVeY2zQajUazTqzEp/83wFeFEO8HYsAfmNu/j2HXPIph2Xw7gJRySgjxV8DD5n4fllJOreD9NRqNRrNExEZurbx//36pffoajWY1+MXRCfpbA+zojtT7UNYcIcQjZr21At2GQaPRNAV/8vVH+czPj9f7MOqODvoajabhyRckk7E06Vy+3odSd3TQ12g0Dc9MIkNBQrawceXs9UIHfY1G0/BMxjMA5PKFRfZsfHTQ12g0Dc9kTAV9nenroK/RaBqeyXga0PIO6KCv0WiagCkt71jooK/RaBqeCS3vWOigr9FoGp7JmJJ3dKavg75Go2l4ivKOzvR10NdoNA2Pcu9ktaavg75Go2l8Jkz3Tk67d3TQ12g0jY927xTRQV+j0TQ02XyBmUTW/Fln+jroazSahmbazPIBctq9o4O+RqNpbFTfnY6wT7t30EFfo9E0OMq50xsNaPcOOuhrNJoGR/Xd6Y36tXsHHfQ1Gk2DozL9vmhAyzvooK/RaBqcyXgat0vQGfHpNgzooK/RaDYAT47M8v98/2mkXP1MfCqeoT3kw+d2I6WxdGIzo4O+RqOpO3c9fJrP3Hecs7OpVX/tiViGrogPj1sAuhWDDvoajWbJ5AuSbx88w42fup9/f3h4xa936OwcAM+cm1vxa5UzGUvTEfbhNYN+sxdzddDXaDRL4v4j47zqU/fz/q8/xtPn5njgxOSKXi9fkDxzbh6AZ87Plzx3dCzGLXc+wKw5o3Y5TMUzdEb8eFxGuGv2Vgw66Gs0mprJ5Qv84ZcPkMzm+cdbnsPeLdGSGa/L4cREnGQ2D8BTZZn+D588xy+PTXLv4bFlv/5kLEOnLdNv9lYMOuhrNE3MyEySN3+29kx6NpkllS3w+y8Y4qZ9/XRG/EytIAsHOHR2FoBtHaEKeefg6RkAHji+vNFEOpdnPp2jM+zD4zYz/SZ38Oigr9E0MU+cmeVXxyc5PDa/+M4YQR+gNeQFoCPkXXGm/9TZOXxuFzft6+fERJyUmfVLKTk4PAPAL48tL+ir7pqGvGNq+jrT12g0zYpyssTSuZLt86ksr/rU/Tx1tjTztoJ+0Aj67WHfioP+obNzXNgXYd9AKwUJh0eNG9DwVJKpeIbdPRFOTSYYmUku+bXVxCyjkGuEO+3e0Wg0TYsKgIl0vmT7qckET52b4+DwdMn28qDfEfIxn86RyS0vkEopOXR2lr39rezpjwJYRV313n/8kp0A/GoZ2b5qtma3bGr3jkajaVpUsI6XZfrq8UyZXl8M+j7AyPSN/ZaX7Z+fSzGdyLJ3IMq2jhBBr5unzxuji4OnZwj53Lz68i20h7zLC/rmguh2947O9DUaTdOiAmA8Uxb0MyrolwbzufJM3wz6k8uUeA6NGAF+75Yobpfgor4WnjaLuQeHZ9g32IrX7eK5Ozp54Pjkkmfslso7WtOHJQR9IYRbCHFQCPFd8/EXhRAnhBCPmv+uMLcLIcSnhRBHhRCPCyGeY3uNW4UQR8x/t6762Wg0miWRrprpG3LPdFmmrzJ/S9MPGUF/ubr+obNzCAF7+gxp5+L+Fp45P08qm+eps7NcsbUdgOfv7GRkJsnpqcSSXn8ynsHrFkQDHu3eMVlKpv8+4OmybX8upbzC/Peoue1GYLf57zbgDgAhRAfwIeBa4BrgQ0KI9hUcu0ajWSHKsx4r0/QXkneCXjc+jxE6VKY/tUx559DZWbZ3hgn7PYAR/GcSWX76zBjZvOTKbW0APG9nJ7B0XX8ylqYz7EcIgdelffpQY9AXQgwCNwGfq2H31wJflgYPAG1CiH7gBuAeKeWUlHIauAd45TKPW6PRrAJWIbdM3omlneWd2WTWyvKhGPRrzfR/fnicT/34iNX07NDZOS7ZErWe39PXAsDXHjoNwJVb2wDY2R2hu8W/ZOvmZDxjHaOV6eugXxOfBD4AlI+LPmpKOJ8QQvjNbQOAvRnHGXNbte0lCCFuE0IcEEIcGB8fr/HwNBrNclCF3HLLppJ3ZpKVmX5bqBj01c9T8domaH38fz/LJ358mA9883Gm4hlGZpLs3dJqPa8cPPcfmWCgLUhPNACAEILn7ejkV0vU9SfjGTojKuibmb6WdxZGCHEzMCalfKTsqduBPcDVQAfwP1bjgKSUn5VS7pdS7u/u7l6Nl9RoNFWoZtmsVsidTWaJ2jJ9r9tFNOBhugZ5Zyqe4fGRWXb3RPiPX5/h7V98GDCKuIrWoJeBtiAAV5jSjuL5OzsZn09X9Ocp50u/PMkd9x7jh0+e4/xskq6IkY96XTrTh9oy/RcArxFCnATuAq4XQnxFSnnOlHDSwL9g6PQAI8BW2+8PmtuqbddoNHUiU8W9E7Np+vbMulzeAUPimapB3rn/yDhSwt++4XLe+7LdPGbOtrXLO1CUeJS0o3j5Jb2EfG4+9eMjVd9jIpbmQ3cf4mM/fIZ3fuXXjM6l6YkaQd/y6Te5ZdOz2A5SytsxsnqEEC8B/kxK+RYhRL+U8pwQQgC/CTxp/srdwHuEEHdhFG1nzf1+BPy1rXj7CvW6Go2mPlSXd4zHuYIkls7REjAC/Wwyy6VlQb897KvI9H/6zCguIXjJRT3WtvsOT9AW8nLZQCuXD7bicwsOnp6xMnHFxf1RfvLMmFXEVXRF/Lzzup18/J7DPHxyiquHOirO58REHIB/uOVKLugIc2Y6wbU7jCKw1XCtySdnLRr0F+CrQohuQACPAu80t38feBVwFEgAbweQUk4JIf4KeNjc78NSyqkVvL9Go1khVeUd201gJpEtCfoVmX7Ix/m50sVP/uYHzzARy/DLD15PwOtGSsn9R8Z54a4u3KaL5j3X73Y8ppsv72dkJsllA20Vz/3hi3bwbw+e5iPfe5pv//HzcZmvpTgxbgT9fQNtbOsMcdlgsV6gWysbLGlylpTyXinlzebP10spL5NSXiqlfIuUMmZul1LKd0spd5rPH7D9/heklLvMf/+yuqei0WiWSrVMP1YW9MG4QSQyedqcMn2bvCOltPrmfPugoeA+c36esfk01124eJ1uT1+UT/zOFZYt1E7Q5+bPbriIx4Zn+O4T5yqePz4Rx+sWDLQHK57z6MlZgJ6Rq9E0NcqzXm7ZjKfztAQMIUBJN+UdNhUdYV+JT38ynrH643/+v08gpeTnhw0n3otrCPqL8borB7ikP8rHfvCM1ZFTcWIixraOkDWasGM1XNPuHY1G06xYhVwHeWewPQQ4BP3yTD/kI5UtkMwYrzFszpq9eV8/R8di/PzwOPcdHmdPXwu9pgVzJbhcgg/euIeRmST3PDVa8tzJiQTbuyKOv6dbKxvooK/RNDFK3snkCyWdMmPpnGWdVMFeyTzRCveO8XgybjQ3G542WiC/87qd9LT4+cefHeXAyemapJ1aed7OTnxuF0+OzFrbCgXJick4O7rDjr/j0a2VAR30NZqmxh4A7RKPkekbQX/anHilmq1VaPpW/x3jeZXp7+gOc+vzh3j45DSZfGFVpB2F1+1iT38LT9iC/tnZJJlcgaFO56CvF0Y30EFfsyQePjm17N7pmo2HPeir4m2hIIln8kSDXlr8HmaSC8s75f13zkwn6Ir4CPk8/O612wh63QS9bvYPrW6rrb1bWnlyZNaaR6Dsmtu7qmT62r0D6KCvWQLnZ1O84Z9/xQ+erHRNaDYn9ht4wtTkE2ZxNOJ30xb2WrJOVU2/rP/O8FTSqge0hXz82Q0X8UfX7cDvca/qsV820MpcKscZU046aQb9avKOXhjdYCU+fU2ToTK+8s6Lms1LJi/xugXZvLQyfeXRD/s9tAV9FYXcck2/U2X6KuhPJ9g32GY9/44Xbl+TY790wJjJ+8TILFs7QhyfiBPyuelp8TvuL4TA7RK6tXK9D0CzeVDujGSZTU6zecnk8rSZmryaoKWCf8TvoS1UmumHfW7L+qiIBry4hOHyyRckZ2eSbHXwya82F/a24HEJq5h7YiLOUGcYo0mAMx6X0O6deh+AZvNgBf2MDvqNQjYvaTd99xWZvs9DW8hnNV2bSWStG4Qdl0vQHjL675yfS5HNS7Z2hNb82ANeNxf2tvCkuXj7iYk426tIOwqv29X08o4O+pqaURm+zvQbh2y+YAVyFexjNnmnPeS1Vs8q77BpR/XfUc6dre1rH/TBkHieHJklkytwZjrJjipFXIXHreUdHfQ1NWMFfZ3pNwyZXMHK9JVlU03Uivg9tAW9zKWy5AuSuWSW1qBzGbDDzPStoN+x9vIOGMXcqXiGh05MkS/Iqs4dhcelM30d9DU1k9CafsORyRcsn71aMrFYyHXTFvIhpeHRd2q2pmgPe5mOZxmeTuISsKVtfYL+3gGjodp/PXYWgKFFgr7XLbRls94HoNk8pHSm33Bk8wWiQaMQqzL98kIuGCtoLRT0Vf+dM1MJ+luDFcXeteLiviguAT88dB6gRnlHZ/oaTU3oTL/xyOQK+Nwuwj6Po2XTmm2byDCTzDgWcsGYlTsdz3B6KmHN5F0Pgj43u3tamE1maQ95qx6fwuty6TYM9T4AzeZBu3cai3xBUpDg87gI+d1WsI+ncwgBIZ/byvTH5lKksoUFM/1cQfLM+fl1ce7Y2Wv69RfT88HM9LWmr9HUhpJ3EjrTbwjUbFyv20XY7yGeUT79PGGfByGElTmfmjQKtFXdO1ZdILeumT7ApebC6ovp+WAUcrV7R6OpESXvpHSm3xCotspetyDi95Rk+mG/0TJBOXtOma6chTJ9xXrZNRVqdazF9HzAmn3czOigr6kZ7dNvLJS27fe4CPncxRm5mRxhv2HNbAl4EQJOTRp9baq7d2xBf53lncsGWnnVZX38xiV9i+7rcetMX/fe0dSM0vITOtNvCOzyTsTv4eyMsc5tPJ0jYgZ9t0vQGvRycsLI9MvbKis6Qvagv77yTsDr5p9+96qa9vW4dKavM31NzagMv3yJOs3mJJsvBv2Qz2ObnJUj7Cvmg21BL2dnjU6WVeWdiBH0fW4XvS0rXx1rrfC6XdqnX+8D0Gwe7A3XVA9zzeZFBX2fxyjkxtK2Qq6/2AZZTdCC6kE/7HPjc7sYaA/iclifdqOgffo66GuWgHLt5AvSKgJqNi/pEnnHXVbItWX6toXQq7l3hBC0h73r7txZKroNg9b0NUvA7tpJZQqrviiGZn1Rwc8o5HpIZvPkC7Ii6Cs7Zovfg3uBLP7W5w9xQcfiDpp6otsw6KCvWQKJbHEN1WQ2TyvOWZ9mc1BeyAWjFUPMVsiFYqbfGlr4er/rJbvW6EhXD8O909yZvpZ3NDWTzBQI+YzsXts2Nz9Zm08/ZGr4c6kc6VyhrJBrZPrV9PzNhNcldBuGeh+AZvOQzOSsSTjK6aHZvGRshVyV2Y/NGbZNeyG3PWxm+g0Q9HUbBh30NTUipSSZzVtBX9s2Nz92eSdkZvajc2mAEnlHBfvGCPp6cpYO+pqayOQLFGRxun0y09xfnEag1LJpZPZj8yrTryzkNkLQ9+rJWTroa2pDefS1vNM4WEHfbZd3KjP9Wgu5mwGP20VeF3I1msVRhVs13V4Xcjc/lrzjscs7jZ3pe9y6kFtz0BdCuIUQB4UQ3y3b/mkhRMz22C+E+LoQ4qgQ4kEhxJDtudvN7c8KIW5YlTPQrAuq346abq81/c1PxpQ5SjL9eSPTtxdye6MBnrejk2uGOtb/IFcZr0tbNpfi038f8DQQVRuEEPuB9rL93gFMSyl3CSHeBHwM+B0hxCXAm4C9wBbgx0KIC6WUOnpsApS802nJO/qybXZUpu9zuxBmEq8yfbu84/O4+Nptz13341sLPG5BviCRUiLExm0XsZbUlOkLIQaBm4DP2ba5gb8FPlC2+2uBL5k/fxN4mTD+uq8F7pJSpqWUJ4CjwDUrO3zNeqEy+46wH9DyTiNg+fQ9wvLlFzP9xpy3qdbubeZibq3yzicxgrtdDHsPcLeU8lzZvgPAMICUMgfMAp327SZnzG2aTYDK7NtCRn91vZDK5idry/TdLkHQ62YqngFKM/1GwmO2kWhm2+aiQV8IcTMwJqV8xLZtC/AG4O9X+4CEELcJIQ4IIQ6Mj4+v9strlonK7INeN0GvW8s7DUAmX0AIrH46Ssd3uwR+T2N6PDw6068p038B8BohxEngLuB64BCwCzhqbg8JIY6a+48AWwGEEB6gFZi0bzcZNLeVIKX8rJRyv5Ryf3d393LOSbMGKE0/5DOCvpZ3Nj+ZfAGv22Vp20rSCfvcDat3e91mpt/EDp5Fg76U8nYp5aCUcgijEPtTKWW7lLJPSjlkbk9IKVW3pbuBW82fX2/uL83tbzLdPduB3cBDq3w+mjXCyvR9bgI66DcE2ZzE7y6GAGXbbFRpB4zWykBTO3jW4up+HvhXM/OfwrhRIKU8JIT4d+ApIAe8Wzt3Ng9Wpu/1EPK5rceazUsmn8drk3EiprzTqEVcMNw7QFN79Zd0daWU9wL3OmyP2H5OYej9Tr//UeCjSzpCzYZAZfYBn4ugT2f6jUA2J/HZMn1L3mngoF+Ud5o302/Mao1m1Ulm8rhdAp/bZcg7ZZn+5+4/zqGzs3U6Os1yyOQLeD1F7T7cVPJO82b6OuhraiKRyRP0GgW+UFmmn8sX+Mj3nubrDw8v8AqajYYq5CrClrzTuCuieS15R2f6Gs2CJLN5Al4jGATLMv2ZZBaAszPJuhybZnlkc4USeUcVchtZ3rEyfR30NZqFSWXz1qpZ5ZbNmYQxoWdkJlWXY9Msj0y+gK+kkNsE8o7K9LW8o9EsTCKTI6gy/TL3znRCZ/qbkWyZvBNqAveOOl+d6Ws0i5DMFghWyfSnzan7s8kssbTus79ZKHfvNEWm79KTs3TQ19REsjzTz+Yx5tzBjJnpA5zT2f6mIZ0vlPj0lXsn7GvcQq7VhqGJJ2fpoK+piWQ2X8z0fW6khLTZsGva1PQBzuigvyxOTcZJ59Z37kN5ITfcFPKOzvR10NfURCKTL5F3oDhLV7l3YHPo+p+7/ziPnJqu92FYpLJ5bvjkfXzjwJl1fV+jkGvz6TeFvKMbrumgr6mJlOnTB1vQN3X9mUSGzrAPt0tsiqD/d//7Wb7+8Ol6H4ZFLJ0jlS0wbvayXy/KC7nbOkJ43YKhrvC6Hsd6YmX6TezeadxbumZVSdotm+b/qr3ydDxLR9hHwOvm7Aa3bebyBVLZAufn1jfALoQaMa13a4tyeeeCzjBPffiVJTeCRsOj3Ts66GtqI+GQ6avVtKYTGdpDPtpDMLLBM/24GWDH5jbOzUkF+/VuYpcpK+QCDR3woejeaeaGa419hTWrQqEgSecKxRm5vnJ5J0tbyMtAe3DDyztx01I6upGCvhns13thmkxZpt8MWD597d7RaKqjgnuorJBryTtmpr+lLcD52RT5Gr5Q3zgwzKd/cmSNjrg6iYwR9KcTWWukUm8SlryzvnMcsnlZMiO3GfBo944O+prFsS+gYv8/mTG8+jOJLG1hL1vaguQKkrH5xbPobxw4wz/de5RMbn2/fLF0MdCvd+G0Gql6yjvuxlwhqxpe7d7RQV+zOCoYOWn6iUyeTL5gZvpBoDbb5umpBKlsgSdG1rcdc9w2Y3ijSDyJOsg7+YIkX5D43I07EcsJj3bv6KCvWZxqmX4ik7cmZrWHvAyYQX+xxmvpXJ5RczTw4InJNTnmatjbRJzfIEHfKuSuo9ykCpn2fvrNgEe3VtbuHc3ilGf6Ia/xsUlm81YLhraQj/7WAFCa6c8msng9wmrbCzAyncTs4MBDJ6Z410uc33dsPsW3fz1CJOChM+xjqCvMnr7ois6lNNPfGPJO0qwzrKe8kzGDfrMVcnVrZR30NTWgZAeV4Qd8xhcnlbVn+j5aAl6iAY8V9KWU/PY//5J9A618/HeusF5veNp4/sLeCAdOTpPLFyz/tJ27Hhrm4/ccth4LAf/9P663RhTLwR70N4ptsx7yTtaspTRbIdftEgih5R2NZkFUoVFl+j63C5cwnDCqrXJ7yAvAQHvICvoHh2c4OhbjsTMzJa93eioBwOuvGiSWzvH0uXnH9z01maAvGuBXt1/PP9xyJVLCE2WvtVSUT78r4ttw8s56uolUpt/ovnwnvC5XU8s7zXfFNUumaNk0BobGkokekpmCtYBKW8gHwEBbwNL0v3NwBDCCt30yzJmpBD6Pi5v3bQGq6/rD0wm2dYTobw3y8ot7cbsEh87Orehc4ukcQhizTzdKIbcePv1szgh6zRj0PW6hLZsazUIkyjR9wFgcPZtnOq40fSPT39IWZGQ6QS5f4LuPnyPsc5MrSE5Nxq3fHZ5OMNgWZEtbkAs6Qzx4YsrxfYenEgx2BK3329kdXnHQj6VzhH0e+loDjG0UTd9WyC2s06QhS9NvMnkHjFm5enKWRrMAKigpLR8g6HORzOSYTmRo8XusjHFLW5C5VI4fHRplMp7h91+4HYCjYzHrd4enkmztCAFw7fYOHj45VRHs0rk85+dSbG0PWdv2bmnl0NmVWTzj6Rxhv5velgDn51LWmgD1xJ7hp9dp3oKaH+FrMp8+GKMb3YZBo1kA5S6xO3BCXo/p3snQFvZa25VX/46fHyUa8PAOM+gfGS0G/dNTCbaaGfw12zuZSWQ5PFaq65+dSSEl1s0BYO+WKKNzaSZiy8/Q4+k8Yb+H3qifRCa/IVb6sls11YzhtSbbxJq+Ie/U/2ZfL5rvimuWTDJjBIgSecfnJpktMJ3I0m7q+WBo+gBPjsxx075+2kI+BtqCHB03gv5cKstsMmtl8Ndu7wDgweOlEs+wWezd2l506lyyxbBrrkTiiaVzRPyGvAMbw7Zpt2qul67f3PKOq24Lo2dyBX77jl/yq2PrOz/FTvNdcc2SSWbz+Dwu3K6iFBD0GvLOTCJjFXGhmOkDvObyAQB29kQseccK5mYGP9geZEtrgIfKdP3h6dL9APb2twKsSOJJZAxNv6dFBf36F3PtQX+9HDzKstmMmb63jpn+dCLDI6em+fXp+i3i03xXXLMoU/EM77vrILOmHdO+Pq4i5DPkHSPTL8o7PS0BPC5Bf2vAyuJ3dUc4Nh6jUJAMTxl2TpXpCyG4ZnsHD52cKtHXh6eSeN2C3mjA2tYa8jLYHlxhpl+Ud2BjBP1EVmf664nH7aqbT1/NE4nXUVZsviuuWZRHTk3znUfPcu/hMcBcH7cs6Ae9bpJmGwa7vON2CV52cQ/veOF2XObIYHdvhFS2wMhM0sr0t9ky+Mu3tjE+n2bM1gBteDrBYHuoZHQBhq7/1AqCvlXIjW4ceSeVydMaNG6c6xb0c805IxcM9069fPrq+uqgr9lQxNJGhv/EGUNGSWSKq2YpAl43sXSO+VTOsmsqPvPW/fzBi3ZYj3f1RADDwTM8naAl4KHV9juXDbSWvB8YXv7B9sqZt3u3tHJiIl5TATaWzlW4c4yg7yHs99Di92yQTD9HZ9i4ca6bvJNvXp++1+2qm09fBft5HfQ1G4lYyvhAPm4G4VQ2by2gogj53FZr4rZgadAvZ1e3LehPJUpsmGAUaF2Cko6bhsOndD+ASweMYu7T56pn+/F0jg//11Nc9n/9iB88eb703MxCLkBP1L8hgn4yU6DDDPrrlelnm1reqZ9PX0l59cz0de8dTQWq5/yTZ2fJF6Rjph/0uVHfm/awr/wlSmgP++gM+8xMP8nO7tKFt0M+Dzu7IzxpBv1Y2mjvUH5zACPTBzg0MsvVQx0Vz993eJzbv/WEtWyjfX5ALl8gnSsQNq2nvdHABgn6OSvor1enzYxVyG1Cn76rfj79RFoF/fot4FPzbV4I4RZCHBRCfNd8/HkhxGNCiMeFEN8UQkTM7X4hxNeFEEeFEA8KIYZsr3G7uf1ZIcQNq342mlVByTuJTJ5j4zFD03eQdxR29041dvVEODI2z/BUokTPV1w20Gpl+kWHT6W809PipyvicyzmnplO8PYvPozf6+Ib73weLX6P1RAOin13wn7j2Puigbpr+lJKktk8nREz6K+TT79Zu2xCfX368czmknfeBzxte/x+KeXlUsp9wGngPeb2dwDTUspdwCeAjwEIIS4B3gTsBV4J/JMQorlWcNgkKHkHDIknmXGWdxTtoYXlHTCC/hMjs6RzBUfZZu9AK2PzacbmUjaPfuV+Qggu2dLqGPSfPT9PviD529dfztVDHbSFvVbrZygOqcOWvBNgbD61bq0PnEjnChQk6y7vZJq0yyYY7p1sveSdzeLeEUIMAjcBn1PbpJRz5nMCCALqr/ha4Evmz98EXmbu81rgLillWkp5AjgKXLMaJ6FZXWLpPP2tAcI+N4+fmSGZdZB3vPagX1umr4qHTsHcKuaOzFqtl51uDmA4eI6MzVcstXhq0rhZDHWGrOMqyfTLgn5v1E82L0v2WW9U4Vb9DddL3mnmGbleV/0arsU3kXvnk8AHgJK/lBDiX4DzwB7g783NA8AwgJQyB8wCnfbtJmfMbSUIIW4TQhwQQhwYHx+v+UQ0q0csnSUa8HLpQKuV6TtZNhXl7h0ndve0WD87yTZ7t0QRZjF3eCpB2OeuOoLYuyVKNi85Uta64dRknBa/x8qa20I+puPFgK4cPxGbvAP1XUFLZfYRv4eA17VuC6k0fSG3TvKOur720fR6s+gVF0LcDIxJKR8pf05K+XZgC4bs8zurcUBSys9KKfdLKfd3d3evxktqlkgsnSMS8HD51jaeOjdHLJ2r0PTVY49LWG6YhVC2TYBBh0w/7PewoyvMkyOznJk2nDvGALGSC3uNG4i9SAtwairBts7i73WEvFa/fygWz1Qht8cM+vXstmlfijLoda+7vONxNV8h15B36pXp56z/69Xsr5bb/AuA1wghTgJ3AdcLIb6inpRS5s3tv21uGgG2AgghPEArMGnfbjJobtNsMGIpw8t+2UArmVyBxAKZflvIVzU42+mN+on4PfS0+CvqAwpVzLV34XRiqDOM2yUqg/5kgqHOojOorUzeiTnIO1DfWbn2pSjVLOf1IJOX+Nyumq5do2HIO/XS9I3rW5DruyaynUWDvpTydinloJRyCKMQ+1PgrUKIXWBp+q8BnjF/5W7gVvPn1wM/lcYt7W7gTaa7ZzuwG3hoNU9GszrE0jla/B4uH2yztjlZNqG2Ii4YBdiL+1vYUWbXtHPpQCujc2mOjcccdX+Fz+Pigo5QhR3zzHSCCzqLv9ce8jGfyllShupgafn0W+ov75Rk+j73uso7zSjtgNmGoW6aflHWqZfEs1yfvgC+JISImj8/Bvyx+dzngX8VQhwFpjBuFEgpDwkh/h14CsgB7zZHCZoNhprAtLUjSFvIcMCUZ+fFoL94EVfx8TdeseDzqpibK0hH3d/Ozp4IR2xB/9xsimxelgZ9s+XzTCJLd4vfKp6FTE3f53HRGfbVNdNXck7IknfWybKZKzSlRx+MuQl1c+/YbuqxdI6eOhzDkoK+lPJe4F7z4Quq7JMC3lDluY8CH13Ke2rWn1jK0PSFEFw20Mr9RyYqNX1L3qkt04fqbhzFXjPog7PDx86ungg/e2aMbL6A1+3ipLky1wU2eUfdkGYSGbpb/NakM3sNojXkZa6ORTWV2Qe8Zqa/ju6dZnTugNFauV6Zvv2mXq8JWs151TVVKRQk8Uze0r33DRqBuJplcymZ/mJEzGIuLH6D2NUdMZdhNGya6v9yeQewirnxdA6XKHUeBb1uUuu4Nm05yWxxgRrVxG49yOSaWd6po6afyRPwGn/3+XR2kb3Xhua86mtILl9gbL7+U/uXi9IcW6yg3wbg0FrZzPTDtWf6tXCpme07NVuzY2/iBoZd0+9x0dtSbMWsRiFTpm1TrY9rL14GveuXXTthX6AmtI6ZfiZfaMrZuGAul1jH1sqqlqQz/Qbhaw8Pc/3f/XzduiWuNpaXPWAE/Rfu6uKWa7dxzfbOkv3Cfg8el7C87qvF7167jXdet9MaaVRjpxn0j42roG+0d3DZLIjKrz9jOnhUh0076ympOKGG+6qQu54N15o2018j986hs7Ps+V8/sGaUO5HI5OlpMVxj9ZqgpRuurTJHR+eJpXPMJLL0tW6+LhPKUaCCY9jv4a9/67KK/cJ+D9961/NLJl2tBtfu6OTaHZ2L7hfxe+hvDdgy/USJng+V8k4ik7eKuIqAt9gttB6o5CDoda+7vNO0mr7bRa4gkVKuqmX1yGiMVLbAM+fnq8qT8XTOWsuhXv13mvOqOzAVz6xKdq4aeM0m11+v+8+DI/zbg6dX9Boq02+pYcLVvsG2igLverLLXIZRSsmpqXiJng9G9uz3uKxM395WWRHwuus6Kktk8nhcAp/Hta7yTjYvm9e9Y44GV7u9svqcLWQBTmbzdNc509dBH2Ooe+On7uMT9xxe8WupCz5Th34u//bgaT7/38dX9Brl8s5GRgX983MpUtmC1XPHTnvIZ2n6cVPTtxP0ukhl66PvQumqZEFzclYtMzX/9OuP8v6vP7rs9800s7xjjnBWW+KZMRO90VnnoJ/JFcjmJZ1hHy6h5Z268qtjk4zOpXl2dH7xnRdhTAX9OmT6M8nMilsFK3mnltYK9WZXT4RkNs8DxycB2NZZOfGrPeyz5J1YOlfRAqL+hdxi2+qg142UkMoWFh1BPXRyakVBO5Mr0LIJbuxrgRrhZAsFgqzeSFV1dD1XJeir+k3Y7yHs8zBfJ6twc97qy/jBk+cAGDG7O1bjm4+c4TM/P1b1+UJBWuu81kPemU1mzSUMl//exaZkGz8gqBW5fvy0sZavc6bvLRZyMzmr2ZoisMaSylQ8w3u/drDq58G+VoFyRC12PJlcgbMzyaoZZS1km9i9o/oNLTfTl1I6tuNW17jaZD/7eg6RgEfLO/Uily/wo0OjAIzMJBccWn/jwDBfPzBc9fmJeNrSCWcTdcj0zfc8v4JgsKmCvungue/ZcTwuwUBbpc3T3l45kc5Xune8bjK5Avk1mqH50IlJ7n7sLL8+Pe34vL2vkQr+i83KPTOdoCCNILLcG3xz+/SVvLM8We/7T5znqo/cU1ELWkzTV730Qz5jjeb4Os2+Lqc5r7qNB45PMRXPcM32DhKZfMmiG+WMz6ethklO2Ls1ziTXV9NPZfOkzc6JK+klU+7e2ch0Rvy0h7zMp3MMtAetL7OdNlunzZiTZdMMuGtVzFXOoPEqslsqWyrvAIs6eNRENFh+s7hmnpFblHeWd6M/dHaW6US2ojvrYpp+wpbph/1a3qkb33viHCGfm1uu2QZgra3qxOhcasEhmT3DXm95x/5+1TTFWoilc/g9rk2TBapsv9yuqegI+5hJZMjkStfHVQRrlFSWiwr61Sbs2dcfrlXeOWW2nAA4P7u8Gk42LzfNNV5tPK6VZfpKwp2Ml/7t1eh+Pp2zRsx2VGYf9Hpo8Wt5py4Y0s55rt/Tw05THz5TRdePpXPEM/kF+2CPml/ssM+94IhhLbAH/ZXKO5tB2lFYQb+KL7ot5KMgi3+TcLmm76ktu14u4zEz068yF8C+QI36f7EJWidtmf5yR3Xppvbpm5n+MjV9dS2n4qWj+Zlk1rI6O30HlUpgZPpuPSO3Hjx0wpB2brqsny1txoSJs1UyfTWMXqgP9uhsCiGM2aLrnenbbzIrknfMBVQ2C7vMyWHlHn2Fav18ZtoIlBU+fTO7TufWOtOvEvSzeYLm6KPWUcepybjVo2gl8o6vWX36StNfZiuGYqZfDPqFgmQmkWFPv/F5dLouKtNXmr7TaGA9aOqg/70nzhH0unnJRT10hH0EvK6q8o5dv6t2sc7PpeiK+OkI++om77hdYmWZfmpzZfoX9hqZfrU+/WpWrhrBVdP0VQ+c1cbS9BfI9ENeJe94rG0LcWoqwUV9LbQGvcu+1s3ehgGW795xyvTn0zkKEi7qM4K+Y6Zv0/RbdNBff6SU/OjQKNfv6SHocyOE4f6oZtu0a7LVhmWjc2n6ogHagt66yTvbu8Irlnc2QxFX8YKdXfzDLVdy3YXOncnbwyroO2f6VtBfc03fOegnMrmKQu5C8k6+IBmeMlpO9EUDyx7VNXMbBnXe2WVo+rl8wdLy7UFf6fkX9UUB59F2vNy9k67PkonNedWBuWSOiViaK7e1WdsG2kNVM337cK1aAWZ0LkVvNEBbyLfuM3LV+13U17JieaeWFgwbBZdLcPO+LbirrPValHeM61q5ApjxFViLoC+ltDT9sfmU4xfcPhGrFnnn7EzSWiymtzVgTQZcCoWCJFeQTRv0laa/nDYMU/EM6jJOxorfceXW648GiAY8jolX0rZgTtjvIVeQluNuPWnOqw6cmzOCQF9rsUvkQFtwgaBfm7zTG/UTDRo2wrXyfjsxl8wiBFzY07KiPkKbTdNfjLZF5J1AjTbJ5TCbzJLNS/pbA6SyhYrPTS5fIJMvWBm+5d5ZwL9tXzegL+pf1g0+Y2a4zSvvLD/Tt4/YpmzuHTWybwt56W8NOmf6mTw+twuv22XNhq6HxNOcV52irbHfFvQH24NMxTOOk2PsF9sp009lDY+/knekZEUzYxfi7EyyYvLXTDJLNOCl3yxIl3uIa2WzafqLEQ14cLtE9ULuGvr0lbRzSb8x5C+XeFRGr4J9wEHe+e7jZ3lyZNZ6fGrKsGsOmfLO+Hx6ydZDFeyadUau8ukvR9NX1zQa8JTIO8qj3xby0tsacCzkJjI5q8ursg7bY8l6ST3NedUpFlr6WouzONWMTicHz+hcik5TH3a6O6sga8g7xbVZ14Jb7nyAv/nh0yXbZpNZM8swgv652YVbSlRjs1k2F0MIQXvIa2Ve6zk5SwWIvVuiJY8V9qUSwSjC+z0u62YgpeT2/3iCj3zvKet3Tk0m8Hlc9EUD9LYGKEiYiC1NSlRWxabN9Ffg3lG1vT190ZK/+6wpr7YGffRF/Y5zZeLpvBXs1efQHktu+vR/80/3Hl3yMS2V5rzqGJm+EFgLGgAMmKs1jcxUXrCxuRTbTZucU6FNefR7WwO0Bo2gvxYOnlg6x8nJRMmsTDBuMK3BYtBf1rDfnMDUSEEfil59WN9CrtLzL9lSW6avflY3g9lklvl0jodOTDFpvtbJibi1WIxawGap1zpj6shNq+m7lu/TVzfui/paSjN9m7zT1xpkIpaukI8SmZx1rS15x5yVm8rmeercHM+cW3nTx8VozqsOnJ9N0h3xl3zwVaZf7uCRUjI6l7ZsgU7yjjVysGf6axD0j5mLhpQHkNmkEfTVAg3LcfCo82okTR+KxVyXwFqfVLGWM3KLmb6xBGR50VUlD+Vr9qqgr+oQBQn3PGX0hzo9lbAayy33Wqtg1LT99FfQWnlsPk1r0Etfa4BkNm9dq5lklojfg9dtjMKkrBzZGYv4lGb6yruv1IX16M7btEH/3GyK/rIGXT0tftwuwchMaRYdS+dIZvNs74pYj8tRGl5v1E9rsHSZvtVErRRV/oFSQb8l4CXi9yyrFYM6r81k2awFVcwtXx8XwG9KHGuxOPr4fBq/x8VgexCf22Vl/gp1o7G3UQ763CSyKugbn0Ofx8UPnjyPlJKTk3G2dRjJhwr6S52g1fSFXMu9s3R5Z3w+TXeL35J6lX1zOpGxRvh9rYZ6UD4CS2Ry1pwM1e01Ztq/lYFkPVx/zXnVMbKj/rL1XT3mXbo801fOnf7WAEGv2zHTH51L4fe4aA16rYs/twZ37aPmmrCzyWyJDq2CPhiOpOVk+ktZNWsz0aGCvsN5CSHWrKe+ChBCCLpb/BVN15JOmb5N3hmeMj6Hr7tygF8em+DYuLEc31CXkel3hn143WLZ8k7TFnIt987yMv2eFr+1/rKSeGYTWWuEX20EFk/nrTYgEb+xr5J31KhuWgf9teP8bKrErqkYaK+0bapheU/Ub06frgwQ5+fS9LUGEEJYwXelhdzJWKUzQ2X6UMz2pZRWIRdY9qSdzbRq1lJoCxt/l/K+O4q1Whx9PJa2lsbrbvFXZvqWb7v49w55PTZ5J0FLwMMbr95KNi/5wi9OAsXmci6XoKclsOS++kV5pzm//lamvwzLppXpR1SmbwTpGdv3r980h5QHfUPTV/KO8VlUCaRKNNdjUmdTXvX5lFEg63cI+oMOs3KVft4bDRDxV8/0e1uM1/N5XEbTtbImaH/674/WbONMZfO85G/v5Qu/OFGy/dh4jLApB6jjiplzAlac6W+iVbOWgmrFUO28Ap61WTJxfD5Nd8QI+j0t/gobbcKSd4pfQ7u8MzydZGt7iCsG2+iN+vnmI2eA0sViepfg1VeWQCvTb3J5Z6mtlaWUjFuZvnFdp0wHz0wiQ5sp67aHvPg8rgrZLZ4pZvrKxaMSLZVozqdyy+7+WStNedWLdk3nTP/8XKqk8l7U6wPW9OlyRudS9NperzXoLXHv/OzZMb716xGrILcYR8dizKdz3H9kwtqWyRU4NZngmu0dAIybjiH1PupD198aYGw+teQPz2ZaQGUpqEJutVrFWq2epbJCMDL98vbKqo4QtGX6RiFXDfkTDLYHcbkEr9zbRyZXwO0SbLHVovpaaxvV/fDJc7zwYz/j3GzS0vSbNdP3LrO1sqrtdTvJO8ksrebnTAjheDNOZvJWpu9yCcI+t/WdU/UbWPtiblNe9eLErMqVlgbaghRkaXFsdC5N2Ocm4nfujiel5Pxsil6b/bM15CsZqp00e6D/5Jmxmo5RyTgHT89YM3tPTcbJFyTP29kJFDN99T7RYFFTXI5/u1HlHZXph3zO5xX0ule9kJvNF5hKZKyg39MSYDqRtbJsKK6QZdf0Q+YNSErJ8FTSWtP3hkv7AGMCoT1Y90Zrk3ceOD7FyEyS27/1hM70lzk5S8mpPS1GqwWvWzAZzyClZCaRpc38/oEhsdrNFFJK4jbLJlCSQI5MJy1TwVoXc5vyqp93mI2rsLz6NolndD5Fj1mciTgsczaXzJHOFUpGDm1BL7O21bNOm776+54dr2n692FzkfZYOscz5+eA4o3gmu2duERxQticbTag/byWOkGrYeWdsJJ3qmj6yyzkpnP5kgVN7KgeLfZMH2DCpusnTUkp5Kss5E7FMySzebZ2GJ/Ha4Y66Aj7GCpbLKYvGiCeyS86nf/ERByPS3Dvs+N87aHTQBMXclXDNQf3Tr4gq7ZPUUmWKs63h3xMxdPEM3lyBWklF2BM+rQnjqlsASlLE4+ImUBm8wXOz6W42Jy5vda6flNedXUH7on6K56zvPq2Yu74XNqaxBXyVS5+cN4m/yjK5Z2TkwmCXjfz6RwPn5xa9BgPj8asIH7gpLG+6jHTubO7J0JXxG9lHmo42BosdQ8s1cpnWTarZMSblcXkneUWcv/mB89wwyfvcwy46trYNX37djB67AhRtI2CcQNKZPKWm0Nl+h63izvfdhV/edPFJe+jEo3FajgnJuLcsLePa4Y6rDWhvZ7m9Okv1Fr5z7/xGLfc+YDj743bgj4Yq7JNxTNWZq7kHcDoizRbbLKnEkW7mSASMIL++dkUBQmXDhhBf1oH/dXn/FySrogPv6cy89viMEFrdD5lBdKIg6Y/6hD020LF9spSSk5Nxnn15f343C5+VoPEc2Rsnufv7KS/NWDdJI6OxdjSatQVeqJFjXi2aqa/9KAf9rlxVelYuVlpW6yQa5sQpfj7nxzhw//1lOP+YAzB73pomFS2wKOnZyqeLw8QKsGwT6pLZo1Vs+xzB5S8M2xqvIPtRQnyqgs6uLC3peR9arnBZ3IFzkwn2Nkd5v99/T5rglqzZvpul7N75+jYPN9+dIRHh2coOGT7Y5a8Y1zLroifyXimOBvXJu/0RgOkcwXru6lWzbJn+mGfEUvUDf5ScxLfWts2m/Kqn6ti1wQjAHRFfNaFMGbjGt0zAcdCrsr0++yZfshrZeATsQyJTJ5L+qNcu6NjUV0/lc1zeirB7p4W9g91cODkNFJKjo7H2GkuD9jTEqjQ9FWm3xH24XO7luzgiaUaq8Omoi3opSPsY2uVJRUDXndFi9v7j0zwnUdHqr7mVx88TTKbRwg4cKpy5FYe9NX/9mJuwrZUoiLo8yBlUcqzB30n+mqYlXt6KkFBwlBXmKGuMH/xqovxe1xWMbLZEELgdYuK1sp33HscKY2lJJ2K4+PzaXxuV8n3bCqesSVdxb+nqheqxCuRVaPoUk0/ls5bqsKlA0bQL2+muNrUHPSFEG4hxEEhxHfNx18VQjwrhHhSCPEFIYTX3C6EEJ8WQhwVQjwuhHiO7TVuFUIcMf/duvqnUxvnZ1P0Rat/mfYNtnHv4TGy+QJzqRypbMHKqMJ+D/FMviQTGHWQi1qDXjK5AqlsUfe9oCvMy/b0cHw8zskJZy0YjC+8lHBhbwv7L2jn/FyKM9NJjo3FrTVhe1r8VtCfTWbxuV1WABFC0Nu69La7jdZsTeFxu7jvAy/lzddsc3w+6HVVZPpzqSyT8YzV88ZOOpfni788yYt2d3FRbwuPnJqu2Ed58rsi/pL/x8sz/fL+/mYWfsSU91oCXhait4b+O+qzpnpHve15Qzz2oVeUBKlmw+NylQT9M9MJvvPoiNUR9YTD93NsPmXp+WAG/VjGyszb7PJOWV1NScIhv13TdxNLZy3nzq6eCB6X2FCZ/vsAe2vHrwJ7gMuAIPAH5vYbgd3mv9uAOwCEEB3Ah4BrgWuADwkh2ldy8Mvl3GzKsYiruOWabYzOpfnxU6OWLVJlaqoYmLBpwOOxNNGAx+qWCEX75EwiW+yB3hHi+j29APx0gWz/yJhRxN3dG2H/kPEnuvuxsySz+ZKgPxlLky9IZpMZokFviUzQHw1yzqFx3EI0atAHQ9qpttCKUyFXFceP2CbDKb7z6FnG59Pc9uIdXD3UUeKwUozPl34mvG4jsy6RdxwyfTX8Pzw6z9Z255FJybH73EQDngXlnRNlQR8o+aw2Ix63KDFU3HnfcQD+6jf3As5Bf3w+TZfNodcZ9jGfzlmGCru8o/7Wx8eN10lY6+OWavrxdJ6R6SQ9LX4CXjdtIe/G0PSFEIPATcDn1DYp5felCfAQMGg+9Vrgy+ZTDwBtQoh+4AbgHinllJRyGrgHeOUqnovFbDLLz54Zc8zSEpkcs8ms1XfeiZfu6WGgLchXHjxltWCwZ/pQ2nRtMp6xMjlFselahlOTcVzCKMpt6wyxqyeycNAfjeFxCYY6w+zpixLxe/j3A8MA7Ow2gn63acucjKVLZuMqLuyL8NDJKd74mV/xvcfP1eQYarQFVGrFyac/ZzqZjoyWdj2UUvK5+4+zp6+FF+7qYv9Qe4nDSmH36Ct6WvwVmX75Sl5qofYTE/FFpR3FYpPxTkzGaQ95mzqzL8frdlmF3PH5NHc9PMzrnjPAlVvb8XtcjiNxNTFL0WHOylU3iKgt6HeEfbSHvBwzg76V6VfIOzlGZpLWtW4NejeMZfOTwAeAishhyjpvBX5obhoAhm27nDG3Vdte/nq3CSEOCCEOjI+P13h4pZyYiPP2Lz7MQYcC20J2TYXbJbjl2m384ugkvzo2CVBSyIXSpmtTsUyFPmq1V05kOTmZYEtb0PJFv2xPDw+emKxqszs8GmN7Vxifx4XbJbhyW5s1WlCZvnKFjM2nrbbKdj5448XcfuMezs4kefe//Zo3f9bZkWCn0RZQqZWg100mV7Cy9Vy+uMrV4dHSTP/nh8c5PBrjD1+0AyEEV11gjMSUw0rhFPS7bZIcGJp+ecatGnLlCrLmoN8bDXBsPFZ1stGJ8ThDXWHH55oVj0tYDdfuvP84mXyBd163E5eZbJ10sOKWX1PVdO34RIyg111xLXd0RzhuOu5Upm93xkV8HjK5Aicm4gyYo7r2svk9a8GiQV8IcTMwJqV8pMou/wTcJ6W8fzUOSEr5WSnlfinl/u7u7mW9hipujc5XZj/FFsgLf6HeuH8rXrfgX8w2COoO77TizWQ8XTXozySznJpKlPirX3xhN9m85ODpSi0YDHnH7tK4esiYgdsW8lofNFU/GJ83M/2yoB/xe/ij63by8z9/KX903Q4OnJqu6MxZzmZbFH21KF9IxX4zPlyW6X/n0bN0hH28+vItgGHx7YsGOFCm6xt9d0oTC6Ppmt27XZnp2x9XKzyX8+rLt3BsPM7t33rCcfWlk5NxtnfqoG/H63aRzUt+9swYd95/nNddOcgOcxQ91BWqkHey+QKT8Uxppm+2Yjg+HrdswXZ2doetTF+10Q6VWTbBkJuVVbwt5NsQmv4LgNcIIU4CdwHXCyG+AiCE+BDQDfypbf8RYKvt8aC5rdr2Vacr4sMlcJyp6LRMohPdLX5eeWk/8UzemokLziveTMUzdJbJO/ZM/9RknAts/VK2mV9mJ809mTGdO70Ra9t+M5vc1R2xdPsemxvE3mGzHLdLcP1FPQAly+45sdkWRV8tVDFVBf25ZHGmbLmm/9CJKa7d3mGN2oQQXDXUziNlcy/sfXcUPS0BxmNpKzAnMpWF3IDtca2Z/hv3b+V9L9vNNx45w0e/93RJ4E9m8pybTZXo+RpD0z86FuO9XzvIxX1RS8sHw+U0PJUsqdOoRdC7S4K+kYCdm03R6iCd7eiOMGHKr06Zvj3BUpNC221W77Vi0aAvpbxdSjkopRwC3gT8VEr5FiHEH2Do9G+WUtrHlXcDbzNdPM8FZqWU54AfAa8QQrSbBdxXmNtWHY/bRVfE2b1i2SsXCfoAb7nWcHvYXTkRS9M3AkShII2gX5bpK4399FSCmUS2JOir13M6vmPjhnNnd08x079iWxselyi5EVgWwLk0s4lsiZ5Yzt6BVoSAx89UD/pSyubV9D2lC6nMmU3xLt/aylQ8Y82iHZlJMjKTtEZeiv0XtHN2NmUthJHI5Iilc46afjYvrS+1Ucgt/XuHSoJ+bZk+wJ+8fDe3Pu8CPvffJ7jj58es7Uqm0PJOKR6X4NHhGXweF3feur/EP7+9M0wmXyhZNlVZbXtsozf7d758pA3F+tvx8ZgVL+yFe7uUqm7wRiG3/pl+Nf4Z6AV+JYR4VAjxf5rbvw8cB44CdwLvApBSTgF/BTxs/vuwuW1N6I0GrCKsnXOzSdpD3prcC9ds7+Di/miJNBMqa4k6k8xSkFTIO8ot8tiZGaDYDhfA73HTEfY5Bn3l3LnQFuBDPg+feetV/PF1u0peoy3k5dxcivl0rqKQW34sO7rCPLFApp/KGpq26vPdTAQqMn0jKO+/wAjuSuJ5+ITxcVUN7xRqPyXxTMxXZoX2x8rOaVg2S7+CIW9lIKgFIQQfevVeXrm3j0/ec8QqBjo5dzSGvON1C/75rVdZ0opC3SDtEk/5vAswRvPKEeb0/VMr7R0bj5PI5Ah6Syc+lgR9m7yTNq3ea8WS0jop5b3AvebPjr9runneXeW5LwBfWNIRLpPeaKCkc53C6KNf25dJCMFX3nENLlF5oZS8M2WunKP6a9t/tzXotbLr8p4p1RplHR6N4XWLiszsZRf3Vuzb0+K3JvFUk3cU+wbb+OWxiarPFztsNp+Vz1onN2MMWFWmf5Vplz0yGuP5O7t46OQULX6P1SNFcXF/CyGfm0dOTvGay7cwHiu1+SqUJHdiIs6FvS0lXRcVAfMm0Bn2VW0QVw2XS/Dul+7ih4fO870nzvG7115gBS6d6Zfynut3EfS6K0ZtADvMv9XJyTgvxqgrls/GBePv3R7yMhHLOAb9bR0hPC5hZPq2tsoKZ3nHiCPTiYxjQ8jVoGFn5Pa1+h29y4t59MvpjPithl1QvFBKo1NaX2e4so9Pm63/zrayolxflT7oR0znTi1tb3taApalcKFMH4zZfqNz6Yr2vopG7bBZC+WLoytNf1d3hGjAU5LpP+eC9gq/v8ft4oqtbfzq+CRfeeCU1b6hpyzo7x1oZaAtyF9++0lOTcZJZh3cO2agH6yxiFvOpQNRdvVE+PavjXLZiYk4PS3+pnRlLcTN+7Y4JlJg3KzDPndJpq9qgeXJnRrhqyVS7XjdLrZ1hjg2HnO8watr0h7yWs+pgvB0fO10/YYN+r1mK9vyYdJCLRhqQVnq1OpZauUcpyntqgFTb9RfUbDraw043pSOjM2X6PkL0dPityZyLJbpX2ZO8a5WzC122Gw+eUdJLOWafmvIy4W9LRwZjTEVz3BkLFYh7Sj2D3VweDTG//zPJ4ln8nzglRdxUVmfnIjfw5d+/2qy+QJv+fyDABXuHXUDWoq0Y0cIwW9dOcCBU9OcnkxwckLbNZeKEIILOsMlQf/eZ8e4bKC1ol+X+t5XS7p2dkc4Nh4nns5VXGuVYNlrN622+T1rReMGfTOw21crSmWNlrV90eUHfbX4gdL0VdDvijgEfTMQX9BR+aXrjQaYiGVK+qtbPXdsev5ClOqLC0+82bslumAxt1EXUKmFgCXvFDV9IQwf9e7eFg6PzVtN75zkAIC3PHcbH7xxD99/74u45/0v5l0v2eXYuG5XTwtf+L391ueyfEau2yXoafFzcV9tN34nfvNKY/rLtw+OcGIibskVmtrZ3hW2Jmidmozz+JlZbt7XX7GfGuE7FXLBCPqnJuPMpbKVQd/M7u01BSXvrKWDp3GDvoNXf3jK0PjLpZalYm+6ppZLa3fI9NUHwe7cUagbj11uOTERR8pi1X8xyotKix3zzu5I9UxfB33SOZXpG9ZVl0twYW+EmUSWHzxxDp/Hxb7BVsfX6GkJ8M7rdnLJlmhJOwwnrrqgg79/85W4XaKkM6viB+97Ebe9eOeyz2egLchzd3Rw18OnmYxndKa/DIa6QgxPJ8nmC3z38XMA3OQQ9IuZvnPStaM7TDYvefb8fMUcGKXxD7RXBv21dPA0bNB36j64WkUttfgBGBOzjFV0Kv+U6oPg9H5qJGKXeFRjtlqdFj1lrZwXY99Aa9VM/8SEURDuamm+qfpBh0xfWWDVJLnvP3GeKwbbVq1nzSv29vHr//Ub3LC3UlfujPhXvKrV664ctHRo7dxZOkOdYfIFyZnpJN99/BxXbmtztNDWIu+A0SO/PNP3uF38f2+4nFufN2Rts9q36Ex/6fQ59BlXnuWVzk60Z/pOfXcUKnA4jSyKN6Wi/HRiwmzM5jAycKJnCZk+GMXcsfm0Yy3hvx47x77B1jVzDGxkKgq5qeJkNyW1ZfIFrt6+uv0BW8ua5K0mN17WZy3OooP+0lF/s588PcrT5+a4ed8Wx/1UYbd60C/+7Z0WJ/rtqwbZZvu+B7xuAl7XmvbfadigHw168JetSH9iIkF7yFuyws1yCPuLq2c59d1RKHmn3K4JtqBflul3RXyLttNVFFtDuGty+1xmShNPlGX7x8djPDEyy2sud/5gNzqqyG5370TNa9Ad8Vtf6Gp6/kakJeDlNy7pxe0SK5YzmxE1Or/z/uMIATddVintADx/Zxcv29Pj+B0HY7SvJnGFarRDt4d8a9pps2GDvhDC6D5oK+SenIivStYT9nlsPv3qQf/5uzr5jUt6HQuzbSGvw00pXvXD44SSd2rJ8gEu6Y/iElRM0rr7sbMIgdVPptlQGbFaHH0ulSUaNLIyIQQX9rTgEljN1TYL//OmS7jzbVc1fRvl5dAZ9tHi9zA6l+bqCzqqOv529UT4/O9dveDfWE3SqnUZ0raQT2f6y6W3JVAh76xGUStsWxx9Mp6u6Luj2NMX5c637Xf8QFg3pdnS47tgCUE/4vcQ8rkd+35UO+7yYq6UkrsfPcu12zsci4rNgBCipKf+XDJrZfoAr768nzdctbXmEdhGoa81YK3foFkaQgi2m8H65suds/xaUbp+uW27Gu1r3FO/sYO+zQtvNZ5ahW6DStMvFCTTiWxF352ajy9aDPqJTI7RuTTbu5Y2FO9p8dMarN1xc9lgK4+dmbGcKofOznF8Is5rr6joct1UBH1uUlk1IzdX0svorc8b4mOv31evQ9PUiaHOMC4BN166sqC/9Ex/bXvqN3TQt69Iv5qNpyKmpj+bzJIvyGWvNdoXDViavrW61hJvSm/Yv5WbqhSZnHjVpf1MxDK8818fIZ3Lc/djZ/G6BTde2rek9200Ah4XyWze6qUf3WRZvWb1+cMX7eBvXrevop3GUlGZfq2aftsa99RvaFO2WpF+LpmrWCd0JYT9HpLZvNU4q3xqdq30tQY4f8i4KS3Vrql490t3Lb6TjZdf0stf/9Zl/MW3n+Cd//oIz5yf58W7u5t+VSW1epaq1USXMHrSNCaXDbZa5oeVcMmWKD63q+auqe0hLzPJLFLKNXF3NfQn275o9IlVzfSNP5ua7OXUd6fW48vkCswksku2a66EW8yW0X/x7ScA+OCNe9b8PTc6Qa+bVCZv9d3Rmb5mtehvDfLQX76sZsNFW9BHviCZX6MRZ9ME/ZMTcboiq9N4Ss2sU5LMSuSd0uOr3a65Um65dhsel+A/Hx3hNy7RxT5VyFV9dxZan0CjWSpLGUlbE7TiWR30l4p9gtbJicSSi6TVUEH/tJnpO/Xdqen4WouLqZycXJpdczV449VbeePVWxffsQkI+tzE0jmrK2q0CbuNajYG9lYM29Zg5N/QhVy1QtXorCHvrFZQVT3nlbzj1HenFqz+QLOpJds1NatLwGu4d9QCKjrT19SL9rDZXnmNHDwNHfQDXmN1qeMTccbn06vWeEpZr05NJar23akFtfTaiYn4suyamtXDCPpa3tHUH9UxV406V5uGDvpgSDwPHp8EWLUWs2FbIbda351a8HlcdEV8PGguw6cz/foR9LpIlhRytbyjqQ/FhVR0pr8seqIBzpoToFYt0zeDfjpXWHYRV9EbDVhtEXRjrPphL+S6RO0TaTSa1Ua5fNZqVm7DB/2+aDETXy1N377W5UqDfl80QL4ggfWxa2qcUT79uWSWloDXcQEUjWY98LhdRAMeLe8sF1Us7YsGau59sRh222e1vju1ovrqr6ddU1NJ0OsmkyswncjqiVmautMW8q1ZIbfhP90q6A+tYpE06HXjElCQLLvvjkLZStfbrqkpRfXUH59P64lZmrrznpfuWrMFjRo+6Kugupp6uRCCsM/DfDq3cnnHzPR1Ebe+qFHg6HyK3pbm7Daq2Tis5fyZppF3VjuTVsXc5fbdURRvSlrPrycBjxH0x+bSWt7RNDQNH/R390a48dI+Xr7KrQZUMXe5fXcU27uM9q17t6y8sZNm+QTMTF932NQ0Og2f0gS8bu54y1Wr/rqqmLtSeWdrR4hffPB6K+PX1IegbaEbPTFL08g0fNBfK0Kmj3u5fXfsNONi5BuNkqCvM31NA9Pw8s5aoTT95fbd0Wwsgr7iV0Fr+ppGRgf9ZRLxu1fUd0ezsQjoTF/TJOiUZpn81nMG2dMfrfdhaFaJgNb0NU1CzWmqEMIthDgohPiu+fg9QoijQggphOiy7SeEEJ82n3tcCPEc23O3CiGOmP9uXd1TWV+uu7Cbd163s96HoVklSjV9nQtpGpelaBPvA562Pf4F8HLgVNl+NwK7zX+3AXcACCE6gA8B1wLXAB8SQrQv77A1mtVFu3c0zUJNQV8IMQjcBHxObZNSHpRSnnTY/bXAl6XBA0CbEKIfuAG4R0o5JaWcBu4BXrnSE9BoVgN7XyYd9DWNTK2Z/ieBDwCFGvYdAIZtj8+Y26ptL0EIcZsQ4oAQ4sD4+HiNh6fRrAy/x+be0fKOpoFZNOgLIW4GxqSUj6zD8SCl/KyUcr+Ucn93d/d6vKVGgxDCaqSne+lrGplaMv0XAK8RQpwE7gKuF0J8ZYH9RwB7t6BBc1u17RrNhiDoc+te+pqGZ9GgL6W8XUo5KKUcAt4E/FRK+ZYFfuVu4G2mi+e5wKyU8hzwI+AVQoh2s4D7CnObRrMhCHhcemKWpuFZ9swiIcR7hRBnMDL2x4UQqsj7feA4cBS4E3gXgJRyCvgr4GHz34fNbRrNhiDgc+uJWZqGZ0lpjZTyXuBe8+dPA5922EcC767y+18AvrDUg9Ro1oOg102LLuJqGhz9CddoTN71kl0lPXg0mkZEB32NxuSmff31PgSNZs3RaY1Go9E0ETroazQaTROhg75Go9E0ETroazQaTROhg75Go9E0ETroazQaTROhg75Go9E0ETroazQaTRMhjK4JGxMhxDiVK3MthS5gYpUOZ7PQjOcMzXne+pybh6We9wVSSsfe9Bs66K8UIcQBKeX+eh/HetKM5wzNed76nJuH1TxvLe9oNBpNE6GDvkaj0TQRjR70P1vvA6gDzXjO0Jznrc+5eVi1825oTV+j0Wg0pTR6pq/RaDQaGzroazQaTRPRkEFfCPFKIcSzQoijQogP1vt41gIhxFYhxM+EEE8JIQ4JId5nbu8QQtwjhDhi/t9e72NdC4QQbiHEQSHEd83H24UQD5rX/OtCCF+9j3E1EUK0CSG+KYR4RgjxtBDiec1wrYUQ7zc/308KIb4mhAg04rUWQnxBCDEmhHjSts3x+gqDT5vn/7gQ4jlLea+GC/pCCDfwj8CNwCXAm4UQl9T3qNaEHPB/SCkvAZ4LvNs8zw8CP5FS7gZ+Yj5uRN4HPG17/DHgE1LKXcA08I66HNXa8Sngh1LKPcDlGOfe0NdaCDEAvBfYL6W8FHADb6Ixr/UXgVeWbat2fW8Edpv/bgPuWMobNVzQB64Bjkopj0spM8BdwGvrfEyrjpTynJTy1+bP8xhBYADjXL9k7vYl4DfrcoBriBBiELgJ+Jz5WADXA980d2mo8xZCtAIvBj4PIKXMSClnaIJrjbGka1AI4QFCwDka8FpLKe8Dpso2V7u+rwW+LA0eANqEEDWv9dmIQX8AGLY9PmNua1iEEEPAlcCDQK+U8pz51Hmgt17HtYZ8EvgAUDAfdwIzUsqc+bjRrvl2YBz4F1PS+pwQIkyDX2sp5Qjwd8BpjGA/CzxCY19rO9Wu74piXCMG/aZCCBEB/gP4EynlnP05afhxG8qTK4S4GRiTUj5S72NZRzzAc4A7pJRXAnHKpJwGvdbtGFntdmALEKZSAmkKVvP6NmLQHwG22h4PmtsaDiGEFyPgf1VK+S1z86ga6pn/j9Xr+NaIFwCvEUKcxJDursfQu9tMCQAa75qfAc5IKR80H38T4ybQ6Nf65cAJKeW4lDILfAvj+jfytbZT7fquKMY1YtB/GNhtVvh9GIWfu+t8TKuOqWN/HnhaSvlx21N3A7eaP98KfGe9j20tkVLeLqUclFIOYVzbn0opfxf4GfB6c7eGOm8p5XlgWAhxkbnpZcBTNPi1xpB1niuECJmfd3XeDXuty6h2fe8G3ma6eJ4LzNpkoMWRUjbcP+BVwGHgGPCX9T6eNTrHF2IM9x4HHjX/vQpD3/4JcAT4MdBR72Ndw7/BS4Dvmj/vAB4CjgLfAPz1Pr5VPtcrgAPm9f5PoL0ZrjXwfwPPAE8C/wr4G/FaA1/DqFtkMUZ276h2fQGB4VA8BjyB4W6q+b10GwaNRqNpIhpR3tFoNBpNFXTQ12g0miZCB32NRqNpInTQ12g0miZCB32NRqNpInTQ12g0miZCB32NRqNpIv5/NWUgEKJWm1sAAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"Y = np.hstack((np.random.normal(4200, 50, size=26),\n", | |
" np.random.normal(4300, 50, size=20),\n", | |
" np.random.normal(4200, 50, size=20),\n", | |
" np.random.normal(4500, 50, size=20),\n", | |
" np.random.normal(4200, 50, size=14)))\n", | |
"\n", | |
"plt.plot(Y)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## SPM\n", | |
"\n", | |
"In SPM we model each condition of the task with a binary regressor (a value of 1 for each volume recorded during the task, 0 otherwise), that gets convolved with a haemodynamic response function (in order to model the slugishness of the recorded signal).\n", | |
"\n", | |
"Another regressor (with the value 1 at every volume) is added, the constant term, to model the fact that the recorded signal is in an (arbitrary) positive range. \n", | |
"\n", | |
"This results in the following design matrix:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 114, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"<matplotlib.image.AxesImage at 0x1f52472e0f0>" | |
] | |
}, | |
"execution_count": 114, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAGMAAAD7CAYAAAB+Diq7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAKkklEQVR4nO2dXYhc5RnHf/+dySTupjF+JaTZ0CTko4htUZZUEYrEFsRK04JI2hJsm5KbWmMp1LQ39sILhaL1SgjakgshFRWUVioljZHeBGNisUlYjVs/siTG1iS6+dyPpxdz1K3uzJ7dOefMM2eeHxxmzsc875P89/++Z86c9zkyMwIf9LQ7geBTQgxHhBiOCDEcEWI4IsRwREtiSLpF0qCkI5K2ZZVUt6LZfs+QVAFeB74FHAVeBr5vZoeyS6+7qLbw2XXAETMbApC0E9gANBSjVum1S+Zc2kKTxbJo9anMY54YHuX0B2Oaal8rYiwF3p20fhT4+mcPkrQF2AIwr7qAG5bf2UKTxXL3s3/OPObWDW823Jf7AG5m281swMwGapXevJvraFoRYxhYNmm9P9kWzJJWxHgZWC1phaQasBF4Lpu0upNZjxlmNibpLuAFoAL8wcwOZpZZF9LKAI6ZPQ88n1EuXU98A3dEiOGIEMMRIYYjQgxHhBiOCDEcEWI4IsRwRIjhiBDDESGGI0IMR4QYjggxHBFiOCLEcESI4YgQwxEhhiNCDEeEGI4IMRwRYjgixHBEiOGIEMMRIYYjQgxHhBiOCDEcEWI4IsRwRIjhiBDDEdOKIWmZpN2SDkk6KGlrsv1ySX+T9Ebyeln+6ZabNM4YA35pZlcD1wM/k3Q1sA3YZWargV3JetAC04phZsfMbH/y/iPgMPVSFRuAHclhO4Dv5pRj1zCjqceSlgPXAnuBxWZ2LNl1HFjc4DP/VzskaExqMSTNB54G7jGzD6VPC8OYmUmaslaSmW0HtgNceskSQ1MWlGmNkpSDTXU2JWkOdSGeMLNnks3vSVqS7F8CnMgnxe5hWmeoboHHgcNm9tCkXc8BdwIPJK/PThfLesRE79xZptqYnjPnM4/ZDtJ0UzcCm4DXJL2abPsNdRGelLQZeBu4I5cMu4hpxTCzfwCNOvqbZ9KYzenhXH/fTD6Sir7BcjgjvoE7oqWqOjNlbJ74YG32TfYNZh6yLYQzHFGoMybmGSNrR4tssqMIZzgixHBEod3U3LmjrFl5bPoDZx45h5jFE85wRKHOWFA9z/pF2Z+H7uGrmcdsB+EMRxTqjL7KBW7oeyPzuOGMIHOKPZtinJXVkSKb7CjCGY4o1BlV9XBVpRzfCfIgnOGIEMMRhXZTPYi5mlNkkx1FOMMRIYYjQgxHFDpmGMa4TRTZZEcRznBEwc6ACzZWZJMdRTjDEcXeHYJx1uLukEaEMxxRrDPMOFuSuRR5EM5wRIjhiIIHcHHecphGVhLCGY4o2Blw3ipFNtlRhDMcMZOpxxVgHzBsZrdJWgHsBK4AXgE2mdnFZjEMhTOaMBNnbKVeHeFjHgQeNrNVwElgc5aJdSNp54H3A98GHkvWBawHnkoOSVWuwoBRq2S+lIW0zvg98CvqYzDUu6ZTZp9cgj1KvZ7I55C0RdI+SftO/Td+y2hGmkn5twEnzOwVSTfNtIHJ5SrWfmWejVqhJ3AdRdpJ+d+RdCswD1gAPAIslFRN3NEPDOeXZneQpsTRr82s38yWAxuBv5vZD4HdwO3JYanKVQTNaaXPuBfYKel+4AD1+iJNuUiV4bEo2NaIGYlhZi8CLybvh4B12afUvRQ6mo6Mz+Wl02uKbLKjiMshjijUGR9dmMeet1ZlHnclZzOP2Q7CGY4o9hvY+R5scH4OgcMZQcYU6ozKBViY/czj0hDOcETBzphgwdC5IpvsKMIZjggxHFFoN6WL49SGT2Ye1+aU4zeScIYjiv2TGh/HTp7KPu6iK7OP2QbCGY4odhrZxDgTI2cyj9sTzgiyptgxw8DGYoJlI8IZjggxHBFiOCLEcETx1xHyeBpZSQhnOKJYZ0ioGpXYGhHOcESxl9AlVAtnNCKc4Yjix4xardAmO4lwhiMKdwbVcvxEmgfhDEcU7AxQNfupwmWpYBXOcESI4YhU3ZSkhdSrI1xDvVf4CTAI/AlYDrwF3GFmzW+KigG8KWmd8QjwVzP7MvA16jVEtgG7zGw1sCtZD1ogTYWES4FvAD8CSCrnXJS0AbgpOWwH9Vmw904TDSrRMzYizf/MCuB94I+SDkh6TFIfsNjMPn426HFg8VQfnlw75OJ4OWYY5UUaMarAdcCjZnYtcIbPdElmZjQ4wzSz7WY2YGYDtUpvfdzolKVg0ohxFDhqZnuT9aeoi/OepCUAyeuJfFLsHqYdM8zsuKR3Ja01s0HgZuBQstwJPEDa2iECy2HMmOgtxxPO0p5n/hx4QlINGAJ+TN1VT0raDLwN3JFPit1DKjHM7FVgYIpdN2eazSw519/X7hQyIc4zHRFiOKIU9019sLYcl1jCGY4oxZ/UyNpyPK0mnOGIgp2Rz2WGNSuPTX9QBxDOcEQpxoz1iwbbnUImhDMcUQpn3NCXTxGrsxPFXoAMZziiFM5YWR3JJe6/LoYzupYQwxGl6KauqpTjl75whiNK4Yy5KsfUtHCGI0IMR4QYjijFmDFu5XjkXDjDEaVwxgUrR3W3cIYjSuGMsxY3JAQZE2I4oiTdVDlmgoczHFEKZ5y3ctQ9DGc4oiTOyL4eSTsIZzgibbmKXwA/pT69+DXqc/qWADuBK4BXgE3JhP3C6RpnSFoK3A0MmNk1QAXYCDwIPGxmq4CTwOY8E+0G0o4ZVeASSaNAL3AMWA/8INm/A/gt8GjWCaZhtFucYWbDwO+Ad6iLcJp6t3TK7JPLpUeBpVN9PspVpCdNN3UZsIF6DZEvAn3ALWkb+Fy5iqAhabqpbwL/NrP3ASQ9A9wILJRUTdzRDwznl2ZzRq0UZ+ipTm3fAa6X1CtJfFquYjdwe3JMunIVQVPS1A7ZK+kpYD8wBhwAtgN/AXZKuj/Z9nieiTZjeOyyXOJeUcnnhupGpC1XcR9w32c2DwHrMs+oiylFZ/vS6TW5xP3e5ftziduIuBziiFI4Y89bq3KJG87oYkrhDBucn0/g6/IJ24hwhiNK4YyF+cw8LpxwhiNCDEeUoptaMHSu3SlkQjjDEaVwRm24+ZMiOoVwhiNK4Qw7eardKWRCOMMRpXDGxMiZdqeQCeEMR5TCGTYWEyyDjAkxHBFiOCLEcEQpBvB2PDksD8IZjiiFM1SNSmxBxpTDGbVwRpAxJXFGrd0pZEI4wxGlcAbVcvwzwhmOCDEcUQp/q9ol88CD4iiFM2IADzJHVmB9P0nvA2+nPPxK4D85ppMlM8n1S2Z21VQ7ChVjJkjaZ2YD7c4jDVnlGt2UI0IMR3gWY3u7E5gBmeTqdszoRjw7o+sIMRzhTgxJt0galHRE0rZ259MMScsk7ZZ0SNJBSVtbCmhmbhbqlUHfBFYCNeCfwNXtzqtJvkuA65L3XwBebyVfb85YBxwxs6GkRu5O6vURXWJmx8xsf/L+I+AwDQpnpsGbGEuBdyetN6wK6g1Jy4Frgb2zjeFNjI5E0nzgaeAeM/twtnG8iTEMLJu03taqoGmQNIe6EE+Y2TOtxPImxsvAakkrJNWol/l+rs05NSSpZvo4cNjMHmo1nisxkhq5dwEvUB8MnzSzg+3Nqik3ApuA9ZJeTZZbZxssLoc4wpUzup0QwxEhhiNCDEeEGI4IMRwRYjjifzVYxGqMbW98AAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"X = np.array([np.convolve([0] * 20 + [1] * 20 + [0] * 60, hrf(np.array(range(100))))[:100],\n", | |
" np.convolve([0] * 60 + [1] * 20 + [0] * 20, hrf(np.array(range(100))))[:100],\n", | |
" [1] * 100]).T\n", | |
"model = GeneralLinearModel(X)\n", | |
"\n", | |
"plt.imshow(model.X, aspect=0.1, interpolation='none')" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The model is then fit to the data, and we can test for the effect of _task_." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 115, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"A > rest: 99.13 (z = 5.55)\n", | |
"B > rest: 261.99 (z = 11.14)\n", | |
"B > A: 162.87 (z = 7.02)\n" | |
] | |
} | |
], | |
"source": [ | |
"model.fit(Y)\n", | |
"con1 = model.contrast(np.array([1, 0, 0]))\n", | |
"con2 = model.contrast(np.array([0, 1, 0]))\n", | |
"con3 = model.contrast(np.array([-1, 1, 0]))\n", | |
"\n", | |
"print(f\"A > rest: {con1.effect[0,0]:.2f} (z = {con1.z_score()[0]:.2f})\")\n", | |
"print(f\"B > rest: {con2.effect[0,0]:.2f} (z = {con2.z_score()[0]:.2f})\")\n", | |
"print(f\"B > A: {con3.effect[0,0]:.2f} (z = {con3.z_score()[0]:.2f})\")" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## FSL\n", | |
"\n", | |
"In FSL the signal is demeaned before fitting the model:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 116, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"[<matplotlib.lines.Line2D at 0x1f5247796d8>]" | |
] | |
}, | |
"execution_count": 116, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABQOklEQVR4nO29eXhkd3nn+/nVvkiq0i61el/sttvttk3bmCUkwRCMYTAQQiAJMAkZ38yQ52aGmZsLydyZSZ5hJjOTCXeykWuWYJgk4LAECA5LDMRstmmv3e2lu92r1Iv2kmrffvePc36nzqlFKqkklVTn93mefqQ6VV11SqfqPd/zfd/f+wopJRqNRqNxF55274BGo9FoNh4d/DUajcaF6OCv0Wg0LkQHf41Go3EhOvhrNBqNC/G1eweaYWBgQO7evbvdu6HRaDRbiieeeGJaSjlY774tEfx3797NsWPH2r0bGo1Gs6UQQlxodJ+2fTQajcaF6OCv0Wg0LkQHf41Go3EhOvhrNBqNC9HBX6PRaFyIDv4ajUbjQloO/kKIkBDicSHEM0KIk0KI3zO37xFCPCaEOCOE+LwQImBuD5q3z5j37251HzQajUazMtZC+eeA10opjwC3AHcLIe4E/hvwUSnlfmAOeL/5+PcDc+b2j5qP02g0mg3h6UvzHB9PtHs32k7LwV8aJM2bfvOfBF4LfMHc/gDwVvP3e83bmPffJYQQre6HRqPRNMPvfe0kf/CN59u9G21nTTx/IYRXCPE0MAl8G3gJmJdSFs2HjANj5u9jwCUA8/4E0F/nOe8TQhwTQhybmppai93UaDQaphZz5Arldu9G21mT4C+lLEkpbwG2A3cAB9fgOe+XUh6VUh4dHKzbmkKj0WhWzEwyT6GsJxiuabWPlHIe+C7wCiAuhFC9g7YDE+bvE8AOAPP+GDCzlvuh0Wg09Ujni2QKJYolrfzXotpnUAgRN38PA68Hnsc4CbzDfNj7gK+Yv3/VvI15/3ekHiSs0Wg2gJlkHoBiSYectejqOQo8IITwYpxMHpRS/r0Q4jngc0KI/ww8BXzSfPwngc8KIc4As8C71mAfNBqNZllmUkbwL5S18m85+EspnwVurbP9LIb/X709C/xCq6+r0Wg0K2U2lQO08ge9wlej0biIacv20cpfB3+NRuMalOevq3108NdoNC6iYvto5a+Dv0ajcQ262qeCDv4ajcY1TOtqHwsd/DUajWvQ1T4VdPDXaDSuwbJ9yhK3ry3VwV+j0bgCKaUV/ME4AbgZHfw1Go0rSOaK5Etl+qIBQFs/OvhrNBpXoFT/cE8I0ElfHfw1Go0rUH19hnuCgFb+OvhrNBpXMJM0Kn1GTOXv9oVeOvhrNBpXUFH+yvbRyl+j0Wg2BRPzGX7ny8fJFUtr/tyzVcFfK3+NRqPZJHztmcv89WMXefri/Jo/93QyR1fQRzToBaCgPX+NRqNZHVJKHn7+Gm/78x/yv/7xdMvPd/LyAgAvXF1s+bmqmUnm6YsG8HuNsFfU1T4ajUazcp4dn+cX/uLHvP+BYzx1cZ4fvjTd8nM+dzkBwAtXFxzbZ5I53nX/j7kwk1r1c8+m8vR3BfB5BKCrfXTw12g0q+K3Pvc052dS/Je3HeZ1Nwwzl8ov/5+WIJ0vcnbaCO7PXXEq/++fnubRs7P8/bNXVv3808kc/dGgpfwL2vPXaDRuJ5kr8sufeJTz080r66uJLG+9ZYxfevlOhnqCzKVbC/7PX1lEStjZF+HU1UVKtmqcpy7OAfDo2ZlVP/9sKk9/NIDPayp/Xe2j0WjczvnpFD88M8NTl+aaeny+WCZTKBEL+wHoiwSYSxcotxBQleXz9tvGyBRKXJxNW/c9dWkegJ+cn11VJVC5LG22j1b+oIO/RqOhEgiTOWdgLZclv/j//ZjvvjDp2J7IFACIRYzg3xsNUCpLFrPFVe/DycsLxCN+XntwCIDnrxi+f7ZQ4rnLCxwY6iJbKPPMpcSKn3shW6BYlmbCV3v+oIO/RqOhUvaYzjmD92K2yGPnZnn8/KxjuxX8lfKPGj9nW7B+Tl5e4NC2Hq4b7sYj4AUz+J+8nKBYltz3mr0IAT9aRWJZLfAa6Ari09U+gA7+Go0Gw8YBSFUF/2TeuD2fLji2Vwf/3ojRKXN2lUnfQqnMi9cWObQtRsjvZc9AlOfNcs+nzJr/n75+kEPbevjxSyv3/VVTN3u1j67z12g0rkfZPqm80/ZRJ4P5KkW/UKP8Wwv+L00lyRfLHNrWA8ANoz2W7fPUpXm294YZ6g7xir39PHVxnmxhZb6/6uvjqPPXwV+j0bidXAPlr25XV/LMZ4zb1cp/teWeJyeMQG8P/uNzGRayBZ6+OM8tO+IAvHLfAPlSmScuNJeYVjhtH1Xto20fjUbjcioJ3+rgbyjsGtsn3UD5r9LzP3l5gZDfw56BLgAOjnQD8IPT00zMZ7h1Zy8At+/pw+sRK7Z+lO3TGwngt6p9tPLXaDQuRwX/dJXtk8w18vyN7T1m8I8EvAR9nqaV/4mJBB/5+nNkzNc7eTnBwZEevKYff3DUuAL4m8cvAljKvyvo4+btsRUnfWdSOXpCPgI+T0X561JPjUbjdlTCt1b5m8E/4wzqiUyBrqDP8s+FEPRFA017/h/73kt8/PvnuO+zx4xSzisLluUDsC0Woifk4/unp/F7heO+V+zt59nxRM2+LsVMKs9AlzHERQV/3dJZo9G4noryrwr+5u1soexIsiYyBcvyUfRGAk2t8i2WyvzgzDR7B6P84Mw0v3j/oyxmixzaFrMeI4Sw1P+NZgWQ4pX7BiiW5bLWz1eenuBPHj7N1565zIWZFP1dhjWlbB+3K39fu3dAo9G0n7zpf6dy9W0fMJK+o7EwYAT/nqrg36zyf2Y8QSJT4CNvu4l0rsRvf/FZAIe6B7hhpJvHz81yq2n5KG7f08u2WIg/fvg0dx0cwmNaRXbKZcmHvnicjO2E9aabRwFsto+7lb8O/hqNZlnbBwzfvxL888TCzvDRGw0wMZ9xbHviwizjcxnuvWXM2vbIqSk8Al69f4C4WSX0+WOXuN5M8ipuMJX/rTvjju1Bn5d/94br+eCDz/DVZy7z1lvHqObaYpZMocS/f9MNvGr/ABdmUhzebjyP1dhNV/u0hhBihxDiu0KI54QQJ4UQv2Vu7xNCfFsIcdr82WtuF0KIPxZCnBFCPCuEuK3VfdBoNK1h2T4Nqn3AWe5Zz/bpi/hrlP+fffcl/t3fPsO1hay17ZHTU9y8PW4F/nfevoMv/stXOqwdgNceHOKewyP89HWDNfv71lvGODwW479/44W6Nf/npowGdTeO9nDDaA933zTKWNw4cemWzgZr4fkXgX8rpbwRuBP4gBDiRuBDwMNSygPAw+ZtgDcCB8x/9wEfW4N90Gg0LWCt8M2XHM3Z7FcCCVvFTyJTIB4OOJ6jNxogkSk4vPRLs2kKJclnfnweMBaLPXNpvm5Ar2aoJ8Sf//LLrJOEHY9H8Dv33MDlRJZP/fBczf2qNfSewWjNfV6PrvaBNQj+UsorUsonzd8XgeeBMeBe4AHzYQ8AbzV/vxf4jDR4FIgLIUZb3Q+NRrN67B0u7T55KlekO2TYO3NVwV81dVOoWv95c/WvlJLxOcMG+qvHLpLJl/jBmWnKEl7TRPBfjlfs6+d1Nwzz5999iWlzBa/i3HSKkN/DcHeo5v8JIfB7ha72WcsnE0LsBm4FHgOGpZRq8sJVYNj8fQy4ZPtv4+a26ue6TwhxTAhxbGpqai13U6PRVJG3BX+7z5/MFS27RNk+2UKJbKFct9oHKqt8p5N5MoUSb755lPl0gS8+Oc4jp6boCfk4sj3GWvDhew6SzBX522Pjju3np1Ps7o/WTQYD+DwerfzX6omEEF3AF4F/LaV0zGCTUkpgRadZKeX9UsqjUsqjg4OtqwSNRtMYZfuAs79PKldkoCtIyO+xmrmpvj71qn2g0krh0pzRj//tt41x8/YYn/rBOR45Nc1PHRi0Omu2yr7BLnb2RTgx4WzzfG46xd46lo/C5xV6he9aPIkQwo8R+P9KSvklc/M1ZeeYP1VD8Algh+2/bze3aTSaNlFooPxTuRLRoJd4OGApenUSiC+j/C+Zw1h29EZ4/6v3cHY6xdWFLK+5bmBN9/2msR6O24J/oVTm4mya3f2Ng7/f69G9fVp9AiGEAD4JPC+l/CPbXV8F3mf+/j7gK7bt7zWrfu4EEjZ7SLNFOD6eYCFbWP6Bmi2BXQUnq2yfaNBHPOK3vPzqds6K6v4+yu/f3hvhnsOjjMYM/30t/H47h7bFuDibthLS43MZimXJnoEllL9H6GqfNXiOVwHvAV4rhHja/HcP8AfA64UQp4HXmbcBHgLOAmeAjwP/ag32QbOBlMqSd/zFj/jsjy+0e1c0a4Td9rGv8k3li3QFffRGAlZb50bBv9cc6GJX/gNdQcIBL36vhw/fcwO//PKd1lqBteLwmJE/OHnFUP9qDvFSto/f63G97dPyIi8p5Q+A+lkVuKvO4yXwgVZfV9M+MoUSuWK5pse7ZuuSL5WNCpiSdIxyTNmU/+nJJNA4+Ad9XrqCPmZTxv2X5tLs6KsE+rcc2cZbjmxb831XK4NPTCR45b6BSpmn2SG0Hj6v0LZPu3dAs/VQnRgzKxyoodm85Itlq55eLfTKFUsUSpKuoI94JGB19mwU/MFQ/6oq6NJshh29kXXf9/6uINtiIU6YMwHOTSfpCfnojdTun0LbPjr4a1aBFfzz7lZOnUShVLaCpfL81ereaMBreP7pPFJK6yRQXe0D0Bcx+vuUypLL8xmH8l9PbhqLceKyYfucm06xZ7ALIx1ZH8P2cffnVwd/zYpRij9TaL6lrmZzUyhVlL8K+qrqJxo0VHSxLEnmiiQyBbpDPmulrJ3eqNHZ80rCSLpuhPIHI/ifm06RzBU5P51m7xLJXlC2j1b+Gs2KsIJ/Xts+nUK+WCbsNwayqISvugLoCvqsVg7z6QILdfr6KJTyvzRrVPrs6NuY4H94LIaU8OSFOSbmM0tW+oCxyMvtyl939dSsGBUctOffOeRLEr/XQ1fQZ7N9Kspfqfz5dKFuUzdFb9RYD6AWeG2U8j80ZiR9v/6sUTW+e5ng7/dqz18Hf82KyWrl33EUSmWCPg+RoNca5Zi0BX/VcXM+k18y+PdFA6TyJV6aSuIRMBqv7a2zHgx1hxjqDvKNk1cBlrd9PHqRl7Z9NCsmrat9Oo580Sj1jAZ8NQnfrmClcmYuXWA+UyDeoJJGrfJ99lKC0VjY6p2/ERwei1mVSMspf93eQQd/zSrQpZ6dR6FUJuDzEA36LLunYvt4rWTwfHo55W9sPz6R2LBKH8Uhc7HXUHeQruDSpoZu76CDv2YVaNun8zCUvxn8q2yfrqDPCvbK869X5gkV5Z/MFdm+QX6/4iZzsddyqh90nT/o4K9ZBZbto4N/x2Cs8PXQFfTWUf4+Aj4jGXx1IUu+WNvOWaH6+8DGJXsVh8020cv5/aDr/EEHf80qqNT5lzC6dWi2OlbCN+CzVvgm80UCPo/l28fCfi7MpKzf69FrD/4bbPuM9IR41+07ePPNy7eQ0HX+utpHswqU4i9LyBXLNbNXNVsPZftUl3ravfPeqJ/z00YJZ/UIR4W9zfNG1fgrhBD8wc/f3NRjjWEu7g7+WvlrVow90VtveLZma1EqS8rSsEIiAaPUU0pp9fJXxMMBLieMxVuNlL/P67EqgTba9lkJRhM7bftoNCvC7vXrip+tjwqCqtqnWJbkimWjl3+govzjET/K5WsU/MFY5RvweRjqDq7rfreCtn207aNZBWlbwE/rpO+WJ2f28vd7BX6vofRTuaLVzllhr+1fKvgr37/R/NzNgG7voIO/ZhVk7cpfB/8tjwqCQZ/Hyt+k8yVSuSKxSMXb77X9vlTw/8WjO8gWN/fnQrd30MFfswrsal97/lufvKX8PVaCN5krkswVGeutVOyohV5CQHeoceh45+07Gt63WfDpRV7a89esnEyhRCTgtX7XbG2U8vd7PUTM4J/OF42Er93zN9V+T8i/qS2dZvB7jPYObi5V1sFfs2Iy+ZK1mEd7/lsfe8K3y6zuSeZKNZ6/mtG7lOWzVfCZaxdKLk766uCvWTGZQiX4a9tn65Oz2T4RU+kns0VreLsiZtb2d0bwN65c3Fzxo4O/ZsXYg79O+G59VHfLgE9YwX4mlaMscSr/SOcof7/HCH1urvjRwV+zYrTt01lYto/XawX7yYUcgGUDQSXhG1tiMPpWQSl/bftoNE0ipTSUvxkIdMJ365O31fmrRP61hSzgVP6xsB8hOkP5K8/fzT39dfDXrIh8qUypLIlHjECgPf+tT96W8A36PPg8gslFQ/nbg7/XI7jn8Civ2jfQlv1cS/we5fm71/bRdf6aFZHNG1+WcMBHxO/Vtk8HYK/zF8JQ/0r5Vw9F+bNfum3D9289UMrfzQu9tPLXrAhl80QCXsIBr7Z9OgB7qScYAb+e8u8k/KbnrxO+Gk2TpPNGu9+w30vI73W0etBsTSoJXyMcRIM+ZlN5wJnw7SR8ZrWPLvXUaJpEKf2Q32u1/9VsbSzbx1T+EZva71Tl79PKXwd/zcpQdf2RgJewX9s+nUDe9L2VFWJX+50a/NV71Z6/RtMkKtiHA4bto4P/1qdgKv+g2c45YuvnY+/t00lUbB+t/FtCCPEpIcSkEOKEbVufEOLbQojT5s9ec7sQQvyxEOKMEOJZIURnlA+4BKX8w6bto1f4bn1Uqaffp5S/EfDDfi/eLd7ArREV20cr/1b5NHB31bYPAQ9LKQ8AD5u3Ad4IHDD/3Qd8bI32QbMB2JW/rvbpDJTyryR8vebPzlT9gDWUXts+LSKlfASYrdp8L/CA+fsDwFtt2z8jDR4F4kKI0bXYD836Y/f8Q/5a5f+FJ8b50ZnpduyaZpXkS2WEwFL5yurp1EofAJ/5Xgva9lkXhqWUV8zfrwLD5u9jwCXb48bNbZotQLra9qlS/v/zWy/yiR+ca8euaVZJvlS2FnhBRfFr5d/ZbEjCVxoTE1b0VxZC3CeEOCaEODY1NbVOe6ZZKfZSz3Ad5T+XznN5PtOOXdOskkJREvRWQoHq79PJwd9q6axLPdeFa8rOMX9OmtsnAPuct+3mNgdSyvullEellEcHBwfXcTc1KyFbKOERxrxXVeqppiFlCyWyhTITOvhvKfKlklXjD5WEb3Vrh05CVfsU9CKvdeGrwPvM398HfMW2/b1m1c+dQMJmD2k2Oel8ibDfixCCsOkNZwuGeppLG6tCF7NFFrKFtu2jZmUUitKqe4fKIq9OVv5+rfzXrNTzb4AfA9cLIcaFEO8H/gB4vRDiNPA68zbAQ8BZ4AzwceBfrcU+aDaGTKFkBf2w32NtA5hLVQL+lfnsxu+cZlUUSmWrrw9UEr0dnfDVnv/adPWUUr67wV131XmsBD6wFq+r2Xgy+RLhgPHFCVcNcZ83lT/A5fkM1490b/wOalZMzkz4KlS1T6cu8IJKS2dd7aPRNEnGtH0A6wogYzZ7m0tXlP+49v1XxeX5DMlccUNfs1AsWzX+4I5qH638dfDXrJC0w/Yxlb/Z438+41T+m50vPjHOw89fa/duOHjX/Y/yp985s6Gvma+yfaJuSPjqxm56mItmZWTzJcvrt4K/ZfsYyn+oO7glgv/H/uklhrqD3HXD8PIP3iAmF7NMmb30N4pCle0z0BWgO+hj72B0Q/djI/Hrls46+GtWRqZQYrA7CFQ8f9Xjfy6VJ+z3smcguiWCfypX5KrcPF/+clmSLZTJFDba9pEO26c75OfJ//B6axVsJ6Lr/LXto1kh6Xyx4vmbP9Uc37l0gd6In7F4mMtboNonmSsyubCxKnspskXj77jRzfJypbKjzh9wrPjtRKz2Dtrz12iaI1soE7ISvrXVPvFIgLHeMFcXsptaVUkpSeWKJM1/mwEV9Dd6QI6R8O3cQF8PIQQ+j9AtnTWaZknni9byf6X8VbCaS+fpjfrZFg9TKktrDuxSfPfFSf7jV04s+7i1Jlcso+xeNay83ai/40Z3Sq2u83cLPq/Q1T4aTbMYi7yqlH++kvCNRwJsi4cBmmrz8LWnL/PAjy9wNbGxAdiu9jdL8Ff22UbbPvmqhK9b8Hs82vbRaJpBJSQbe/550/MPAc2Ve16cTQPw2LmZ9djlhqRswX+z+P7ptto+7gsFPq+2fTSaplAJSaX4/V6B1yNI50uUy5JEpkBvJMBorHnlf2lOBf/qcRDri135X90kyl/ZPRtt++TrJHzdgM/rbuWvSz01TWMf4QhG0ixidvZcyBYoS4hHAkSDPuIRv0P5J3NFSmVJLOy3tmULJa6ZqvvxJYJ/Mlfkbx67SCjgpT9q2EpHtsdaqkZJ5SoBdrPYPurvu+G2j1uVv0ds6qKE9UYHf03TWINcApWGX6GAl2yhZLV26I0YwX1bzFnu+Zt//STpXIkHf+MV1jZ1ZXDdcBenriWZTuYY6ArWvO63Tl7lIw8979j2hd94BUd39636vWxm2ydTMK6kPBtUZ18oSfcmfF28yMt9R1yzapS3r5S/+j2dL1ntnHsjAQDGesOW8r+ayPJPp6Z4dmKesu3Lpvz+d7xsOwA/aaD+L8ykEQJ++KHX8gXz5PH0pfmW3kvKXJg20BXYdLYPGNVIG4WR8HVXqSeohK97lb8O/pqmUcEpYlP+kYAxzUt19Iybyn8sHraU/deeuYyU1Ax6GTeD/z2HRwn5PQ19/0tzaUZ7QozFwxzd3cdwT5CTlxdaei9K+e8d6NpEtk/laiSd35i1B6WypFSWrqz20aWeGk2TpPO1yj9kev6ql79S/tviIWuoy1eemSBqnjDOTCWt/3tpLkPA52FbLMzLdvU2Dv6zabb3Razbh7bFOHk50dJ7SZqe/97BKJMLOWsaWTuxK/+NqvhRyteVto/Ho6t9NJpmsOb3Bpy2T6aO7aNq/b9/apoTEwv82qv3APDSpC34z6bZ0RvG4xG8fE8/L1xdIJGunQB2aTbDjl578O/hpamUZUOtBkv5D0bJl8qOdtTtwh7wW3lvKyGvgr8Llb/fK1xd7eO+I65ZNaoKpcb2KZSYTxfwCOgOGTUEKvh/7J/O4BHwnlfsoi8a4PS1SvC/OJtmh6no79jTh5Twk/NO9Z8rlri2mGVHX9jadmhbD6Wy5IWri6t+L6lckYDPw1jceP3NYP20RfmbuQV32j5a+Ws0TVFd6gnGVUCmYCj/eCRgVaiMmcH/xMQCr9o/wFB3iP1DXU7bZzZtKfpbdsQJeD01i70m5jJISZXyjwG0ZP0kc0W6gj5GYkZ10aYI/vmND/55V9s+7VX+933mGF9+arxtr+++I65ZNUqZhuvYPkZrh0oN/2BX0KogecuRbQBG8J9MIqWxIGwhW7QUfcjv5ZYd8Zp6/0tzRoJ4h83z394bpifkaynpm86XiAa9DHUbq5E3W/DfKNunUDSCnxuVv9/raWud/3dfnOTxc3Nte333HXFN0+SKJT74+ae5ZFbl1FP+EZvyV34/gMcjGI2FCfo83H3TCAD7B7tIZApMJ/PWc9oV/R17+jhxecFR6WI9zmb7CCG4cVtPS8E/mSsSDfgY6lHKv/21/uk22D6uVv5trPPPF8sUStKx3mSjcd8R1zTN2akUX3pqgq8+cxmwKX9/vYRvwVrgpXjtwSHec+cuukPG9gPDXQCcmUzagnol+B/ZEadUljx/pRLUL82lCXg9DJsKXXFoW4wXriysWrmlckWiQR9Bn5e+aGBTKP9svmStgN6oUs98USV83Vfn72tjYzd1fNsZ/PUKX01DVP+b4+OGt57Olwh4PdbwazDsmlyxzGwqx6FtPY7//5/ecshxe/+QCv6L1onEHvwPj8Ws13vZLmP17vhshu1mRZCdQ9t6yBXLnJ1Ocd1w95LvI5Mv4fcKx36nckVi5pXKUHdwUwT/dL5EfzRAIlPYONun5N6Er9/bvvYOKfPKblErf81mJJk1PpjPjs8Dhg8d8js/MqryZ2oxRzzsVP7VjPSE6Ar6TOWfoSfkc/T6Ge4JMtAV5PhERflfrKrxV9w0tnzSt1Aq86ffOc2R3/8Wn/rhOed7yxXpCnrN1w1tCtsnUyjRFzVOSLrOf/0xqn3ao/wzm0D5u++Ia5pGKf/LCWOouDHIxXmxqJK/ZQm90UDNc9gRQrBvMMqZqSSX5tIO1a/uPzzWw4mJSkC/NGesBahm70CUoM/DyYn6vv+z4/P8sz/5AX/4rVMUSmXO2NYXgNHYLWq+l+GezaH8M/mS9TfcqM6eeReXevo9om3tHVRjQR38NZsSe9vj4xPzZAplR6UPYI10BBzVPo3YP9TNmckkF2fT7Kyj6A+PxTg9uUgmX2IxW2A+Xag5SYCh2g6O1k/6ZgslfunjjzGXznP/e17G9cPdNYu4UnnD8wfjimQ6mWt7h8dMoUQ04CXk92xYZ8+8i22fdrZ3UL2l2jlC1H1HXNM0yvYBeHY8QSZfcgR7cC74slf7NGL/UBfXFnJcnKlV/gCHxmKUJTx3ZYFLs2aZZ2/t48Dw/U9eTtS0ZrgwkyaZK/I799zAzx0aoTcSsHoPQWV+b9S0fYZ6QpQlTCfztJN0vkQ44LOa5W0ESvkHXWv7tOeEnzaVvw7+mk2J+mDuH+oygn+h6Aj24Kz8aU75G0nfYlnWtXNU0vfERMIa9GIv87RzaFsPC9ki43POoTEXZlIA7O6PAtAb9TuUf7ZgzO9Vyn+4Z3PU+mcLJcJ+L5GAb8NsH1Xt4kbl72/jIi+l/LOFctuuON13xDVNY9TCezmyPc6z4wlDmfobB/9mlP8BM/gDdRO5o7EQ/dEAxycSddcC2FErfZ+74rR+LswY/08F/3gkwFyqourVSa3LZvtAeyd6SSnNnIqXsNkpdSNwfcK3TYHXfnztg4U2EvcdcU3TJLNFukI+juyIMZ3McX46VeP5h1do++zoi1iBpp7nL4TgprEYJyYSjM9l6DangtWjUjrqTOZemE0RC/uJmf+vLxJgPlOw7CGVZLMnfAEm2xj88yXjaiQc8Jq2z8bW+buxn7/PKyi0qdonZQv+yQ061tXo4K9pSNJcCKWsmLl0oVb5B1Zm+3g9gr0DhiJX/X+qMZK+SU5PLrK9L9JwXGNX0Me2WKg2+M+k2d1fObHEI35KZclC1plkU7ZPf1cQr0e0tdzTvno6bK6a3ghc3dXT0z7ln7Z5/e2q+HHfEdc0TTJXpDvo44bRHnzmIqtGnn/Y761JBjfihtEetveGGz7+prEYpbLksbOzdfMCdvaZ/YLsXJhJs8u0fKByRaKSviqZqmwfr0cw2BVsq+1j75sU0bbPhuDzCsoSx3S5jcKu/BezLgv+Qoi7hRAvCiHOCCE+1K790DQmmTNsn5Dfy/Ujxira6oCtlH91a4el+J17buDTv3p7w/sPbzeuNIplWbciyI5qFqe+wPlimfG5NLtsyr83auzbrOn7K6UVCVbeS7tr/dO2dtntqPZxZcLXfM+FNlT82G09Vyl/IYQX+DPgjcCNwLuFEDe2Y180jUlmi5Y6vnl7HKDW8zdPBvEm/H7FYHeQ/UONWzJsi4Wsla7LKf/9Q11kCiUuJ4yKn4n5DGVJA+VvVPxUJ3wBYpGAZQu1A6X0Qxts+7i5vYO6mm1HrX/akfB1UfAH7gDOSCnPSinzwOeAe9u0L5oGKM8f4GZTjUcaVPsodb0WCCGsPkHLKv9BZ9JXlXk6lL8Z/NW0sVSV5w8Q9nvIbpDarod9PrJqlrcRuDvha4S/9gT/otUqpV39fdoV/MeAS7bb4+Y2CyHEfUKIY0KIY1NTUxu6c2uBlJKrifa3DGgF5flDJfhXK3+f10PA61mR8m8GlWRuxvYBe/A3yjzrBX9l+1jKP2AP/huntuthT/hGNjThKwl4PQ2T6p2MOuG1w/ZJ5UrWLAm3Kf9lkVLeL6U8KqU8Ojg42O7dWTHfOzXFq/7bd5iYzyz/4E2IlNLy/AEOjvTw66/ew88eHKp5bDzit2rl14q337add9+x06oMakR/V5DeiJ+XpirBPxLwMtgVtB7THfLh9QjL9lF11XbPfyOtlnooGyAc8BI2F3ltxFD5QqnsymQvGC2dYe2V/+Rilhv/wzc4VjWS1E46X2So2/iMtiv4t6ul8wSww3Z7u7mtYzhzLUmpLLkyn2lY0riZyRbKlMrSska8HsG/f3P9tMxn3/9y64O8Vuwf6uK/vv1w04+12z47q8pDPR5BPOy3bJ903pjf669qTd1O2ydrm5UQ9nuR0jgG1Vdaa02+WHal5QNGtQ+w5s3dLs6kSedLHJ9IcHR3X93HpHIltsVDBHwe19k+PwEOCCH2CCECwLuAr7ZpXxwkMoU1OROrypFEprDMI9eeR05N8affOd3ScyhrpDu4vD64fqR72Y6e68n+oW5Om+MhL8w6K30U8YjfkfDtqnpfoTbbPpVqH59VTrsR+1MolV2Z7IWK7bPWbZ3V52yp0uFMwejj1BX0ucv2kVIWgd8Evgk8DzwopTzZjn2p5r2fepzf/sKzLT+POvDz6Y0P/n/31AR/8p0zLdkGli8e2vzzfvYPdTGfLjCVNBrG7e6vtYp6IwFHqWc0WJu4LpZl21r82qekKbXfzCrfP/rWi7zzL3686tfNa9tnzRd6zZuC79oSOb+U2TolGvS2rb1D277ZUsqHgIfa9fr1ODuV5JlL82tyJp40V4vOt0H5z2cK5IplEpnCqhOxqqNnV3DtqnjWC5X0/eGZafKlMjvrKP/eaMDqFZS09fJXqKqlbKHUFiWshnuo9g7GtuWDwuPnZ3l6fB4p5aqStvli2ZWre8GW8F1jz18tJryyRPBP50tEAj6iAZ/7FnltRv7hxFUAJuYyS6rm7704yX956Pkln+tqG20f9ZpLffiWo9ICYX0957VABf9/fH4SoIHy99sSvnVsn3W2WnLFEh988GnrBFRNplDC6xH4vWJFts/FmTT5YnnVV5hutn0s5d9CtU+91cHq+9do0aCU0pwn4aU75DLbZ7Py9WevAMaXrnr4h52vPXOFT//ofMP7pZSV4J/e+B7xSnm0Umpa8fw3v/LfFgsRCXh55EWjJLie598bCTCXzlvdM6PBBso/vz62z+lrSb705ATfO1W/bDmdLxHxexFC2GyfpYN/tlDisnmMV9uaIl90se3TovJ/6uIcB//DN7iScFb02T3/eiIyWygjpZHfiQZ9VnvnjcadR70O56dTPHdlgTv2GNn5ibnGJZqTi1nyxXJDfziRKViLZ9ph+yQyxoeplV41yZyx31vB8zfGQ3axmCvi9wpGY7XVVfFIgFyxTKZQMhev1V+stl7Kf2rRsAGnGhyTbKFkXX00a/vYryJWe6wLJenaah+/tzXP//kri+SLZc5NpRzb1Xc+WyizkKkN7CqXEw16iQZ9jqFJG4kO/iZfP26o/l9/9R4AJubrX55Dxc9PN0jU2L+IG237SClJZJb3HJej4vlv/uAPFetnR18Er6c2mPWZK5Dn0gXH/F5FOGB8FdY7+E8u1u8canjARtBXc5KX2xe1oA2WTi4uhbsTvq1V+0wuGn/zmZTz6t4+Ne7KQq2ITNsW9HUHfW2b5uXOo16Hh45f4ZYdcW4363Krp0PZuWYe9EZ9uFVr4GjAu+HVPplCybqMvZpY/QKzZM7Z+XKzo4L/rgYrglXiey6VN6t9qjx/X/NJ1tUwlTSVf4Pgn7ENylE/l7N9zs9UFGcrto9rPX/V2G2Vyl8dy9mq4J/IFKwS6XrWa8pS/qbto4N/+7gwk+Lk5QXedHiUeMRP2O/l8nzjy3N74rAeSoUdGO7ecOVvP9lcbaE/fTJXwOsRVv+RzY4V/Oske8HZ4iGVb5zwzRbbo/yNum8z+DeZ8L0wk6Y75GOgK7DqjqSFkq72We0KX3Usq5X/XDrPwVGjcWG942KtMA94Tc+/1Ja20u486lU8dNyo8nnj4RGEEIz1hhvaPnbl1uhyTamw64a7Njz4q9fzekRryt/s6LlVer5cN2x82fYONgr+hu1zNZF1zO9VVBK+6+z5L6H8K7aPugpZWhFemDXWNAz3hFad3NftHVZf7VNR/s5jOp8uWJ/Hq4na4522KX91hdCOpK87j3oV3zx5lSPbY2w3Z8WOxcMNe/Ionw+WUP4LWXojfoa6Q8yn8xt6VlfBf89AtMVqn9KWsXzAeL/3v+dl/Pxt2+ver1Ygj5tD4bvalPCdTubqfh7s85FDTdo+F2ZS7OqPMNITWvVVnpttn1br/OvZPsVSmcVskcHuIP3RAFfreP7Vyt++bSNx51Gv4qXJJLfu7LVuj/WGG1b72Ef9LRX8h3tCxCN+ynJjZ3Qq2+f6kW4WssVV+4nJXGFLBX+Anzs0UqPoFfGwofxVLidSk/Bd5+Bvev7FsmS2Tvlv1lzuD8ZVW9DnWXJfCqUy43MZdvVHGI6FVj1/2Kj2cWcYsFo6r0L5Symt4D+TrBxPNRMiHvY3vCLLFMxhQgGfVXXWjqSvO4+6jcVsgcVckZFYpSvlWDzMXLpQd3m93cNLLlHtM9wToscMOIkNTPoumMr/BnPy1moTgfaOnp2Az+uhO+Szgn9NwncFq2pXw9RijlHzM1bP+jGUf+XruNwox4m5DKWyZFd/lJGeEDOpPLlV5Ctybq7z96xe+ScyBWv+sV35q0qfeCTAaKz+FZlS+dGAscgLdPBvCyqYj9qC/3ZzelQ99W9P2DVW/jlGekKW2lwv339qMVcTSObNMs/rR4xhKKstAbRP8eoUeiMBm+1THfyNr0J2HZR/Ol8kmSty46hxTOolfTOFkuNqpHqU4/dPT/HDM9PW7Qtmjf9uM/hDpQR5JRgJ362R11lr/C0Mc1Hfu56Qzxn8ze96LOJnOBaqm/BVojIS9Fklx/ZYshGtvEEHf6sW3t6PXrVgHq/j+19byNJv+sf1ztaFUpnpZI7hnqBVXrhe5Z7/5vNP88EHn3ZsS2SMKp19ZuJztbX+9TpfbnV6owHrSqh6kVfA68EjjIU5a830ohEc1HSyeso/ky855iNXzxf4yNef5//+4rNWYLBPLBs2hctqKn5cnfC1unqu/JirE/jBkR5m03lKZh5HXeXHw8aMi9lUvkZQKOUf9lc8f3ss+T8++wS/8+XjK96nleLOo25DBUf7qtAxU/lfrhP8JxdybDcXEtWzhaaTOaSE4ViI2Dor/5OXE46FPmCcaGJhv/V+WrJ9Oi34mzkYqFX+Qoh1m+Y1lTSOwY3blPJ3HpNiqUy+VLaqfMDwg5XtI6VkfC7D+FyGk5cXADg/nSbk9zDUHbSEy2qOtasTvh5V57965X9wtBspK3aPuvLujQQsK7n6iiydLxL2e/F6RMX2sa3yferSPC9cWVjxPq0Udx51GyohMxyrDCMZ6g7h84i6ts+1hSwjPUGigfqtWK/ariTiZnmh+kCsJTPJHHPpApOLzv4hiUyBeNhPOOAlFvavuuInme0szx8qtf5gXHJXs17TvFSg2NkXpSvoqwkG9nbO1r7Y5vgmMgVLGX7DbD54cTbF7v4oQohK8F/hsS6XJcWymxO+qs5/Ncrf+Ftfb+bWlPWjrvLt0+2qT8rpfMm68oxWlXpmCyWmFnMb0hbGnUfdxpVEloGuAEFf5Yvn9QhGYqG65Z6TizmGe0J0NViWrS69h3sqyn89bB81uSpbKDv2I5EpWInm0VhoVbZPuSxJ5UsNK2e2KupkDM75vYr1mualgv9gd5Ch7qBV+aOwgn/Aafukze0qSR3wefiHE0YbkvMzaXaaq5l7wj6CPs+KbR81u1bbPqtT/iG/h119hr2qFnrNpQsIAd0hv6X86wV/day7qmwf9X3diM4A7jzqNq4mMo5KH8VYvLbcM1sokcgUGOoONlyWrUpBh3tChPxegj6PVYGzlpwxZ9aCM4GYyBSsk85ILFS3zng5lAppZorXVqLPpvzrtapeN9tnMYdHQF80wEB3kKlq5Z9vpPyN46AauL391jFemkpx6toiF2fT7DbnGwshzGO9soSvaj7o2hW+ntW3d5hczDHUHaIvWlk5DkYX356QH69HMGxdkTm/g8YgF+O7FfR58HmEZfuogoSNWB/kzqNu40oiy0hPbRdIY5Wv86Cpy/WhnpDRja9O8L+6kMXnEVZS2D4+cLUk0oWapJFS/vb9AkMxKIU70hOqu8JwObbSFK+VEDePSdDnsWq87ayb7ZPM0d8VxOsRSyp/p+df2Rel/N9vNh389I/Oky+WHa2rh3tCK67sUl63W7t6ejwCj1h9tc9gd5D+LuMzpZT/fKby/esJGSM5q7+D9iZ+QgiHkFSCsyxZ99m+rg/+VxeyjjJPxfZ4mGsLWYcqUD6fsn3qK/8sQ91BPGYNcTwccHj+qVyRD37+6Zoe4Evxtj//Yc3wmJemUkTND5A9gVit/KeTOUvhNctW6+jZLKrFQ6P3FfJ516XUc2oxx2CXkVMa6q5dkKVKOkNVto+6Irg0l6Yn5OPAcDe37YzzhSfGAefQGmOVb3PBX+WILOXvq70Kcgs+r8eyv1bC1GKOoe5gpWdUsuL5qxJvlY+ptuNSVfMkDAvZONZ2wTm/zrNAXB38M3mjSVtd26c3TFk6k2gVSyfYcPbmtYWsVXoHEAv7HdU+T16c40tPTfC3x8ab2sdEusDZ6RQ/OD3t2P7SZNKaPaA85XJZspCtfPhGV1kCaCn/jgv+xhe1US4jFPCSWYdST6USwfD9U/mSQzioPEOkyvZRJ4XxuYzVeuSNN41aQXunrYOpYfvUHx5i5+lL89z+kYc5MZGwhI1blT+A3yNWpfwnzWMa8BmLB1V/n/lMgZjNXhyuc1K293ECFfyNGGHvJrzUQKm1wNXBX6nvesp/LG58sexnYiuZ272E7ZPIMtxtC/5Vts95szTz4Rcmm9rH05OLAJydTjFt2gWpXJGJ+Qy37uwl4PNYnv9itoiUWAlf5TmuOvh3mO2jgr/9i2cn7PesW8JXBf8h86e91t/q717H9pFScmk2bS08vPumEcAI2NviFbtyuCfU1DjHY+dnmU7m+Hd/+4yV23FrwhcM5b/Sap9csZL7A+iPBizbJ5HOW+ILzJNyolb52xf02YXkxFyGoHk85rTyXz+ssswGyh+cq3yvLWYJeD3EI37D9qlT5z+5kHM8X7xK+V80F+c8c2m+YYdHO6euVbz9Jy7MAXDWnBy0f6iLoe6gZSOo11GLy1St/0orfjrW9okubfusNuFbLktHDsaOlJKppFP5gzNJX8/zDwd8SGm0Xxify7DDVPk7+iIc2tbDzqqhNaqs8Nri0sf63HQKr0fwwtVFPvrtU4B7E75gnEQLdRKrUsqGJwV79RYYifxZm+ffG3EG/2sLWUfyNp1zKn+7kJyYz3CDuRJ8vdvCuPeoU3+Bl0JdDdiV/9SC8SUWQhAJ1Hr+qVyRxVzRUtxQa/ucn0lbVR3ffXF59X/q2iJhv5eAz8Ox87MAvGRW+qjgrxKIKrdgef6tKv9OC/7L2D6rTfj+5Y/O8/qP/hPnplM19yUyBQolWfH8e2qVv/L2HSt8zXYT43NpMoWSpfwB/uc7j/A/fuGI43VGzHUqy9X6n5tOcXgsxttuHeObJ68BuLbOH4y2zvWC/Ee/fYqf+cPv1bXRaoN/kNmUUZ2TqLJ9RnpCFMvS0fO/2vPvDhnBv1AqcyWR4aYxI/hr5b+OKC/O3tpBEfJ7GewO1ij/YfPL2xX0UihJRzK1UuNfWTAWj/hJ50tW060LMyletb+f4Z4g323C+jk9uciB4S6ObI/xk/OG8j8zmcTrEezuj5oJROPDWFH+ZrVB2EfY71258u/Q4B/yewn5PY0TvnXq/B88dokP/NWTDZ+zUCrzye+fRUp47OxMzf3VgWLItATtSfqK8q/sl/pdXfkpzx+MlgK32brQQvMW3/npFHsHovzHf3YjA+YJyd22T63nP5fK84kfnGN8LlNTmQWVqzZ1LAe6DNtH2a5226f6uJTKkmzBuZo7agpJNWvixtGYsR9a+a8fVxIZY3JXAw94WzzMuG2oy7WFnHUwK324K+q/3slEqYBEpkC5LI367P4orz04zCOnppatxDl9LcmBoW6O7u7j5OUEmXyJM5NJdvVFCPg8DPUErQ+j8ntjtmqD0Tqe43Io26fTFnkB7BvssiyUakJ+b80kr0dfmuHrx680bOL30PErXE5kEQKOmbacnergHw/78XmEw/ZJ16vzNz+Tp64ZOZ8dfbVXp3ZUIFqqtDeTL3E5kWX3QJR4JMB/ffthfOaCRrfi93pqFnl9+kfnrWNyfrp2qFM922culbeUun0xoXIQlABTJ/qow/M3lL9yGXb2RegJ+Uho5b9+XE1k66p+xZHtMY6dn2POvGRTffqBug2Z1Nl9qMr2AaPV8uRijmyhzK6BKHcdHCKVL/H4udmGr59IG//nuuEuju7qpVCSPDM+z5mpJPvMsYVD3UESmYK1AA1qlcdKe74kc0WCPk9HKsIv/MYr+bc/d13d+8J+42rOXt67kDX+pi9N1Xr6Uko+/v2z7B00jucT9YJ/0hkoPB7BYHfQafsUlO1T+XurE8HpOsq/HgGfxxwe0vhYX5g1bKk95uKw1984zInfe4M1dcqN+DzC0dgtmSvy6R+dtzqwnpuuPe6TizmEwFrL0xcNUDSFHTQK/kZgT+dUR09ntU8qV7QW8431humNBrTyX0+uJOrX+Ct+6eU7yRXLfOGJcTL5kjWhByqWiD3pq77QQ3bbx9biwerE2Bfhlfv7Cfg8fGcJ6+eUWelzYLiLl+0yLvMfPTvD+emUNbNWKb6pxZwV/HvCzg/flfnMitrEdmJTN0U44G3ocVujHG2+/0LGOL72xLvix2dnODGxwL/4qb3cvruPc7aKLEW1SlS/OxK+ZqMv+8jMiE3595oFBssxXKem3M65KWfwB2eewY34vB5HY7e/eewiiUyB37/3EH6v4FwD5d8fDVgLBdVCLyUQYuGK5z/YHSQS8FpFGql8bXK/K+SjLI21O2B8Z+ORgPb8WyFbKPFPp6asM2o1VxNZRuONL6cPjvRwdFcvf/XYBUtRLWX7zKTyBLweR1sEq7lbumB14NzdHyUS8PHKff08/MK1hoFZqb4DQ93EIwGuG+7iy09NUCxL9g0awX+wp1I9ksgUCPk9ji/0dSPdXE5kefOf/IAHf3KpqUVMnTbIpVnUIit70lcp/9Om/WLnE98/R380wNtuHePobuPkXK3+pxZzBH3Oz8RQHeVfXX6q9uXcdGpZ1a+oV1Zo55wpPnbbgr/b8XuFlfDNFkp8/PtneeW+fo7u7mNHX4TzdZL4U4tZBm3l3H1R4zuoArxd+Qsh2DfYZZ0YVLyIVNk+AC9eXWCoO0jI7yUebr0zwHJ0dPBP5oq871OP8/Dz12ruyxZKzKTyjC5h+wC85xW7OD+T5ktPGouy7Alf4zUqgWI2macvGnAoOHtb5/MzKXwewba48Zp3HRziwkyas3U+YGCovkjAa80XeNmuPusEopS/qiKZWswyn85br6f4tVft4T+/9SaKJclvf/FZfu6jjyx7AujEQS7NUBnibrN9zKupU1XB//S1Rb7zwiTvfcVuQn4vN43FHBVZClXjb/9MDHaHmLIlfNNVvfyhogyLZemo9FmK4Z4Q43PphvmJc1MpBruDrjy2jTBsH0N8PXjsEpOLOT7ws/sB2NMf5fxMveCfc1zJKfvnrGkRxau+g3sHo9aJQeUS7J6/iiUvXl20Ssx7I/516QZsp6ODf18kgN8ruFannl5VyCyX7Lr7phH6owE+9YNzwPLKXzV6UsTNS8D5TIEL5mIddbn4musGARr6/qcnFzkw1GW1irh9d6XCQw1rsZcOGu2cna8f8Hn4lTt38Y1//VP8/r2HuDib5vlleoUnc8WOTPYuR70h7moma7Xt89VnLuMR8Ct37gQg6PNy81isJulrr/FXDHYHmUnlHYqzWvlH/JW/f6MEdTVvOjxKMlfkN/73E3VHOp6fSbGnX6t+O4btU+bERIL/8tDz3Lm3j1fu6weMK6TzM6maBmuTZmsHhfrOqwBfLcD2DXYxMZ8hky/ZpnjZPX/j8ZcTWUvoxSMB5lNa+a8aj0cw1F2/4VVlde/Sqiro8/LO23dYXp066OrMnawK/sr/U3SHfAhhrPy7MJNil+3LNxYP4xFwpU7raDACzgFbMu7oLqOdw0hPiO6Q8YHpjwbxiIrtU/3BUwghuOuGYQBOTCSWfM/JXLHjOno2QzjgHOVYLBntssN+LxPzGceJ/vFzsxzaFqO/qxIEXra7lxMTCceVlb2vj2KoO4iUlWZg9ha/ilCg8tVsVvm/+sAAf/DzN/P909N88PPPWNOlFOemUw6/X2PYPlOLOe77zDF6IwH+5N23WVdpuweiZAtlx8K5clkyXXVCV8H/SiJLd8hX0zRwrynUzk4n6yp/e4fZivIPsGjW/q8XHR38wbBp6lVAWGWZTZS5/dIdOxHCUNEquHbVUf6zqZx1CajweAQ9Ib+h/GfSjk6MPq+Hga76+zefzjO1mOOAae+AUe431B3kwHBlm9cjGOgKMrmQYz5dcCR7q9kWC9EfDfDs+PLB35Wev8+p/NWJ/cgOo+76tLmKN1cs8dSleW7f3ef4/0d39VEoScfft9oigIqAUFefmXzJUeYJTk+42eAP8M6jO/jde27g68ev8P985YS1fSFbYDqZ135/FT6Ph5emUsym83z8vUcdx0pdJdkX782bi/bsyj/k91pNFu1+v0Ll585OpWyev7PaR6HyO/Zc4XrRUvAXQvyCEOKkEKIshDhadd+HhRBnhBAvCiHeYNt+t7ntjBDiQ628fjM0qoC4skRrh2p29EW46+Awu/sjlipQl20O2yeZt5I/duIRP+emUyxmiw7lr16/Xh92FWjsZXhCCP743bfy4Tfe4HisUeufZcHWTrYeQghuGotxfDnl71LPvzrhqyp91BWX8v2PjyfIF8tWYz2Fqsg6dsGw8QqlMrPpfF3bByrjHTOFWuVvPxnsaDLhq/gXr9nL+1+9h79+7KK1zypxqZW/E1X59T/ecYSbxmKO+3YPGH93e61/veotgD7zir/adgXjby6EUQ2UrlftYw/+lu2jcoXr5/u3qvxPAG8HHrFvFELcCLwLOATcDfy5EMIrhPACfwa8EbgReLf52HXDCP61wfWqeYnWbJD76C8e4bPvf7l1O+jz4vcKK+GbLZRI50s1tg8YHqBSg7v7nV/kRn3Y1ZfWrvIB7tzbb82CVQx1h5g0R781sn0UN2+PcXoyuWTSt5NLPZeikvA1g79Z6XPTWA9Bn8eq+HncTOraczBgXP7vHYzyhLkSezaVR8raQKHWgagywuouj2Bc0al1FmMrUP6K3/jpfXg9gi8/NWG+lg7+9fjVV+3mD3/hCP/syLaa+7bFwgR8HkfSV63MHup2ikYl+uqJr5DfKNo4O5WySsOrWzor7LYPrO8q35aCv5TyeSnli3Xuuhf4nJQyJ6U8B5wB7jD/nZFSnpVS5oHPmY9dN0ZiIZK5Yk0HziuJzJI1/tV0h/yOnj1gHECVwFH+bbXtA87+Pruqgn+jPuynryWJ2ip9lmKoO8jl+QzpfKmm0qCam8ZilMqS5xokffPFMrli2dXBv6L8K43y9g12WUnfn5ybZd9g1OH3K47u6uWJi3N86clxPvjg00BtoBjtCXHz9hh/+M0XeebSfN1qHzDU4UBXwGEBNctgd5CfOjDAV56aoFyWnJtOIUTt58/tvGr/AO942fa693k8gt39EYftoxyD6hO6+t43El97zXLPTL6ER2B17gTniUB9363gn9q8yr8RY8Al2+1xc1uj7TUIIe4TQhwTQhybmppa9Y6o0sxq6+dKIsvIMsne5YgGKt34ZszFPdXVPlDpsilE7UrNkVjIWqFr5/TkIvuHuhwlgo0Y6g5aCiG2hO0DcNi8tG2U9FU2lhs9/3C17WMq/56Qn+uGuzh9bZFSWXLswlyN5aM4uruP+XSBDz74DOen0/yrn9nHq/cPOB7j8Qg+8b6jDHQH+LVP/4S5dL5um+mI38vYCi0fO2+7dYzLiSyPnjMWBm6LhV2/qGul7O6POoL/916cZKAr4JilAJXvfSPbdZ9Z7rmYNUY4Vi/oE8Io71QnAsvzX8dB7st+w4UQ/wiM1Lnrd6WUX1n7XTKQUt4P3A9w9OjRVQ+ztBorJbJW4gWM4H9wpLVl7fZpXpbyr6MGY2HjzzxqzvWtt39XzZ4rilPXkvy0WQq6HHYVspztMxoLMdDVOOnbqU3dmkEdG9VlU3n+PWFjitbfPX2ZJy7MsZgt1iR7FW++eZSZZJ479vRy287ehifvoe4QD/zqHbzjL35sVPvUCcojsVBLrRd+7sYRuoI+vvzkhK70WSV7BqJ878UpsyFbie+8MMkvvGyHo502VJR/Pc8fjKRvplDi7HSqJr8jhKAr4HPYe5WE7/op/2W/4VLK163ieSeAHbbb281tLLF9XRiu0+c8nS8ytZirOXuvFPsQBjXGrZ7toz4Q1cleqDSBu7pQCf6L2QJTiznHyWop7KsNlwv+KunbSPm7O/gbF8I5s9mepfzDfisI//VjFwAaKv9IwMe//Jl9Tb3e3sEuPvm+o/zyJx6rsRQB/vKf34Hft/opW+GAl7tvGuEfTlxFAPfeWutra5Zm90CUfKnM5fkMT16cI1so8+abR2set5zyV+Wex8fnLSfATnfI57B4u4I+fB6xrp7/en3Dvwr8tRDij4BtwAHgcUAAB4QQezCC/ruAX1qnfQBswdXW7VBl71ste4sGfSxmlfI3bZ86CV/1gVDVA479i9XaUmoV7546j6+Ho5dQnQ9WNTePxXjk1JRRYlilQtQlbrWn6QYCXg8eYVf+BYSAroCP68zE+0PHr7ItFmq65cJy3Lqzl0d/5y666vj6y1l4zfD2W8esmb97BpoTE5oKak7y+ZkUf//sFYZ7gnWv+irBv/73b78p5ObSBccENsXv3XuTIwcphDAWeq2j8m+11PNtQohx4BXA14UQ3wSQUp4EHgSeA74BfEBKWZJSFoHfBL4JPA88aD523YgGfXQHfY7gqrL3u1tc7Vht+1T39VGo2vudfbWvZ7d9FCoAN3tyGlqB7QNG0rcsqZv0/dozlxnoCnJrVb94NyCEcEzzWsgW6Qn58XgEO3ojhPwe8qUytzdQ/atFvcZ6cOfefiuoNCsmNBWUVXZ8IsE/vTjFPYdH6x6rfqvUs/73z95WI1rnRP/6G4drSk3jkfXt79Nqtc+XpZTbpZRBKeWwlPINtvs+IqXcJ6W8Xkr5D7btD0kprzPv+0grr98sQz1BR/BfaXBtRNQW/Ov19VGoD0R1mScYVUTRgNdR8VPp/tnc/tlV+nLVPgCHtxsfsuPj847tC9kCD78wyZtvHq3xNN2CfZrXQqZAj5mv8XiE1U+pkd+/GfF4BPfeYtRU7NXKf8UM9wQJ+7088KPz5Etl3nxzfevsyPY4r7lukCM74nXvNxq8Gd9ne2uHpeiN+Ne1s2fHr/AFtZDKpvynUwytQYOraMBreeSzdfr6KG7ZGeenrxts6BOrOZ+Kc9NpRnpCDYfMVBP0eW3Tu5YP/iM9IQa6ghyfcCr/b528Rr5Y5i23uNcbtk/zWsgW6AlV/p7XDRm+f6PjuFn5wM/u40/efate3bsKhBDs6o9wbSHHWDzMbTvjdR/X3xXkM792x5J26V7T+qmn/Oth2D6bVPlvFYZtow7BsH3W4osQDfpI5UtIKZmu09dHMdQd4oFfu6NuJRDUtuI9P5NacT32UHeQ7pCvKcUuhODwWE9N0vcrT0+woy/MrQ3Uixtw2D6ZoiP4v+GmEX72+kHLv90qdIf8dRcxaZpDJWvfdPNoU6XXjVDKv1lRp5X/GjBsKmvVne/cdHpNuhtGgz5KZUmuWK7b16fp/etxBv8LMysvyxvqDjXl9ysOb49zenLRWkQynczxo5dmeMuRbS19wLc64YDXWnOxkK3YPgBvODTCX/7qHevmz2s2Jyo3WK/KZyVUlH9zwV8r/zVgpCdEsSzNIcsFppO5NVH+9uZusw36+jS7f5OLOcplae5fvm5Z6FK85cg23n5b/ZWK9Xj9DcP4PB7e95ePk0gXeOj4FUplyVuO1F1z5xpCvirPP9R6xY1ma/OLt+/g37/pBmuB5GpRpduRJu3meMRPrli2qs/WGlcUc9tX+aqhWWtR+aBW482k8qQa9PVphpGYcXKaTuUse2ql+/fO23cs/yAbh7fH+Niv3Ma//N9P8p5PPUZZSq4f7ub6Fhe+bXVCAa/VimMhW2wqh6LpbHb1R/n1n9q7Bs8TIR7xN92ltdLfJ0840Fo3gnq4QvlbC70Wsms6yk5N4FFjIluxfQCuJXJrVonUDHfdMMzHfuU2XriyyImJBVcnehVhv4dsvmT18tfKX7NWhPxeHvntn+Vdt+9s6vH2+d/rgauC/9WFrNXattkyyqVQyl8tympU7bMcI+u0f81w1w3D/MV7buPorl5+fgW2UaeiEr6qisvu+Ws0rdIT8jddRq0WjK3XQi9XfLKNGapwbSHH+Gya0VjzZZRLoYL/RaX8G1TzLIeaKXB1Icv5mZWVea4Frz04zGsPDm/Y621mVJ2/sn608te0i96o8dlbrxYPrlD+fnNi1rWEYfu0urJXoRK+rdo+A11BvB7BtUR2VWWemrUj5DeqfSpN3XTw17QHu+e/Hrgi+IOR9L22aNgqa+WnW7aPGfzr9fVpBq9HMGiOczyvuy+2FSv4W+2cXXFxrNmEqNLtxDq1dXZN8B/pCXHq6iJz6QJ71yr4ByoJ30Z9fZplOBbizGSSmdTKyzw1a0fY76VQklaLbq38Ne0i5PcS9nvXbaCLa4L/UE+Iy+ZCqrVW/rliuWFfn2YZ6Qlas3V1A672ofrqT5rtNnTw17QTY5WvVv4tMWLrl75WwdXv9VhzVldb6aMY6QlRMlcg6x4s7UMNcZ80B3Vr20fTTuKRwLoNcXdN8FcLvTwCdrQ4xMWOSvqudoGXYtjWy3ujyjw1tSjlf20hi0c034RLo1kP4uuo/F3zyVa1/tviYYK+tSujjAa9zKZWX+mjUFcmG13mqXFSsX1ydK9jn32Nphnec+cu8qXyujy3a4L/iDXQYm1VtVKGq+3ro1D7p8s820s4YFwMX1vM6gVemrbzxsOtNZNbCvfYPuac27Wq8Vesle2jlL8u82wvIV9F+esFXppOxjXBPx7x886j21tuy1qNqvhp1fYZjYWJBrwc2tazFrulWSUq4av7+mg6Hddc1woh+O/vOLLmz6uUf6vVPuGAl+/+Xz9DXxMD2DXrh/L8Qff10XQ2+tPdIhFTKa62r4+doe7Q8g/SrCuO4K+Vv6aDcY3ts16sle2j2RzYK630Ai9NJ6ODf4tYtk+LCV/N5iCklb/GJWjbp0XecGiEfKncUl8fzeYh5K/oIe35azoZ/elukcPbYxze3tpsT83mIeD14BFQllr5azobbftoNDaEEFbSV3v+mk5GB3+NpgqV9NVN3TSdjA7+Gk0VIa38NS5AB3+Npgpt+2jcgA7+Gk0V2vbRuAEd/DWaKkI+r+7lr+l4Wgr+Qoj/IYR4QQjxrBDiy0KIuO2+DwshzgghXhRCvMG2/W5z2xkhxIdaeX2NZj0IBby6l7+m42lV+X8buElKeTNwCvgwgBDiRuBdwCHgbuDPhRBeIYQX+DPgjcCNwLvNx2o0m4aw30O3tnw0HU5Ln3Ap5bdsNx8F3mH+fi/wOSllDjgnhDgD3GHed0ZKeRZACPE587HPtbIfGs1a8t5X7OaaOcBdo+lU1lLe/BrwefP3MYyTgWLc3AZwqWr7y+s9mRDiPuA+gJ07d67hbmo0S/Oq/QPt3gWNZt1ZNvgLIf4RGKlz1+9KKb9iPuZ3gSLwV2u1Y1LK+4H7AY4ePSrX6nk1Go1G00Twl1K+bqn7hRD/HHgzcJeUUgXpCWCH7WHbzW0ssV2j0Wg0G0Sr1T53A78NvEVKmbbd9VXgXUKIoBBiD3AAeBz4CXBACLFHCBHASAp/tZV90Gg0Gs3KadXz/1MgCHxbCAHwqJTyN6SUJ4UQD2IkcovAB6SUJQAhxG8C3wS8wKeklCdb3AeNRqPRrBBRcWo2L0ePHpXHjh1r925oNBrNlkII8YSU8mi9+/QKX41Go3EhOvhrNBqNC9HBX6PRaFzIlvD8hRBTwIUWnmIAmF6j3dkquPE9gzvftxvfM7jzfa/0Pe+SUg7Wu2NLBP9WEUIca5T06FTc+J7Bne/bje8Z3Pm+1/I9a9tHo9FoXIgO/hqNRuNC3BL872/3DrQBN75ncOf7duN7Bne+7zV7z67w/DUajUbjxC3KX6PRaDQ2dPDXaDQaF9LRwd8t84KFEDuEEN8VQjwnhDgphPgtc3ufEOLbQojT5s/edu/rWmOOB31KCPH35u09QojHzGP+ebN7bEchhIgLIb5gzs9+Xgjxik4/1kKIf2N+tk8IIf5GCBHqxGMthPiUEGJSCHHCtq3usRUGf2y+/2eFELet5LU6Nvi7bF5wEfi3UsobgTuBD5jv9UPAw1LKA8DD5u1O47eA5223/xvwUSnlfmAOeH9b9mp9+V/AN6SUB4EjGO+/Y4+1EGIM+D+Bo1LKmzA6Ar+LzjzWn8aYe26n0bF9I0a7/AMYUw8/tpIX6tjgjzEz+IyU8qyUMg+oecEdh5TyipTySfP3RYxgMIbxfh8wH/YA8Na27OA6IYTYDrwJ+IR5WwCvBb5gPqQT33MMeA3wSQApZV5KOU+HH2uM9vNhIYQPiABX6MBjLaV8BJit2tzo2N4LfEYaPArEhRCjzb5WJwf/MWrnBY81eGzHIITYDdwKPAYMSymvmHddBYbbtV/rxP+LMUyobN7uB+allEXzdice8z3AFPCXpt31CSFElA4+1lLKCeAPgYsYQT8BPEHnH2tFo2PbUozr5ODvOoQQXcAXgX8tpVyw32eO2OyYul4hxJuBSSnlE+3elw3GB9wGfExKeSuQosri6cBj3YuhcvcA24AotdaIK1jLY9vJwX+pOcIdhxDCjxH4/0pK+SVz8zV1GWj+nGzX/q0DrwLeIoQ4j2HpvRbDC4+b1gB05jEfB8allI+Zt7+AcTLo5GP9OuCclHJKSlkAvoRx/Dv9WCsaHduWYlwnB3/XzAs2ve5PAs9LKf/IdtdXgfeZv78P+MpG79t6IaX8sJRyu5RyN8ax/Y6U8peB7wLvMB/WUe8ZQEp5FbgkhLje3HQXxrjUjj3WGHbPnUKIiPlZV++5o4+1jUbH9qvAe82qnzuBhM0eWh4pZcf+A+4BTgEvAb/b7v1Zx/f5aoxLwWeBp81/92B44A8Dp4F/BPrava/r9P5/Bvh78/e9wOPAGeBvgWC7928d3u8twDHzeP8d0Nvpxxr4PeAF4ATwWYzZ4R13rIG/wchrFDCu8t7f6NgCAqOi8SXgOEY1VNOvpds7aDQajQvpZNtHo9FoNA3QwV+j0WhciA7+Go1G40J08NdoNBoXooO/RqPRuBAd/DUajcaF6OCv0Wg0LuT/Bzgk9Lqv2jlMAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"Y -= np.mean(Y)\n", | |
"\n", | |
"plt.plot(Y)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"We then model the task in the same way, but we do not add a constant term regressor.\n", | |
"\n", | |
"This leads to the following design matrix:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 117, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"<matplotlib.colorbar.Colorbar at 0x1f5247b0eb8>" | |
] | |
}, | |
"execution_count": 117, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAH8AAAD7CAYAAABDnEvdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQ5klEQVR4nO2dfZBk1VnGf8/M7OyyQ9glWcAVSIBiWaUwFkgpFFUGgejmo8AqI4YYJBHlnxhjEqNELaE0fxA/grGMKBICScUgYkq3KhG0CB+lRZAlUCQstQE3QHYD4Wt3Zdmvme7HP+7tYVimu29vn9N9557zq7q10z2333t3nn7POfec9z5XtsmkycS4TyAzPrL4CZPFT5gsfsJk8RMmi58wQ4kvaYOkLZKekHRlqJPKjAYd6nW+pEngu8DbgW3AA8AltjeHO71MTKaG+OxPA0/Y3gog6RbgIqCr+NOTK33YslVDHHJ07J3dxYHWHg0T4xd+bsYvvtSqtO+Dj+y/w/aGYY43KMOIfyzw/QWvtwE/c/BOkq4ArgBYMXUEZ59w2RCHHB33PXnz0DFefKnF/9zx5kr7Tq59fM3QBxyQYcSvhO3rgesBVq1Ym9RcsoE27XGfRleGEX87cPyC18eV72VKjJl1tWZ/HAwj/gPAOkknUoj+XuB9Qc6qQTQy823PSfot4A5gErjR9qPBzqwBGNOq8arpUH2+7a8DXw90Lo2kTUPFz/TGQCuLny458xPFwGxT+/xMb4xzs58shlZ9tc/ix6SY4asvWfyoiBZDrQ1FJYsfkWLAl8VPkuI6P4ufLO2c+WmSMz9hjGjVuEY2ix+Z3OwnihEHPDnu0+hKfdukBlBM8kxU2voh6UZJz0n6TpffS9Jfl2X0j0g6o1/MLH5kWuVET7+tAjcBvap73wGsK7crgOv6BczNfkRs0XKY/LJ9r6QTeuxyEfBFFzdifFPSaklrbT/T7QNZ/Mi0R3ept1gp/bFAFn8cFAO+yn/iNZI2LXh9fVn2Ho0sfkQ6A76KvGD7zCEON3ApfR7wRaZlVdoCsBH4tXLUfxawq1d/DznzoxJyhk/SV4BzKbqHbcBVwDIA239HUUX9TuAJYA/wwX4xs/iRaYcb7V/S5/cGPjRIzCx+RIqFnfr2rFn8iBgxW+Pp3Sx+RGyCTfLEIIsfFY1ykmdgsvgRMTnzkyYP+BLFKBdzpEpRul3fP3F9z6wR5Js2ksWEm+GLQRY/MjnzE8VWrTO/75lJOl7SXZI2S3pU0kfK998o6T8lPV7+e2T8011aFAO+yUrbOKjytZwDPm77VOAs4EOSTgWuBO60vQ64s3ydeQ1FDV+VbRz0PartZ2x/q/z5ZeAxitqwi4COR+nNwC9GOsclSzHgU6VtHAzU55fVo6cD9wPHLKgUeRY4pstnXuO9mxqNmOGTdDjwL8Dv2P4/6dVvq21LWtSA5DXeu4etNYrwLa+p6VEjZvgkLaMQ/su2v1q+/cNOXbiktcBzsU5yKTNAAefI6Su+ihT/PPCY7c8s+NVG4DLgmvLff+sXyxOivXL5IZ5qdyZe2Rc8ZghsmG0vYfGBc4BLgW9Lerh87w8oRL9V0uXAU8DFUc5wCVM0+0tYfNv/BV2nqc4f5GBeNsHe42YG+UglZrbUM/Mhz/AlS+dSr66MVPy5FeKl9eEPObMleMhALPFmPzMcuYavpL3C7F4/O8pDjpVitF/f0u36tkkNoDPJE2p6t99DLCW9uVyEe6h053hnr3hZ/Mi0y/Ltfls/yodYfo7CgeNU4JJygW0hfwTcavt0imce/W2vmCNt9pcvn+WUk3reOHqokSPEHJ7Ao/0qD7E00FlAWQX8oFfAPOCLzACj/X7mDFUeYnk18B+SPgzMABf0OuBIxT9iah/nHR3+uuwe3ho8ZghsMVdd/GHNGQAuAW6y/ZeSzga+JOk024s6v+fMj0zAZr+K88bllI5dtu+TtAJYQ5dFt5GKPzO5n7NnHg8et7aZT1DxqzzE8mmKKfebJP04sAJ4vlvAnPmRCSV+t4dYSvoTYJPtjcDHgX+Q9FGK794HStOGRRntaJ8WJ03tHuUhx0roYo7FHmJp+48X/LyZYhW2EjnzI5OndzsH0wRHTdbzmjwGNswt8WKOzBDkJd1EaUQBZygmEMu1bJSHHDvO4qdLHvAlip37/IQRrTzaLzCmtfgaQ2PJfX6i5OrdBRjY77lRHnK8uLa3EQI586OTR/slbcweJ1S9mwd8aZOb/ZK2zZ46/zUikEf7iWJn8ZMmX+qVtBH7avzHiEGde7mc+RExop1H+wVtYF+NnzkTgxonfs78qDRlwFfeKLgJ2G773WX9+C3Am4AHgUttH+gVwyi5zK9z6g/SIX2Ewn2zw6eBa22fDOyguFskcxC2Km3joJL4ko4D3gXcUL4WcB5wW7lLJfvVQYyIB9nqioF2W5W2cVA18/8K+D2KMRsUTf1Oe36JrvOs9tch6QpJmyRt2vliWmv5GLCqbWOgiuX6u4HnbD94KAewfb3tM22fueqNk8x6KvhWZ+xqWxX6OXOU+1y8wB7/H3vFq2rCeGFp8bGC4ub/zwKrJU2V2d/3We3JEmjAt8CZ4+0ULe0DkjaWt2h19lkHfBI4x/YOSUf3ilnFcv2Tto+zfQLFnaHfsP2rwF3Ae8rdKtmvpke1wV7FAd+8M0d5VdVx5ljIbwKfs70DwHZPP+Rh2szfB26R9CngIQp/3p4cYIrtc4k9kKN65odw5jgFQNJ/U9zJe7Xt27sdcCDxbd8N3F3+vJXi25jphsHVR/IhnDmmgHXAuRRd8b2SfsL2zm47j4zdreXcu+uUUR6yBozUmWMbcL/tWeB7kr5L8WV4YLGA9V11aAquuPVn3plD0jTF+GvjQfv8K0XWI2kNRTewtVvAkWb+y/tXcM+TJwePexJ7gscMRqDRfkVnjjuAn5e0GWgBn7D9YreY9b5IXup0JnlChevvzGHgY+XWl9GKv28Cbzk8QuD6Zn4u5kiZMc3bV2Gk4k/uh9XhndhqzeLPHKsHOfNjUn0kPxZGnPltjti6d5SHHDPjW7GrQs782OTMT5galzCMVHwdaDG9fUfwuF5W0+9w4Ov80NT0r9Yc8mi/Q6uFd+wMH/foNeFjhqLG4ueFnYQZrS1Lu0V79yvB407UOPNzs58qJk/vzmPwXEKGTFDrPj9nfmRys58yWfyEyeKniZyb/dei+o5+o5BH++mSM7+DhKbSetJG7vNTJff5ryIJTefMrws58yOjXMxRIqHp6ZEeMtOdvKQbm3D36lVy5ij3+yVJltTzrt+RZz5TCfU0AQd8VZw5yv3eQOGcdn+/mDnzYxMu86s4cwD8KYVN3r5+AUec+aCp8NZpNR5Qj9SZQ9IZwPG2vybpE/0OmFAbPHrEQKP9oZw5JE0AnwE+UPUzWfyYhJ3k6efM8QbgNODuwiOTHwE2SrrQ9sIWZZ5K4ktaTeG+eRpFQ/brwBbgn4ATgCeBizsuUD0CpTXgg5B90rwzB4Xo7wXeN38YexcwX8wo6W7gd7sJD9UHfJ8Fbrf9Y8BPUnjwXgncaXsdcGf5OnMwgQZ8pd9hx5njMeDWjjOHpAsP5dT6pqGkVcDPUvYl5UjzgKSLKP1fKLx376awZ+sVDSbTusAIObffz5njoPfP7RevihInAs8DX5D0kKQbJM0Ax9h+ptznWeCYxT680Hv3QKu+DhrRCDjJE5oq4k8BZwDX2T4deIWDmvjSC2bR/8JC793pyZVFv78ktmH/tMVfRO1q2zioIv42YJvtzozRbRRfhh9KWgtQ/tvT6jNZapz5fft8289K+r6k9ba3AOcDm8vtMuAaqnrvChyhz2+vXB48pifClF81YT3/w8CXS/O/rcAHKVqNWyVdDjwFXBznFJc4S1182w8Di80+nR/0bA6RvcfNBI/p/w3QQmVPnnQRzWj2M4dIFn8hEer2X1of/r8xd3+g88ziJ0wWPy67188Gj9leEUC1XLqdOFn8DorS559y0jP9dxqQHcvDtCa5dDthcrMfmfOO3hI85mNTfesf+5MneRInix+Xs2fCm/jfNLl/6Bh5hi9x1K6v+o0Q/6Sp3cFjLqc1fJDc56dNbvZTJosfl6Mmw1fyTClMxVHO/JTJ4sdlucJbvUyEKN91nt5Nlrpf56d1+8w4sKttFejnzCHpY5I2S3pE0p2S3tIrXhY/Mh0L1n5b3zivOnO8AzgVuETSqQft9hBwpu23Utxf8We9Yjai2W85fMfqECO1sJM8884cAJI6zhzztiy271qw/zeB9/cK2Ajx68wAA76hnTkO4nLg33sdsBHi73f4p3eESthROXO85pjS+ynus3hbr/0aIX5tMZUHcxXo58wBgKQLgD8E3ma759JkI8Tf4wgFnIFyP+ClXk9nDgBJpwN/D2yw3ffG2Tzaj81onTn+HDgc+GdJD0va2CtmIzK/roSe5OnnzGH7gkHiNUL8PeH61XnaIWLauZgjaeqrfTPE3xfhMeXtIL4s9Z7bb4T4tcVAbvbjss/h/XyDTRjXV/tmiF9nlnyzL+mjwG9QfI+/TeHJs5bC9vtNwIPApaVB48iJkfkO1efXuNnvO8kj6VjgtymWCk8DJilmlz4NXGv7ZGAHxUJCZiFVJ3jqasW2YL/DJM0CK4FngPN4dXrxZuBq4LrQJ1iF2SiZPzzFJM8Sznzb24G/AJ6mEH0XRTO/s5xyhGJ58djFPp+8/Wq74jYGqjT7R1IUDZwI/CgwA2yoeoDX2a8mhuxK2zio0uxfAHzP9vMAkr4KnAOsljRVZv+iy4ujYtbhL1ocYuKo5rdrVVnVexo4S9JKFY9w6Niv3gW8p9ynmv1qchRz+1W2cVDFe/d+SbcB3wLmKIoErwe+Btwi6VPle5+PeaK92D53ZPCYB0I1ZDUe8FW1X70KuOqgt7dSFBVmupFv2ojPvbtOCR5zdyuQ4cNSz/zMENRX+2aIf8+TJweP+fL+FUHiqF3fdr8R4tcWM7YJnCo0QnxvOTx80H3D17aK8U3gVKER4teaLH5cVod3YuMHwzuxFWTxEyX3+WmTR/uROWLr3uAxJ/eHEK268UIVJG2geK7xJHCD7WsO+v1y4IvATwEvAr9i+8lu8fLtWjHp3KgZwJmjojnD5cCOsrrqWopqq640IvOnt+8IHlMHAjhwQsg+v685Q/n66vLn24C/kaTycbevI2d+ZAIWcyxmznBw9dT8PmWdxS6KAttFaUTme8fO8EFbgTK/ep/fz5kjOI0Qv7bY0Krc7vdz5qhiztDZZ5ukKWAVxcBvURohfnv3K8Fjuj3yzO9HX3MGYCNFVdV9FFVW3+jW30NDxK81gcS3PSepY84wCdzYMWcANtneSFFN9SVJTwAvUXxButII8T0X3pApyDp84Bs1K5gz7AN+uWq8RohfXwwRPAJDkcWPiRlkwDdysvixyat6CZPFj0yER7SGKbwMu7ATmmaIX1cM5CXduGgq/JM2mA3UmuTMT5WBpndHTiPE13T4zNdcmLt0na/zE6bGnjyNEF/T0+GD7s19fmYY7Dzaj85UhP9GqLmDnPmpYhyqIigCWfyYZO/d+GgqvA9fIAPOvKSbKgacMz8ydR3wORdzJE2dB3zqUdwZ/mDS88BTFXdfA7wQ+BQGifkW20cNczBJt5fHrMILtis7m4ZgpOIPgqRNoZ4wGTPmUibfrpUwWfyEqbP4Me5Ti3rv21Kjtn1+Jj51zvxMZLL4CVM78SVtkLRF0hOSrgwU80ZJz0n6Toh4TaFW4lf0nTkUbmKAR8OkQq3EZ4HvTPmMvo7vzFDYvpfiluXMAuomfhXfmUwg6iZ+ZoTUTfwqvjOZQNRN/HnfGUnTFLYiG8d8To2lVuKX3nEd35nHgFttPzpsXElfoTApWi9pm6T83F/y9G7S1CrzM6Mli58wWfyEyeInTBY/YbL4CZPFT5j/B1ikPPqd9DyLAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 2 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"X = np.array([np.convolve([0] * 20 + [1] * 20 + [0] * 60, hrf(np.array(range(100))))[:100],\n", | |
" np.convolve([0] * 60 + [1] * 20 + [0] * 20, hrf(np.array(range(100))))[:100]]).T\n", | |
"\n", | |
"plt.imshow(X, aspect=0.1, interpolation='none')\n", | |
"plt.colorbar()" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Importantly, the model (i.e. all regressors) is then demeaned, too.\n", | |
"\n", | |
"The actual design matrix hence looks like this (pay attention to the range of values):" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 118, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"<matplotlib.colorbar.Colorbar at 0x1f524848080>" | |
] | |
}, | |
"execution_count": 118, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAIgAAAD7CAYAAACxDNw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARyklEQVR4nO2df5CdVXnHP9/dZRNJTEDDjxRQYAi0qW3VyfhjnFEK2CJ1oDO1FNpabGmd6VRr1VqxdqrT2hlsOyp/YKcpIhRtgVLbZkYrbSPItKMpQRitZAIYURPCDyGhxJBk995v/3jfxWW7995z85733vfe+3xm3sneu+8+52z2e5/znPM+5zmyTRB0YmrYHQiaTQgk6EoIJOhKCCToSggk6EoIJOhKJYFIulDSTkkPSboqV6eC5qCjXQeRNA08ALwR2A3cDVxu+/583QuGzUyFn30V8JDtXQCSbgYuAToKZHb6WL/gmLUVmhwcz849zZHWQVWx8bM/vcpPPtVKuveerx++3faFVdqrgyoCOQX43qLXu4FXL71J0tuBtwOsnFnDa0+/okKTg+MrD99Y2caTT7X479tfknTv9PoH11VusAaqCCQJ25uBzQBrV66fqHV9A23aw+5GJaoIZA9w2qLXp5bvBSXGzDltiGkqVQRyN7BB0hkUwrgM+OUsvRojJtaD2J6X9A7gdmAauN72N7P1bAwwpjXiT8srxSC2vwB8IVNfxpI2EyyQoDsGWiGQoBvhQYKOGJib5Bgk6I5xDDFBFwyt0dZHCKROipXU0SYEUiuiRaXnfUMnBFIjRZAaAgk6UKyDhECCLrTDgwSdCA8SdMWI1ojnhYdAaiaGmKAjRhzx9LC7UYnR9n8Np1gom0q6Uui1zUTSSyTdIeleSV+XdFHV3yEEUjOtcrGs19WLcpvJtcCbgI3A5ZI2Lrntj4Bbbb+CIsPvk1X7H0NMjdii5WyfwZRtJgbWlF+vBR6p2mgIpGba+aa5KdtMPgz8m6R3AquAC6o2GkNMjRRB6kzSBayTtH3R9fajaPJy4AbbpwIXATdJqvQ3Dg9SIwtBaiLft72py/dTtplcCVwIYPsrklYC64DHUzuxlPAgNdOykq4EnttmImmWIgjdsuSe7wLnA0j6MWAl8ESV/ocHqZGcK6mdtplI+hNgu+0twHuBv5H0bgoH9jZXrFIYAqmZdr5ZzLLbTGz/8aKv7wdel61BQiC1UjysG+1RPARSI0bMjfhSewikRmxyLpQNhRBIrSjnQtlQCIHUiAkPEvQggtSgI0aRMBR0ptj2MNr/xaPd+8YTG6eCLpi8K6nDIARSM+FBgo7YGnkP0rP3kk4rE2Hvl/RNSe8q33+RpH+X9GD57/H1d3e0KILU6aSrqaTIex54r+2NwGuA3ymTZa8CttreAGwtXwfPo8hJTbmaSs+e2d5r+2vl188AOyjyIy8BFupV3wj8fE19HFmKIFVJV1PpKwaRdDrwCmAbcJLtveW3HgVO6vAzz6vVPmlMzEqqpNXAPwK/Z/t/pR+q3rYlLZu59Lxa7S9Yb1TDp6WhheImZiVV0jEU4vis7c+Vbz8mab3tvZLWUyExdpzpI2m5kfQUiApX8Slgh+2PLfrWFuAK4Ory33/pZctTor16xVF2tTNTzxzKbjMHNsy1x1wgFDmObwW+Iem+8r0/pBDGrZKuBL4DXFpLD0eYYogZc4HY/k/ouBx4fj+NtWenOHjqqn5+JInVO5rpQSBWUoMuLExzR5mBCqS1QuzbkL/J1Tuym8zEBAwxQTUiJ7UP2ivNgXPmBtnkUClmMc19zpLCaPu/hrOwUJZrqT3lIGtJly56sPp3VX+HGGJqJtcQs6jC0HMHWUvasvgga0kbgA8Ar7O9T9KJVdsdqEBWrJjj7DP39r6xf8s12KxO5llMSoWh3wKutb0PwHbl1e0YYmqm7amki94FZJarMHTKknvOBs6W9F+Sviqp8kneA/Uga2YOcd6JO7Pb/TI/md1mDmwxnz7N7VVAJoUZYANwLkWBmbsk/YTt/VUMBjWScYhJqTC0G9hmew74tqQHKARz99E2OlCBrJo+zGtXPZjdbmM9CFkFknKQ9T9T1Cn7tKR1FEPOriqNhgepmVwCSawwdDvwM5LuB1rA+2w/WaXdwc5iaHHmzIFBNjlUcicMJVQYMvCe8spCeJCaiaX2fhrTFCdMN3PNog5smJ+AhKGgAvG4P+jIxCQt52IKsULHDLLJoeMQSNCNCFKDjtgRgwRdEa2YxaRjTMujfpp9f0QMEnQkstr7xMBhzw+yyeHixm4bTiY8SM3ELKYP2piDnqCs9ghSg17EENMHbZuDo/4/1icxiwk6YodAgh7ENLcP2ohDI/4f1i+jPqKGB6kRI9oxi0mnDRxqcNHYOhhxBxIepFYmKUgtNw9vB/bYfnO5P+Nm4MXAPcBbbR/pZsNo4jzIqLuQfgbId1FUWV7go8DHbZ8F7KM4Nz5Ygq2kq6kkCUTSqcDPAdeVrwWcB9xW3pJUiruf4vb9XE3FQLutpKuppHqQTwB/QBFnQjGs7LefezS73E5zoCjFvbBjff+Tk5ULggEr7WooKceBvBl43PY9R9OA7c22N9netPZF08x5JvvVZOy0K4WUCkPlfb8gyZKqVgtILqR7saSLgJXAGuAa4DhJM6UXWW6neQDZgtSUCkPlfS+kiBe35Wg35TiQD9g+1fbpFDvKv2T7V4A7gLeUtyWV4p480gLUxCD1uQpD5WxxocLQUv6UYgKRpbpwFf/8fuBmSR8B7qWo596VI8ywZ37CDqZK9yDrJG1f9HpzeVLGAstVGHr1YgOSXgmcZvvzkt53FL39f/QlENt3AneWX++iUHXQCYPTZyiVKgxJmgI+BrztaG0sx0AjvAOtFdz19NmDbLIBDKzC0AuBlwF3lmf5nAxskXSx7cWeqS+aPQUYB/KtpHatMGT7aWDdwmtJdwK/X0UcMGCBPHN4JV9++Kzsds/kYHab2cgkkMQKQ9kJD1InCwtlucz1qDC05P1zc7Q5WIEcmsI7V9dguLkeJBKGgu40+DlLCgMVyPRhOC5/FcxGs/xZoKNDeJA6MSOfDzJgD9Jmza5nB9nkkGn2k9oUwoPUTXiQoCsjngIzUIForsXsI/uz2/VMQ7PKMq+DDIPwIDUTs5h+mG/hp/blt3viut73DIsRF8hob/sKamewJajaLdoHfpDd7lSDPUgMMUFnTCy194XB8xNUxA5GPgYJD1IzMcQE3QmBBF0JgQSdkGOI6R+NdlTfNzGLCboRHqQfJDQ7O9Amh04IJOhIxCD9IQnNTJgmQyBBNxQJQ30wiTHIiBOP++vGiVcCvSoMSXqPpPslfV3SVkkvrdr9gXsQJikGyRikJlYYuhfYZPugpN8G/hz4pSrthgepm3wepGeFIdt32F7Yh/pVihIRlRiwBwHVkGDc6InCACsMLeFK4F+TW+/ABPn7wSP6msVUqjD0vHalXwU2AW+oaisEUid5F8p6VRgCQNIFwAeBN9g+XLXRJIFIOo6iyvLLKJzmbwA7gVuA04GHgUttd09Zn7QgFQZWYQhA0iuAvwYutP14jkZTg9RrgC/a/lHgpyhqtl8FbLW9Adhavg6WkilILevRLlQY2gHculBhSNLF5W1/AawG/kHSfZIqVx3q+XGWtBZ4PWX1vDKCPiLpEuDc8rYbKaofvr+HNZierIlTzmcxvSoM2b4gX2sFKX+tM4AngE9LulfSdZJWASfZ3lve8yhw0nI/vLhW+5FWcysB1UbGhbJhkBIQzACvBN5pe5uka1gynNi2tPxnpZyqbQZYu/JkM1WDB6kjCSmHSY/+s5iUv9ZuYLfthdrft1EI5jFJ6wHKf7MERWPHuHsQ249K+p6kc2zvBM4H7i+vK4Cr6aNWu6fyf9rbx67IbjNXPyclH+SdwGclzQK7gF+n8D63SroS+A5waT1dHHEmQSC276NYmVvK+X23WEO88Oypq7Lb9LcyxEoNHz5SmLBVq8EiJmeICY6SEEgDeOqc/L/G/LZMQ2EIJOhKCKQfVEuQeuCcuew22ysz/GVj20PQkxDI8Dn7zL29b+qTfSvyeKVRX2ofC4E0mRhiGsB5J+7MbnPHTIZTRWOhLOhJCGT4vHZV/kNobpiunM4ZK6lBb9QebYWMhUDOnDmQ3eYKWtWNRAwS9CKGmKA7IZDhc8J0/oyyGeXJnQ0PEnQnBDJ8VuiY7DancqS1T0hWe3CULKyDpFxJ9noXkFkh6Zby+9sknV71dwiB1I2ddvVgUQGZNwEbgcslbVxy25XAPttnAR8HPlq1+yGQmsnoQXoWkClf31h+fRtwvlQtAWcsYpA5Z1jUWoJzRJd5F8pSCsg8d4/teUlPAy8Gvn+0jY6FQJpMH0FqrwpDQ2EsBHLYNaQcZvroZ6wwlFJAZuGe3ZJmgLXAk8k9WIaIQerEZAtSWVRAptzheBmwtP7HFoptsABvAb5kpxnvxFh4kEO1xCB5yLWSWsYUCwVkpoHrFwrIANttbwE+Bdwk6SHgKQoRVWIsBNJoBltA5hDwi/laDIHUSiQMNYSD1YbZZWnnsGlHwlDQg9HWx3gI5JBrKEqTpQZVDDFBNwzEEDN8Djl//fdsT+lHWx/jIZAmMxFDjKR3A79J8Xn4BkWNsvUUTxRfDNwDvLV8yjhw6vAgzhWDjPgQ03OpXdIpwO9SHFTzMopVvMsocg0+XuYe7KPIRQgWk1oCs8EaSh1iZoAXSJoDjgX2Aufxw2LyNwIfBv4qdwdTmKvFg1SnWChr8F8/gZ4exPYe4C+B71II42mKIWV/WWAeityEU5b7+Ykvxd1OvBpKyhBzPEWm0hnAjwCrgAtTG7C92fYm25tmp4896o6OKrKTrqaSMsRcAHzb9hMAkj4HvA44TtJM6UWWPdxmUMw5/2TMORbfGh5fpJCSD/Jd4DWSji3zGxdKcd9BkXMAfZTiniyKZzEpV1NJqdW+TdJtwNeAeYqjNzcDnwdulvSR8r1P1dnRbuyZPz67zSO5HGKDh48UUktxfwj40JK3d1FkWgedGIONU2OxknrX02dnt3mglakozSR4kKACo62P8RDIlx8+K7vNZw6vzGJH7dEeY8ZCII3FNHoRLIWxEIh3rs5v9FD1HSGi2YtgKYyFQBpNCGT4HPdAfpuPVK+CWRACCToSMUjQi5jFNIA13342u83pwzn+sMn7bish6UXALcDpwMPApbb3Lbnn5RT5OmuAFvBntm/pZTs2b9dJ3s3b3bgK2Gp7A7CVJSejlxwEfs32j1Oka3xC0nG9DI+FB5nds6/3TX2iI5k2hA9mhLkEOLf8+kbgTuD9i2+w/cCirx+R9DhwArC/m+GxEEiT6WMdpEoBmZNsL5yq9ChwUtc+Sa8CZoFv9TI8FgLxvv35jbYyeZB0gXQtICPpP4CTl/nWB5/fnC113mwhaT1wE3CF7Z7+bSwE0lhsaOUZY2xf0Ol7kh6TtN723lIAj3e4bw1FHs8HbX81pd2xEEj7wA+y23R74B6kCguVha6mQ3ZfWZXon4C/tX1bquGYxdTNYGYxVwNvlPQgRQ7x1QCSNkm6rrznUuD1wNsk3VdeL+9leCw8iOfne9/Ut9FMNgaQb2r7SYpc4aXvb6fYEYntzwCf6df2WAikuRh6x4GNJgRSJyZbkDosQiB1E09zg66EQBpAtXr1y5Pl7zqYh3V1Mh4CaSoG4nH/8NFM/hOnmMvklcKDBJ3Jt9Q+LMZCIJrN70E0n2d3f8LzsEYzFgJpNA3euZ/CWAhEs7P5jT4bMQiMiUAaix2zmEYwU8OvkWttJTxI0BnjXJlpQyIEUidRq70ZaCZ/ndRMhZbjcX/QGQMOD9IAmhqkOhKGgh6MepCqiseq9teY9ATwncTb11HhSPEMNl9q+4QqjUn6YtlmCt+3nVzBelAMVCD9IGl7j5OoG2Fz3IltD0FXQiBBV5oskNSNy8O2OdY0NgYJmkGTPUjQAEIgQVcaJxBJF0raKekhScuVUjoam9dLelzS/+SwN0k0SiCSpoFrgTcBG4HLJW3MYPoG+jhGLfghjRIIxfkzD9neVZ7BezNF/a1K2L4LeKqqnUmkaQI5BfjeotcdT9MMBkPTBBI0jKYJZA9w2qLXQz1NM2ieQO4GNkg6o6ypdRlF/a1gSDRKIOUZvO8Abgd2ALfa/mZVu5L+HvgKcI6k3ZKurGpzUoil9qArjfIgQfMIgQRdCYEEXQmBBF0JgQRdCYEEXQmBBF35PxGEeUkX2K7lAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 2 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"X-= np.mean(X, axis=0)\n", | |
"model = GeneralLinearModel(X)\n", | |
"\n", | |
"plt.imshow(X, aspect=0.1, interpolation='none')\n", | |
"plt.colorbar()" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The (demeaned) model is then fit to the (demenaed) data, and we can test for the effect of _task_." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 119, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"A > rest: 99.13 (z = 5.58)\n", | |
"B > rest: 261.99 (z = 11.20)\n", | |
"B > A: 162.87 (z = 7.06)\n" | |
] | |
} | |
], | |
"source": [ | |
"model.fit(Y)\n", | |
"con1 = model.contrast(np.array([1, 0]))\n", | |
"con2 = model.contrast(np.array([0, 1]))\n", | |
"con3 = model.contrast(np.array([-1, 1]))\n", | |
"\n", | |
"print(f\"A > rest: {con1.effect[0,0]:.2f} (z = {con1.z_score()[0]:.2f})\")\n", | |
"print(f\"B > rest: {con2.effect[0,0]:.2f} (z = {con2.z_score()[0]:.2f})\")\n", | |
"print(f\"B > A: {con3.effect[0,0]:.2f} (z = {con3.z_score()[0]:.2f})\")" | |
] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3", | |
"language": "python", | |
"name": "python3" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 3 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython3", | |
"version": "3.6.8" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 2 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment