Neuroimage/Sign-consistency.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook implements the algortihm proposed in:\n",
    "\n",
    "Vanessa Gómez-Verdejo, Emilio Parrado-Hernández, Jussi Tohka and for the Alzheimer's Disease Neuroimaging Initiative, \"Sign-consistency based voxel importance for machine learning in brain imaging\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Create simulated data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def generate_data(N, Drel, Dirrel, sigma = 1):\n",
    "    # Case 1: all data condense to 2 points\n",
    "    # Case 2: completely separable, all variables are relevant\n",
    "    # Case 3: completely separable, some irrelevant variables start to take role\n",
    "    # Case 4: the irrelevant variables have same intensity as the \n",
    "    # Create labels\n",
    "    Y = np.concatenate((np.ones((N/2,)),-1*np.ones((N/2,))), axis=0)\n",
    "\n",
    "    # Relevant features\n",
    "    Xrel = Y[:, np.newaxis] * np.ones((1,Drel)) + np.random.normal(0, sigma, (N, Drel))\n",
    "\n",
    "    # Irrelevant features\n",
    "    Xirrel = np.random.normal(0, sigma, (N, Dirrel))\n",
    "    X = np.concatenate((Xrel,Xirrel), axis=1).T\n",
    "    GT = np.concatenate((np.ones((Drel,)),np.zeros((Dirrel,))), axis=0)\n",
    "\n",
    "    return X, Y, GT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Number of data\n",
    "Ndata = 200\n",
    "Drel = 1000\n",
    "Dirrel = 10000\n",
    "\n",
    "X, Y, GT = generate_data(Ndata, Drel, Dirrel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voxel importance with ensembles of linear SVMs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Train an ensemble of bagged SVMs**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "\n",
    "num_bagging_iter = 10000\n",
    "perc = 0.5\n",
    "conformal_learning = False\n",
    "\n",
    "def compute_bagging_clf(X, Y): \n",
    "    \n",
    "    w_sign_tot = np.zeros((X.shape[0],))\n",
    "    \n",
    "    ntr = X.shape[1]\n",
    "    pos_tr = np.arange(ntr)\n",
    "    p_1 = np.where(Y == 1)[0]\n",
    "    p_0 = np.where(Y == -1)[0]\n",
    "    n_1 = p_1.shape[0]\n",
    "    n_0 = p_0.shape[0]\n",
    "       \n",
    "    # Define the classifier to be used for the voxel selection\n",
    "    K = np.dot(X.T,X)\n",
    "    clf = SVC(C=100, kernel='precomputed')\n",
    " \n",
    "    for i in range(num_bagging_iter):\n",
    "        # Select training samples (at random) and balancing by clases\n",
    "        np.random.seed(i)\n",
    "        # In case of unbalanced problems, we balance the subsampling\n",
    "        samples_sel_1 = np.random.permutation(p_1)[:int(np.round(0.5*perc*ntr))]\n",
    "        samples_sel_0 = np.random.permutation(p_0)[:int(np.round(0.5*perc*ntr))]\n",
    "        samples_sel = np.hstack((samples_sel_0, samples_sel_1))\n",
    "        # Create a vector with the training positions\n",
    "        pos_tr_sel = pos_tr[samples_sel]\n",
    "       \n",
    "        Y_sel =Y[pos_tr_sel]\n",
    "        \n",
    "        # Train the classifier\n",
    "        clf.fit(K[pos_tr_sel[:, np.newaxis],pos_tr_sel], Y_sel) \n",
    "        w = np.dot(X[:,pos_tr_sel[clf.support_]],clf.dual_coef_.ravel())\n",
    "       \n",
    "        # Create a vector counting the number of positive weigths (indirectly we are also counting the negatives)\n",
    "        w_sign = np.zeros(w.shape)\n",
    "        w_sign [w>0] = 1\n",
    "        # Add it to the total sum of positive weigths\n",
    "        w_sign_tot = w_sign_tot + w_sign\n",
    "    return w_sign_tot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nPositives = compute_bagging_clf(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Voxel selection: hypothesis test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = num_bagging_iter\n",
    "# Estimate the probability of being positive\n",
    "hatp = nPositives/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHWRJREFUeJzt3XuwHOWZ3/HvTxISGCOBYJEKCyRzs2WIzZIY2LXXKwMG\nhGNgtwwBBxswu0UCxg7ExEAqkShXJZarCJjyOphdhRW+gLnsghbjwBJx2MUgwCvLYHEPV4nlYC6C\nghQY0JM/+j06c45nzukz09PTM/37VJ06PT3vdL/d09NPv9399qOIwMzM6mtarytgZma95UBgZlZz\nDgRmZjXnQGBmVnMOBGZmNedAYGZWc20HAkkrJQ1LeqBh3LclPSxpvaQbJM1ueO8CSY+n94/otOJm\nZlaMTloEVwJHjht3G7BfRBwAPA5cACDpI8AJwGJgKfA9Sepg3mZmVpC2A0FE3AW8Om7c7RGxJb1c\nCyxIw8cA10TEuxHxNFmQOKjdeZuZWXG6eY3gy8AtafgDwHMN721K48zMrMe6Eggk/WfgnYi4emRU\nk2J+toWZWQXMKHqCkk4BjgYObRi9Edi94fUC4PkWn3eAMDNrQ0S0de210xaBaDjal3QU8J+AYyLi\n7YZyq4ETJc2U9EFgb+C+VhONCP9FsGzZsp7XoSp/VVkXaQtlpEFb53VRhT+vi9G/TrTdIpD0Y2AJ\nsLOkZ4FlwIXATODv001BayPizIh4SNK1wEPAO8CZ0WnNzcysEJ3cNfSFiNgtImZFxB4RcWVE7BMR\nCyPiwPR3ZkP5/x4Re0fE4oi4baJpH3HE5zniiM9zwgmn8cYbb7RbxdqaP38RkpDE/PmLel6Hbtej\n6OWtwvozK5OqdmCeXSO4DoDttjuH++77Gfvvv3+Pa9UbQ0NDLFmyZMqfy1pjI9+rOm42tmNsHTqv\nx0TroujlnWh6VVi3c+fO59VXh7e+njdvIS+88HTp9WjH/PmLGB5+Biim3u3+RgaRJKLNawQVDQRZ\nnWbP3p+f//yagQgERf8AJlKFnVXRgSD/vAY/EJS5botWhfU3qDoJBH7WUEmyIJBdZBwJCDa4fHpp\n8DR+p9Onb1/aqc8yDEwg8A/POtG4/RTBgb9/5N13NH6nW7b8P0bvHuv/73hgAkHdfnj5A98sB8gc\nGrcfq5e67TuaGZhAUDf5N963c5ZrrQqtrTLvQrI6m1XL7cyBwCZVhSOmsUfs9Thyq0IArp/RA6ci\ntrN+OYBxIOiJ9o46ij6P7R1NHmO/qzIVHYDzft/lbRfFH31XbZvulwOYjm4flbQS+NfAcER8NI3b\nCfgJsBB4GjghIl5L711Glo/gTeDUiFjfZJpjbh/dZpvNvPzyJmDi2y7z3pZW5m2cE9Wvndv/Wk9j\n4lscO70Vsp1b/oq+xXGi6bVTv8btYFTzdTZ+vt1aZ3mn0e66be/73pbsKLmY30ve30G7v9POt9X8\n21me76Ds26h7dftos+Q05wO3R8SHgDWMJqdZCuwVEfsAZwCX55lBFgSKi6btHGXlbd71SzNwsOU7\nyhx/pGatdH6NqR1VOB1ZJx0FgmiSnAY4FliVhlel1yPjr0qfuxeYI2leJ/MvS97mXb80AztT9buQ\nij3Ha1YH3bhGsGtEDANExAvArmm8k9MMhN4cIRah6Gss9VONa1tWvDIvFpeYnKbbR61VPyqug1lT\n3rnUoa9Ad/uXtNfaqsN673eFJ6YBhiXNi4hhSfOBF9P43MlpYDkAb731YvO3JzWywcLwcDeOQro9\n/d5qfiG1nHlNm/a+1Gtz7PDvGv0Omh9j1NPoTneybbOIbXjW1kDcTw++6wd5LpYPDQ0xNDRUzAwL\nSIawCHiw4fUK4Btp+HzgW2n4aOCnafgQslwFzaYXEAERs2fvF42vs+o2N75cq8/knV7rz7SeRrNy\nk9ev3XrkX8aprovJlnd0eFbj1daYN29hrnUxlXU70TrLV25sHTufXnvbYzuK2M4mml475YpYFxNN\nL89yzJu3sOk2N9nnpl6/zn/DRf3W834mor39eEenhpQlp7kb2FfSs5JOA74FfEbSo8Bh6TURcQvw\nlKQngO8DZ7aY7ATaac727j7w8vRqGadyqqBXp9Ma69g94+8Yq4a820XVfyNj6+c7iopX+cdQv/76\nBkZ/xO3dOz92J1DMPfHNpjGVe93Hlst3r3a7y5h3nU20HHnnm3caRcyrCuWmurwTnQob/923049g\nsnvdy1gX401lu21nOfJs0xMpph9B899w3n4Ek/Vrybsc0WY/gm5cIyjJrIocvRRRj8G+5mDjt5HR\nH/WWLaM7ikH57su8xlQdnf2GG6/tZMrdFvr4ERNjT0tUox5VV5VTAFO/46e/9dM20rnB6qxX9LY6\nOr3GnAZ5PzP+tGpRt+b2cSCwqXPwNIOp7ECL3lZHpzc2p0G+z4xvaRV1a27NA4H7A5jVkfs2jFXz\nQNC/vWTNBk/3ThlW866u6ujji8VWPVW5gG/9qXudBHt9Mbbqat4iaFTPzETFqtu5/7pd9LZB5UCw\n1fjOUS8U8CPvbEfh5mzVFR342r2rq84BaaJ1NojrpTt3/nUlEEg6R9KvJT0g6UeSZkpaJGmtpEcl\nXS2p4qeliviRdzaNcm/DG8QfTb9p966uKrTEerX9TLTOqrBeitadO/8KDwSSdgPOBg6MLGvZDOAk\nsmcQXRxZwprNwOlFz3swVOEHZTZV3n76WbdODU0Htk9H/duRPWX008AN6f1VwJ90ad59zj8os/7W\nf63rwgNBRDwPXAw8S5Z85jVgHbA5IrakYhuB3Yqet5lZ7/XfwVzh5+kl7UiWlnIhWRC4jixh/XgT\nrKXlQCf5CMwGjW/NtfGG0v/lHU+pGxdsDweejIhXACT9LfCHwI6SpqVWwQRJaWBkwbbd9np++9vf\ndKGKZv3GiXhsvCXp//L0/6K2p9SNawTPAodI2lbZIcxhwAbgDuD4VOYU4KYuzNvMzKaoG9cI7gOu\nB34J/Irs8OUKsmxl50p6DJgLrCx63mZmNnVduZc/Ii7id9spTwEHd2N+ZmbWPvcsNjOrOQcCM7Oa\ncyAwM6s5BwIzs5pzIDAzqzkHAjOzmnMgMDOrOQcCM7Oa61ZimjmSrpP0sKQNkg6WtJOk21Jimlsl\nzenGvM3MbGq61SL4DnBLRCwGPgY8QvaIidtTYpo1wAVdmreZmU1BNzKU7QD8UURcCRAR70bEa2SP\npl6Viq0Cjit63mZmNnXdaBHsCbwk6UpJ6yRdIel9wLyIGAaIiBeA3+vCvM3MbIq68dC5GcCBwFkR\n8QtJl5CdFppCup7lgBPTmJm1NpT+L+94SoooNp2apHnAPRGxZ3r9SbJAsBewJCKGJc0H7kjXEMZ/\nPkZixuzZ+/P66xsYm5Cj2fBE7xVdrsx5DUq5Ktap6uWqWKeql6tincotFxFtZS3qRj6CYeA5Sfum\nUSOJaVYDp6ZxTkxjZlYRXclHAHwV+JGkbYAngdOA6cC1kr5MlsXs+Ak+b2ZmJelWYppfAR9v8tbh\n3ZifmZm1zz2LzcxqzoHAzKzmHAjMzGrOgcDMrOYcCMzMas6BwMys5hwIzMxqzoHAzKzmuhYIJE1L\nTx9dnV4vkrQ2Jaa5WlK3ejWbmdkUdLNF8DXgoYbXK4CLU2KazcDpXZy3mZnl1K1UlQuAo4G/ahh9\nKHBDGl4F/Ek35m1mZlPTrRbBJcB5pOejStoZeDUitqT3NwK7dWneZmY2BYWfp5f0WWA4ItZLWjIy\nOv01miARwnLAiWnMzFobSv+Xdzylblyw/QRwjKSjge2AHYBLgTmSpqVWwQLg+daTWA7Atttez29/\n+5suVNHMrN8tSf+Xp/8XtT2lbiSmuTAi9kgZyk4E1kTEycAdjOYgcGIaM7OKKLMfwfnAuZIeA+YC\nK0uct5mZtdDVe/kj4k7gzjT8FHBwN+dnZmZT557FZmY150BgZlZzDgRmZjXnQGBmVnMOBGZmNedA\nYGZWcw4EZmY150BgZlZzhQcCSQskrZH0kKQHJX01jd9J0m0pMc2tkuYUPW8zM5u6brQI3gXOjYiP\nAH8AnCXpw2SPmLg9JaZZA1zQhXmbmdkUdeOhcy9ExPo0/AbwMNnTRo8lS0hD+n9c0fM2M7Op6+o1\nAkmLgAOAtcC8iBiGLFgAv9fNeZuZWT5de+icpPcD1wNfi4g3JE2QiGa85YAT05iZtTaU/i/veEqK\nmML+Oe9EpRnAzcDPIuI7adzDwJKIGJY0H7gjIhY3+WyMJC+bPXt/Xn99A6PJzNRieKL3ii5X5rwG\npVwV61T1clWsU9XLVbFO5ZaLiPGZIHPp1qmh/wU8NBIEktXAqWnYiWnMzCqiGzmLPwH8W+BBSb8k\nC1cXAiuAayV9GXiW0WxlZmbWQ4UHgoj4OTC9xduHFz0/MzPrjHsWm5nVnAOBmVnNORCYmdWcA4GZ\nWc05EJiZ1ZwDgZlZzTkQmJnVnAOBmVnNlR4IJB0l6RFJj0n6RtnzNzOzsUoNBJKmAd8FjgT2A05K\nSWvMzKxHym4RHAQ8HhHPRMQ7wDVkCWvMzKxHupaPoIUPAM81vN5IFhzG+R8AvP32S2XUycys1soO\nBM2eld0kIcJ/BODtt5t9rNVwmeWqWKeql6tinaperop1qnq5KtapzHLtKTsQbAT2aHi9AHi+sUC7\niRXMzKw9ZV8juB/YW9JCSTOBE8kS1piZWY+U2iKIiPckfQW4jSwIrYyIh8usg5mZjdWVnMVmZtY/\netazeLKOZZJmSrpG0uOS7pG0R7PpDIIc6+IcSRskrZf095J270U9y5C3w6Gkz0vaIunAMutXpjzr\nQtIJadt4UNIPy65jWXL8RnaXtEbSuvQ7WdqLenabpJWShiU9MEGZy9J+c72kA3JNOCJK/yMLQE8A\nC4FtgPXAh8eV+ffA99LwvwGu6UVdK7Iu/hjYNg3/uzqvi1Tu/cCdwN3Agb2udw+3i72BfwJmp9e7\n9LrePVwX3wfOSMOLgad6Xe8urYtPAgcAD7R4fynw0zR8MLA2z3R71SLI07HsWGBVGr4eOKzE+pVp\n0nUREXdGxFvp5Vqy/hiDKG+Hw28CK4C3m7w3KPKsiz8H/iIiXgeIiEHteJNnXWwBZqfhHYFNJdav\nNBFxF/DqBEWOBa5KZe8F5kiaN9l0exUImnUsG79z21omIt4DNkuaW071SpVnXTQ6HfhZV2vUO5Ou\ni9TUXRARt5RZsR7Is13sC3xI0l2S7pZ0ZGm1K1eedXER8EVJzwE3A2eXVLeqGb+uNpHjwLHsfgQj\n8nQsa9ZrYhCvbOfsZAeSTgb+JdmpokE04bqQJOAS4JRJPjMI8mwXM8hOD32KrH/OP0rab6SFMEDy\nrIuTgCsj4hJJhwA/JHueWd3k3p806lWLYNKOZWRRbXcASdPJzoNO1CTqV3nWBZIOBy4APpeax4No\nsnWxA9mPe0jSU8AhwE0DesE4z3axEbgpIrZExNPAo8A+5VSvVHnWxenAtQARsRbYVtIu5VSvUjaS\n9ptJ0/3JeL0KBHk6lv0do0d+xwNrSqxfmSZdF5J+H7gcOCYiXu5BHcsy4bqIiNcjYteI2DMiPkh2\nveRzEbGuR/Xtpjy/kRuBQwHSTm8f4MlSa1mOPOviGeBwAEmLgVkDfM1EtG4Jrwa+BJBaRpsjYniy\nCfbk1FC06Fgm6SLg/oi4GVgJ/EDS48DLZF/+wMm5Lr4NbA9cl06PPBMRx/Wu1t2Rc12M+QgDemoo\nz7qIiFslHSFpA/Au8PVBbDXn3C6+DvylpHPILhyf0nqK/UvSj4ElwM6SngWWATOBiIgrIuIWSUdL\negJ4Ezgt13TTbUZmZlZTTlVpZlZzDgRmZjU3aSBo1qVZ0k6SbpP0qKRbJc1peK9p92ZJp6Tu4Y9K\n+lLxi2JmZu3I0yK4kizHcKPzgdsj4kNkd/NcAJCe77FXROwDnEF2pwuSdgL+K/Bxsm7PyxqDh5mZ\n9c6kgaBFl+bGxz+sYrS7d6vuzUcCt0XEaxGxmezq/1GdV9/MzDrV7jWCXUfuTY2IF4Bd0/hWXcHb\n6vZsZmbdV3Q/glaPhZjKYxR8P6uZWRuizVS/7bYIhkeeaCdpPvBiGt+qe3OuxyiM6PWjXrv5t2zZ\nsp7XwcvX/l/aQhk5jul1ffz9edlG/jqRNxCM79K8Gjg1DZ8K3NQwvln35luBz0iaky4cfyaNMzOz\nHpv01FCLLs3fInvcwZeBZ8meBUS06N4cEa9K+ibwC7JDqYsiu2hsZmY9NmkgiIgvtHjr8Bblv9Ji\n/F8Df523YoNqyZIlva5CVw368g26Qf7+BnnZOlW5Zw1JiqrVyWxE9sy/ke1THZ+bNSuKJKLki8Vm\nZjYgHAjMzGrOgcDMrOYcCMzMas6BwMys5hwIzMxqzoHAzKzmHAjMzGrOgcDMrOYcCMzMaq6jQCDp\nHEm/lvSApB9JmilpkaS1KTfx1ZJmpLIzJV2T8hnfI2mPyaZvZmbd13YgkLQbcDZwYER8lOwBdicB\nK4CLI8tnvBk4PX3kdOCVyPIZXwp8u5OKm5lZMTo9NTQd2D4d9W9Hlmzm08AN6f1VwHFpuDHP8fXA\nYR3O28zMCtB2IIiI54GLyfIRbAJeA9aRJaPZkoqN5CyGhrzFEfEesFnS3Hbnb2Zmxejk1NCOZEf5\nC4HdgO2BpU2Kjj6zd9wkaJG32MzMytNJ8vrDgScj4hUASX8L/CGwo6RpqVXQmJt4JJ/x85KmA7Mj\n4tVmE16+fPnW4SVLljihhJnZOENDQwwNDRUyrbYT00g6CFgJfBx4G7gSuB/4FPA3EfETSf8T+FVE\nXC7pTGD/iDhT0onAcRFxYpPpOjGNVZYT01hVdZKYpqMMZZKWAScC7wC/BP6MrBVwDbBTGndyRLwj\naRbwA+D3gZeBEyPi6SbTdCCwynIgsKrqWSDoBgcCqzIHAqsqp6o0M7O2ORCYmdWcA4GZWc05EJiZ\n1ZwDgZlZzTkQmJnVnAOBmVnNORCYmdWcA4GZWc05EJiZ1ZwDgZlZzTkQmJnVXKfJ6+dIuk7Sw5I2\nSDpY0k6SbkvJ62+VNKeh/GUpef16SQd0Xn0zM+tUpy2C7wC3RMRi4GPAI8D5wO0pef0a4AIASUuB\nvVLy+jOAyzuct5mZFaCTxDQ7AOsjYq9x4x8B/jgihiXNB+6IiMWSLk/DP0nlHgaWRMTwuM/7MdRW\nWX4MtVVVrx5DvSfwkqQrJa2TdIWk9wHzRnbuEfECsGsqvzV5fbKJ0cT2ZmbWI53kLJ4BHAicFRG/\nkHQJ2WmhVodIzSJV07LOWWxmNrGq5CyeB9wTEXum158kCwR7kU75THJqaOsppHHT9akhqyyfGrKq\n6smpobQDf07SvmnUYcAGYDVwahp3KnBTGl4NfAlA0iHA5vFBwMzMytdp8vqPAX8FbAM8CZwGTAeu\nBXYHngWOj4jNqfx3gaOAN4HTImJdk2m6RWCV5RaBVZWT15uVxIHAqsrJ683MrG0OBGZmNedAYGZW\ncw4EZmY150BgZlZzDgRmZjXnQGBmVnMOBGZmNedAYGZWcw4EZmY150BgZlZzHQcCSdNSYprV6fUi\nSWtTzuKrJc1I42dKuiblLL5H0h6dztvMzDpXRIvga8BDDa9XABennMWbgdPT+NOBV1LO4kuBbxcw\nbzMz61BHgUDSAuBoskdRjzgUuCENrwKOS8PHptcA15PlLzAzsx7rtEVwCXAe6bm8knYGXo2ILen9\njYzmJd6aszgi3gM2S5rb4fzNzKxDbecslvRZYDgi1ktaMjKa381NHA3vjZkEzllsZtaWquQs/m/A\nycC7wHbADsCNwBHA/IjYklJSLouIpZL+dxq+V9J04J8jYtcm03ViGqssJ6axqupVzuILI2KPlLz+\nRGBNRJwM3AEcn4qdwticxaek4eOBNe3O28zMitONfgTnA+dKegyYC6xM41cCu0h6HPgPqZyZmfWY\ncxabTYFPDVlVOWexmZm1zYHAzKzmHAjMzGrOgcDMrOYcCMzMas6BwMys5hwIzMxqzoHAzKzmHAjM\nzGrOgcDMrOYcCMzMaq7tQCBpgaQ1kh6S9KCkr6bxO0m6LeUsvlXSnIbPXJZyFq+XdEARC2BmZp3p\npEXwLnBuRHwE+APgLEkfJnuq6O0pZ/Ea4AIASUuBvVLO4jOAyzuquZmZFaKTfAQvRMT6NPwG8DCw\ngLG5iVel16T/V6Xy9wJzJM1rd/5mZlaMQq4RSFoEHACsBeZFxDBkwQIYyUK2NWdxsonRfMZmZtYj\nHQcCSe8Hrge+lloGrR7Q3uw52X6Yu5lZj7WdvB5A0gyyIPCDiBhJSTksaV5EDEuaD7yYxm8Edm/4\n+ALg+WbTdfJ6M7OJVSJ5PYCkq4CXIuLchnErgFciYoWk84EdI+J8SUcDZ0XEZ1NS+0sj4pAm03SG\nMqssZyizquokQ1nbgUDSJ4B/AB4k+2UEcCFwH3At2dH/s8DxEbE5fea7wFHAm8BpEbGuyXQdCKyy\nHAisqnoSCLrFgcCqzIHAqso5i83MrG0OBGZmNedAYGZWcw4EZmY150BgZlZzDgRmZjXnQGBmVnMO\nBGZmNedAYGZWcw4EZmY150BgZlZzpQcCSUdJekTSY5K+Ufb8zcxsrFIDgaRpwHeBI4H9gJNSnuPa\nKOr54VU16Ms36Ab5+xvkZetU2S2Cg4DHI+KZiHgHuIbRnMa1MOgb46Av36Ab5O9vkJetU2UHgvF5\nizfivMVmZj3VUarKNlQub/F5553HI488AsCNN97I9OnTe1kdM7PSlZqYJqWoXB4RR6XX5wMRESsa\nyjjTh5lZG/oiQ5mk6cCjwGHAP5OltTwpIh4urRJmZjZGqaeGIuI9SV8BbiO7PrHSQcDMrLcql7PY\nzMzK1bOexXk7lkn6vKQtkg4ss36dyrN8kk6QtEHSg5J+WHYdOzHZ8knaXdIaSeskrZe0tBf1bIek\nlZKGJT0wQZnLJD2elu2AMuvXqcmWT9IXJP0qLdtdkv5F2XXsRJ7vL5X7uKR3Jf1pWXXrVM5tc4mk\nX0r6taQ7ck04Ikr/IwtATwALgW2A9cCHm5R7P3AncDdwYC/q2q3lA/YG/gmYnV7v0ut6F7x83wfO\nSMOLgad6Xe8pLN8ngQOAB1q8vxT4aRo+GFjb6zoXvHyHAHPS8FGDtnypzDTg/wA3A3/a6zoX+N3N\nATYAH0ivc+1XetUiyNux7JvACuDtMitXgDzL9+fAX0TE6wAR8VLJdexEnuXbAsxOwzsCm0qsX0ci\n4i7g1QmKHAtclcreC8yRNK+MuhVhsuWLiLUR8Vp6uZY+6+uT4/sDOBu4Hnix+zUqTo5l+wJwQ0Rs\nSuVz7Vd6FQgm7ViWmtsLIuKWMitWkDwd5/YFPpSa3ndLOrK02nUuz/JdBHxR0nNkR11nl1S3Moxf\n/k302c5yCv4M+FmvK1EkSbsBxwGX07xvUz/bF5gr6Q5J90v6Yp4Pld2hbMSEHcskCbgEOGWSz1RV\nno5zM8hOD30K2AP4R0n7jbQQKi7P8p0EXBkRl6T+Iz8ke77UIKhcx8hukPRp4DSy0xGD5FLgGxER\n2a6mr/Ytk5kBHAgcCmwP3CPpnoh4YrIP9cJGsp3fiAXA8w2vdyDbaQyloDAfuEnSMRGxrrxqtm2y\n5Rspc09EbAGelvQosA/ZdYOqy7N8p5M9XJCIWCtpW0m79NkpsFY2Ars3vG62/H1N0keBK4CjImKy\n0yz95l8B16R9yy7AUknvRMTqHterCBuB30TEW8Bbkv4B+BjZNb2WenVq6H5gb0kLJc0ETgS2fgkR\n8XpE7BoRe0bEB8nOU36uT4IATLJ8yY1kURtJu5AFgSdLrWX78izfM8DhAJIWA7P6LAiI1keKq4Ev\nwdbe8psjYrisihWk5fJJ2gO4AfhiRPzfUmtVnJbLl/YrI/uW64Ez+ywITLRt3gT8kaTpkt5HdjPD\npH21etIiiBYdyyRdBNwfETeP/wh91HzLs3wRcaukIyRtAN4Fvt4vR145v7+vA38p6RyyC8entJ5i\ntUj6MbAE2FnSs8AyYCbZ41CuiIhbJB0t6QngTbLTJ31jsuUD/gswF/heOmp+JyIO6lV9pyrH8jXq\nq1N6ObbNRyTdCjwAvAdcEREPTTrddIuRmZnVlFNVmpnVnAOBmVnNORCYmdWcA4GZWc05EJiZ1ZwD\ngZlZzTkQmJnVnAOBmVnN/X9+QHTwsx0iSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe7022bee90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distribution irrelevant voxels (upper) vs. relevant ones (lower)\n"
     ]
    }
   ],
   "source": [
    "# As it is a sythetic problem, we know the ground truth, let's plot the distributions over each sign\n",
    "pos_GT1 = np.where(GT==1)[0]\n",
    "pos_GT0 = np.where(GT==0)[0]\n",
    "\n",
    "plt.figure(1)\n",
    "plt.subplot(2,1,1)\n",
    "plt.hist(hatp[pos_GT0],100)\n",
    "plt.subplot(2,1,2)\n",
    "plt.hist(hatp[pos_GT1],100)\n",
    "plt.show()\n",
    "print 'Distribution irrelevant voxels (upper) vs. relevant ones (lower)'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "# Statistical test\n",
    "def computeZTest(hatp, N, p0=0.5, gamma =0.5):\n",
    "    s = ((1-gamma)/gamma)*hatp*(1-hatp)\n",
    "    # To avoid divide by zero  \n",
    "    s[np.where(s==0)]=np.finfo(float).eps\n",
    "    z = np.divide(hatp-p0, np.sqrt(s))\n",
    "    p = norm.cdf(z)\n",
    "    return z, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "z, p = computeZTest(hatp, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Irrelevant voxels\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X20XFWZ5/HvL4EEBXlVkpEAUYgQmXGQXg3p1mkjMkBw\nCXQvQexWEbGXNvgy0HYbmLUm16VrpuMaRF22gy9pJviGgK0gjYJMDCONwZcYQYIYhpcQ6FwGIbLE\naQy5z/xx9s2tW1Tde+55qTq36vdZq9atOrXr1L6nznlq1z5nP1sRgZmZDb45/a6AmZn1hgO+mdmQ\ncMA3MxsSDvhmZkPCAd/MbEg44JuZDYlpA76kNZJGJd3Vsuzjku6VtEnSNyTt2/LcJZK2pOdPrqvi\nZmY2M3la+FcCp7QtuwU4JiKOBbYAlwBIeiVwNrAUWAF8VpKqq66ZmRU1bcCPiNuBp9qW3RoRY+nh\nBmBRun86cHVEPBcRD5F9GRxfXXXNzKyoKvrw3wXclO4fAjzS8tyjaZmZmfVZqYAv6T8DOyPia+OL\nOhRz7gYzswbYo+gLJZ0LnAac2LJ4G3Boy+NFwGNdXu8vAjOzAiKi0LnRvC180dJ6l3Qq8LfA6RHx\nbEu5G4BzJM2T9DLgSOBH3VYaEb5FsGrVqr7XoSm3pmyLtIcy/gN1mLdFE27eFhO3MqZt4Uv6KrAc\nOEjSVmAVcCkwD/heughnQ0RcEBGbJV0DbAZ2AhdE2RqamVkl8lyl8+cR8dKImB8Rh0XElRGxJCIO\nj4jj0u2ClvL/LSKOjIilEXHLVOs++eQ3c/LJb+bd7/4Au3btquL/GSoLFy5GEpJYuHDxwNeh6vdq\nwvYz6yX1qwGe9eFfC8DcuW9jx44n2GefffpSl35bv349y5cvn/Hrsl9X45+fSv/cK6LqOky1Lap+\nr6nW14Rte+CBC3nqqVEAFiw4nO3bH+p5HYpauHAxo6MPA9XUvegxMogkEQX78Psc8LP33nPPfXjy\nye0DEfCr3tGn0oSg1Ms6DFvAb0IdiprNdW+6MgHfuXQqlgX77GTfeOC3weVuocHT+pnOnbv37vuD\n8BnPuoDvA8zKaN1/quAv+Nkjb+xo/UzHxn7HxNVas/8znnUBf9gOsPxfcPP9RZhD6/5jw2XYYkcn\nsy7gD5v8O+mzOct114RfT02ogw2D4WwgOeDbbk1oATWhDr3U+gU3bMGnv8o3kFrNls/RAb9WxVoR\nVfczu9Wcx/xJB2yvTO5iKh98ZhJ4erdfzK80GDYxuFb9Odal6AQoB0i6RdJ9km6WtF/Lc59OE6Bs\nknRs3oocccQxAzioplgroup+5mFrNbdqDw7dTXxWs7l/fyaBZ/J+sb3G42Xytm2tU5HjdLYE1yYq\nOgHKSuDWiDgKWMfEBCgrgCMiYgnwHuCKvBV5/PGtVPkBFglyeXe+ZnyZDLt8rcb24GDdVNvFkdcw\nN0b6odAEKMAZwNp0f216PL78qvS6O4H9JC2opqr1y7vzDcdO2vSTWt1bjWbWWdE+/IMjYhQgIrYD\nB6flngBlYPSnxVeFqs+BDJ9ife7e7s1X9UnbHk+AUmcrtNoTTVZEsROpw3Ctfb3jM4r9ehqG7T7b\nFZ0AZVTSgogYlbQQeDwtzz0BSmYEgF27fl+wGuM7JoyOVt2qmFh3Pevvv9a8P718rzlzXphGMNLx\n8YTJn0Hn9sRwmgiu0+2bVRwj83d/4c62JG5Nlyf31vr161m/fn0l75creZqkxcC3I+LfpcergScj\nYrWklcD+EbFS0mnAhRHxRknLgE9GxLIu65yUPG3nzmfIk2ypPSlTlcmvuq978jpmknSrSAKpov/j\nTLfF5OXPX8fE/b3IAkemdcfMu52n27bd/o985Z5fx6rfK+//NdMkYZ0+gyL7Wd7PNO/+U8W2KHss\ndQuGU22zqZQ/rrq/V1V1yvua2pKnKZsA5Q7gFZK2SjoP+DvgP0q6D3hDekxE3AQ8KOl+4HPABV1W\nO40iP0P7cx11b/Xrf8z7E7+f3WCtdaxP/ss8e2l+zvo0/RiZXL/huDiitxqTHrm9hd+rVklbnaZc\nX7kWfvdWcv565G+tVtvCz78tipYr38Ivv83y7D95W9Dt3WWt3VZ5W6tFf0nWvc3qbOEX/Qx618Kf\nOI6L/Oro3I3auxZ+4UnMe2d+Q1ojZesx+OcErH0fmfi8x8YmAsKgfPa9PAfUHOXOibSee8n0dl+Y\nBakVmjICsjddBtVoyk/3vF0Ng2I27SPlDc6gtjqOl4l1tubUz/ua9u7Qqi55nQUB32bOX5Jm/U2r\nMbHOyTn1872m/ZdTVZe8DknAb/qoUTOr2uD8AqnOkAT82Ttq1Gzw1NfV59G+U5sFJ22teZpyIt1m\np9YLGKrdjyafFPU+2m5IWvitnDKhvGHrmx+2k882qIYw4LcPItpe8mAuf4a/mYN5bELVX3BFvkCa\ncuVVv0y1zQbxC7mez7tUwJd0kaRfSLpL0lckzZO0WNIGZZOjfE1Sw7uNyh7M5c/w9/bk0iAeHLNN\nkX2uKVde9Wv/mWqbDeIvzno+78IBX9JLgfcDx0XEq8jOB7wVWA1cliZH2QGcX0VFB08TDxyz6Xj/\nmc3KdunMBfZOrfgXkGXGfD3wjfT8WuBPS77HgPKBYza7zb5fy4UDfkQ8BlwGbCWb6OQ3wEZgR0SM\npWLbgJeWraSZWfPMvkZb4f51SfuTTWl4OFmwvxZY0aHoFFtjBCiTD99s0PiSV2u3Pv0dKb2mMidU\nTwIeiIgnASR9E/hjYH9Jc1IrP9cEKHPn/nfGxnaWqIrZoKjvGnWbrZanvyPp70cKr6lMH/5WYJmk\nvZQ1Sd4A3AN8HzgrlTkXuL7Ee5iZWUXK9OH/CLgO+Bnwc7LmyOeBlcDFkn4FHAisqaCeZmZW0qyY\nAKWXE130v1wT69T0ck2sU9PLNbFOTS/XnDrVNsWhmZkNBgd8M7Mh4YBvZjYkHPDNzIaEA76Z2ZBw\nwDczGxIO+GZmQ8IB38xsSJSdAGU/SddKulfSPZJOkHSApFvSBCg3S9qvqsqamVlxZVv4nwJuioil\nwL8HfkmWWuHWNAHKOuCSku9hZmYVKDPj1YuA/xARVwJExHMR8RuylMlrU7G1wJmla2lmZqWVaeG/\nHHhC0pWSNkr6vKQXAgsiYhQgIrYDL6miomZmVk6ZfPh7AMcBF0bETyRdTtadM4NsbCOAJ0AxM+tu\nffo7UnpNhbNlSloA/DAiXp4ev5Ys4B8BLI+IUUkLge+nPv721ztbZt/fa1DKNbFOTS/XxDo1vVxz\n6tTzbJmp2+YRSa9Ii8YnQLkBeGdadi6eAMXMrBHKdOkAfAD4iqQ9gQeA84C5wDWS3kU2K9ZZU7ze\nzMx6pFTAj4ifA3/Y4amTyqzXzMyq55G2ZmZDwgHfzGxIOOCbmQ0JB3wzsyHhgG9mNiQc8M3MhoQD\nvpnZkHDANzMbEqUDvqQ5KVvmDenxYkkb0gQoX5NUdjSvmZlVoIoW/geBzS2PVwOXpQlQdgDnV/Ae\nZmZWUtkpDhcBpwFfbFl8IvCNdH8t8Kdl3sPMzKpRtoV/OfA3pLydkg4CnoqIsfT8NuClJd/DzMwq\nULh/XdIbgdGI2CRp+fjidGs1RcL9EcAToJiZdbc+/R0pvaYyE6D8V+BtwHPAC4AXAd8CTgYWRsSY\npGXAqohY0eH1ngCl7+81KOWaWKeml2tinZperjl16scEKJdGxGFpxqtzgHUR8Tbg+0zkwPcEKGZm\nDVHHdfgrgYsl/Qo4EFhTw3uYmdkMVXKNfETcBtyW7j8InFDFes3MrDoeaWtmNiQc8M3MhoQDvpnZ\nkHDANzMbEg74ZmZDwgHfzGxIOOCbmQ0JB3wzsyFROOBLWiRpnaTNku6W9IG0/ABJt6QJUG6WtF91\n1TUzs6LKtPCfAy6OiFcCfwRcKOlostQKt6YJUNYBl5SvppmZlVUmedr2iNiU7v8WuBdYBJxBNvEJ\n6e+ZZStpZmblVdKHL2kxcCywAVgQEaOQfSkAL6niPczMrJzSydMk7QNcB3wwIn6b5bnPawTwBChm\nZt2tT39HSq+p8AQoAJL2AG4EvhMRn0rL7gWWR8SopIXA9yNiaYfXegKUvr/XoJRrYp2aXq6JdWp6\nuebUqecToCT/AGweD/bJDcA7031PgGJm1hBl5rR9DfAXwN2Sfkb29XMpsBq4RtK7gK1MzH5lZmZ9\nVDjgR8Q/A3O7PH1S0fWamVk9PNLWzGxIOOCbmQ0JB3wzsyHhgG9mNiQc8M3MhoQDvpnZkHDANzMb\nEg74ZmZDoraAL+lUSb+U9CtJH67rfczMLJ9aAr6kOcBngFOAY4C3pslRzMysT+pq4R8PbImIhyNi\nJ3A12cQoZmbWJ6Xz4XdxCPBIy+NtZF8CbT4BOB++mVkv1BXwO+Vq7pB4/68BGBvr9LJu94ehXBPr\n1PRyTaxT08s1sU5NL9eUOhVTV8DfBhzW8ngR8FhrgaIJ/M3MrJi6+vB/DBwp6XBJ84BzyCZGMTOz\nPqmlhR8RuyS9D7iF7EtlTUTcW8d7mZlZPqXmtDUzs9mj9pG20w3AkjRP0tWStkj6oaTDOq1nEOTY\nFhdJukfSJknfk3RoP+rZC3kH5kl6s6QxScf1sn69lGdbSDo77Rt3S/pyr+vYKzmOkUMlrZO0MR0n\nK/pRz7pJWiNpVNJdU5T5dIqbmyQdm2vFEVHbjewL5X7gcGBPYBNwdFuZvwI+m+6/Bbi6zjr165Zz\nW7wO2Cvdf+8wb4tUbh/gNuAO4Lh+17uP+8WRwE+BfdPjF/e73n3cFp8D3pPuLwUe7He9a9oWrwWO\nBe7q8vwK4J/S/ROADXnWW3cLP88ArDOAten+dcAbaq5Tv0y7LSLitoj41/RwA9l4hkGUd2DeR4HV\nwLO9rFyP5dkWfwn8fUQ8DRART/S4jr2SZ1uMAfum+/sDj/awfj0TEbcDT01R5AzgqlT2TmA/SQum\nW2/dAb/TAKz2ILa7TETsAnZIOrDmevVDnm3R6nzgO7XWqH+m3RbpJ+qiiLiplxXrgzz7xSuAoyTd\nLukOSaf0rHa9lWdbfAR4u6RHgBuB9/eobk3Tvq0eJUcDsa7r8MflGYDVaXTBIJ5JzjkYDSS9DfgD\nsi6eQTTltpAk4HLg3GleMwjy7Bd7kHXr/AnZ+JYfSDpmvMU/QPJsi7cCV0bE5ZKWAV8my9c1bHLH\nk1Z1t/CnHYBF9i11KICkuWT9lFP9lJmt8mwLJJ0EXAK8Kf2sHUTTbYsXkR3E6yU9CCwDrh/QE7d5\n9ottwPURMRYRDwH3AUt6U72eyrMtzgeuAYiIDcBekl7cm+o1yjZS3Ew6xpN2dQf8PAOwvs1ES+4s\nYF3NdeqXabeFpFcDVwCnR8Sv+1DHXplyW0TE0xFxcES8PCJeRnY+400RsbFP9a1TnmPkW8CJACm4\nLQEe6GkteyPPtngYOAlA0lJg/gCf0xDdf9neALwDIP3S2RERo9OtsNYunegyAEvSR4AfR8SNwBrg\nS5K2AL8m+5AHTs5t8XFgb+Da1K3xcESc2b9a1yPntpj0Ega0SyfPtoiImyWdLOke4DngQ4P4Kzjn\nfvEh4AuSLiI7gXtu9zXOXpK+CiwHDpK0FVgFzAMiIj4fETdJOk3S/cAzwHm51psu6zEzswHnKQ7N\nzIaEA76Z2ZDIHfAlzZH0M0k3pMf/U9IDadlGSa9qKTvzIb9mZlarmZy0/SBwDxOj3AL464j4Zmuh\nlNviiIhYIukEsqtOllVRWTMzKy5XC1/SIuA04Is5Xl9oyK+ZmdUrb5fO5cDf8PyRXB9L3TaXSdoz\nLSs05NfMzOo1bZeOpDcCoxGxSdLylqdWRsRoCvRfAD4MfIycQ34l+XpQM7MCouAUsXla+K8BTpf0\nAPA14ERJV42P6krD/68ky3QHMxjy2+8UpE25rVq1qu91aMrN28Lbwtti6lsZ0wb8iLg0Ig6LiJeT\njYJdFxHvkLQQdie6OhP4RXpJoSG/ZmZWrzKpFb6S8nqIbKKC9wJEwSG/ZmZWrxkF/Ii4jWwGIiKi\n60QlEfG+kvUaKsuXL+93FRrD22KCt8UEb4tq9C2XjqTo13ubmc1WkogaT9rW5pBDjuaQQ47m+ONP\nZOfOQU39bmbWDLlb+JLmAD8BtkXE6ZIWk805eQCwEXh7RDyX8lhfRTZj0xPAWyJia4f1BdwLwB57\nHMdTTz3OPvvsU/4/MjMbYL1q4X8Q2NzyeDVwWUQcBewgm4mG9PfJiFgCfJIsx3sXRwNHk32XmJlZ\nncqkVjgR+Ea6v5bs0kzIUiusTfevA7qe3DUzs94plFpB0kHAUxExlp5vnV1+d2qFiNgF7JB0YGU1\nNjOzQqYN+K2pFZhIm9BprsVoeW7SKsgxm7qZmdUrz3X446kVTgNeALyIrG9+P0lzUiu/NX3CeGqF\nxyTNBfaNrvNvjgCwa9fv+cEPfsCKFSuK/ydmZgNo/fr1rF+/vpJ1zeg6fEmvI8uBf7qkrwP/GBFf\nl/Q/gJ9HxBWSLgD+bURcIOkc4MyIeN7E5NlVOtl777nnPjz55HZfpWNmNo1+XYe/ErhY0q+AA4E1\nafka4MWStgD/KZUzM7M+6+tIW7fwzcxmZtaOtDUzs+ktXLgYSWTJiYtzwDcza7jR0YfJekRqzocv\nab6kOyX9TNLdklal5VdKeiAt3yjpVS2v+bSkLWn6w2NL1dDMzCox7WWZEfGspNdHxO/SZZb/LOm7\n6ekPRcQ/tpaXtAI4IiKWSDoBuAJYVnnNzcxsRnJ16UTE79Ld+WRfEuMjbDt1KJ1BljyNiLiT7Hr9\nBSXraWZmJeXNpTNH0s+A7cD3IuLH6amPpW6by9Jk5tCSWiF5lIm0C2Zm1id5W/hjEfFqshG1x0t6\nJbAyIpYCfwgcBHw4Fe/U6ndqBTOzPpvRVToR8TTZFIenjk9MHhE7gSuB41Ox8dQK41rTLrQZAUZ2\np1YwM7PJJtIqjDCejqaoaQdepYnKd0bEbyS9ALgZ+DtgY0RsV3Zh6CeA/xcRl6acOxdGxBslLQM+\nGRHPO2nrgVdmZvlkYXYiP2XRgVd5kqf9G2BtmvFqDvD1iLhJ0v9KXwYCNgHvBUjPnSbpfuAZ4Lwi\nFTMzs2o5tYKZWcNV1cL3SFszsyHhgG9mNiQc8M3MhkSZXDqLJW2QdJ+kr0naIy2fJ+nqlEvnh5IO\nq/ufMDOz6U0b8CPiWeD1aeDVscCKlCNnNXBZRBwF7ADOTy85H3gyIpaQTYX48VpqbmZmM1I0l04A\nrwe+kZavBc5M989IjwGuA95QSU3NzKyUQrl0gP8D7EgTmEM2unY8X87uXDoRsQvYIenASmttZmYz\nlmfgFSmwv1rSvsA3gaWdiqW/7deHtl5A2mYEYHdqhRUrVuSpjpnZ0JicWqGcGQ+8kvRfgN8Bfwss\njIixlEJhVUSsSLnyV0XEnSl//r9ExMEd1uOBV2ZmOfRs4JWkF0vaL91/AXASsBn4PnBWKnYucH26\nf0N6THp+XZGKmZlZtcrk0rkXuFrSR4GfAWtS+TXAlyRtAX4NnFNDvc3MbIacS8fMrOGcS8fMzGbE\nAd/MbEjkOWm7SNI6SZtTaoX3p+WrJG2TtDHdTm15zSUptcK9kk6u8x8wM7N88py0fQ64OCI2SdoH\n+Kmk76XnPhERn2gtLGkpcDbZtfqLgFslLYl+nSwwMzMgXy6d7RGxKd3/LXAvE6NqO504OAO4OiKe\ni4iHgC1MzHdrZmZ9MqM+fEmLyRKo3ZkWXShpk6Qvjl+rT0tqheRRJr4gzMysT3IH/NSdcx3wwdTS\n/yxwREQcS5Zj57Lxoh1ePkVqhZHdqRXMzGyyyakVRkqtK9d1+CnX/Y3AdyLiUx2ePxz4dkS8StJK\nICJidXpud6qFttf4Onwzsxx6fR3+PwCbW4O9pIUtz/8Z8It0/wbgnDQRysuAI4EfFamcmZlVZ9qr\ndCS9BvgL4O6UIjmAS4E/l3QsMAY8BLwHICI2S7qGLN/OTuACX6FjZtZ/Tq1gZtZwTq1gZmYz4oBv\nZjYkiqRW+EBafoCkWyTdJ+nmluvwkfTplFphU+rnNzOzPsvTwh9PrfBK4I/IBlsdDawEbo2Io8gm\nObkEQNIKsuvzl5CdyL2ilpqbmdmMFE2tsIgshcLaVGxtekz6e1Uqfyewn6QFFdfbzMxmqGhqhQ3A\ngogYhexLARift9apFczMGqhMaoVu13POILWCmZn1Sp70yOOpFa4DvhQR45OVj0paEBGjadTt42n5\nNuDQlpcvAh7rvOYRgN25dFasWDHT+puZDbTJuXTKyZtL5yrgiYi4uGXZauDJiFid8ufsHxErJZ0G\nXBgRb5S0DPhkRCzrsE4PvDIzy6GqgVfTBvyUWuF/A3endxxPrfAj4Bqy1vxW4KyI2JFe8xngVOAZ\n4LyI2NhhvQ74ZmY59Czg18UB38wsH6dWMDOzGXHANzMbEg74ZmZDIk8unTWSRiXd1bJslaRtkjam\n26ktz12S8ujcK+nkuipuZmYzk6eFfyVwSofln4iI49LtuwCSlgJnA0uBFcBnlZ1tMDOzPsuTS+d2\n4KkOT3UK5GcAV0fEcxHxELAFOL5UDc3MrBJl+vAvTOmPv9iSGtl5dMzMGqpowP8sWQrkY4HtwGVp\n+Qzz6IwAI7tTK5iZ2WSTUyuMlFpX3tQKhwPfjohXTfVcSrEQEbE6PfddYFVKk9z+Og+8ssZauHAx\no6MPA7BgweFs3/5QfytkQ63XA69ES+s9JUsb92fAL9L9G4BzJM2T9DLgSLIUDGaVWrhwMZJ23xYu\nXFxqHe2vz4J9lklkPPCbzXZ5Lsv8KnAH8ApJWyWdB3xc0l2SNgGvAy4CiIjNZPl1NgM3ARdEv3I3\nWKNMFVyLaA3IMwnKrfWYHNS3T/oC6ZUqvrjM8nIuHeuJDj9JZ7yO1m6WTOs68q2zvR6d7z//ubqO\nk8n16f5e7f+7u5mGi3Pp2NCZ3KrvrrXVPHfu3pW23PvVIi/6i8aslQO+DZzW4Dg29jtaA2WV627v\nCmoN/u6qsSYqmlrhAEm3SLpP0s0t1+Ej6dMptcImScfWVXEbTMVb5/Nr7H+fat3P0qnV3a8Wub9o\nbCpFUyusBG6NiKOAdcAlAJJWkF2fvwR4D3BFhXW1gqo+YVrnuou3zp/NWa6IKtY9v+LPYH7HL8a8\nv0BsOBVNrXAGsDbdX5sejy+/Kr3uTmA/SQuqqaoVVfUlht2vdCm27tb1Da7WXwLbp/h/834xTKxv\n8hdjvvf1F8BwKtqHf3BEjAJExHbg4LS8cGqFI444xjtiQ7V3ExQ5eTrZ/ELra7683UpT/WLo3EVU\njdb3Lbb+On8tWv32qHh9BVIrwK5dv+fxx7fuLjo6OsgtvX6bvzsgzZnzwtQ65HmPWy/7mwjI47p9\nPvM7BLvWyx7HPUu+9c02rf9XFf9Tp+1ZpYn1T3WZZ7fLYUdH95pUP18qWo/JqRXKKdrCHx3vqkmj\nbh9Py7eRTWo+bhHwWPfVjAAjzJ07r2A1Zo86TqZ1O8FZrCtg8uOpux2mX/fsb603QZ3nJSavvzWg\n5/9FV/4Xg01v+fLl6d4IZYN+odQKZCkU3pnuvxO4vmX5OwAkLQN2jHf9DLs6TqZ1O8FZzYFXd7Cx\nZpnfJcCX//x95dDMtW+zqkw70jalVlgOHASMAquAbwHXkrXmtwJnRcSOVP4zwKnAM8B5EbGxy3on\njbTdufMZejGysV86jags+/9ONWq0dX1FR5fOnnJNrFPTy9XxXnuRNRSe313Yvr5BO76rliNeFPoW\naExqBQf8KgP+xIE3oUnBZjYEr0Ev19865Ukf0e2cUt2akCm1roA/q0faVn3FQJH11f1ztdj/6L50\nmx26XeJbRTdl+7GZ9zzXVJcaz/qrlCKiLzcgIAIi9txz72h9DPN3f+pz5rywNXLFggWHx7jJryHK\nmmp9CxYcnqMOk1/X+ppO5Tr9v1Ovf37b+uiyvqneaxDLNbFOTS/Xzzq178f51tftuGo9Xoocc62v\nny4OFIk53WLHVOVybAuK3Ep16Uh6CPgNMAbsjIjjJR0AfB04HHgIODsiftPhtUFqfXbq0ul8P3s8\nXuf2Lo0y/8t06+v2XKefXt3q116u+3PtXTIzXcewlWtinZperol1ytftWf0xl+9Yn+65bqZ6TfcM\nsM3s0hkDlkfEqyNifLLyjmkXqjPznCnV/wyrM28L+AoZs3Z1H3PTv281793fQYdlA746rKM97cKZ\nJd+jzcyDYd4+ucmm+qAdkM16a6pjrs4vg6nOh3XOZzR1w7K/59fKBvwAbpb0Y0nvTssWxOS0Cy8p\n+R45TQ7QrR/AVPIOKilSDzPrhbINsKKt+O6DGJs6CK1saoU/jojtkl4C3CLpPvrW7J08XH9srL2/\nblzdw9WrHl5vZvWqI9VH3XGmmFIBP7XgiYj/K+lbwPGktAsRMdqWdqGDESDLpdM7g5rHxcyao8qG\n3/r0d6Tkeih+lY6kFwJzIuK3kvYGbgE+ArwBeDIiVkv6MHBARKzs8PooepXOYJdrYp2aXq6JdWp6\nuSbWqenlmlOnolfplGnhLwC+mQVu9gC+EhG3SPoJcI2kd5HSLpR4DzMzq0jhgB8RDwLPm8IwIp4E\nTipTKTMzq96sTq1gZmb5OeCbmQ2J2gK+pFMl/VLSr9LJWzMz66NaAr6kOcBngFOAY4C3Sjq6jvcy\nM7N86mrhHw9siYiHI2IncDVZygUzM+uTugL+IcAjLY+3pWVmZtYnZVMrdNNpUMDzRnjtu++bAHjm\nmV6OtDUzG051BfxtwGEtjxcBj7UXevrpG9uWKMf9YSjXxDo1vVwT69T0ck2sU9PLNaVOxdQyp62k\nucB9ZGkW/gX4EfDWiLi38jczM7NcamnhR8QuSe8jy68zB1jjYG9m1l+1tPDNzKx5ah9pO90ALEnz\nJF0taYukH0o6rNN6BkGObXGRpHskbZL0PUmH9qOevZB3YJ6kN0sak3RcL+vXS3m2haSz075xt6Qv\n97qOvZLjGDlU0jpJG9NxsqIf9aybpDWSRiXdNUWZT6e4uUnS8/KadVR09vM8N7IvlPvJJjTfE9gE\nHN1W5q+onLpsAAAC70lEQVSAz6b7bwGurrNO/brl3BavA/ZK9987zNsildsHuA24Aziu3/Xu435x\nJPBTYN/0+MX9rncft8XngPek+0uBB/td75q2xWvJklPe1eX5FcA/pfsnABvyrLfuFn6eAVitc+Be\nR3aidxBNuy0i4raI+Nf0cAODO3Yh78C8jwKryWaTGFR5tsVfAn8fEU8DRMQTPa5jr+TZFmPAvun+\n/sCjPaxfz0TE7cBTUxQ5A7gqlb0T2E/SgunWW3fAzzMAa3eZiNgF7JB0YM316oeZDkY7H/hOrTXq\nn2m3RfqJuigibuplxfogz37xCuAoSbdLukPSKT2rXW/l2RYfAd4u6RHgRuD9Papb07Rvq0fJ0UCs\n6zr8cXkGYHW62HQQzyTnGowGIOltwB+QdfEMoim3hbLJQC8Hzp3mNYMgz36xB1m3zp+QjW/5gaRj\nxlv8AyTPtngrcGVEXC5pGfBlsnxdwyZ3PGlVdws/zwCsR4BDYff1+/tGxFQ/ZWarXIPRJJ0EXAK8\nKf2sHUTTbYsXkR3E6yU9CCwDrh/QE7d59ottwPURMRYRD5GNcVnSm+r1VJ5tcT5wDUBEbAD2kvTi\n3lSvUbaR4mbSMZ60qzvg/xg4UtLhkuYB5wA3tJX5NhMtubOAdTXXqV+m3RaSXg1cAZweEb/uQx17\nZcptERFPR8TBEfHyiHgZ2fmMN0XExj7Vt055jpFvAScCpOC2BHigp7XsjTzb4mHSjHqSlgLzB/ic\nhuj+y/YG4B0A6ZfOjogYnW6FtXbpRJcBWJI+Avw4Im4E1gBfkrQF+DXZhzxwcm6LjwN7A9embo2H\nI+LM/tW6Hjm3xaSXMKBdOnm2RUTcLOlkSfcAzwEfGsRfwTn3iw8BX5B0EdkJ3HO7r3H2kvRVYDlw\nkKStwCpgHhAR8fmIuEnSaZLuB54Bzsu13nRZj5mZDThPcWhmNiQc8M3MhoQDvpnZkHDANzMbEg74\nZmZDwgHfzGxIOOCbmQ0JB3wzsyHx/wHQnl3AqGyORgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe6fc233cd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(2)\n",
    "plt.subplot(2,1,1)\n",
    "plt.hist(hatp[pos_GT0],100)\n",
    "plt.subplot(2,1,2)\n",
    "plt.hist(p[pos_GT0],100)\n",
    "\n",
    "print 'Irrelevant voxels'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relevant voxels\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGfRJREFUeJzt3X/sVfV9x/HnCyhda/UrSIBYQCetretWGZk/tnXdt9Uq\n2GyapixoVGR0M9E6rakpNtnQNFlGE6MztmVuxGK7SltclbRuEAesMwVnRykWxR9zLYL6dU6+NjNp\ng/LeH+dz4fLtvdzrPZdzzvee1yP5xnM+99xzPx/v5b7uOfd+zlsRgZmZ1deEsjtgZmblchCYmdWc\ng8DMrOYcBGZmNecgMDOrOQeBmVnNdQwCSasljUja2dQ2RdJGSU9J2iBpqOm2OyU9I2mHpHlN7Usk\nPZ3uc2X/h2JmZr3o5ojgHuDCMW3LgYcj4n3AJuBmAEkLgbkR8V7gamBVap8C/BVwFnAOsKI5PMzM\nrDwdgyAiHgH2j2m+GFiTltek9Ub7vel+jwJDkmaQBcnGiHgtIkaBjcCC/N03M7O8ev2OYHpEjABE\nxEvA9NT+buD5pu32prax7ftSm5mZlazfXxarxXq0aCe1m5lZySb1eL8RSTMiYkTSTODl1L4XmN20\n3SzghdQ+PKZ9c6sdS3JAmJn1ICJafejuqNsjAnHkp/r1wFVp+Srgwab2KwEknQuMplNIG4CPSRpK\nXxx/LLW1FBED+7dixYrS++Dx9f6XXqE0DmjL7o+fP4+t8ZdHxyMCSd8g+zR/kqQ9wArgb4BvS/pT\nYA+wKP2jeEjSRZKeBV4Hlqb2/ZK+APww/Qu6NbIvjc3MrGQdgyAiLmtz0/lttv90m/avAl/ttmNm\nZlYMzywu2PDwcNldOKYGfXyDbpCfv0EeW17Ke26p3yRF1fpk1iA1fggHoNznZs36RRJxjL8sNjOz\nAeUgMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjWXKwgkfUbSTyTtlPSP\nkiZLOlXStlSb+D5Jk9K2kyWtTfWMt0qa058hmJlZHj0HgaSTgeuA+RHxQbIL2F0KrARui6ye8Siw\nLN1lGfBqZPWM7wC+mKfjZmbWH3lPDU0Ejkuf+t9BVoTmI8D96fY1wCVpubnO8TrgvJyPbWZmfdBz\nEETEC8BtZPUI9gGvAdvJitEcTJs1ahZDU93iiHgTGJU0tdfHNzOz/shzauhEsk/5pwAnA8cBC1ts\nevhSjWN2gesWm5mVrteaxZAVpnkuIl4FkPQd4PeAEyVNSEcFjZrFcLie8QuSJgInRMT+Vju+5ZZb\nDi0PDw/7OuJmZmNs2bKFLVu29GVfPdcjkHQ2sBo4C/glcA/wGPBh4J8i4puSvgL8OCJWSboG+M2I\nuEbSYuCSiFjcYr+uR2CV5XoEVlV56hHkKkwjaQWwGDgA/Aj4FNlRwFpgSmq7PCIOSHo78DXgt4H/\nBRZHxE9b7NNBYJXlILCqKi0IjgUHgVWZg8CqyhXKzMysZw4CM7OacxCYmdWcg8DMrOYcBGZmNecg\nMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzm8havH5L0bUlPStol6RxJUyRtTMXr\nN0gaatr+zlS8foekefm7b2ZmeeU9Ivhb4KGIOAM4E9gNLAceTsXrNwE3A0haCMxNxeuvBlblfGwz\nM+uDPIVpjgd2RMTcMe27gT+MiBFJM4HNEXGGpFVp+ZtpuyeB4YgYGXN/X4baKsuXobaqKusy1KcB\nr0i6R9J2SXdLeicwo/HmHhEvAdPT9oeK1yf7OFzY3szMSpKnZvEkYD5wbUT8UNLtZKeF2n1EapVU\nLbd1zWIzs6OrSs3iGcDWiDgtrX+ILAjmkk75dDg1dOgU0pj9+tSQVZZPDVlVlXJqKL2BPy/p9NR0\nHrALWA9cldquAh5My+uBKwEknQuMjg0BMzMrXt7i9WcC/wC8DXgOWApMBL4FzAb2AIsiYjRtfxew\nAHgdWBoR21vs00cEVlk+IrCqcvF6s4I4CKyqXLzezMx65iAwM6s5B4GZWc05CMzMas5BYGZWcw4C\nM7OacxCYmdWcg8DMrOYcBGZmNecgMDOrOQeBmVnN5Q4CSRNSYZr1af1USdtSzeL7JE1K7ZMlrU01\ni7dKmpP3sc3MLL9+HBFcDzzRtL4SuC3VLB4FlqX2ZcCrqWbxHcAX+/DYZmaWU64gkDQLuIjsUtQN\nHwXuT8trgEvS8sVpHWAdWf0CMzMrWd4jgtuBm0jX5ZV0ErA/Ig6m2/dyuC7xoZrFEfEmMCppas7H\nNzOznHquWSzp48BIROyQNNxo5ldrE0fTbUfsAtcsNjPrSVVqFv81cDnwBvAO4HjgAeACYGZEHEwl\nKVdExEJJ/5KWH5U0EXgxIqa32K8L01hluTCNVVVZNYs/HxFzUvH6xcCmiLgc2AwsSpst4ciaxUvS\n8iJgU6+PbWZm/XMs5hEsB26U9DQwFVid2lcD0yQ9A9yQtjMzs5K5ZrHZW+BTQ1ZVrllsZmY9cxCY\nmdWcg8DMrOYcBGZmNecgMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzmeg4CSbMk\nbZL0hKTHJf1Fap8iaWOqWbxB0lDTfe5MNYt3SJrXjwGYmVk+eY4I3gBujIjfAH4XuFbS+8muKvpw\nqlm8CbgZQNJCYG6qWXw1sCpXz83MrC/y1CN4KSJ2pOX/A54EZnFkbeI1aZ3033vT9o8CQ5Jm9Pr4\nZmbWH335jkDSqcA8YBswIyJGIAsLoFGF7FDN4mQfh+sZm5lZSXIHgaR3AeuA69ORQbsLtLe6TrYv\n5m5mVrKei9cDSJpEFgJfi4hGScoRSTMiYkTSTODl1L4XmN1091nAC6326+L1ZmZHV4ni9QCS7gVe\niYgbm9pWAq9GxEpJy4ETI2K5pIuAayPi46mo/R0RcW6LfbpCmVWWK5RZVeWpUNZzEEj6feD7wONk\n/zIC+DzwH8C3yD797wEWRcRous9dwALgdWBpRGxvsV8HgVWWg8CqqpQgOFYcBFZlDgKrKtcsNjOz\nnjkIzMxqzkFgZlZzDgIzs5pzEJiZ1ZyDwMys5hwEZmY15yAwM6s5B4GZWc05CMzMas5BYGZWcw4C\nM7OaKzwIJC2QtFvS05I+V/Tjl61f1w+vqkEf36Ab5OdvkMeWV6FBIGkCcBdwIfAB4NJU8L42Bv3F\nOOjjG3SD/PwN8tjyKvqI4GzgmYj4WUQcANZyuLi9mZmVIFepyh6MLWC/lywcSnPffffx4osvAnDD\nDTcwYYK/NjGzeim0MI2kTwIXRMSfp/XLgbMi4vqmbVzpw8ysB70Wpin6iGAvMKdp/VcK2Pc6EDMz\n603R50EeA94j6RRJk4HFwPqC+2BmZk0KPSKIiDclfRrYSBZCqyPiySL7YGZmR6pc8XozMytWaT+R\n6XZimaRPSjooaX6R/curm/FJ+hNJuyQ9LunrRfcxj07jkzRb0iZJ2yXtkLSwjH72QtJqSSOSdh5l\nmzslPZPGNq/I/uXVaXySLpP04zS2RyT9VtF9zKOb5y9td5akNyR9oqi+5dXla3NY0o8k/UTS5q52\nHBGF/5EF0LPAKcDbgB3A+1ts9y7g34AfAPPL6OuxGh/wHuA/gRPS+rSy+93n8f0dcHVaPgP477L7\n/RbG9yFgHrCzze0Lge+l5XOAbWX3uc/jOxcYSssLBm18aZsJwL8C3wU+UXaf+/jcDQG7gHen9a7e\nV8o6Iuh2YtkXgJXAL4vsXB90M74/A74UET8HiIhXCu5jHt2M7yBwQlo+EdhXYP9yiYhHgP1H2eRi\n4N607aPAkKQZRfStHzqNLyK2RcRraXUb2fyfcaOL5w/gOmAd8PKx71H/dDG2y4D7I2Jf2r6r95Wy\ngqDVxLIjXmzpcHtWRDxUZMf6pOP4gNOB96VD7x9IurCw3uXXzfhuBa6Q9DzZp67rCupbEcaOfx/j\n7M3yLfgU8M9ld6KfJJ0MXAKsAgbt5+qnA1MlbZb0mKQrurlT0fMIGlr9zz/0rbUkAbcDSzrcp6qO\nOr5kEtnpoQ+Tza34d0kfaBwhVFw347sUuCcibpd0LvB1sutLDYJuxj/uSfoIsJTsdMQguQP4XERE\n9lYzrt5bOpkEzAc+ChwHbJW0NSKe7XSnMnSaWHY82ZvGlhQKM4EHJf1xRGwvrps96zhxLm2zNSIO\nAj+V9BTwXrLvDaqum/EtI7u4IBGxTdKvSZo2zk6BtbMXmN203mr845qkDwJ3AwsiotNplvHmd4C1\n6b1lGrBQ0oGIGIQ5TXuB/4mIXwC/kPR94Eyy7/TaKuvU0FEnlkXEzyNiekScFhG/Tnae8o/GSQhA\ndxPnHiBLbSRNIwuB5wrtZe+6Gd/PgPMBJJ0BvH2chYBo/0lxPXAlQDraGY2IkaI61idtxydpDnA/\ncEVE/FehveqftuNL7yuN95Z1wDXjLASO9tp8EPgDSRMlvZPsxwwd52qVckQQbSaWSboVeCwivjv2\nLoyjw7duxhcRGyRdIGkX8Abw2fHyyavL5++zwN9L+gzZF8dL2u+xWiR9AxgGTpK0B1gBTAYiIu6O\niIckXSTpWeB1stMn40an8QF/CUwFvpw+NR+IiFIvDvlWdDG+ZuPqlF4Xr83dkjYAO4E3gbsj4omO\n+00/MTIzs5ryNZfNzGrOQWBmVnMdg6DVlGZJUyRtlPSUpA2Shppuazn1XtKSdDmCpyRd2f+hmJlZ\nL7o5IriH9DPAJsuBhyPifcAm4GaAdD2ZuRHxXuBqsgkbSJoC/BVwFtm32Cuaw8PMzMrTMQjaTGm+\nGFiTltdw+PIC7abeXwhsjIjXImKU7NcmC/J338zM8ur1O4Lpjd9NR8RLwPTU3u7SA3Wakm9mNq70\nex7B2N/6i/ZzAFr+blWuWWxm1pPosdRvr0cEI42rLUqayeEr+LWbet/NJQkOKftSr8fyb8WKFaX3\nwePr/S+9Qml8jim7P37+PLbGXx7dBsHYKc3rgavS8lVk05ob7a2m3m8APiZpKH1x/LHUZmZmJet4\naqjNlOa/Ab4t6U+BPcAigGgz9T4i9kv6AvBDso9St0b2pbGZmZWsYxBExGVtbjq/zfafbtP+VeCr\n3XZsUA0PD5fdhWNq0Mc36Ab5+RvkseVVuWsNSYqq9cmsIbsGW+P1qdznZs36RRJR8JfFZmY2IBwE\nZmY15yAwM6s5B4GZWc05CMzMas5BYGZWcw4CM7OacxCYmdWcg8DMrOYcBGZmNZcrCCR9RtJPJO2U\n9I+SJks6VdK2VJv4PkmT0raTJa1N9Yy3SprTaf9mZnbs9RwEkk4GrgPmR8QHyS5gdymwErgtsnrG\no8CydJdlwKuR1TO+A/hino6bmVl/5D01NBE4Ln3qfwdZsZmPAPen29cAl6Tl5jrH64Dzcj62mZn1\nQc9BEBEvALeR1SPYB7wGbCcrRnMwbdaoWQxNdYsj4k1gVNLUXh/fzMz6I8+poRPJPuWfApwMHAcs\nbLHp4Wv2jtkFbeoWm5lZcfIUrz8feC4iXgWQ9B3g94ATJU1IRwXNtYkb9YxfkDQROCEi9rfa8S23\n3HJoeXh42AUlzMzG2LJlC1u2bOnLvnouTCPpbGA1cBbwS+Ae4DHgw8A/RcQ3JX0F+HFErJJ0DfCb\nEXGNpMXAJRGxuMV+XZjGKsuFaayq8hSmyVWhTNIKYDFwAPgR8Cmyo4C1wJTUdnlEHJD0duBrwG8D\n/wssjoiftting8Aqy0FgVVVaEBwLDgKrMgeBVZVLVZqZWc8cBGZmNecgMDOrOQeBmVnNOQjMzGrO\nQWBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjXnIDAzqzkHgZlZzeUtXj8k6duSnpS0S9I5kqZI2piK\n12+QNNS0/Z2peP0OSfPyd9/MzPLKe0Twt8BDEXEGcCawG1gOPJyK128CbgaQtBCYm4rXXw2syvnY\nZmbWB3kK0xwP7IiIuWPadwN/GBEjkmYCmyPiDEmr0vI303ZPAsMRMTLm/r4MtVWWL0NtVVXWZahP\nA16RdI+k7ZLulvROYEbjzT0iXgKmp+0PFa9P9nG4sL2ZmZUkT83iScB84NqI+KGk28lOC7X7iNQq\nqVpu65rFZmZHV5WaxTOArRFxWlr/EFkQzCWd8ulwaujQKaQx+/WpIassnxqyqirl1FB6A39e0ump\n6TxgF7AeuCq1XQU8mJbXA1cCSDoXGB0bAmZmVry8xevPBP4BeBvwHLAUmAh8C5gN7AEWRcRo2v4u\nYAHwOrA0Ira32KePCKyyfERgVeXi9WYFcRBYVbl4vZmZ9cxBYGZWcw4CM7OacxCYmdWcg8DMrOYc\nBGZmNecgMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmssdBJImpMI069P6qZK2pZrF90malNonS1qb\nahZvlTQn72ObmVl+/TgiuB54oml9JXBbqlk8CixL7cuAV1PN4juAL/bhsc3MLKdcQSBpFnAR2aWo\nGz4K3J+W1wCXpOWL0zrAOrL6BWZmVrK8RwS3AzeRrssr6SRgf0QcTLfv5XBd4kM1iyPiTWBU0tSc\nj29mZjn1XLNY0seBkYjYIWm40cyv1iaOptuO2AWuWWxm1pOq1Cz+a+By4A3gHcDxwAPABcDMiDiY\nSlKuiIiFkv4lLT8qaSLwYkRMb7FfF6axynJhGquqsmoWfz4i5qTi9YuBTRFxObAZWJQ2W8KRNYuX\npOVFwKZeH9vMzPrnWMwjWA7cKOlpYCqwOrWvBqZJega4IW1nZmYlc81is7fAp4asqlyz2MzMeuYg\nMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjXnIDAzqzkHgZlZzfUcBJJm\nSdok6QlJj0v6i9Q+RdLGVLN4g6ShpvvcmWoW75A0rx8DMDOzfPIcEbwB3BgRvwH8LnCtpPeTXVX0\n4VSzeBNwM4CkhcDcVLP4amBVrp6bmVlf5KlH8FJE7EjL/wc8CcziyNrEa9I66b/3pu0fBYYkzej1\n8c3MrD/68h2BpFOBecA2YEZEjEAWFkCjCtmhmsXJPg7XMzYzs5LkDgJJ7wLWAdenI4N2F2hvdZ1s\nX8zdzKxkPRevB5A0iSwEvhYRjZKUI5JmRMSIpJnAy6l9LzC76e6zgBda7dfF683Mjq4SxesBJN0L\nvBIRNza1rQRejYiVkpYDJ0bEckkXAddGxMdTUfs7IuLcFvt0hTKrLFcos6rKU6Gs5yCQ9PvA94HH\nyf5lBPB54D+Ab5F9+t8DLIqI0XSfu4AFwOvA0ojY3mK/DgKrLAeBVVUpQXCsOAisyhwEVlWuWWxm\nZj1zEJiZ1ZyDwMys5hwEZmY15yAwM6s5B4GZWc05CMzMas5BYGZWcw4CM7OacxCYmdWcg8DMrOYK\nDwJJCyTtlvS0pM8V/fhmZnakQoNA0gTgLuBC4APApanOcW306/rhVTXo4xt0g/z8DfLY8ir6iOBs\n4JmI+FlEHADWcrimcS0M+otx0Mc36Ab5+RvkseVVdBCMrVu8F9ctNjMrVa5SlT2oXN3im266id27\ndwPwwAMPMHHixDK7Y2ZWuEIL06QSlbdExIK0vhyIiFjZtI0rfZiZ9WBcVCiTNBF4CjgPeJGsrOWl\nEfFkYZ0wM7MjFHpqKCLelPRpYCPZ9xOrHQJmZuWqXM1iMzMrVmkzi7udWCbpk5IOSppfZP/y6mZ8\nkv5E0i5Jj0v6etF9zKPT+CTNlrRJ0nZJOyQtLKOfvZC0WtKIpJ1H2eZOSc+ksc0rsn95dRqfpMsk\n/TiN7RFJv1V0H/Po5vlL250l6Q1Jnyiqb3l1+doclvQjST+RtLmrHUdE4X9kAfQscArwNmAH8P4W\n270L+DfgB8D8Mvp6rMYHvAf4T+CEtD6t7H73eXx/B1ydls8A/rvsfr+F8X0ImAfsbHP7QuB7afkc\nYFvZfe7z+M4FhtLygkEbX9pmAvCvwHeBT5Td5z4+d0PALuDdab2r95Wyjgi6nVj2BWAl8MsiO9cH\n3Yzvz4AvRcTPASLilYL7mEc34zsInJCWTwT2Fdi/XCLiEWD/UTa5GLg3bfsoMCRpRhF964dO44uI\nbRHxWlrdxjib69PF8wdwHbAOePnY96h/uhjbZcD9EbEvbd/V+0pZQdBxYlk63J4VEQ8V2bE+6Wbi\n3OnA+9Kh9w8kXVhY7/LrZny3AldIep7sU9d1BfWtCGPHv49x9mb5FnwK+OeyO9FPkk4GLgFW0Xpu\n03h2OjBV0mZJj0m6ops7FT2hrOGoE8skCbgdWNLhPlXVzcS5SWSnhz4MzAH+XdIHGkcIFdfN+C4F\n7omI29P8ka+TXV9qEFRuYuSxIOkjwFKy0xGD5A7gcxER2VvNuHpv6WQSMB/4KHAcsFXS1oh4ttOd\nyrCX7M2vYRbwQtP68WRvGltSKMwEHpT0xxGxvbhu9qzT+BrbbI2Ig8BPJT0FvJfse4Oq62Z8y8gu\nLkhEbJP0a5KmjbNTYO3sBWY3rbca/7gm6YPA3cCCiOh0mmW8+R1gbXpvmQYslHQgItaX3K9+2Av8\nT0T8AviFpO8DZ5J9p9dWWaeGHgPeI+kUSZOBxcChJyEifh4R0yPitIj4dbLzlH80TkIAOowveYAs\ntZE0jSwEniu0l73rZnw/A84HkHQG8PZxFgKi/SfF9cCVcGi2/GhEjBTVsT5pOz5Jc4D7gSsi4r8K\n7VX/tB1fel9pvLesA64ZZyFwtNfmg8AfSJoo6Z1kP2boOFerlCOCaDOxTNKtwGMR8d2xd2EcHb51\nM76I2CDpAkm7gDeAz46XT15dPn+fBf5e0mfIvjhe0n6P1SLpG8AwcJKkPcAKYDLZ5VDujoiHJF0k\n6VngdbLTJ+NGp/EBfwlMBb6cPjUfiIizy+rvW9XF+JqNq1N6Xbw2d0vaAOwE3gTujognOu43/cTI\nzMxqyqUqzcxqzkFgZlZzDgIzs5pzEJiZ1ZyDwMys5hwEZmY15yAwM6s5B4GZWc39P2D6Cty3eUjw\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe6d4d80750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(3)\n",
    "plt.subplot(2,1,1)\n",
    "plt.hist(hatp[pos_GT1],100)\n",
    "plt.subplot(2,1,2)\n",
    "plt.hist(p[pos_GT1],100)\n",
    "\n",
    "print 'Relevant voxels'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Voxel selection**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Irrelevant voxels selected:\n",
      "1077\n",
      "Total Irrelevant voxels:10000\n",
      "Relevant voxels selected:\n",
      "1000\n",
      "Total Relevant voxels:1000\n"
     ]
    }
   ],
   "source": [
    "# Set the confidence level\n",
    "alp = 0.05\n",
    "\n",
    "# It's a two tailed distribution, we select the voxels of the lower and upper tails\n",
    "print 'Irrelevant voxels selected:'\n",
    "print np.where(np.any([p[pos_GT0]<0.5*alp, p[pos_GT0]>(1-0.5*alp)], axis=0))[0].shape[0]\n",
    "print 'Total Irrelevant voxels:' + str(pos_GT0.shape[0]) \n",
    "print 'Relevant voxels selected:'\n",
    "print np.where(np.any([p[pos_GT1]<0.5*alp, p[pos_GT1]>(1-0.5*alp)], axis=0))[0].shape[0]\n",
    "print 'Total Relevant voxels:' + str(pos_GT1.shape[0]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Transductive refinement of the voxel importance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create information for conformal analysis**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Training data\n",
    "Ndata = 200\n",
    "Drel = 1000\n",
    "Dirrel = 10000\n",
    "\n",
    "X, Y, GT = generate_data(Ndata, Drel, Dirrel)\n",
    "\n",
    "# Create test data\n",
    "Ndata = 100\n",
    "X_test, Y_test, GT_test = generate_data(Ndata, Drel, Dirrel)\n",
    "\n",
    "pos_GT1 = np.where(GT==1)[0]\n",
    "pos_GT0 = np.where(GT==0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "conformal_iterations = 10\n",
    "n_data_test = X_test.shape[1]\n",
    "# The number of conformal samples is set to 1% of the training data\n",
    "n_data_conf = int(0.01* X.shape[1])\n",
    "data = []\n",
    "labels = []\n",
    "for i in range(conformal_iterations):\n",
    "    temp = []\n",
    "    for j in range(n_data_conf):\n",
    "        temp.append(random.choice([-1,1]))\n",
    "    labels.append(temp)\n",
    "    data.append(random.sample(range(n_data_test), n_data_conf))\n",
    "conformal={}\n",
    "conformal['data'] = np.array(data)\n",
    "conformal['labels'] = np.array(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Train an ensemble of bagged SVMs with conformal analysis**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "\n",
    "num_bagging_iter = 10000\n",
    "perc = 0.5\n",
    "\n",
    "def compute_bagging_clf_conf(X, Y, X_test, n_conf): \n",
    "    \n",
    "    w_sign_tot = np.zeros((X.shape[0],))\n",
    "    \n",
    "    ntr = X.shape[1]\n",
    "    pos_tr = np.arange(ntr)\n",
    "    p_1 = np.where(Y == 1)[0]\n",
    "    p_0 = np.where(Y == -1)[0]\n",
    "    n_1 = p_1.shape[0]\n",
    "    n_0 = p_0.shape[0]\n",
    "       \n",
    "    Xall = np.hstack((X, X_test))\n",
    "    K = np.dot(Xall.T,Xall)\n",
    "    clf = SVC(C=100, kernel='precomputed')\n",
    "        \n",
    "    for i in range(num_bagging_iter):\n",
    "        # Select training samples (at random) and balancing by clases\n",
    "        np.random.seed(i)\n",
    "        # In case of unbalanced problems, we balance the subsampling\n",
    "        samples_sel_1 = np.random.permutation(p_1)[:int(np.round(0.5*perc*ntr))]\n",
    "        samples_sel_0 = np.random.permutation(p_0)[:int(np.round(0.5*perc*ntr))]\n",
    "        samples_sel = np.hstack((samples_sel_0, samples_sel_1))\n",
    "        # Create a vector with the training positions\n",
    "        pos_tr_sel = pos_tr[samples_sel]\n",
    "       \n",
    "        # For the conformal learning\n",
    "        Y_sel = np.hstack((Y[pos_tr_sel], conformal['labels'][n_conf]))\n",
    "        pos_tr_sel = np.hstack((pos_tr_sel, (ntr-1)+ conformal['data'][n_conf])).astype('int')\n",
    "        \n",
    "        # Train the classifier\n",
    "        clf.fit(K[pos_tr_sel[:, np.newaxis],pos_tr_sel], Y_sel) \n",
    "        w = np.dot(Xall[:,pos_tr_sel[clf.support_]],clf.dual_coef_.ravel())\n",
    "       \n",
    "        # Create a vector counting the number of positive weigths (indirectly we are also counting the negatives)\n",
    "        w_sign = np.zeros(w.shape)\n",
    "        w_sign [w>0] = 1\n",
    "        # Add it to the total sum of positive weigths\n",
    "        w_sign_tot = w_sign_tot + w_sign\n",
    "    return w_sign_tot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_w_count = np.zeros((conformal_iterations, X.shape[0]))\n",
    "for n_conf in range(conformal_iterations):\n",
    "    all_w_count[n_conf,:] = compute_bagging_clf_conf(X, Y, X_test, n_conf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Voxel selection applying the conformal with the intersection of the selected subsets**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "nPositives =all_w_count\n",
    "\n",
    "N = num_bagging_iter\n",
    "hatp = nPositives/N\n",
    "z, p = computeZTest(hatp, N)\n",
    "\n",
    "p1= np.where(p>0.5)\n",
    "paux = copy.copy(p)\n",
    "paux[p1] = 1-paux[p1]\n",
    "# If p<0.5, we get the maximum; and if p>0.5, we get the minimum\n",
    "pos = np.argmax(paux, axis=0)\n",
    "p = np.array([pp[pos[i]] for i, pp in enumerate(p.T)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Irrelevant voxels selected:\n",
      "29\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFfNJREFUeJzt3X+MZeV93/H3B7AXY5sNtsOOtPzYmGBCozaLI8XbNio3\nseV4XSX4j9DYkoshpKVV0qSJ0hbyR9moaiVbipwgt9qmodaS1AbsyGbjkkDQMmmtGOwEtkY24HUc\nYCHeaWy8sQ0Uxd5v/5izO3dmZ+ae2bk/5t7zfkmjPffMc888c/bezzznuc/znFQVkqTZd9akKyBJ\nGg8DX5I6wsCXpI4w8CWpIwx8SeoIA1+SOmJg4Cd5U5JHkzzS/Ps3SX4xyQVJ7k/yZJL7kmzve85t\nSY4kOZxk92h/BUlSGwMDv6q+VFVXVdWbgR8GXgA+AdwMPFBVVwCHgFsAkuwFLquqy4GbgP2jqrwk\nqb2Ndum8DfiLqjoKXAMcaPYfaB7T/HsHQFU9DGxPsmMIdZUkbcJGA/9ngI802zuqagGgqo4BFzb7\ndwJH+57zXLNPkjRBrQM/ySuAnwI+1uxaa02GrLLP9RskacLO2UDZvcCfV9XXmscLSXZU1UKSOeD/\nNvufBS7ue95FwF+tPFgS/whI0hmoqtUa1gNtpEvnPcBH+x4fBK5vtq8H7unbfx1Akj3A8ZNdPytV\nlV9V3HrrrROvw1b58lx4LjwX639tRqsWfpJXsfiB7T/v2/1+4O4kPws8A1zbhPi9Sd6Z5Mssjui5\nYVM1lCQNRavAr6qXgO9dse95Fv8IrFb+FzZfNUnSMDnTdgvo9XqTrsKW4blY4rlY4rkYjmy2T+iM\nf3BSk/rZkjStklBj+NBWkjTFDHxJ6ggDX5I6wsCXpI4w8CWpIwx8SeoIA1+SOsLAl6SOMPAlqSMM\nfEnqCANfkjrCwJekjjDwJakjDHwJmJvbRZJTX3NzuyZdJWnoXB5ZYnHJWeh/PWbTt5OTRsHlkSVJ\nAxn4ktQRBr4kdYSBL0kdYeBLUkcY+JLUEa0CP8n2JB9L8niSLyR5S5ILktyf5Mkk9yXZ3lf+tiRH\nkhxOsnt01ZcktdW2hf9bwL1VdSXwQ8ATwM3AA1V1BXAIuAUgyV7gsqq6HLgJ2D/0WkuSNmzgxKsk\nrwUOV9VlK/Y/AVxdVQtJ5oAHq+rKJPub7buaco8DvapaWPF8J15py3DilabFqCdevRH4WpIPJ3kk\nyW8nOQ/YcTLEq+oYcGFTfidwtO/5zzX7JEkTdE7LMm8Gfr6q/izJB1nszlmr+bPaX55Vy+7bt+/U\ndq/Xo9frtaiOJHXH/Pw88/PzQzlWmy6dHcBnquqNzeMfZTHwL6PpqhnQpXOq62fFce3S0ZZhl46m\nxUi7dJqgPprkTc2utwJfAA4C1zf7rgfuabYPAtc1FdsDHF8Z9pKk8Wu1WmaSHwJ+B3gF8BXgBuBs\n4G7gYuAZ4NqqOt6U/xDwDuAF4IaqemSVY9rC15ZhC1/TYjMtfJdHljDwNT1cHlmSNJCBL0kdYeBL\nUkcY+JLUEQa+JHWEgS+tahtJSMLc3K5JV0YaCodlSqw+LHPpsUM0tXU4LFOSNJCBL0kdYeBLUkcY\n+JLUEQa+JHWEgS9JHWHgS1JHGPiS1BEGviR1hIGvzpqb23Vq+QSpC1xaQZ21fDkFl1bQdHBpBUnS\nQAa+JHWEgS9JHWHgS1JHGPiS1BGtAj/JU0n+T5JHk3y22XdBkvuTPJnkviTb+8rfluRIksNJdo+q\n8pKk9tq28E8Avaq6qqp+pNl3M/BAVV0BHAJuAUiyF7isqi4HbgL2D7nOkqQz0Dbws0rZa4ADzfaB\n5vHJ/XcAVNXDwPYkOzZZT0nSJrUN/ALuS/K5JD/X7NtRVQsAVXUMuLDZvxM42vfc55p9kqQJOqdl\nuX9QVceSfC9wf5InWT4tsd9qM8BWLbtv375T271ej16v17I6ktQN8/PzzM/PD+VYG15aIcmtwLeB\nn2OxX38hyRzwYFVdmWR/s31XU/4J4OqTVwN9x3FpBU1U+6UVzgVeBmDHjks5duypMdVQOt1Il1ZI\ncl6S1zTbrwbeDjwGHASub4pdD9zTbB8ErmvK7wGOrwx7abq8zGL4FwsLT0+6MtIZa9OlswP4RJJq\nyv+Pqro/yZ8Bdyf5WeAZ4FqAqro3yTuTfBl4AbhhRHWXJG2Aq2WqszayWqYrZ2qrcLVMSdJABr4k\ndYSBL0kdYeBLUkcY+JLUEQa+JHWEgS9JHWHgS1JHGPiS1BEGviR1hIGvTpmb20WSZlkFqVtcS0ed\nsvb6Oa6lo+ngWjqSpIEMfEnqCANfkjrCwJekjjDwJakjDHxpQ7adGtaZhLm5XZOukNSawzLVKcMY\nlrmynK9jjZPDMqU19E+0crKVus4Wvmba8hY9bKTlbgtfW5EtfEnSQAa+JHVE68BPclaSR5IcbB7v\nSvJQkieTfDTJOc3+Vya5M8mRJJ9JcsmoKi9Jam8jLfxfAr7Y9/j9wG9U1RXAceDGZv+NwPNVdTnw\nm8AHhlFRSdLmtAr8JBcB7wR+p2/3jwO/32wfAN7VbF/TPAb4OPDWzVdTkrRZbVv4HwT+Dc3whCSv\nB75RVSea7z8L7Gy2dwJHAarqu8DxJK8bWo0lSWfknEEFkvxjYKGqDifpndzdfPXrH7e27BAsH8d2\nyr59+05t93o9er3easUkqbPm5+eZn58fyrEGjsNP8p+A9wLfAV4FvBb4JPB2YK6qTiTZA9xaVXuT\n/FGz/XCSs4GvVtWFqxzXcfgaOcfha9aMdBx+Vf1aVV1SVW8E3g0cqqr3Ag8C1zbF3gfc02wfbB7T\nfP/QmVRMkjRcmxmHfzPwK0m+BLwOuL3ZfzvwhiRHgH/dlJMkTZhLK2im2aWjWePSCpKkgQx8SeoI\nA1+SOsLAl6SOMPAlqSMMfGlTtnl/W00Nh2Vqpo1jWGb/93xNa9QclilJGsjAl6SOMPAlqSMMfEnq\nCANfGhpH7GhrM/A1c+bmdp0K3vF6mcURO8XCwtNj/tnSYA7L1MxZPhRzvMMyHaKpUXNYpiRpIANf\nkjrCwJekjjDwJakjDHxJ6ggDX5I6wsDX1Osfdz/+sffS9HAcvqZe+yWQ1/ue4/A1HRyHL0kaaGDg\nJ9mW5OEkjyZ5LMmtzf5dSR5K8mSSjyY5p9n/yiR3JjmS5DNJLhn1LyFJGmxg4FfVy8CPVdVVwG5g\nb5K3AO8HfqOqrgCOAzc2T7kReL6qLgd+E/jASGouSdqQVl06VfVis7kNOIfFjsofA36/2X8AeFez\nfU3zGODjwFuHUlNJ0qa0CvwkZyV5FDgG/DHwF8DxqjrRFHkW2Nls7wSOAlTVd4HjSV431FpLkjbs\nnDaFmmC/Ksn5wCeAK1cr1vy78tPjlcMaTtm3b9+p7V6vR6/Xa1MdSeqM+fl55ufnh3KsDQ/LTPLv\ngReBfwvMVdWJJHuAW6tqb5I/arYfTnI28NWqunCV4zgsU0PhsEx1yUiHZSZ5Q5LtzfargLcBXwQe\nBK5tir0PuKfZPtg8pvn+oTOpmCRpuAa28JP8XRY/hD2r+bqrqv5jku8D7gQuAB4F3ltVf5tkG/C7\nwFXA14F3V9VTqxzXFr6Gwha+umQzLXxn2mrqGfjqEmfaSlvOtmXr+3hTc20FtvA19bZqC39lOV/v\nGgZb+NKWt83WvibOFr6m3rS08O3f1zDYwpckDWTgS1JHGPiS1BEGviR1hIEvSR1h4EtSRxj4ktQR\nBr4kdYSBr6kxN7fL2arSJjjTVlNj+Yzapdmq0zfT9lzgZQB27LiUY8eeOv2XldawmZm2rW5xKGmY\nXuZk+C8snNH7VjojBr6m1LamZS+pLfvwNaVOtpLtFpTaMvAlqSMMfEnqCANfkjrCwJekjjDwJakj\nDHxJ6oiBgZ/koiSHknwxyWNJfrHZf0GS+5M8meS+JNv7nnNbkiNJDifZPcpfQJLUTpsW/neAX6mq\nvwP8feDnk/wAcDPwQFVdARwCbgFIshe4rKouB24C9o+k5pp5/WvnOMlK2ryBgV9Vx6rqcLP9beBx\n4CLgGuBAU+xA85jm3zua8g8D25PsGHK91QELC0+zNLnKCVbSZm2oDz/JLmA38BCwo6oWYPGPAnBh\nU2wncLTvac81+yRJE9R6LZ0krwE+DvxSVX07yVpNrtWuvVctu2/fvlPbvV6PXq/XtjqS1Anz8/PM\nz88P5VitlkdOcg7wKeAPq+q3mn2PA72qWkgyBzxYVVcm2d9s39WUewK4+uTVQN8xXR5Zp5mb29V0\n5Zw0zCWLh3GM4ZfzfaCN2MzyyG27dP478MWTYd84CFzfbF8P3NO3/7qmYnuA4yvDXlrL8n57ScM0\nsIWf5B8C/wt4jKV34q8BnwXuBi4GngGurarjzXM+BLwDeAG4oaoeWeW4tvB1mpU3OZn9Fv7SzVAA\nzjrrPE6ceBHw5iha3WZa+N7xSltK9wJ//XK+R7TSOLp0JElTzsCXpI4w8CWpIwx8SeoIA1+SOsLA\nl6SOMPAlqSMMfE1c/zLIkkbHwNfEuZzCWrad+kM4N7dr0pXRDGi9WqakcXuZk38EFxa8+tHm2cKX\npI4w8CWpIwx8SeoIA1+SOsLA19j1D8N0KOaZ6T+HjuBRW66Hr7FbvuY9jG/t+XH+rNHeCnHlfQN8\nL3WH6+FLkgYy8CWpIwx8jYXLJ0iTZx++xmK096rtQh/+8pudL7IPv4s204fv0grSVFhaZmGRV0ra\nOLt0NDJ240hbi106Gpm1u3Hs0hl2Od9L3THSYZlJbk+ykOTzffsuSHJ/kieT3Jdke9/3bktyJMnh\nJLvPpFKSpOFr06XzYeAnVuy7GXigqq4ADgG3ACTZC1xWVZcDNwH7h1hXSdImDAz8qvo08I0Vu68B\nDjTbB5rHJ/ff0TzvYWB7kh3Dqaq2OpdMmJSlG6Wcffarl/0fuOyC+p3pKJ0Lq2oBoKqOJbmw2b8T\nONpX7rlm38KZV1HTYunOVScZ+uOxNILnxInlff3eOEX9hj0sc7VX15qfJu3bt+/Udq/Xo9frDbk6\nkjTd5ufnmZ+fH8qxWo3SSXIp8AdV9feax48DvapaSDIHPFhVVybZ32zf1ZR7Arj65NXAimM6SmfG\nTG5RNEfprFfO99lsGcfiaWF56/0gcH2zfT1wT9/+65pK7QGOrxb2mh2OtZemx8AWfpKPAD3g9Sz2\nxd8KfBL4GHAx8AxwbVUdb8p/CHgH8AJwQ1U9ssZxbeHPgK2xZIItfFv43bGZFr4Tr7QpBv5WL7e0\nBs+OHZdy7NhTaLq5Hr6kNZwcwVPNKKpFK4fQOnyzG2zha1Ns4U9XuZPvudU+YPf9OB1s4UuSBjLw\ntWGOzJlF2+ze6QC7dLRhW28VTLt0htGls1o5bT126WikXCNHmg3e8UoDuUaONBsMfKkztnmF1nF2\n6WhVfjA7i5bG5KubbOHrlLm5Xcsm5yz/IFDStLOFr1OW+uptAUqzyMCXtMLaY/L7u/ocrz99DPwZ\n5Dop2pzV19+B5VeBK7+nrc/AnxH9Ib+8a2btRbNW3v9UOt02XyMzxMCfEW373/vLnTjxIv1/GKTT\n9Y/sWfka2eaV5JRxlE4nOP5ao7B083TwhunTwBZ+Jzj+WuPgAmxbnYE/xZwcpa2l/8PeY4b/FmTg\nb3HrDYNz3Ly2rtXD3z8Ak2Xgb0FrjbhxGJym0/IPftu2/tcaUeYfjDPnevhb0Hrrzfefs+ldl971\n8Gev3JkeY+km62eddV4zcqzf4PdB17ge/hRYORnqzFosjonWrFlq/bcfJuxw0DM1ksBP8o4kTyT5\nUpJ/N4qfsVWt1ee+cjJU/4t7ZR/n2tYbEy11xcouotUnFvqH4HRD79JJchbwJeCtwF8BnwPeXVVP\nrCg3tV06/atKrn8ZunS5unw/bI3L8K6V24p1mpVyk63T6rduXP7+27HjUo4de4ppt9W6dH4EOFJV\nT1fV3wJ3AtesVvC88y7gvPMu4KKLLudb3/rWCKqyMW1bB+1nqzr+XRq9bWtcHa/9YXF/l+rKJUbO\npLt11FcWwxqCPYrA3wkc7Xv8bLPvNC+99BVeeukrPP/8t5YF/nrrvax1Mtv2ka9Xbr0RMY55l7aq\ntg2r1T8vWNloW6u7db0/DG1H07XNtpU5Nawh2KNYWmG1RFy1lueffx0AL774N+zevYe//uujpz3l\nxInll28LC+eeCt3Tu1OWyvU/r/8565VbXvXVliNYrZyk2bW0fMTKLGqbHWvl1HrZ1l/u9OOfuVEE\n/rPAJX2PL2KxL/803/zmp05tLw97WP4Lrv7Lnt53vrLcWidpvXLrndi1yrU93kZ+7rh+VtfKbcU6\nzUq5rVinSZ6LRevnVNu8GVS2nVF8aHs28CSLH9p+Ffgs8J6qenyoP0iStCFDb+FX1XeT/AJwP4uf\nEdxu2EvS5E1spq0kabxGPtN20CSsJK9McmeSI0k+k+SS1Y4zC1qci19O8oUkh5P8cZKLJ1HPcWg7\nOS/JTyc5keTN46zfOLU5F0n+SfPaeCzJ7427juPS4j1ycZJDSR5p3id7J1HPUUtye5KFJJ9fp8xt\nTW4eTrK71YGramRfLP5B+TJwKfAK4DDwAyvK/EvgvzTbPwPcOco6Teqr5bm4Gji32f4XXT4XTbnX\nAH8C/Cnw5knXe4Kvi+8H/hw4v3n8hknXe4Ln4r8CNzXbVwJ/Oel6j+hc/CiwG/j8Gt/fC/zPZvst\nwENtjjvqFn6bSVjXAAea7Y+z+GHvLBp4LqrqT6rq/zUPH2KN+QszoO3kvP8AvJ/l05VnTZtz8c+A\n/1xV3wSoqq+NuY7j0uZcnADOb7a/B3hujPUbm6r6NPCNdYpcA9zRlH0Y2J5kx6Djjjrw20zCOlWm\nqr4LHE/yuhHXaxJaT0hr3Aj84UhrNDkDz0VziXpRVd07zopNQJvXxZuAK5J8OsmfJvmJsdVuvNqc\ni18H/mmSo8CngH81prptNSvP1XO0aCCO+p62bSZhrTbYdBY/SW49IS3Je4EfZrGLZxatey6yOPvk\ng8D7BjxnFrR5XZzDYrfOP2Jxjsv/TvKDJ1v8M6TNuXgP8OGq+mCSPcDvAT848pptPa3zpN+oW/ht\nJmEdBS6GU2P4z6+q9S5lplWrCWlJ3gbcAvxkc1k7iwadi9ey+CaeT/KXwB7gnhn94LbN6+JZ4J6q\nOlFVT7E4z+Xy8VRvrNqcixuBuwGq6iHg3CRvGE/1tpRnaXKzseYE136jDvzPAd+f5NIkrwTeDRxc\nUeYPWGrJXQscGnGdJmXguUhyFbAf+Kmq+voE6jgu656LqvpmVV1YVW+squ9j8fOMn6yqRyZU31Fq\n8x75JPDjAE24XQ58Zay1HI825+Jp4G0ASa4Ets3wZxph7Svbg8B1AM2VzvGqWhh0wJF26dQak7CS\n/Drwuar6FHA78LtJjgBfZ/E/eea0PBcfAF4NfKzp1ni6qt41uVqPRstzsewpzGiXTptzUVX3JXl7\nki8A3wF+dRavglu+Ln4V+G9JfpnFD3Dft/YRp1eSjwA94PVJngFuBV4JVFX9dlXdm+SdSb4MvADc\n0Oq4zbAeSdKM8xaHktQRBr4kdYSBL0kdYeBLUkcY+JLUEQa+JHWEgS9JHWHgS1JH/H8+wuGbaVF5\n7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe6d4c2f4d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(4)\n",
    "plt.hist(p[pos_GT0],100)\n",
    "\n",
    "alp = 0.05\n",
    "# It's a two tailed distribution, we select the voxels of the lower and upper tails\n",
    "print 'Irrelevant voxels selected:'\n",
    "print  np.where(np.any([p[pos_GT0]<0.5*alp, p[pos_GT0]>(1-0.5*alp)], axis=0))[0].shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relevant voxels selected:\n",
      "995\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEF9JREFUeJzt3X+sZGV9x/H3B1aqFllRwhJ3QVLAdW1TEdO6xrZOSw0/\nrEJRWjAqILUk/mw1xrVpu5fYGDGhIsGWWqhCWwXEJhB/FEKXqyEKRtcFRH6s8cfugq6xsLb6B1H3\n2z/uWRhu73LnztyduXef9yu54TnPPOfMd04m53Oec+YsqSokSe06YNIFSJImyyCQpMYZBJLUOINA\nkhpnEEhS4wwCSWrcvEGQ5MokO5Pc1dd3aJKbk9yf5KYkK/teuzTJ1iRbkhzf139Okge6dd64+B9F\nkjSMQWYEHwdOmtW3AbilqtYCm4D3ASQ5BTimqo4DLgAu7/oPBf4W+C3gJcDG/vCQJE3OvEFQVbcB\nj8zqPg24qmtf1S3v6b+6W+8OYGWSVcwEyc1V9ZOq2gXcDJw8evmSpFENe4/g8KraCVBVPwQO7/pX\nA9v7xu3o+mb3P9j1SZImbLFvFmeO5Zqjn65fkjRhK4Zcb2eSVVW1M8kRwI+6/h3AkX3j1gAPdf29\nWf23zrXhJAaEJA2hquY66Z7XoDOC8MSz+huBc7v2ucANff1vBEiyHtjVXUK6CXhFkpXdjeNXdH1z\nqqpl+7dx48aJ12D9k6+jxfqXc+3Ltf5Vq5474CH8yc07I0jySWbO5p+dZBuwEfgg8OkkbwK2AWd2\nB/DPJzk1ybeBnwHndf2PJHk/8DVmLgldWDM3jSVJQ9q58/s8fpV9qMkAMEAQVNXr9vLSH+5l/Nv2\n0v8J4BODFiZJGg+fLF5kvV5v0iWMxPonaznXv5xrh+Vf/yhStbTuzSappVaTJC1FyZ4fZgKE2sc3\niyVJ+ymDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxB\nIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS\n1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGjdSECT5yyTfTHJXkn9P\nclCSo5PcnuT+JJ9KsqIbe1CSa5JsTfKVJEctzkeQJI1i6CBI8hzg7cAJVfWbwArgbOAi4OKqWgvs\nAs7vVjkfeLiqjgMuAT40SuGSpMUx6qWhA4Ff7c76nwY8BPw+8Jnu9auA07v2ad0ywPXAiSO+tyRp\nEQwdBFX1EHAxsA14EPgJsBnYVVW7u2E7gNVdezWwvVv3l8CuJM8a9v0lSYtjxbArJnkmM2f5z2Um\nBD4NnDLH0NqzyuxN9L32BFNTU4+1e70evV5v2DIlab80PT3dtaZG3laq5jwWz79i8lrgpKp6c7f8\nBuClwGuBI6pqd5L1wMaqOiXJf3btO5IcCPygqg6fY7s1bE2S1JKk/3w6VNXsE+6BjHKPYBuwPslT\nM1PNicA9wK3Amd2Yc4AbuvaN3TLd65tGeG9J0iIZekYAkGQjcBbwc+AbwJ8Ba4BrgEO7vtdX1c+T\n/Arwr8CLgP8Gzqqq782xTWcEkjSAxZoRjBQE+4JBIEmDWQqXhiRJ+wGDQJIaZxBIUuMMAklqnEEg\nSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLU\nOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0z\nCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjRgqCJCuTfDrJvUnuSfKSJIcmuTnJ/UluSrKyb/yl\nSbYm2ZLk+NHLlySNatQZwUeAz1fVOuCFwH3ABuCWqloLbALeB5DkFOCYqjoOuAC4fMT3liQtglTV\ncCsmzwC2VNUxs/rvA15eVTuTHAHcWlXrklzeta/txt0L9Kpq56z1a9iaJKklSYA9x8tQVRlmO6PM\nCH4N+HGSjyfZnORjSZ4OrNpzcK+qHwKHd+NXA9v71n+w65MkTdAoQbACOAH4aFWdAPyMmctCezud\nnyupPPWXpAlbMcK6O4DtVfW1bvkzzATBziSr+i4N/ahv/JF9668BHpprw1NTU4+1e70evV5vhDIl\naf8zPT3dtaZG3tbQ9wgAknwReHNVPZBkI/D07qWHq+qiJBuAZ1bVhiSnAm+tqlcmWQ9cUlXr59im\n9wgkaQCLdY9g1CB4IXAF8BTgO8B5wIHAdcyc/W8DzqyqXd34y4CTmbmMdF5VbZ5jmwaBJA1gSQTB\nvmAQSNJglsKvhiRJ+wGDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQ\nSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEk\nNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLj\nDAJJatzIQZDkgCSbk9zYLR+d5PYk9yf5VJIVXf9BSa5JsjXJV5IcNep7S5JGtxgzgncC3+pbvgi4\nuKrWAruA87v+84GHq+o44BLgQ4vw3pKkEY0UBEnWAKcCV/R1/wHwma59FXB61z6tWwa4HjhxlPeW\nJC2OUWcEHwbeAxRAkmcDj1TV7u71HcDqrr0a2A5QVb8EdiV51ojvL0ka0dBBkOSVwM6q2gJkT3df\ne4/qe+0Jm+h7TZI0IStGWPdlwKuTnAo8DXgGM9f+VyY5oJsVrAEe6sbvAI4EHkpyIHBIVT0y14an\npqYea/d6PXq93ghlStL+Z3p6umtNjbytVI1+Up7k5cC7q+rVSa4F/qOqrk3yj8CdVXV5krcAv1FV\nb0lyFnB6VZ01x7ZqMWqSpP1d0n9hJVTV7CsvA9kXzxFsAN6V5AHgWcCVXf+VwGFJtgJ/0Y2TJE3Y\noswIFpMzAkkazFKeEUiSlhGDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CS\nGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlx\nBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQ\nSFLjDAJJatzQQZBkTZJNSb6V5O4k7+j6D01yc5L7k9yUZGXfOpcm2ZpkS5LjF+MDSJJGM8qM4BfA\nu6rqBcBLgbcmeT6wAbilqtYCm4D3ASQ5BTimqo4DLgAuH6lySdKiGDoIquqHVbWla/8UuBdYA5wG\nXNUNu6pbpvvv1d34O4CVSVYN+/6SpMWxKPcIkhwNHA/cDqyqqp0wExbA4d2w1cD2vtUe7PokSRO0\nYtQNJDkYuB54Z1X9NEntbegcfXOOnZqaeqzd6/Xo9XojVilJ+5fp6emuNTXytlK1t+P2ACsnK4DP\nAl+oqo90ffcCvarameQI4NaqWpfk8q59bTfuPuDle2YPfdusUWqSpFYk4fHz6VBVc51wz2vUS0P/\nAnxrTwh0bgTO7drnAjf09b8RIMl6YNfsEJAkjd/QM4IkLwO+BNzNTCQV8FfAV4HrgCOBbcCZVbWr\nW+cy4GTgZ8B5VbV5ju06I5CkASzWjGCkS0P7gkEgSYNZKpeGJEnLnEEgSY0zCCSpcQaBJDXOIJCk\nxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqc\nQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkE\nktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXFjD4IkJye5L8kDSd477veXJD3RWIMgyQHAZcBJwK8D\nZyd5/jhr2Nemp6cnXcJIrH+ylnP9y7l2WP71j2LcM4LfBrZW1fer6ufANcBpY65hn1ruXybrn6zl\nXP9yrh2Wf/2jGHcQrAa29y3v6PokSRMy7iDIHH015hr2iUcffZQkXHjhhSThjDP+dNIlSdJAUjW+\n43CS9cBUVZ3cLW8Aqqou6huzXwSDJI1bVc11sj2vcQfBgcD9wInAD4CvAmdX1b1jK0KS9AQrxvlm\nVfXLJG8DbmbmstSVhoAkTdZYZwSSpKVnYk8Wz/dgWZILktyV5BtJvrTUnjcY9MG4JK9NsjvJCeOs\nbz4D7P9zkvwoyebu702TqHMug+z7JH+S5J4kdyf5t3HX+GQG2Pd/333vNye5P8nDk6hzbwao/8gk\nm7r6tyQ5ZRJ17s0A9R+V5JYkd3af4zmTqHMuSa5MsjPJXU8y5tIkW7t9f/xAG66qsf8xE0DfBp4L\nPAXYAjx/1piD+9qvAr4wiVqHrX/PZwC+CHwZOGHSdS9w/58DXDrpWoes/Vjg68Ah3fJhk657od+d\nvvFvA66YdN0L3P//BFzQtdcB35103Qus/zrg9V27B1w96br7avsd4Hjgrr28fgrwua79EuD2QbY7\nqRnBvA+WVdVP+xYPBnaPsb75DPpg3PuBi4BHx1ncAAatf6hfIOxjg9T+ZuCjVfU/AFX14zHX+GQW\n+lDl2cCnxlLZYAapfzdwSNd+JvDgGOubzyD1vwDYBFBV03O8PjFVdRvwyJMMOQ24uht7B7Ayyar5\ntjupIBjowbIkb0nybeCDwDvGVNsg5q2/m5KtqarPj7OwAQ36YN8Z3fTyuiRrxlPavAap/XnA2iS3\nJflykpPGVt38Bn6oMslRwNF0B6UlYpD6LwTekGQ78Fng7WOqbRCD1L8FeA1AkjOAg5McOp7yRjb7\n8z3IAA/tTioIBnqwrKr+oaqOBd4L/M0+r2pwT1p/kgAfBt49zzqTMsj+vxE4uqqOB/4LuGqfVzWY\nQWpfwczlod8DXgdckeSQ/7fWZCzkocqzgOurm+cvEYPUfzbw8ao6EnglsJTu0QxS/3uAXpKvA7/L\nzMH0F/u6sEUy1EO7kwqCHcBRfctrgIeeZPy1wOn7tKKFma/+ZzDzj+pNJ/kusB64YQndMJ53/1fV\nI93UGeCfgRePqbb5DPLd2QHcUFW7q+p7zDy7ctx4ypvXQr77Z7G0LgvBYPWfz8x1dqrqduCpSQ4b\nT3nzGuS7/4Oqek1VvRj4667vf8dX4kh2AEf2Lc93bJ0xoRseB/L4DZuDmJmKrZs15ti+9quAr076\nRs1C6p81/lbgRZOue4H7/4i+9h8DX5503Quo/STgE137MOD7wKGTrn0h3x1gLfCdSdc75P7/HHBO\n114H7Jh03Qus/9k8/tP6v2PmX0OYeO199R0N3L2X107l8ZvF6xnwZvEkP8zJzJypbQU2dH0XAn/U\ntS8BvglsZubSxF4PtEux/lljN7GEfjU04P7/QLf/v9Ht/+dNuuaF7HvgYuAe4E7gzEnXPET9G4EP\nTLrWIb8764DbuoPsZuDESde8wPpfAzwA3Ad8DHjKpGvuq/2TzJzhPwpsA84DLgD+vG/MZV3Y3Tno\ncccHyiSpcf6vKiWpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmN+z/g3E7AV/Yc+QAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe6d4923610>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(5)\n",
    "plt.hist(p[pos_GT1],100)\n",
    "\n",
    "alp = 0.05\n",
    "# It's a two tailed distribution, we select the voxels of the lower and upper tails\n",
    "print 'Relevant voxels selected:'\n",
    "print np.where(np.any([p[pos_GT1]<0.5*alp, p[pos_GT1]>(1-0.5*alp)], axis=0))[0].shape[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}