{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "A notebook for thoughts about Cyril Pernet's article on GLM misconceptions" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import numpy.linalg as npl\n", "import matplotlib.pyplot as plt\n", "import matplotlib.mlab as mlab" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import statistical distributions from scipy:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.stats import t as t_dist, f as f_dist, gamma" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "A routine to scale the design matrix for display:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def scale_design_mtx(X):\n", " \"\"\"utility to scale the design matrix for display\n", "\n", " This scales the columns to their own range so we can see the variations \n", " across the column for all the columns, regardless of the scaling of the \n", " column.\n", " \"\"\"\n", " mi, ma = X.min(axis=0), X.max(axis=0)\n", " col_neq = (ma - mi) > 1.e-8\n", " Xs = np.ones_like(X)\n", " mi = mi[col_neq]\n", " ma = ma[col_neq]\n", " Xs[:,col_neq] = (X[:,col_neq] - mi)/(ma - mi)\n", " return Xs" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Display the design matrix nicely:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def show_design(X, design_title, **kwargs):\n", " \"\"\" Show the design matrix nicely \"\"\"\n", " plt.figure()\n", " plt.gray() # Gray colormap\n", " plt.imshow(scale_design_mtx(X), interpolation='nearest', **kwargs)\n", " plt.title(design_title)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Very simple t statistic from contrast and Ordinary Least Squares fit:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def fit_ols(Y, X):\n", " \"\"\" betas, fitted data, and residuals from OLS linear fit.\n", " \n", " This is OLS estimation; we assume the errors to have independent\n", " and identical normal distributions around zero for each $i$ in \n", " $\\Epsilon_i$ (i.i.d)\n", " \"\"\"\n", " Y = np.asarray(Y)\n", " X = np.asarray(X)\n", " # Get the estimated betas\n", " betah = npl.pinv(X).dot(Y)\n", " fitted = X.dot(betah)\n", " resid = Y - fitted\n", " return betah, fitted, resid" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "def t_stat(Y, X, C):\n", " \"\"\" betas, t statistic and significance test given data, design matrix, contrast\n", " \n", " Ordinary least squares estimation - see `fit_ols` function.\n", " \"\"\"\n", " Y = np.asarray(Y)\n", " X = np.asarray(X)\n", " C = np.atleast_2d(C)\n", " # Calculate the parameters\n", " B, fitted, resid = fit_ols(Y, X)\n", " # Residual sum of squares\n", " RSS = (resid**2).sum(axis=0)\n", " # Degrees of freedom - number of observations - number of fitted parameters\n", " df = X.shape[0] - npl.matrix_rank(X)\n", " # Mean residual sum of squares\n", " MRSS = RSS / df\n", " # Standard error of contrast estimate C.dot(B)\n", " SE = np.sqrt(MRSS * C.dot(npl.pinv(X.T.dot(X)).dot(C.T)))\n", " t = C.dot(B)/SE\n", " ltp = t_dist(df).cdf(t) # lower tail p\n", " p = 1 - ltp # upper tail p\n", " return B, t, df, p" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make the simulated data from C.P.s figure 2 (was figure 1):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "a_block = np.array([-0.1, 0, 0.1])\n", "baseline = 10\n", "activation = 11\n", "on_off = np.hstack((a_block + baseline, a_block + activation))\n", "n_on_off = 3\n", "e = np.tile(a_block, (2 * n_on_off,))\n", "y = np.tile(on_off, (n_on_off,))\n", "plt.plot(y)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ "[]" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UlVW6B/DvUTBL0/AHB/Bgkkb8FCnNYomeVLRi5GI6\nJeqVBdrMmhnvHctphmzdldSM4dhcw7FpzZ3U6Oaycs2ETCp3JGS0jKwwbdIBJBjQQwcF0QgMgX3/\n2B4EOZzf79nvft/ns1Zr6YEDT+5zvme/+33e/RoYYwyEEEI0bYjoAgghhCiPwp4QQnSAwp4QQnSA\nwp4QQnSAwp4QQnSAwp4QQnTAadhnZ2fDaDQiPj6+97G9e/ciNjYWQ4cORUVFhd3nNTQ04KGHHkJs\nbCzi4uKwbds231VNCCHELU7DPisrC8XFxf0ei4+Px3vvvYfZs2cP+rzAwEBs3boVX331FcrLy/Hq\nq6/izJkz3ldMCCHEbQHOviE5ORl1dXX9HouKinL6g0NCQhASEgIAGDlyJKKjo2GxWBAdHe1ZpYQQ\nQjzmlzX7uro6nDhxAjNnzvTHryOEEHITpzN7b7W1tWHp0qXIz8/HyJEjB3zdYDAoXQIhhGiSO7vd\nKDqzv3btGpYsWYKVK1ciPT190O9jjKnqv+eff154DVSTtuqimqgmX//nLq/DfrBfyhjD6tWrERMT\ng3Xr1nn7awghhHjBadhnZGQgKSkJlZWVCA8Px86dO1FYWIjw8HCUl5cjNTUVjzzyCADAYrEgNTUV\nAPDRRx/hrbfewuHDh5GYmIjExMQBXT2EEEL8w+ma/Z49e+w+bm9ZJiwsDPv37wcAzJo1Cz09PV6W\nJ4bZbBZdwgBUk+vUWBfV5BqqSTkG5snijy8LMBg8Wn8ihBA9czc7absEQgjRAQp7QgjRAQp7QgjR\nAQp7QgjRAQp7QgjRAQp7QgjRAQp7QgjRAQp7QgjRAQp7Qgjxg54eoLNT3O+nsNcRiwWYMwf4+99F\nV0I8sXcvMGuW6CqIJzo6gGXLgN/+VlwNFPY6ceoU8OCDwL/+BXz5pehqiDsYA156CVi/HvjkE7Gz\nQ+I+qxV46CEgIAD4xS/E1UFhrwMHDwLz5/NZxerVQGOj6IqIq65dA9as4bP6jz8GjEbgm29EV0Vc\ndfo08MADwIIFwO7dwPDh4mqhsNe4114DsrOBwkLgiSeA0FAKe1m0tgIPPwxcuAAcOQJMmEDjJ5OS\nEsBsBnJzgRdeAETflI/CXqO6u4Gnnwby84EPPwSSkvjjFBZyqK3lYxYfD7z3HmC7oyeNnxxefx1Y\nsYIfka1aJboaTvF70BL/++47YPly4Ntv+aF/UNCNr4WE0DKA2pWXA4sXA889B6xd2/9rNH7q1tPD\nx23vXn40ds89oiu6gcJeYywWYNEiYOpU/oIbNqz/12lmqG579wI/+xmwaxdw/aZv/dD4qVdHB5CZ\nyd+D5eXAuHGiK+qPlnE05NQpfjJoyRJg586BQQ8AwcFAczPQ1eX/+sjg+nbc/O1v9oMeoLBXq74d\nNyUl6gt6gMJeM2wdN1u2ABs2DH4yKCAAGDsWaGryb31kcDd33EybNvj3Utirj5o6bhyhZRwN+MMf\ngBdf5B03thOxjtgCIyxM+dqIY62t/Ejsttv4Gq/tROxgKOzVpaSEnx97+WX1nIgdDM3sJWbruNm2\nrX/HjTMUGOpg67iJi+Mf1M6CHqCxUxM1dtw4QjN7Sdk6bq5cGdhx4wwFhni2jpsNG4D/+A/Xn2c0\n8r777m5g6FDl6iODU3PHjSM0s5eQxQLMng2MGQP83/+5F/QAhb1oe/fyjqnXX3cv6AF+0n30aODi\nRWVqI47Z9rg5epR/YMsS9ICTsM/OzobRaER8fHzvY3v37kVsbCyGDh2KioqKQZ9bXFyMqKgo3H33\n3di8ebPvKta5kyedd9w4Q2Evhq3j5umngUOHBu+4cYbGTwwZOm4ccRj2WVlZKC4u7vdYfHw83nvv\nPcyePXvQ53V3d2Pt2rUoLi7G6dOnsWfPHpw5c8Y3FetY3z1uHHXcOEMX5vifrePm3Xf5jNBRx40z\nISEU9v4mS8eNIw7X7JOTk1FXV9fvsaioKKc/9Pjx45gyZQomTZoEAFi2bBn27duH6OhojwvVO1vH\nzb59rp+IHQzNDP2rb8fN0aOunYh1JDSUPqz9SaaOG0cUOUF7/vx5hIeH9/7dZDLhk08+GfT7N27c\n2Ptns9kMs9msRFlS6u7m26IePMg7biZP9v5nUtj7T20tX66ZPx/YutU3J1Vp/Pxnxw5+FL13L78X\nhEhlZWUoKyvz+PmKhL3BzfWFvmFPbvCm48YR28yQMfE78WmZpx03zoSGAmfP+u7nkYHU2HFz80Q4\nNzfXrecr0o0zYcIENDQ09P69oaEBJpNJiV+lWd523DgyfDhfUmhp8d3PJP1503HjDM3slSVzx40j\nXoU9Y8zu49OnT0d1dTXq6urQ2dmJd955B2lpad78Kl2x7XHz2GOed9w4Q4GhDF913DhCY6ecpia5\nO24ccRj2GRkZSEpKQmVlJcLDw7Fz504UFhYiPDwc5eXlSE1NxSOPPAIAsFgsSL3+yg4ICMD27dux\ncOFCxMTE4IknnqCTsy46eBCYN4933Dz3nHLLLBQYvufLjhtHaOyUoYWOG0cMbLDpub8KMBgGPULQ\nG1vHzZ//7H3HjTMrV/IXtczdBWrSt+Nmzx7vO24c+e47PuNsb6dzLr4iY8eNu9lJV9CqgKd73HiD\nerV9x5M9brwxYgQQGAhcvqzs79GLHTv4HjfvvitP0HuC9sYRTKmOG2dCQ4Fz5/zzu7RMqY4bZ2wf\n1nfc4b/fqTVq7LhREs3sBVKy48YZWvf1npIdN87QhVXe0WrHjSMU9oKcOgU8+KCyHTeOUNh7jjEg\nL4/fVUqpjhtnaPw819QEzJ2rzY4bRyjsBSgu5ldUbt6sbMeNIxQWnrl2DXjySeCdd5zfVUpJNH6e\nsXXcpKRos+PGEVqz97PXXgNeeMH1u0ophcLCfb7e48YbNH7u++ADICNDro4bX6KZvZ/YOm7y8/3X\ncePIqFG8prY2sXXIwt8dN85Q2Ltnxw7eCKH1jhtHaGbvB999x1u7Ll/2b8eNIwbDjcC4+27R1ahb\neTk/t7JhA7B2rehqOAp71+it48YRmtkrzGLhu+UFBfm/48YZCgznbB03f/qTeoIeoLFzhR47bhyh\nsFeQreNm8WIxHTfO0E1MBqeGjhtHKOwd02vHjSO0jKOg9HTg178G/v3fRVdiHwXG4D74gK/zfvwx\nMGGC6GoGCgriM9eODuDWW0VXoz5r1/Kum//+b9pSwoZm9gppa+NBumKF6EoGR2E/uC+/BB5+WJ1B\nD/AAoyOzwX35JbB6NQV9XxT2Cjl7lt9VaoiK/4Up7AdXXQ1ERoquwjEaP/u6u4G6Ot/c1U1LVBxF\ncquuVn+XC4XF4Gj85FVfD4wfT8tbN6OwV0hVFc0MZVZVRWEvKxneeyJQ2CuEZobyunoVsFqBO+8U\nXYljNH72yfDeE4HCXiEyvODGjeNbK3d2iq5EXWpqgEmTeNuemlHY2yfDe08ECnuFyHAoOWQIEBxM\nHR03k2EJB6CwH4wM7z0RKOwV0NrKlwJCQkRX4hy17w0kQycOQGE/GJrZ20dhr4DqamDKFDl6fCkw\nBpIlLOjWkgN1dvI7sEVEiK5EfSjsFSDTYSSF/UCyLOMEBwMtLUBXl+hK1KO2FjCZ1Lc1iRpQ2CtA\nlpkhQGFvjyzLOAEBwNixfB8Ywsn03vM3CnsFyPSCo7Dvr62Nn3NR6zYJN6Px60+m956/OQz77Oxs\nGI1GxMfH9z7W0tKClJQUREZGYsGCBWhtbbX73JdeegmxsbGIj4/H8uXL8f333/u2chWTZWYIUFjc\nrLpa/dtc9EXj159MS6j+5vAlnZWVheLi4n6P5eXlISUlBVVVVZg3bx7y8vIGPK+urg5/+tOfUFFR\ngS+//BLd3d14++23fVu5SjEmz5ovQGFxM5k+qAEav5vRzH5wDsM+OTkZQTfdbaOoqAiZmZkAgMzM\nTBQWFg543qhRoxAYGIj29nZ0dXWhvb0dE2Q5LvbSxYu8C2fsWNGVuIbCoj/ZwoLGrz/Zxs+f3L5G\n0Gq1wmg0AgCMRiOsVuuA7xkzZgzWr1+PiRMn4tZbb8XChQsxf/78QX/mxo0be/9sNpthNpvdLUs1\nbDNDGdouAcBo5Cf4uruBoUNFVyNeVRUwe7boKlwXGgr84x+iq1CHjg45trnwVFlZGcrKyjx+vlcX\nhBsMBhjspFpNTQ1eeeUV1NXVYfTo0fjhD3+I3bt3Y8Ugm7v3DXvZybSEA/AWtdGjgeZm3sqnd9XV\nwJo1oqtwXWgov5MW4dtcRERod9Jy80Q4NzfXree7fRrKaDTim+uXXDY2NiLYTkJ89tlnSEpKwtix\nYxEQEIDHHnsMx44dc/dXSUnGw0haCrhBtvGjC6tukG3s/M3tsE9LS0NBQQEAoKCgAOnp6QO+Jyoq\nCuXl5ejo6ABjDCUlJYiJifG+WgnIdoIPoLC3uXSJb3NxfZVSCqGhtN2FDXXiOOYw7DMyMpCUlITK\nykqEh4dj165dyMnJwaFDhxAZGYnS0lLk5OQAACwWC1Kv35U5ISEBq1atwvTp0zF16lQAwI9+9COF\n/1fUQbZlHIDC3ka28y3AjbBnTHQl4tHM3jEDY2JfJgaDAYJL8BnGgNtvB86f5+vgssjJAUaNAjZs\nEF2JWLt3A3/9KyBbl3BQEL8NpiwdYEqZMwd4/nlg7lzRlfiHu9kpyaUjcmhsBEaMkCvoAZrZ28h4\nVAbQ+NnQMo5jFPY+RGEhNxnPtwA0fgDw7bf8RjxhYaIrUS8Kex+Sdc2QwoKj8ZPX2bNybXMhAv3T\n+JCsM0O6gYl821z0RWFPSziuoLD3IdnDQiPnyT1y4cKNLYNlQ2Ev71GZP1HY+5CsL7iRI/lVh1eu\niK5EHFnHDqALqwC5x89fKOx9pLsb+PprfjtCGel9dijrURlAF1YBtIzjCgp7H2lo4EsAI0aIrsQz\neg97Wc+3ADR2AM3sXUFh7yOyv9j0Hhgyj5/ex66lhd9onDbyc4zC3kdknhkCFBgyL+OMGsWXEdva\nRFcihozbXIhAYe8jMocFoO+wZ4z3acs6fgaDvsdP5qMyf6Kw9xHZX3B6DguLhe9pNGqU6Eo8p+fx\nk/295y8U9j4i+zKOni+skv2oDNB32FMnjmso7H3g2jXgX/8C7rpLdCWe03NYyP5BDdD4yf5h7Q8U\n9j5QV8c3YLrlFtGVeI7CQnQV3tHr+DGmjfHzBwp7H9DCzHDMGKC9nd+0WW+0sIyj12W4piYgMJC/\nfoljFPY+oIWwMBj0Gxha+LDW68yeZvWuo7D3Aa284PQYGN3dQG0t3x5XZnocO0A77z1/oLD3AS3M\nDAF9BkZ9PTBuHHDbbaIr8Y4exw6gThx3UNj7gBaWcQB9BoZWPqjHjeO7lnZ2iq7Ev2hm7zoKey9d\nvcrXuSdNEl2J9/S4Va5WwmLIEL43jN7OuWhl/PyBwt5LX38N3Hknv/GF7PS4Va5WjsoA/R2Z9fTI\nvc2Fv1HYe4nCQm5aWcYB9Dd+Fgvf4uL220VXIgeHYZ+dnQ2j0Yj4+Pjex1paWpCSkoLIyEgsWLAA\nra2tdp/b2tqKpUuXIjo6GjExMSgvL/dt5SqhpcNIvYUFQOMnMy2NnT84DPusrCwUFxf3eywvLw8p\nKSmoqqrCvHnzkJeXZ/e5P//5z/Hoo4/izJkzOHXqFKKjo31XtYrQzFBe167xbhyZt7noS2/XSVAn\njnschn1ycjKCgoL6PVZUVITMzEwAQGZmJgoLCwc87/Llyzh69Ciys7MBAAEBARg9erSvalYVLS3j\nBAcDzc1AV5foSvyjthYwmYBhw0RX4ht6+7Cmmb173D6taLVaYTQaAQBGoxFWq3XA99TW1mL8+PHI\nysrCyZMncd999yE/Px+3DdLMvHHjxt4/m81mmM1md8sSRksvuIAAfmvFpia+14/WaWnsAH2GfVKS\n6Cr8p6ysDGVlZR4/36seEoPBAIOd28N0dXWhoqIC27dvx4wZM7Bu3Trk5eXhhRdesPtz+oa9TNra\n+C3RwsNFV+I7tsDQQ9hr6agM0F/Y620Z5+aJcG5urlvPd7sbx2g04pvrC4ONjY0ItnPjR5PJBJPJ\nhBkzZgAAli5dioqKCnd/leqdPcsvsx+ioZ4mPQWGls63APoau+5uvtus7Ntc+JPbMZWWloaCggIA\nQEFBAdLT0wd8T0hICMLDw1FVVQUAKCkpQWxsrJelqo/WlgEAfZ3k09r4GY3AhQs8CLWuvh4YPx64\n9VbRlcjDYdhnZGQgKSkJlZWVCA8Px65du5CTk4NDhw4hMjISpaWlyMnJAQBYLBakpqb2Pvf3v/89\nVqxYgYSEBJw6dQobNmxQ9v9EAK3NDAF9zQ61towzbBgwejRw8aLoSpSntyUcX3C4Zr9nzx67j5eU\nlAx4LCwsDPv37+/9e0JCAj799FMvy1O3qipg1izRVfhWaChw+rToKpR39SpgtfKrn7XE9mF9vYdC\ns7R2VOYPGlpt9j8tvuD0MrOvqeH7GWlhm4u+9DJ+WnzvKY3C3gu0jCMvrS3h2OjlnAst47iPwt5D\nra38Fn4hIaIr8S29hL0WP6gBfY2fFj+slURh76HqamDKFH47Py2x7XzJmOhKlKXVsNBD2Hd2AufO\nARERoiuRC4W9h7Q6Mxw+nN+1qaVFdCXK0uoyjh7CXmvbXPgLhb2HtBoWgD5uYqLVD2s9hL1Wj8qU\nRmHvIS2/4LR+E5O2Nn7OZcIE0ZX4HoU9GQyFvYe0OjMEtB8Y1dXa2+bCxjZ2Wj7nQp04ntHgy115\njGl7GUcPYa/VsBgxAggMBC5fFl2Jcmhm7xkKew9cvMi7cMaOFV2JMvQQ9loOCxo/Yg+FvQdsM0Ot\ntV3aaD0stHxUBmj7BHtHhza3ufAHCnsPaD0stB72Wl7GAbR9gr2mhvfXDx0quhL5UNh7QOuHkXoI\nexo/OWl97JREYe8BPcwMtRoWly7xHS+1vCuklsePOnE8R2HvAa3PLkaN4jcdb2sTXYnvaf18C6Dt\nsNf6e09JFPZuYkz7LziDQbvrvlofO4DCnthHYe+mxkbeyzx6tOhKlKXVwND6yXVAu2MH0DKONyjs\n3aSXmYVWA0MP46fVsfv2W+DKFSAsTHQlcqKwd5MeZoaAdgND6yfXASAoiPejd3SIrsS3zp7V7jYX\n/kD/bG7SQ1gA2gx7rW9zYWMwaPPCKlrC8Q6FvZv0sAwAaDPsL1zgF+NodZuLvrR4gl0v7z2lUNi7\nSQ8zQ0CbYa+XozJAu+Onh/eeUijs3dDTA3z9Nb8dodZpMSz08kENaHf89PJhrQSHYZ+dnQ2j0Yj4\n+Pjex1paWpCSkoLIyEgsWLAAra2tgz6/u7sbiYmJWLRoke8qFqihgS8BjBghuhLlaXHNV08zQy2G\nvZ7GTwkOwz4rKwvFxcX9HsvLy0NKSgqqqqowb9485OXlDfr8/Px8xMTEwKCRyxX1NDMcP57vid7Z\nKboS36FlHHm1tPDXYnCw6Erk5TDsk5OTERQU1O+xoqIiZGZmAgAyMzNRWFho97nnzp3DgQMHsGbN\nGjCN3DZHT2ExZAh/Y1mtoivxHT19WGst7PWwzYXSAtx9gtVqhfH6LlJGoxHWQdLgqaeewpYtW3Dl\nyhWnP3Pjxo29fzabzTCbze6W5Rd6O4y0BUZ4uOhKvMcY79PWy/hpMez1MnaDKSsrQ1lZmcfPdzvs\n+zIYDHaXaN5//30EBwcjMTHRpeL6hr2aVVUBKv0cUoSWAsNiAW6/nW/ypgdaGjuAwh4YOBHOzc11\n6/lud+MYjUZ8c72Bt7GxEcF2FtGOHTuGoqIiREREICMjA6WlpVi1apW7v0p19LSMA2grMPS0hAPw\ncy4tLXz3Ui2gThzvuR32aWlpKCgoAAAUFBQgPT19wPds2rQJDQ0NqK2txdtvv425c+fizTff9L5a\nga5dA+rrgbvuEl2J/2gp7PU2MwwI4J1jTU2iK/ENvY2fEhyGfUZGBpKSklBZWYnw8HDs2rULOTk5\nOHToECIjI1FaWoqcnBwAgMViQWpqqt2fo4VunLo6Hn633CK6Ev/RWtjrbWaolfHTw7bi/uBwzX7P\nnj12Hy8pKRnwWFhYGPbv3z/g8Tlz5mDOnDkelqceegwLLfXaV1UBGlhJdItWwr6pCQgMBMaMEV2J\n3OgKWhfpcWahpf1V9Dp+Wgh7PY6dEijsXaS3E3yAdsKiuxuordXHNhd9aWX8KOx9g8LeRXpdxmlq\n4nsCyay+Hhg3DrjtNtGV+JdWwp46cXyDwt5FepxdDBvG+9IvXhRdiXf0OHaAdsJer+PnaxT2Lrh6\nlb9pJk0SXYn/aSEw9HhUBmhj7AAKe1+hsHfB118Dd97Je5f1RguBocfzLYA2uql6evS1zYWSKOxd\noOeZhRbCXq/jFxrKN7KTeR9Ci4UvJd5+u+hK5Edh7wK9zgwB7YS9Hpdxhg/nJ6VbWkRX4jm9flAr\ngcLeBXoNC0D+pQA9bnPRl+wf1tSJ4zsU9i7Q8+xC9guramuBCRN4Z5EeyR72en7v+RqFvQtoGUd0\nFZ7T81EZoI3x0+t7z9co7J347jvg0iVt3MDDE7KHhZ4/qAFtjJ+eP6x9icLeibNn+XrvEJ3+S9nC\nQtaODr3PDGUO++5uvtvs5MmiK9EGnUaY6/Q+Mxw5Ehg6FHDh7pKqRMs48oZ9fT2/Ccutt4quRBso\n7J3Qe1gAcgeG3j+sZe6moiUc36Kwd0LvywCAvGF/9Sq/qOjOO0VXIo7M3VT03vMtCnsn9D4zBOQN\n+5oavp+RHre5sJF17AAKe1+jsHeClnHkXQqgD2q+1UB3N9DWJroS99Eyjm9R2DvQ2gp0dPCw0zNZ\nlwJoZggYDPLO7mn8fIvC3oHqan53Iw3cL90rMocFzQzlHL/OTuDcOSAiQnQl2kFh7wCFBSdjWAC0\njGMj4/jV1gImk363uVAChb0DdBjJyRgWAI2fjYzjR2PnexT2DtAJIk7GsGhr4+dcTCbRlYgn4/jR\nUbXvOQ377OxsGI1GxMfH9z7W0tKClJQUREZGYsGCBWhtbR3wvIaGBjz00EOIjY1FXFwctm3b5tvK\n/YBmF9yYMUB7Oz9ZLYvqan6ZvV63uehLxm4qWoLzPadvhaysLBQXF/d7LC8vDykpKaiqqsK8efOQ\nl5c34HmBgYHYunUrvvrqK5SXl+PVV1/FmTNnfFe5whijsLcxGHhgyNSRQ2N3g4zdVDR+vuc07JOT\nkxEUFNTvsaKiImRmZgIAMjMzUVhYOOB5ISEhmDZtGgBg5MiRiI6OhsVi8UXNfnHxIg+5sWNFV6IO\nss0OaRngBlrGIQDg0bWFVqsVRqMRAGA0GmG1Wh1+f11dHU6cOIGZM2fa/frGjRt7/2w2m2E2mz0p\ny6dsMwu9t13ayDY7rKoCZs8WXYU6yBb2HR1AUxMwcaLoStSlrKwMZWVlHj/f6wvJDQYDDA4Ssa2t\nDUuXLkV+fj5Gjhxp93v6hr1a0GFkf7IFRnU1sHq16CrUYdw4vmtpZ6ccrYy2bS6GDhVdibrcPBHO\nzc116/kenb4yGo345vo0r7GxEcHBwXa/79q1a1iyZAlWrlyJ9PR0T36VMNSJ05+MYU/jxw0ZAgQH\ny3NkRmOnDI/CPi0tDQUFBQCAgoICu0HOGMPq1asRExODdevWeVelADSz70+msL90ie94eX2lkUCu\n8aNOHGU4DfuMjAwkJSWhsrIS4eHh2LVrF3JycnDo0CFERkaitLQUOTk5AACLxYLU1FQAwEcffYS3\n3noLhw8fRmJiIhITEwd09agZhX1/MoUFnW8ZSMbxI77ldM1+z549dh8vKSkZ8FhYWBj2798PAJg1\naxZ6enq8LE8MarscSLawoGWA/mQbvxUrRFehPXTJiR2NjcCIEcDo0aIrUQ+ZwoKWAQaSqXWWxk8Z\nFPZ20Kx+oOBgoLkZ6OoSXYlzNH4DydI6++23vHMoLEx0JdpDYW8HdeIMFBDAt01oahJdiXO0jDOQ\nLEdmZ8/ybcVpmwvfo39SO2hmaJ8Ms0PGaBnAHlnCnsZOORT2dlDY2ydDYFy4wC/GoW0u+pNh7AB6\n7ymJwt4OWsaxT4bAoCUc+4xG/kHY3S26Esdo/JRDYX+Tnh7g66/5uiHpT4awp2UA+4YN491lFy+K\nrsQxGj/lUNjfpKGBLwHcdpvoStRHhrCnZYDB0fjpG4X9TegwcnCyhAWNn31qH7+WFuDaNd7mS3yP\nwv4mX3xBM4vBqD0surqAL7+k8RuM2i+sOnmStrlQktdbHGvJzp3Ali3AX/8quhJ1UnPYX7kCPPEE\ncOedwNSpoqtRJzW3zh49CmRkAC+/LLoS7aKZPfhJ2Q0bgE2bgCNHgPvvF12ROtluTciY6Er6a2gA\nZs3iQb9/PxAYKLoidVLrh/VbbwFLlgD/+7/AypWiq9Eu3c/sOzqAzEzAYgHKy/mNHoh9t97K/7t0\niV9Nqwaffw78278BTz0FPP00LQE4EhrKZ9BqwRiQmwsUFACHDwOxsaIr0jZdh31TEw+KiAigpAQY\nPlx0Repnmx2qIez37QPWrAH+53+AxYtFV6N+aprZf/89H7uqKj7JonsPKE+3yzhnzgAPPACkpAC7\nd1PQu0oNgcEYsHUr8NOfAgcOUNC7Sg1jB/AN9VJS+FH14cMU9P6iy5n9Bx/cOBm0apXoauQiOjC6\nuoD//E++HHHsGF+nJ66xjR1j4pa7qquB1FT+Af3SS7ThmT/pLux37OAnY999F+hz717iIpFhb+u4\n6ekBPvyQ7jfgrhEj+Mnry5eBO+7w/+8/cgR4/HHgxReBJ5/0/+/XO918rvb0AM8+y2cTR45Q0HtK\nVNjf3HFDQe8ZUeP31lvA0qW844aCXgxdhH1HB7BsGT/0Ly8H7rlHdEXyEhEWn38OPPgg75p67TW+\ntz7xjL8vrGIM2LgR+K//4uvzKSn++92kP82/bajjxrf8HRbUceNb/rywijpu1EXTM/vTp6njxtf8\nFRbUcaNTs4l2AAAM8klEQVQMfx2ZUceN+mh2Zk8dN8rwR1hQx41y/DF+1HGjTpochh07gOXLeccN\nBb1vjR7Ndyb87jtlfv6VK8CiRUBNDe+4oaD3LaXD/sgRIDkZeOYZYPNmCno1cTgU2dnZMBqNiI+P\n732spaUFKSkpiIyMxIIFC9Da2mr3ucXFxYiKisLdd9+NzZs3+7bqQVDHjfIMBuUCgzpulKdk2FPH\njbo5DPusrCwUFxf3eywvLw8pKSmoqqrCvHnzkJeXN+B53d3dWLt2LYqLi3H69Gns2bMHZ86c8W3l\nN7F13Bw5Anz8MXXcKEmJwKCOG/9QYuyo40YODsM+OTkZQUFB/R4rKipCZmYmACAzMxOFhYUDnnf8\n+HFMmTIFkyZNQmBgIJYtW4Z9+/b5sOz+mpqAuXN5QHzwATB+vGK/isD3gbFvH/Dww8Dvfw+sX0+b\nmSnJ12P3/fd8qfTgQd5xQ5uZqZfb8yer1Qrj9VPrRqMRVqt1wPecP38e4eHhvX83mUz45JNPBv2Z\nGzdu7P2z2WyG2Y31l9OngR/8gG+NmptLQeEPvgoMxoBXXuEn0Q8cAGbM8P5nEseCgvhRcEcH38HU\nG83N/CRscDCf0dOtPJVVVlaGsrIyj5/v1cGywWCAwU662nvMkb5h7w5bx82WLfzwn/iHL8KeOm7E\nMBhuXCtx112e/xzquPG/myfCubm5bj3f7SEyGo345nqjdWNjI4Lt3DBywoQJaGho6P17Q0MDTCaT\nu7/Kob4dNxT0/uXthVXUcSOWtx/W1HEjJ7eHKS0tDQUFBQCAgoICpKenD/ie6dOno7q6GnV1dejs\n7MQ777yDtLQ076sFddyogTcXVlHHjXjejB913MjLYdhnZGQgKSkJlZWVCA8Px65du5CTk4NDhw4h\nMjISpaWlyMnJAQBYLBakpqYCAAICArB9+3YsXLgQMTExeOKJJxAdHe11sdRxow6ezgyp40YdPBk/\n6riRn4ExsXcUNRgMcKWEvnvc7NxJWx+IZLUCcXHAhQuuP4f2uFGPF18Erl4FfvMb176/7x43RUW0\n9YFauJqdNlKsttEeN+oyfjzfE72z0/n30h436uPOzJ72uNEO1Yd9SQlfl3/+eeCFF6i1Ug2GDOGB\nb6frtp+uLuBnP+NHYseOUWulWrga9tXVfNntwQd5IwS1VspN1WG/YwewYgV13KiRs8Cgjhv1ciXs\nqeNGe1Q5hNRxo36OAqNvx83771PHjdo4C3vquNEm1fVDdHTwWbzFwi+/HjdOdEXEnsEC4/PP+Yn0\np54Cnn6alt3UaPx4oKWFL7P17YhijF+F/sYbfH2etj7QFlWFfVMTkJbGr+yju0qpm70Lq6jjRg4B\nAcDYsfycy4QJ/DFbx01lJZ9khYSIrZH4nmqWcWwdNwsWUMeNDPpemEMdN/LpO359O27KyijotUoV\nYW/ruNm4kTpuZGFbxqGOGznZxo86bvRDFcs4to4bOhErj9BQ3mmzaBE/of7hh3QiViahocCf/8y3\nJn7xRToRqweqmNlTx418wsL40ht13MgpLAz4y1+o40ZPpNkugagLY8Dx48D999Oym4ysVqC9nW8/\nQuTkbnZS2BNCiIQ0uTcOIYQQ71DYE0KIDlDYE0KIDlDYE0KIDlDYE0KIDlDYE0KIDlDYE0KIDlDY\nE0KIDlDYE0KIDlDYE0KIDngc9vn5+YiPj0dcXBzy8/MHfP3ixYt4+OGHMW3aNMTFxeGNN97wpk6/\nKisrE13CAFST69RYF9XkGqpJOR6F/T/+8Q+8/vrr+PTTT3Hy5Em8//77qKmp6fc927dvR2JiIr74\n4guUlZVh/fr16Orq8knRSlPj4FJNrlNjXVSTa6gm5XgU9v/85z8xc+ZMDB8+HEOHDsWcOXPwl7/8\npd/3hIaG4sqVKwCAK1euYOzYsQgIUMX2+YQQojsehX1cXByOHj2KlpYWtLe3Y//+/Th37ly/73ny\nySfx1VdfISwsDAkJCXaXegghhPiHx1sc79y5E3/4wx8wYsQIxMbG4pZbbsHWrVt7v/7rX/8aFy9e\nxCuvvIKamhqkpKTg5MmTuP322/sXQJuhE0KIR9yJb4/XVbKzs5GdnQ0A2LBhAyZOnNjv68eOHcNz\nzz0HAJg8eTIiIiJQWVmJ6dOne1wsIYQQz3jcjdPU1AQAqK+vx3vvvYfly5f3+3pUVBRKSkoAAFar\nFZWVlbjrrru8KJUQQoinPF7GmT17NpqbmxEYGIitW7fioYcewh//+EcAwI9//GNcvHgRWVlZqK+v\nR09PD5599tkBHwiEEEL8hAly8OBBds8997ApU6awvLw8UWX0U19fz8xmM4uJiWGxsbEsPz9fdEm9\nurq62LRp09gPfvAD0aUwxhi7dOkSW7JkCYuKimLR0dHs448/Fl0S27RpE4uJiWFxcXEsIyODXb16\n1e81ZGVlseDgYBYXF9f7WHNzM5s/fz67++67WUpKCrt06ZIq6vrFL37BoqKi2NSpU9nixYtZa2ur\n8JpsXn75ZWYwGFhzc7Mqatq2bRuLiopisbGx7Je//KXwmj755BM2Y8YMNm3aNDZ9+nR2/Phxpz9H\nSNh3dXWxyZMns9raWtbZ2ckSEhLY6dOnRZTST2NjIztx4gRjjLFvv/2WRUZGqqIuxhj73e9+x5Yv\nX84WLVokuhTGGGOrVq1iO3bsYIwxdu3aNb8Hxc1qa2tZREREb8A//vjj7I033vB7HUeOHGEVFRX9\n3pjPPPMM27x5M2OMsby8PParX/1KFXX97W9/Y93d3Ywxxn71q1/5vS57NTHGJ10LFy5kkyZN8nvY\n26uptLSUzZ8/n3V2djLGGGtqahJe05w5c1hxcTFjjLEDBw4ws9ns9OcI2S7h+PHjmDJlCiZNmoTA\nwEAsW7YM+/btE1FKPyEhIZg2bRoAYOTIkYiOjobFYhFcFXDu3DkcOHAAa9asUcUJ7cuXL+Po0aO9\nJ+gDAgIwevRooTWNGjUKgYGBaG9vR1dXF9rb2zFhwgS/15GcnIygoKB+jxUVFSEzMxMAkJmZicLC\nQlXUlZKSgiFDeATMnDlzQPu0iJoA4Omnn8Zvf/tbv9ZiY6+m1157Dc8++ywCAwMBAOPHjxdeU2ho\nKC5fvgwAaG1tdem1LiTsz58/j/Dw8N6/m0wmnD9/XkQpg6qrq8OJEycwc+ZM0aXgqaeewpYtW3rf\nmKLV1tZi/PjxyMrKwr333osnn3wS7e3tQmsaM2YM1q9fj4kTJyIsLAx33HEH5s+fL7QmG6vVCqPR\nCAAwGo2wWq2CKxpo586dePTRR0WXgX379sFkMmHq1KmiS+lVXV2NI0eO4IEHHoDZbMZnn30muiTk\n5eX1vt6feeYZvPTSS06fIyQ91N5b39bWhqVLlyI/Px8jR44UWsv777+P4OBgJCYmqmJWDwBdXV2o\nqKjAT3/6U1RUVGDEiBHIy8sTWlNNTQ1eeeUV1NXVwWKxoK2tDbt37xZakz0Gg0F1r//f/OY3GDZs\nmPAGivb2dmzatAm5ubm9j6nhNd/V1YVLly6hvLwcW7ZsweOPPy66JKxevRrbtm1DfX09tm7d2nuU\n7YiQsJ8wYQIaGhp6/97Q0ACTySSilAGuXbuGJUuWYOXKlUhPTxddDo4dO4aioiJEREQgIyMDpaWl\nWLVqldCaTCYTTCYTZsyYAQBYunQpKioqhNb02WefISkpqXdbjsceewzHjh0TWpON0WjEN998AwBo\nbGxEcHCw4IpueOONN3DgwAFVfDDW1NSgrq4OCQkJiIiIwLlz53Dffff1tnmLYjKZ8NhjjwEAZsyY\ngSFDhqC5uVloTcePH8fixYsB8Pff8ePHnT5HSNhPnz4d1dXVqKurQ2dnJ9555x2kpaWJKKUfxhhW\nr16NmJgYrFu3TnQ5AIBNmzahoaEBtbW1ePvttzF37ly8+eabQmsKCQlBeHg4qqqqAAAlJSWIjY0V\nWlNUVBTKy8vR0dEBxhhKSkoQExMjtCabtLQ0FBQUAAAKCgpUMYkAgOLiYmzZsgX79u3D8OHDRZeD\n+Ph4WK1W1NbWora2FiaTCRUVFcI/HNPT01FaWgoAqKqqQmdnJ8aOHSu0pilTpuDvf/87AKC0tBSR\nkZHOn6TE2WNXHDhwgEVGRrLJkyezTZs2iSqjn6NHjzKDwcASEhLYtGnT2LRp09jBgwdFl9WrrKxM\nNd04X3zxBZs+fbqwtj17Nm/e3Nt6uWrVqt7uCX9atmwZCw0NZYGBgcxkMrGdO3ey5uZmNm/ePKGt\nlzfXtWPHDjZlyhQ2ceLE3tf6T37yEyE1DRs2rPffqq+IiAi/d+PYq6mzs5OtXLmSxcXFsXvvvZcd\nPnxYSE19X1Offvopu//++1lCQgJ74IEHWEVFhdOf4/FFVYQQQuShjvYOQgghiqKwJ4QQHaCwJ4QQ\nHaCwJ4QQHaCwJ4QQHaCwJ4QQHfh/VW914fVCxAQAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "x_on = np.hstack((np.zeros(len(a_block)), np.ones(len(a_block))))\n", "x_off = 1 - x_on\n", "X_over_part = np.column_stack((x_on, x_off, np.ones_like(x_on)))\n", "X_over = np.tile(X_over_part, (n_on_off, 1))\n", "show_design(X_over, 'over parametrized')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAH0AAAEICAYAAAByGPvjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEdtJREFUeJzt3X1MlfX/x/HXAQ6gAnIAORw53AoIiAOUpdQYEALJkuHK\nBBLxZlNzbmRqNUaCaQJubS5dizVt6KQop+IdJDOPkIV0gzl3XLAE5T5FCBDZ4eb9/aNf168jhIo3\n54LP+7GxyXWum8+5npwPR+CcS0FEBCYUM1MPgD1/HF1AHF1AHF1AHF1AHF1AHN2EEhIScPjw4ae6\nz5ycHKSlpY25jsVTPSID8PeJ/+OPPx4a9OzZs0/92AqF4qHryOKRPjg4OKH3/7iICKb8mdhDo1+/\nfh1RUVFQqVQICgrCqVOnAACXL1+GRqMxGvzx48cRHBwMABgeHkZeXh58fHzg5OSE5cuXo7OzEwDQ\n0NAAMzMzHDx4EB4eHli0aNGI4+p0Omi1WuTm5mLGjBnw8vJCUVGRdPuZM2cQGhqK6dOnw93dHTt2\n7JBu+6/9L1u2DBqNBvb29oiMjIRer5e2WbVqFTZu3IiEhATY2toiIiICbW1tyMjIgEqlQkBAAK5c\nuSKt39LSgtdeew3Ozs7w9vbGvn37AABlZWXIzc1FcXExbG1tERoaCgCIiopCVlYWXnrpJdjY2ODG\njRuIiorCgQMHAADBwcGwtbWVPszMzFBRUQEAqKqqwosvvgiVSoWQkBBcvHhRGkd9fT0iIyNhZ2eH\nuLg43Llz52FJARqDwWCgWbNmUW5uLg0MDNB3331Htra2VFtbS0REs2bNovLycmn9119/nfLz84mI\naO/evRQeHk7Nzc1kMBho/fr1lJKSQkRE9fX1pFAoKD09nfr6+qi/v3/EsS9cuEAWFha0ZcsWMhgM\ndPHiRZo2bRr9/vvvRESk0+no2rVrRER09epVUqvVdOLEiTH3/8UXX1Bvby8ZDAZ6++23KSQkRDpe\neno6OTk50a+//kr9/f308ssvk4eHBx0+fJiGh4cpKyuLoqOjiYhoaGiI5s2bRzt37qSBgQG6ceMG\neXt707fffktERDk5OZSWlmZ0fyIjI8nDw4P0ej0NDQ3RwMAARUVF0YEDB0bc94KCAgoICKCenh5q\namoiR0dHKi0tJSKi8vJycnR0pDt37hAR0cKFC6VzVFFRQba2tiOO/aAxo1dUVJCLi4vRspSUFMrJ\nySEioqysLFqzZg0REXV3d9O0adPo1q1bREQUEBBA58+fl7ZraWkhpVJJQ0NDUpT6+vr/PPY/0fv6\n+qRlb7zxBu3cuXPU9TMyMmjz5s1ERI+0/87OTlIoFNTd3U1ERKtWraJ169ZJt+/bt48CAwOlz69e\nvUr29vZERFRVVUXu7u5G+9u9ezetXr2aiIiys7NpxYoVRrdHRUVRdnb2iGUPRq+srCRnZ2eqq6sj\nIqK8vLwREePj46mwsJBu3rw54hylpqaOOPaDxpzeW1pa4ObmZrTMw8MDzc3NAICUlBQcO3YMBoMB\nx44dw/z586X1GxoasHTpUqhUKqhUKgQGBsLCwgLt7e3Svh7c94NUKhWmTJlidOyWlhYAf397iY6O\nhrOzM+zt7VFQUICOjg6j7f+9/+HhYbz//vvw8fHB9OnT4eXlBQBG06Gzs7P0b2tra6PPp0yZgt7e\nXgDAzZs30dLSIt03lUqF3Nxc/Pnnn2Pen4fd38bGRixfvhyHDh2Cj4+PdKxvvvnG6FiXLl1CW1ub\nNIYHz9HDjBl95syZaGxsNPq+ffPmTWi1WgBAYGAgPDw8UFpaiqKiIqSmpkrrubu7o6ysDJ2dndJH\nX18fNBqNtM7Dnmn+s82/j+3q6goASE1NRVJSEpqamtDV1YUNGzZgeHjYaPt/7//IkSM4efIkzp8/\nj7/++gv19fUAMK4nVG5ubvDy8jK6b93d3Th9+jQAwMxs9NM61v29f/8+kpKSsHnzZsTHx0vL3d3d\nkZaWZnSsnp4evPvuu9BoNKOeo4ed1zGjL1y4EFOnTsWePXswMDAAnU6H06dPIzk5WVonNTUVe/fu\nRWVlJZYtWyYt37BhAzIzM3Hr1i0AwO3bt3Hy5MkxBzOa7OxsDAwMoLKyEmfOnJGO0dvbC5VKBUtL\nS1RXV6OoqGjMO9vb2wsrKys4ODjg3r17yMzMNLr9ceK/8MILsLW1xZ49e3D//n0MDQ3h2rVr+Pnn\nnwEAarUaDQ0NI/Y51jHWrFmDgIAAbN261Wj5ihUrcOrUKZw7dw5DQ0Po7++HTqdDc3MzPDw8EBYW\nJp2j77//XvrCG8uY0ZVKJU6dOoXS0lLMmDEDmzZtwuHDh+Hn5yetk5KSgoqKCsTExMDBwUFanpGR\ngcTERMTFxcHOzg7h4eGorq6Wbn+U/0+6uLhApVJh5syZSEtLQ0FBgXTsTz/9FNu3b4ednR127tyJ\n5cuXG2374P5XrlwJDw8PuLq6IigoCOHh4UbrKBSKMT//9z7Nzc1x+vRpXLlyBd7e3pgxYwbWrVuH\n7u5uAJC+MB0dHREWFvZI97m4uBgnTpwwegZ/6dIlaLValJSUYPfu3XB2doa7uzs+/vhjaVYrKirC\n5cuX4eDggA8//BDp6ekPPa8KGs/89hzodDqkpaWhsbHR1EOZdGTxwxn2fMk6+qN8C2CPT7bTO3t2\nZP1IZ8/GhPstm1yn/Ik0YQrxSKf/+63Wgx/Z2dlPZflEI7voZWVl8Pf3h6+vL/Lz8009nElJVtGH\nhoawadMmlJWVQa/X48svv8T169dNPaxJR1bRq6ur4ePjA09PTyiVSiQnJ6OkpOSZHS8qKuqpLJ9o\nZBW9ubnZ6DdRWq1W+o3esyBqdFk9e5frM/MH6XQ66HQ6Uw9j3GQV3dXV1ehn7Y2NjdKvceUkKirK\n6FH/7z/VmghkNb2HhYWhrq4ODQ0NMBgMKC4uRmJioqmHNenI6pFuYWGB/fv3Iz4+HkNDQ1i7di0C\nAgJMPaxJZ8L97H083/ef9V1UKBQT6oc0spre2fPB0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE\n0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QXE0QUkq5c1ParHfTXJ474q\nZiK9WmU8+JEuII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uII4u\nII4uII4uII4uII4uII4uII4uII4uII4uII4uII4uINm9ls3T0xN2dnYwNzeHUqlEdXW1qYc06cgu\nukKhgE6ng4ODg6mHMmnJcnqf7K8aNTXZRVcoFFi0aBHCwsLw+eefm3o4k5LspvdLly5Bo9Hg9u3b\niI2Nhb+/PyIiIozWycnJkf794MVun4eJflVlWV92c8eOHbCxscGWLVukZeO5xOWzflMCvuzmE+jr\n60NPTw8A4N69ezh37hzmzp1r4lFNPrKa3tvb27F06VIAwODgIN58803ExcWZeFSTj6yn99Hw9P7k\nZDW9s+dDVtP7o+J3i3oy/EgXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcX\nEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcXEEcX\nEEcXEEcXEEcX0IR8oyG5vaXYRMOPdAFxdAFxdAFxdAFxdAFxdAFxdAFxdAFxdAFxdAFxdAFxdAFx\ndAFxdAGZJPqaNWugVquNrrl29+5dxMbGws/PD3Fxcejq6jLF0IRgkuirV69GWVmZ0bK8vDzExsai\ntrYWMTExyMvLM8XQxEAmUl9fT0FBQdLns2fPpra2NiIiam1tpdmzZ4+63XiGDOCxPsaz/4lENt/T\n29vboVarAQBqtRrt7e0mHtHkJcs/l1IoFGP+iRNfVfkJmWqKGW16b21tJSKilpYWnt6fIdlM74mJ\niSgsLAQAFBYWIikpycQjmsRM8ZWWnJxMGo2GlEolabVaOnjwIHV0dFBMTAz5+vpSbGwsdXZ2jrrt\neIYMfqQb4asqj2I8+59Ip1E20zt7fji6gDi6gDi6gDi6gDi6gDi6gDi6gDi6gDi6gDi6gDi6gDi6\ngGT5lzMPw28c9GT4kS4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4g\nji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4gji4g\nji4gji4gji6gCfnuUnK7hstEw490Acnmqso5OTnQarUIDQ1FaGjoiAvwsqdHNldVVigUeOedd1BT\nU4Oamhq88sorphiaEEwSPSIiAiqVasTyyf69VC5k9T193759CA4Oxtq1a9HV1WXq4Uxasnn2/tZb\nb2H79u0AgA8++ABbtmzBgQMHRl2Xr6r8hEx1vc8Hr6r8qLeNZ8jga60akc303traKv37+PHjRs/s\n2dNlkuk9JSUFFy9exJ07d+Dm5oYdO3ZAp9PhypUrUCgU8PLyQkFBgSmGJgS+qvIoxrP/iXQaZTO9\ns+eHowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuIowuI\nowuIowuIowto0kT/r5cZjeflR09zX3LE0Z/xvuRo0kRnj46jC2hCvqxJjibSaZTN69Mf1UQ6uXLF\n07uAOLqAZB/97t27iI2NhZ+fH+Li4nD06FH4+/vD19cX+fn50npqtRrW1tawsrKCjY0NfH194eDg\nAEdHR8ydOxc6nQ7Tp0+X3rJs/vz5UKvVcHJygq+vL4KDg1FTU4PGxkZER0fDy8sL5ubm0tucbdu2\nDdHR0ZgzZw6CgoIQGRlptB2AEcfYtWuXqU7b2Ez2HhiPaNu2bZSfn09ERLt37yZ7e3uqr68ng8FA\nwcHBpNfraXBwkCwsLKimpoZKSkrI1taW9Ho9ffbZZxQUFERBQUF04cIFWrJkibTfiooK+uSTT8jG\nxoaIiKqqqmjBggXU2tpKNTU1dOHCBVq8eDH5+fmRXq+XlhMRHT16lKZNm0Z6vV7ajohGHEOuZP9I\nP3nyJNLT0wEAc+fOhcFggKenJ5RKJZKTk1FSUoLq6moolUq4ubnh7NmzePXVV1FSUoL169ejr68P\ng4ODAIyfBEZERKCqqkp6a7MFCxagq6sLCoUCISEhAABzc3MEBASgpaUFLi4u0vLy8nJp+T/btbe3\njziGXMk+ent7O9RqNQDg/v37UkAA0Gq1aG5uRnNzMywsLLBo0SJ89dVX6OnpQXNzMwDAxcUFAwMD\nUCgU+OGHHxAcHIyEhATo9Xq0tbVBqVQa7a+pqQnA3/81rKysxNmzZ5Gfnw+9Xi+tV1dXh8bGRixY\nsMBou9GOIUey+C9bbGws2traRiz/6KOPjD43Mxv9a1ShUGDJkiU4cuQI4uPj8dNPP8HS0tLo9nnz\n5qGxsRFTp05FaWkpkpKS4ObmNuq+AMDPzw/e3t7Yvn07rKyskJSUhNraWvT29uKXX37Be++9Bxsb\nm4ceo7a2dlzn5FmSRfTy8vL/vE2tVqOtrQ0uLi6wsrKChcX/D7mxsRFarRaurq7o6OgAAHh5eaGl\npQX9/f0AgLa2NlhaWsLW1lbabvHixdi4cSMcHBxQV1cnLW9qaoKrqysGBgawatUqrFy5EklJSQCA\njRs3or29HStXrkRgYCC8vLxGbDfaMe7evQsHB4cnPENPl+yn98TERBQWFgIAfvvtN1haWqKhoQEG\ngwHFxcVITExEYGAgamtr0dDQgOjoaFy/fh1Lly5FVVUV7OzsYGFhgfb2dun7bXV1NYgICQkJ6Ozs\nBABUVVXB3t4ezs7OWLt2LTw9PZGRkSGtPzw8jG3btiEwMBBZWVk4dOiQ0XZqtXrUY8gtOAD5P3vv\n6OigmJgY8vX1pdjYWPr666/Jz8+PPD09afbs2UREtGvXLnJyciIrKyuytLQkNzc3mjVrFk2ZMoVs\nbGxIqVSSvb09zZw5k7RaLXl7e1NsbCxpNBoyMzMjc3Nz0mq1lJmZSVu3biWFQkGurq5kbW1N1tbW\n5O/vT6mpqaRQKCg4OJhCQkLIycmJNBoNubq6UmZmJhER7d+/n+bMmUPBwcEUHh5OP/74oylP3X+a\ncD97Z09O9tM7e/o4uoA4uoA4uoA4uoA4uoD+B94/HAXVmqseAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [ "X_well_part = np.column_stack((x_on, np.ones_like(x_on)))\n", "X_well = np.tile(X_well_part, (n_on_off, 1))\n", "show_design(X_well, 'well parametrized')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAHoAAAEICAYAAACQxOCaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEIRJREFUeJzt3XtMU/cbx/FPkeKcFArCCtJykyITHUPNnFHEifUWQY3G\naabOS2Rx0czoLtHpgMwpuCVL1Cy6i6BxEeYSxSwbY8usZlOGi6ghaPCGQrnIuGyIMKA8vz8Wz0+0\n3JzQlud5JSZwetpzefd8ORbaoyIighjwXOy9AqJ/SGgmJDQTEpoJCc2EhGaiX0KvXLkS27dvBwCY\nzWYYDIb+WKxD2rVrF9auXftUH7Mn+9T1qS6xEyqVCiqVqj8WZTdmsxnLly9HaWlpl/Nt2bKln9ao\no34buvv7dZm2trZ+XV5PWK1Wuy27y9Dp6elISEhQvjcajVi8eLHyvcFgwOXLlwEAV69ehclkwrBh\nwxAREYFjx4492Qq5uGDv3r0YMWIEfH198e677ypPkhs3bmDatGnw8fGBr68vli1bhr/++ku5b3Bw\nMHbv3o0XXngBGo0GVqsVqampCAsLg4eHByIjI3HixAll/oyMDEyaNAmbNm2Cl5cXwsLCcPbsWaSn\npyMwMBA6nQ6HDx9W5v/nn3/w9ttvIygoCH5+fli3bh2am5vR2NiI2bNno7y8HBqNBh4eHqioqEBy\ncjIWLVqE5cuXw9PTExkZGUhOTsby5csBAOvXr4dGo1H+qdVqpKSkAADKy8uxcOFCPPfccwgNDcXe\nvXuV9WhqasLKlSvh7e2NyMhInD9/vvsdS124efMmabVaIiKyWCwUFBREBoOBiIhu3LhBXl5eRER0\n79490uv1lJGRQVarlQoKCsjHx4eKioqIiGjlypW0bds2IiI6deoU6fX6TpepUqlo2rRpVFdXR3fu\n3KHw8HD68ssviYjo+vXr9PPPP1NLSwtVV1fTlClTaOPGjcp9g4KCKDo6msrKyqi5uZmIiI4dO0YV\nFRVERJSVlUVDhw6lyspKIiJKT08nV1dXysjIoPb2dtq2bRsFBATQ+vXrqaWlhXJzc0mj0VBjYyMR\nEW3cuJHmzZtHdXV11NDQQPHx8bRlyxYiIjKbzY9tV1JSEqnVasrOziYioqamJkpOTqZly5Y9tt0F\nBQXk6+tLFy9eJKvVSmPHjqUPP/yQWltb6ebNmxQaGko//vgjERG99957NGXKFKqrq6PS0lKKjIxU\nunSmy9BERAaDgS5cuEBHjx6lxMREmjBhAl29epUOHjxI8+bNIyKizMxMiomJ6XC/xMRESklJIaLe\nh36wQUREn332GcXFxdmc9/jx4xQdHa18HxwcTOnp6V1uz4svvqjs+PT0dDIajcptly9fJpVKRXfv\n3lWmDRs2jC5dukTt7e00dOhQunHjhnLb2bNnKSQkpNPtSkpKotjY2MemPRr67t27FBQURFlZWURE\nlJeXR4GBgR3m2blzJ61atYqIqEN0IqLPP/+8y31KRNTtyVhsbCzMZjOuX7+O2NhYaLVanD59GufO\nnUNsbCwA4Pbt2/j999/h5eWl3K+trQ0rVqzofkix4eEzyMDAQJSXlwMAqqqq8NZbb+HXX39FQ0MD\n2tvb4e3t3el9AeDw4cP49NNPUVJSAgC4d+8eampqlNt1Op3y9ZAhQwAAvr6+Habdu3cP1dXVuH//\nPsaNG6fcRkRob2/vclv0en2Xt7e2tmLRokVYtmyZ8mPx9u3bKC8v77A/rVYrpkyZAuDfYf3RfdSd\nHoU+efIkSkpK8P7770Or1eLIkSPIy8vDhg0blAXFxsYiNze308fpzVn3nTt38PzzzytfBwQEAAC2\nbt2KQYMGobCwEFqtFidOnFDWwdZybt++jcTERPzyyy+YOHEiVCoVoqOjn+jE0MfHB0OGDEFRURH8\n/f17tH22/rfx6PcbNmyAVqvFjh07lGmBgYEICQlBcXGxzXXx9/d/bB91p9uz7tjYWJw6dQrNzc0Y\nPnw4Jk+ejJycHNTW1iI6OhoAMHfuXBQXF+PIkSNobW1Fa2srzp8/j6tXrwL495nfm537ySefoL6+\nHqWlpdizZw9effVVAP8ejUOHDoWHhwcsFgs+/vjjLh+nsbERKpUKPj4+aG9vR3p6OgoLC3u8Hg9z\ncXHB2rVrsXHjRlRXVwMALBaL8uTW6XSoqanB33//rdzH1jY/PO3AgQM4c+YMjhw50mGel156CRqN\nBrt370ZTUxOsVisKCwvxxx9/AAAWL16MXbt2ob6+HmVlZR1O1Dpd/+5mMBqN0Gg0iImJAQB4eHhg\nxIgRmDRpkvLsdHd3R25uLjIzMxEQEAB/f39s2bIFLS0tAB5/Znd3dM+bNw/jxo1DdHQ05s6di9Wr\nVwMAkpKScOHCBXh6eiI+Ph4LFy7s8rFGjRqFzZs3Y+LEifDz80NhYSEmT57cYT26O+IelpaWhrCw\nMLz88svw9PSEyWRSjrqIiAgsXboUoaGh8Pb2RkVFRaeP/2BaZmYmbt26heHDhytn3qmpqXBxccF3\n332HixcvIjQ0FL6+vkhMTFSeRElJSQgKCkJISAhmzZqFFStWdLtPVfQk41gfcnFxwfXr1xEaGmrv\nVRlQ5LVuJhwu9EB/qdReHG7oFn3D4Y5o0Tf65bdXfa2/hntnHvzYHtEP/m//8L+kpCSb05058ANO\nETonJwcREREwGo1IS0uz9+o4JYcPbbVasX79euTk5KCoqAhHjx7FlStX7L1aTsfhQ+fn5yMsLAzB\nwcFQq9VYsmQJsrOz+2RZU6dO7ZPHdQQOH9pisXT4TY1er4fFYumTZQ3k0A5/1m2vF1DMZjPMZrNd\nlt0XHD50QEBAhz+4Ky0t7fZ3vE/D1KlTOxzhD/7Ex1k5/NA9fvx4XLt2DSUlJWhpaUFWVlaHv2MT\nPePwR7Srqyv27duHmTNnwmq1Ys2aNcov3EXPDYjXup/k53hvN1ulUjn1CycOP3SLp0NCMyGhmZDQ\nTEhoJiQ0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJCc2EhGZCQjMhoZmQ0ExIaCYkNBMSmgkJ\nzYSEZsLh35LTU878Lor+IEc0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJCc2EhGZCQjMhoZmQ\n0ExIaCYkNBMSmgkJzYSEZkJCMyGhmZDQTEhoJiQ0ExKaCQnNhFO89yo4OBgeHh4YNGgQ1Go18vPz\n7b1KTscpQqtUKpjNZnh7e9t7VZyW0wzd8m7J/8YpQqtUKkyfPh3jx4/HF198Ye/VcUpOMXT/9ttv\n8Pf3R3V1NUwmEyIiIhATE9NhnuTkZOXrRy8g+iQG2tVmne6ShSkpKXB3d8fmzZuVaf1xOUG5ZGEf\nu3//PhoaGgAAjY2NyM3NxZgxY+y8Vs7H4YfuqqoqLFiwAADQ1taG1157DTNmzLDzWjkfpxu6bZGh\nu3sOP3SLp8Phh+6e6u01pJ356HwSckQzIaGZkNBMSGgmJDQTEpoJCc2EhGZCQjMhoZmQ0ExIaCYk\nNBMSmgkJzYSEZkJCMyGhmZDQTEhoJiQ0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJCc2EhGZC\nQjMhoZmQ0ExIaCYGzIfVcPvwmd6SI5oJCc2EhGZCQjMhoZmQ0ExIaCYkNBMSmgkJzYSEZkJCMyGh\nmZDQTDhM6NWrV0On03W4plVtbS1MJhPCw8MxY8YM1NfX23ENnZvDhF61ahVycnI6TEtNTYXJZEJx\ncTHi4uKQmppqp7UbAMiB3Lp1i0aPHq18P3LkSKqsrCQiooqKCho5cqTN+/XHZjjYruo1hzmibamq\nqoJOpwMA6HQ6VFVV2XmNnJfT/CmRSqXq8tpWcrXZbth7SHmYraG7oqKCiIjKy8tl6P4PHHroTkhI\nwKFDhwAAhw4dwvz58+28Rk7M3s+0B5YsWUL+/v6kVqtJr9fTwYMHqaamhuLi4shoNJLJZKK6ujqb\n9+2PzXCgXfVE5GqzDrSMvuTQQ7d4eiQ0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJCc2EhGbC\naf7CpDtd/fWJLc78m6gnIUc0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJCc2EhGZCQjMhoZmQ\n0ExIaCYkNBMSmgkJzYSEZkJCMyGhmZDQTEhoJiQ0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJ\nCc2EhGZCQjMxYD6ViNunDPWWHNFMOExoW1ebTU5Ohl6vR3R0NKKjox+7SKnoOYcJbetqsyqVCps2\nbUJBQQEKCgowa9YsO62d83OY0DExMfDy8npsuvzsfTocJnRn9u7di6ioKKxZs0YuFP4fOPRZ97p1\n6/DBBx8AALZv347Nmzfjq6++sjmvXG22G3a8XOJjHr3abE9v64/NcLBd1WsOPXRXVFQoXx8/frzD\nGbnoHYcZupcuXYrTp0/jzz//hMFgQEpKCsxmMy5evAiVSoWQkBAcOHDA3qvptORqsw60jL7k0EO3\neHokNBMSmgkJzYSEZkJCMyGhmZDQTEhoJiQ0ExKaCQnNhIRmQkIzIaGZkNBMSGgmJDQTEpoJCc2E\nhGZCQjMxoEN39paa3k4fCCR0D6YPBAM6tPg/Cc3EgHlLTn9w5l3lMG+y+y+cOUB/kaGbCQnNhFOG\nrq2thclkQnh4OGbMmIFvv/0WERERMBqNSEtLU+bT6XR45plnMHjwYLi7u8NoNMLb2xvDhg3DmDFj\nYDab4enpqXy81bhx46DT6eDj4wOj0YioqCgUFBQoj/fo/Dt27LDH5j8Z+33YwpN75513KC0tjYiI\ndu7cSVqtlm7dukUtLS0UFRVFRUVF1NbWRq6urlRQUEDZ2dmk0WioqKiI9u/fT6NHj6bRo0fTqVOn\nKD4+XnncM2fO0J49e8jd3Z2IiPLy8mjChAnK7Y/O70yc8og+efIkXn/9dQDAmDFj0NLSguDgYKjV\naixZsgTZ2dnIz8+HWq2GwWDA999/j7lz5yI7OxtvvPEG7t+/j7a2NgAdT+RiYmKQl5enfAzWhAkT\nUF9fj6qqKmUectITP6cMXVVVBZ1OBwBoampSogGAXq+HxWKBxWKBq6srpk+fjszMTDQ0NMBisQAA\n/Pz80NraCpVKhbNnzyIqKgpz5sxBUVERKisroVarOzxeWVkZANic31k47H+vTCYTKisrH5v+0Ucf\ndfjexcX2c1WlUiE+Ph5ff/01Zs6cifPnz8PNza3D7WPHjkVpaSmeffZZ/PDDD5g/fz4MBoPNxwJg\nc/7i4uL/spn9xmFD//TTT53eptPpUFlZCT8/PwwePBiurv/fjNLSUuj1egQEBKCmpgYAEBISgvLy\ncjQ3NwMAKisr4ebmBo1Go9xv9uzZePPNN+Ht7Y1r164p08vKyhAQEAAANuevra2Ft7f309noPuSU\nQ3dCQgIOHToEALh06RLc3NxQUlKClpYWZGVlISEhAaNGjUJxcTFKSkrwyiuv4MqVK1iwYAHy8vLg\n4eEBV1dXVFVVKT9z8/PzQUSYM2cO6urqAAB5eXnQarXKjwlb8ztDZADOedZdU1NDcXFxZDQayWQy\n0TfffEPh4eEUHBxMI0eOJCKiHTt2kI+PDw0ePJjc3NzIYDDQiBEjaMiQIeTu7k5qtZq0Wi0NHz6c\n9Ho9hYaGkslkIn9/f3JxcaFBgwaRXq+nrVu30v79+4mIaN++fRQZGUlRUVE0ceJEOnfunD13Q68M\niNe6RfeccugWvSehmZDQTEhoJiQ0ExKaif8BZk60XOkWeokAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "X_times2 = X_well.copy()\n", "X_times2[:, 0] *= 2\n", "show_design(X_times2, 'times 2')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAD8AAAEICAYAAADoVyx+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACutJREFUeJztnW1sU1UYx/93bTcgvXvLtGvaYjvXvelSp5hFI2HICgwd\nLyaaIUYyiB80RhMJHzCKzEzDouEL+6JGA4QpLqLOCC4j6DS+1JGsmmBrlugmXdkWHA6Gb2vr4wez\nK4Ot9J3B8/ySJr33Pjv3+e2cnnvbe889ChERmJJ1tRO4mog8V0SeKyKfCk6dOgVVVXEtHTkTlrfb\n7fj000+15cWLF2NychKKoqQksUQ4cuQI7rnnHhQUFMBsNuOxxx7DhQsX5oxPWF5RlHlXy+fPn8fO\nnTsxMjICv9+PYDCI7du3z/0HlACPPPIIZWVl0cKFC8loNNIrr7xCg4ODpCgKRSIRIiJatmwZPffc\nc3T33XeT0WikxsZGOnPmDD388MOUm5tLd955Jw0NDWll+v1+qq+vp8LCQiovL6fOzk5t25EjR6iq\nqopUVSWLxUKvvvpqTHm+//77VF1dPef2hOSJiOx2Ox0/flxbnk3e6XTSzz//TOfOnaOqqioqLS2l\n48ePUzgcpkcffZSam5uJiOjChQtktVpp3759FIlEyOv1UlFREfn9fiIiKi4upi+//JKIiCYmJqi/\nvz+mHJ9++mnauHHjnNvT1tsrioLm5mY4HA7k5uaioaEBZWVluPfee6HT6fDggw/C6/UCAD7++GM4\nHA5s3rwZWVlZuO222/DAAw+gs7MTAJCdnY0ffvgB58+fR15eHmpqaq64/2PHjuHAgQN48cUX54xJ\n66HOZDJp7xcsWIAbb7xxxvJ0Z/TLL7/g22+/RUFBgfZ6++23MTY2BgA4fPgwjh49Crvdjrq6Ong8\nnqj79Xg82LRpEw4fPozS0tI54/SJisXbq0eLX7x4MZYtW4aenp5Zty9ZsgQffvghIpEI9u7di4ce\neginTp2aNdbr9WLdunXYt28fli9fHjWnhGveZDLhp59+ihpDFx0NKMqR4b777sPAwAAOHjyIUCiE\nUCiEEydO4Mcff0QoFEJHRwfOnTsHnU4HVVWh0+lmLefkyZNYvXo12tvbsWbNmis6JCy/Y8cOtLa2\noqCgAHv27AFwee1evKwoypzbVVVFT08PDh06BIvFArPZjB07dmBqagoAcPDgQTgcDuTl5eH1119H\nR0fHrDnt2bMH4+Pj2LJlC1RVhaqqqK6untNBoWhVcp0j5/ZcEXmuJHycvxKZ+naXTH89r2qe/vuu\ncdnrhRdemHV9siQs393djYqKCjidTrS1tSWdyNUgIflIJIInn3wS3d3d8Pl8eOedd+D3+1OdW9pJ\nSL6vrw+lpaWw2+0wGAxoampCV1dXqnPTqKurS0u5CckHg0HYbDZt2Wq1IhgMpiypS0mXfEK9/dX6\nna63txe9vb0pKy8heYvFgkAgoC0HAgFYrdaUJTUXdXV1M1pBS0tLcgXG9HvQJYRCISopKaHBwUH6\n+++/yeVykc/nmxEDIO5XvCSYvkZCNa/X69He3o5Vq1YhEolg69atqKysTK4WrgJp+0qbSL8QbyrJ\n/nw+r87wMo3Ic0XkuSLyXBF5rog8V0SeKyLPFZHnishzReS5IvJcEXmupO2eHCC5+2UyAeuaF3mu\niDxXRJ4rIs8VkeeKyHNF5Lki8lwRea6IPFdEnisizxWR50rC1+rsdjtyc3Oh0+lgMBjQ19eXyrwy\nQlIPCert7UVhYWEq88koSTX7+X4V9kok9Ty8+vp6LFmyBG+88UYqc8oYCTf7r776CmazGWfOnIHb\n7UZFRQWWLl06I2bXrl3a+0sHASdCqkdRp2Q4aUtLC4xGI7Zt2/Z/wRl4UuJVGU76xx9/YHJyEgDw\n+++/o6enJ+ozqOYrCTX7sbExbNiwAQAQDoexadMmrFy5MqWJZYK0jqK+Lpv99UJa78aKdwx9ps8b\nWNe8yHNF5Lki8lwRea6IPFdEnisizxWR54rIc0XkuSLyXBF5rog8V0SeKyLPFZHnishzhbW8PB6K\nKyLPFZHnishzReS5IvJcEXmuiPxcbNmyBSaTacaYubNnz8LtdqOsrAwrV67ExMRE2pNMF1Hlm5ub\n0d3dPWPd7t274Xa7MTAwgBUrVmD37t1pTTCtXGnu1sHBQbr11lu15fLychodHSUiopGRESovL5/1\n72IoOmmS3Ufcn/mxsTGYTCYAgMlkwtjYWIqrI3Mk9TOWoihRx87N91HUCTX7kZERIiI6ffo0r2a/\ndu1a7N+/HwCwf/9+rF+/PnU1kWmi/WeamprIbDaTwWAgq9VKb731Fo2Pj9OKFSvI6XSS2+2m3377\nLS21EgvJ7kNGUXNF5Lki8lwRea6IPFdEnisizxV5QhJXRJ4rIs8VkeeKyHNF5Lki8lwRea6IPFdE\nnisizxWR54rIc0XkuSLyXBF5rog8V0SeKyLPFXk8FFfiHkW9a9cuWK1W1NTUoKam5rKBxtcScY+i\nVhQFzzzzDLxeL7xeL1avXp3WBNNJVPmlS5eioKDgsvXz/bMcKwl95vfu3QuXy4WtW7de0w8PiLu3\nf/zxx7Fz504AwPPPP49t27bhzTffnDX2uhtFHeu2GIpOmmT3EXezHxkZ0d5/8MEH1+Qc1NNEbfYb\nN27E559/jl9//RU2mw0tLS3o7e3Fd999B0VR4HA48Nprr2Uq15Qjo6i5IvJcEXmuiDxXRJ4rIs8V\nkeeKyHNF5Lki8lwRea6IfCaZ63JTtMtQKb1EdREizxXW8mm9XJUJkkk/bXdjpfs6XSpg3exFPhVc\nOtHHe++9h4qKCjidTrS1tWlxJpMJCxYsQE5ODoxGI5xOJ1wuF9atWweTyYSSkhLk5eVpt7rdcccd\n2u1wTz31lBbv9XoB/HcYvDi+tbU19qSTvDNEY/v27dTW1kZERC+//DLl5+fT4OAgTU1NkcvlIp/P\nR+FwmPR6PXm9Xurq6iJVVcnn85HH46HKykrq7+8nh8NBjY2NWrlffPEF9ff300033UQNDQ1EROTx\neKi2tpaIiD777LMZ8fGQspr/6KOPsHnzZgBAdXU1pqamYLfbYTAY0NTUhK6uLvT19cFgMMBms+Ho\n0aO4//770dXVhdraWvzzzz+IRCLTFaKVO3073OTkpFZ+bW0tJiYmtPk0KMHONWXyF0/08eeffyIc\nDmvbrFYrgsEggsEg9Ho96uvrcejQIUxOTiIYDGox0/f7fP3113C5XFizZg18Ph8AIBQKwWazzShz\neHgYiqLMGh8LcR3q3G43RkdHL1v/0ksvzVjOypr9f6ooChobG9HR0YFVq1bhxIkTyM7OnrF94cKF\nOHnyJBYtWoRPPvkE69evR09PD4DLa1hRFNx+++0IBAIz4gcGBmLyiUv+2LFjc24zmUwYHR1FcXEx\ncnJyoNf/X3QgEIDVaoXFYsH4+DgAwOFw4PTp0/jrr78AAMPDwyguLkZWVhYWLVoEAGhoaMATTzyB\niYkJGAwGBAIBrczh4WFYLBaoqqqtm44/e/YsCgsLr+iTsmZ/8UQf33//PbKzszE0NISpqSm8++67\nWLt2LaqqqjAwMIChoSEsX74cfr8fGzZsgMfjQX5+PoqKihAOh7Ua7uvrAxEhPz8fqqriwIEDAKDF\nT8+jc2l8LOIAUtfbXzrRR2dnJ5WVlZHdbteme2ltbaWioiLKycmh7OxsstlsdPPNN5PFYqGqqioy\nm82k0+lIr9eTzWajkpISqqys1CYVMRqNdMMNN5DFYqFnn32WiIja29vplltuIZfLRXfddRd98803\nMeectnP7awE5w+OKyHNF5LnyL7ZCgOxqFn0AAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "B_over = npl.pinv(X_over).dot(y)\n", "B_over" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": [ "array([ 4., 3., 7.])" ] } ], "prompt_number": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "# Print betas, t statistic, df, p value for different designs, contrasts\n", "print(t_stat(y, X_over, [1, 0, 0]))\n", "print(t_stat(y, X_well, [1, 0]))\n", "print(t_stat(y, X_times2, [1, 0]))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(array([ 4., 3., 7.]), array([[ 185.90320062]]), 16, array([[ 0.]]))\n", "(array([ 1., 10.]), array([[ 24.49489743]]), 16, array([[ 2.05391260e-14]]))\n", "(array([ 0.5, 10. ]), array([[ 24.49489743]]), 16, array([[ 2.05391260e-14]]))\n" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The t statistic on the over-parametrized design does not relflect the estimability of the contrast:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def is_estimable(X, C):\n", " \"\"\" Is a contrast C estimable on a design X\n", " \n", " To be estimable, the contrast needs to be orthogonal to the null-space of X.\n", " The null space of X is all the vectors $k$ such that Xk = 0.\n", " If $k$ is a vector in the null-space of X, then:\n", " \n", " Y = X B + e\n", " \n", " and\n", " \n", " Y = X (B + k) + e\n", " \n", " give the same fit (X B == X (B + k)).\n", " \n", " The null space of X is the orthogonal complement of the row-space of X. So C\n", " has to be in the row space of X.\n", " \n", " A less mathematical way of seeing this is that the information we will need to form\n", " our betas is made of linear combinations of the rows of X, because the betas are\n", " constructed from Xt Y. \n", " \"\"\"\n", " C = np.atleast_2d(C)\n", " rankX = npl.matrix_rank(X)\n", " return rankX == npl.matrix_rank(np.vstack((C, X)))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "is_estimable(X_over, [1, 0, 0])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 15, "text": [ "False" ] } ], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "is_estimable(X_over, [1, -1, 0])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "True" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our case, the null space of $\\mathbf{X}$ is:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import sympy\n", "X = sympy.Matrix(X_over)\n", "sympy.pretty_print(X.nullspace())" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\u23a1\u23a1-1.0\u23a4\u23a4\n", "\u23a2\u23a2 \u23a5\u23a5\n", "\u23a2\u23a2-1.0\u23a5\u23a5\n", "\u23a2\u23a2 \u23a5\u23a5\n", "\u23a3\u23a3 1 \u23a6\u23a6\n" ] } ], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Any contrast vector not orthogonal to this null space is not estimable. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way of seeing this is to get a minimal set of vectors that \"span\" the rows of X. These can be found using the svd :" ] }, { "cell_type": "code", "collapsed": false, "input": [ "[u, s, vt] = npl.svd(X_over, full_matrices=False)\n", "tol = s.max() * max(X_over.shape) * np.finfo(s.dtype).eps\n", "nz = np.where(s > tol)[0]\n", "print nz\n", "vt[np.abs(vt) < tol] = 0\n", "print vt.T[:,nz]\n", "print \"\\nNormalize the columns \\n\"\n", "print vt.T[:,nz].dot(np.diagflat([1./vt[0,0], 1./vt[1,1]]))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[0 1]\n", "[[-0.40824829 -0.70710678]\n", " [-0.40824829 0.70710678]\n", " [-0.81649658 0. ]]\n", "\n", "Normalize the columns \n", "\n", "[[ 1. -1.]\n", " [ 1. 1.]\n", " [ 2. 0.]]\n" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, any contrast that is made of these two is estimable. In particular, what do you think the first column would test for?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare to an F test:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def f_stat(Y, X, col):\n", " \"\"\" betas, F statistic and significance test given data, design and column to test\n", " \n", " We do the F statistic long hand, by fitting the model with and without the column to test.\n", " \n", " This is OLS estimation; we assume the errors to have independent\n", " and identical normal distributions around zero for each $i$ in \n", " $\\Epsilon_i$ (i.i.d).\n", " \"\"\"\n", " Y = np.asarray(Y)\n", " full_X = np.asarray(X)\n", " # Delete column to make reduced design\n", " reduced_X = np.delete(full_X, col, 1)\n", " # fit both designs\n", " _, _, full_resid = fit_ols(Y, full_X)\n", " _, _, reduced_resid = fit_ols(Y, reduced_X)\n", " n_obs = Y.shape[0]\n", " # Number of parameters used by each design\n", " n_p_full = npl.matrix_rank(full_X)\n", " n_p_reduced = npl.matrix_rank(reduced_X)\n", " # Extra sum of squares\n", " full_SS = (full_resid ** 2).sum()\n", " reduced_SS = (reduced_resid ** 2).sum()\n", " extra_SS = (reduced_SS - full_SS) / (n_p_full - n_p_reduced)\n", " # F statistic\n", " F = extra_SS / (full_SS / (n_obs - n_p_full))\n", " ltp = f_dist((n_p_full - n_p_reduced), n_obs - n_p_full).cdf(F) # lower tail p\n", " p = 1 - ltp # upper tail p\n", " return F, p" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 19 }, { "cell_type": "code", "collapsed": false, "input": [ "f_stat(y, X_well, 0)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 20, "text": [ "(600.00000000000409, 4.1189274213593308e-14)" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "What happens with the over-parametrized design?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "f_stat(y, X_over, 0)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stderr", "text": [ "-c:24: RuntimeWarning: divide by zero encountered in double_scalars\n" ] }, { "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "(inf, nan)" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Was that what you were expecting?" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Correlated regressors" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Remind ourselves about the t statistic on the over-parametrized design\n", "print(t_stat(y, X_over, [0, 0, 1]))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(array([ 4., 3., 7.]), array([[ 514.39284598]]), 16, array([[ 0.]]))\n" ] } ], "prompt_number": 22 }, { "cell_type": "code", "collapsed": false, "input": [ "# Add a tiny bit of noise to the over-parametrized design\n", "X_over_tweaked = X_over.copy()\n", "X_over_tweaked[0, 2] = 1e-14\n", "print(t_stat(y, X_over_tweaked, [0, 0, 1]))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(array([ 10.8875, 9.9 , 0.1125]), array([[ 1.2456822]]), 15, array([[ 0.11599406]]))\n" ] } ], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The general case of correlated regressors" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# first a quick function to yield two correlated regressor with correlation c\n", "def correlated(c, n=20, seed=None):\n", " np.random.seed(seed) # To get predictable random numbers\n", " # generate two uncorrelated regressors\n", " Y = np.random.normal(0,1,(n,2))\n", " #make them correlated with mixing matrix M:\n", " #we want a correlation c between Y[0] and Y[1], with unit variance\n", " #this is the mixing matrix needed:\n", " c1 = np.sqrt( (1 + np.sqrt(1-c*c))/2 ) \n", " c2 = .5 * c / c1\n", " M = np.asarray([[c1, c2],[c2, c1]])\n", " Yc = Y.dot(M) # Make Yc = MY\n", " return Yc" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 24 }, { "cell_type": "code", "collapsed": false, "input": [ "n = 20\n", "Xc = correlated(.8, n)\n", "x1c = Xc[:,0]; x2c = Xc[:,1]; \n", "print np.corrcoef(Xc.T)[0,1]\n", "\n", "y_corr = x1c + 1.5*np.random.normal(size=n) + 10\n", "#x1c = np.random.normal(size=n)\n", "#x2c = np.random.normal(size=n)\n", "X_corr = np.column_stack((x1c, x2c, np.ones_like(x2c)))\n", "\n", "C = np.array([1, 0, 0]) # The contrast\n", "t_stat(y_corr, X_corr, C)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.830592959914\n" ] }, { "metadata": {}, "output_type": "pyout", "prompt_number": 25, "text": [ "(array([ 1.01047766, 0.14863412, 9.42780916]),\n", " array([[ 1.9797293]]),\n", " 17,\n", " array([[ 0.03207994]]))" ] } ], "prompt_number": 25 }, { "cell_type": "code", "collapsed": false, "input": [ "XtX = X_corr.T.dot(X_corr)\n", "XtX" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "array([[ 18.34278453, 14.98641369, -2.36167903],\n", " [ 14.98641369, 17.77134663, -4.45422709],\n", " [ -2.36167903, -4.45422709, 20. ]])" ] } ], "prompt_number": 26 }, { "cell_type": "code", "collapsed": false, "input": [ "npl.pinv(XtX)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 27, "text": [ "array([[ 0.17851102, -0.15384086, -0.01318282],\n", " [-0.15384086, 0.19217721, 0.02463391],\n", " [-0.01318282, 0.02463391, 0.05392957]])" ] } ], "prompt_number": 27 }, { "cell_type": "code", "collapsed": false, "input": [ "C.dot(npl.pinv(XtX)).dot(C)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 28, "text": [ "0.17851102422647533" ] } ], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [ "#Xc = correlated(.0, n, seed=42)\n", "#y_corr = Xc[:,0] + 1.*np.random.normal(size=n) + 10\n", "\n", "\n", "n=20\n", "x1c = np.random.normal(size=n)\n", "x2c = np.random.normal(size=n)\n", "y_corr = x1c + 1.*np.random.normal(size=n) + 10\n", "# if you make the data exactly the same, for the \n", "# experiment to work, you would have to make the correlated regressor \n", "# different (use the correlated function with the seed not fixed). \n", "# Otherwise, you are not going to change your t since \n", "# the amount of the first regressor that is in the data is just going\n", "# to be scaled, and the scaling doesnt change the t-value. In other\n", "# words you are not changing the amount of data variance that can be \n", "# explained by the part of x1 that is orthogonal to x2, since x1 and \n", "# x2 directions remains exactly the same.\n", "\n", " \n", "x2c_corr = np.linspace(0.1, 0.999, 15)\n", "\n", "dtype = np.dtype(dict(names=['t', 'mss', 'b0', 'b1', 'b2', 'xvar', 'df'],\n", " formats=['f8'] * 7))\n", "res = np.zeros(len(x2c_corr), dtype=dtype)\n", "\n", "for i, x2c_corr in enumerate(x2c_corr):\n", " new_x1c = (1 - x2c_corr) * x1c + x2c_corr * x2c\n", " #Xc = correlated(x2c_corr, n, seed=42)\n", " #new_x1c = Xc[:,0]; x2c = Xc[:,1];\n", " #print np.corrcoef(new_x1c,x2c)[0,1]\n", " X_corr = np.column_stack((new_x1c, x2c, np.ones_like(x2c)))\n", " y_corr = new_x1c + 1.*np.random.normal(size=n) + 10\n", " \n", " betas, res[i]['t'], res[i]['df'], p = t_stat(y_corr, X_corr, C)\n", " res[i]['b0'], res[i]['b1'], res[i]['b2'] = betas[:]\n", " fitted = X_corr.dot(betas)\n", " res[i]['mss'] = ((y_corr - fitted)**2).sum() / res[i]['df']\n", " res[i]['xvar'] = C.dot(npl.pinv(X_corr.T.dot(X_corr)).dot(C))\n", "print(mlab.rec2txt(res))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " t mss b0 b1 b2 xvar df\n", " 4.824 1.360 1.144 -0.018 9.782 0.041 17.000\n", " 5.043 0.683 0.913 0.149 10.290 0.048 17.000\n", " 4.343 0.901 0.978 -0.043 10.037 0.056 17.000\n", " 4.319 1.373 1.309 0.067 9.772 0.067 17.000\n", " 4.767 0.896 1.283 -0.565 10.077 0.081 17.000\n", " 4.142 0.828 1.192 -0.370 9.882 0.100 17.000\n", " 4.805 0.449 1.144 -0.208 10.143 0.126 17.000\n", " 2.846 1.016 1.165 -0.427 9.295 0.165 17.000\n", " 4.372 0.777 1.826 -0.910 10.194 0.224 17.000\n", " 0.404 1.138 0.245 0.992 10.018 0.323 17.000\n", " 0.674 1.115 0.505 0.055 9.971 0.504 17.000\n", " 0.133 1.243 0.141 0.791 10.225 0.893 17.000\n", " 1.310 0.863 1.720 -0.847 10.129 1.999 17.000\n", " -0.884 0.946 -2.411 3.074 9.777 7.873 17.000\n", " 0.773 1.399 167.393 -165.998 10.192 33482.838 17.000\n" ] } ], "prompt_number": 29 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "HRF and temporal derivatives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now consider a design convolved with an HRF:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def spm_hrf(t):\n", " \"\"\" Return SPM hrf sampled at times `t`\n", " \"\"\"\n", " # gamma.pdf only defined for t > 0\n", " hrf = np.zeros_like(t, dtype=np.float)\n", " hrf[t > 0] = gamma.pdf(t[t > 0], 6, 0, 1) - gamma.pdf(t[t > 0], 16, 0, 1) / 6.\n", " return hrf / np.sum(hrf)\n", "\n", "def spm_hrf_d(t):\n", " \"\"\" Return temporal derivative of SPM HRF sampled at times `t`\n", " \"\"\"\n", " # This is what spm does!\n", " return spm_hrf(t) - spm_hrf(t - 1)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "t = np.arange(24)\n", "plt.plot(t, spm_hrf(t), label='HRF')\n", "plt.plot(t, spm_hrf_d(t), label='TD')\n", "plt.legend()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The design without HRF convolution" ] }, { "cell_type": "code", "collapsed": false, "input": [ "block_length = 24\n", "on_off = np.hstack((np.zeros(block_length), np.ones(block_length)))\n", "off_on = 1 - on_off\n", "X_hrf_over = np.column_stack((on_off, off_on, np.ones_like(on_off)))\n", "X_hrf_over = np.tile(X_hrf_over, (10, 1))\n", "show_design(X_hrf_over, 'before convolution', aspect=0.01)\n", "S = npl.svd(X_hrf_over, compute_uv=False)\n", "print('Smallest singular value (if close to 0 matrix may be rank deficient)')\n", "print(S.min())\n", "# Is -1 -1 1 still in the null space? Yes.\n", "np.allclose(X_hrf_over.dot([-1, -1, 1]), 0)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Design with convolution:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "hrf = spm_hrf(t)\n", "x0 = np.convolve(hrf, X_hrf_over[:, 0], mode='same')\n", "x1 = np.convolve(hrf, X_hrf_over[:, 1], mode='same')\n", "X_conv_over = np.column_stack((x0, x1, np.ones_like(x0)))\n", "show_design(X_conv_over, 'Convolved design', aspect=0.01)\n", "S = npl.svd(X_conv_over, compute_uv=False)\n", "print('Smallest singular value (if close to 0 matrix may be rank deficient)')\n", "print(S.min())" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the matrix with convolved columns is no longer rank deficient, and the vector:\n", "\n", "$$\n", "\\left[\\begin{matrix}-1\\\\-1\\\\1\\end{matrix}\\right]\n", "$$\n", "\n", "is no longer in the null space (because there is no null space for the convolved matrix)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Is -1 -1 1 still in the null space?\n", "np.allclose(X_conv_over.dot([-1, -1, 1]), 0)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, let's make some data based on the design:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x0_hrf = X_conv_over[:, 0]\n", "X_conv_well = np.column_stack((x0_hrf, np.ones_like(x0)))\n", "y = (x0_hrf - x0_hrf.mean()) * 1 + 10" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "npl.pinv(X_conv_well).dot(y)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We don't reconstruct 10 correctly because the ``x0`` regressor does not have mean 0:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x0_hrf.mean()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "X_conv_orth = np.column_stack((x0_hrf - x0_hrf.mean(), np.ones_like(x0)))\n", "npl.pinv(X_conv_orth).dot(y)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a temporal derivative:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "dhrf = spm_hrf_d(np.arange(24))\n", "x0_dhrf = np.convolve(dhrf, x0, mode='same')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The derivative is not orthogonal to the HRF, after convolution:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Basis functions, before convolution with block - nearly orthogonal\n", "hrf.dot(dhrf)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# After convolution with block - not orthogonal\n", "x0_hrf.dot(x0_dhrf)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can orthogonalize using the same mechamism as we use for the fit:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def orth_x_wrt_y(x, y):\n", " \"\"\" Orthonalize vector `x` with respect to matrix `y`\n", " \"\"\"\n", " y_in_x_beta = npl.pinv(y).dot(x)\n", " y_in_x = y.dot(y_in_x_beta)\n", " return x - y_in_x" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "x0_dhrf_orth = orth_x_wrt_y(x0_dhrf, X_conv_orth)\n", "x0_dhrf_orth.T.dot(X_conv_orth)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "x0_hrf.dot(x0_dhrf_orth)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we know, adding an orthogonal regressor to the whole the rest of the design cannot affect the fit for other regressors:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X_dconv_orth = np.column_stack((X_conv_orth[:, 0], x0_dhrf_orth, np.ones_like(x0)))\n", "npl.pinv(X_dconv_orth).dot(y)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With or without noise:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "ye = y + np.random.normal(size=(y.shape))\n", "print(npl.pinv(X_conv_orth).dot(ye))\n", "print(npl.pinv(X_dconv_orth).dot(ye))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But - if the new regressor is not orthogonal to the whole of the rest of the design, the fit can change, even for regressors that *are* orthogonal to the new regressor:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Design in which hrf regressor not orthogonal to constant\n", "partial_orthed = orth_x_wrt_y(x0_dhrf, x0_hrf[:, None])\n", "X_colinear = np.column_stack((x0_hrf, np.ones_like(x0_hrf)))\n", "X_mixed = np.column_stack((x0_hrf, partial_orthed, np.ones_like(x0_hrf)))\n", "X_mixed.T.dot(X_mixed)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the second temporal derivative regressor is orthogonal to the first HRF regressor, but not orthogonal to the constant. Estimating:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print(npl.pinv(X_colinear).dot(ye))\n", "print(npl.pinv(X_mixed).dot(ye))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimation in the presence of noise caused the parameter for the first HRF regressor to change, even though the added regressor is orthogonal to the first HRF regressor." ] } ], "metadata": {} } ] }