{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed and Component models\n",
    "\n",
    "Example of a fit of fixed and component pcm-models to data from M1. The models and data is taken from [Ejaz et al. (2015)](http://www.diedrichsenlab.org/pubs/Ejaz_NN_2015.pdf). *Hand usage predicts the structure of representations in sensorimotor cortex*, Nature Neuroscience.  \n",
    "\n",
    "We will fit the following 5 models \n",
    "\n",
    "* null: G=np.eye, all finger patterns are equally far away from each other, Note that in many situations the no-information null model, G = np.zeros, maybe more appropriate   \n",
    "* Muscle: Fixed model with G = covariance of muscle activities \n",
    "* Natural: Fixed model with G = covariance of natural movements\n",
    "* Muscle+Natural: Component model of muscle and natural covariance \n",
    "* Noiseceil: Noise ceiling model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary libraries\n",
    "import PcmPy as pcm\n",
    "import numpy as np\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the activity Data (`Data`), condition vector (`cond_vec`), partition vector (`part_vec`), and model matrices for Muscle and Natural stats Models (`M`): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('data_demo_finger7T.p','rb')\n",
    "Data,cond_vec,part_vec,modelM = pickle.load(f)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are build a list of datasets (one per subject) from the Data and condition vectors "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = list()\n",
    "for i in range(len(Data)):\n",
    "    obs_des = {'cond_vec': cond_vec[i],\n",
    "               'part_vec': part_vec[i]}\n",
    "    Y.append(pcm.dataset.Dataset(Data[i],obs_descriptors = obs_des))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect the data\n",
    "Before fitting the models, it is very useful to first visualize the different data sets to see if there are outliers. One powerful way is to estimate a cross-validated estimate of the second moment matrix. This matrix is just another form of representing a cross-validated representational dissimilarity matrix (RDM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate and plot the second moment matrices across all data sets\n",
    "N=len(Y)\n",
    "G_hat = np.zeros((N,5,5))\n",
    "for i in range(N):\n",
    "    G_hat[i,:,:],_ = pcm.est_G_crossval(Y[i].measurements,\n",
    "                            Y[i].obs_descriptors['cond_vec'],\n",
    "                            Y[i].obs_descriptors['part_vec'],\n",
    "                            X=pcm.matrix.indicator(Y[i].obs_descriptors['part_vec']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAFhCAYAAAAsiOM3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAH41JREFUeJzt3X9s1PWe7/HXzLQzQ+l0sBAKPS3C3V1FxZ8IpLBBExtJ5A/5Sz0xhiUuGtMaGvYPJTmxYc1NzYnXmKtEifcAm6um6B+EBF0MqQeInrJw2uOGHx7O0btxh5S2cpZMabEzMPPZP7gdqHR6OuX9nV99PpL5o8On7/l8p69+fTkznfE555wAAAAM+Au9AQAAUD4oFgAAwAzFAgAAmKFYAAAAMxQLAABghmIBAADMUCwAAIAZigUAADBTkc8bS6fT6uvrUyQSkc/ny+dNYxqcc7p06ZLq6+vl99t1UHJQWrzIARkoLZwLkEsG8los+vr61NjYmM+bhIFYLKaGhgazeeSgNFnmgAyUJs4FmEoG8losIpGIJOmT3y1WVbXtszBbPv5H03mSVPf7pPlMSXIB+2aeCto/q5W6Mqqeg/8z83OzMjZv8Xv/JP+skOnsym+qTedJUvgv3rzrfSrowcyQB9lKjupP/+efTXMwNuuO32xRoMo2A74jc0znSVL4Ytp8piSlKu1/XiML7WemE6P6f//bNgPS9Ryc+f0CRYz/m7C2+x9M50lS7b9Wmc+UpMgnJ+yHrrjHdNzVVEJf9f6vKWUgr8Vi7KGuqmq/ZkdsQxQIh03nSVJFhTcvQXEV9r/4vkrvXi5j/RDl2Dz/rJD8VbY/t0DIPgeBoEcfp+NBsZAHxWKMZQ7GZgWqQvbFwpMMeFMs5EGxCJRIBm6cF6n2q8b4vwnW5xZJClTaz5SkCl+lB0O92etUMsCLNwEAgBmKBQAAMEOxAAAAZqZVLHbs2KHFixcrHA5r1apVOn78uPW+UOTIACRyADKAm+VcLPbu3autW7eqvb1dvb29uv/++7Vu3ToNDg56sT8UITIAiRyADGBiOReLt956S5s3b9amTZt099136/3331dVVZV27drlxf5QhMgAJHIAMoCJ5VQsksmkenp61NzcfH2A36/m5mZ1d3ebbw7FhwxAIgcgA8gup/exuHDhglKplOrq6sZdX1dXpz/+8Y83rU8kEkokEpmvh4aGprlNFItcMyCRg3LEuQCcC5CNp38V0tHRoWg0mrnw1q0zEzkAGYBEDmaKnIrFvHnzFAgENDAwMO76gYEBLViw4Kb127ZtUzwez1xisdit7RYFl2sGJHJQjjgXgHMBssmpWASDQS1fvlxdXV2Z69LptLq6utTU1HTT+lAopJqamnEXlLZcMyCRg3LEuQCcC5BNzp8VsnXrVm3cuFEPP/ywVq5cqbffflsjIyPatGmTF/tDESIDkMgByAAmlnOxePrpp/Xjjz/qtddeU39/vx544AEdPHjwphfwoHyRAUjkAGQAE5vWp5u2traqtbXVei8oIWQAEjkAGcDN+KwQAABghmIBAADMUCwAAIAZigUAADBDsQAAAGam9Vcht2rLx/+oQDhsOjN997DpPEk6H5htPlOSon+2nxm6lDKf6XPOfOaNKr+pViBkm4PL9WnTeZLku+pN//bZb1XOg616MXOM78gc+YwzEL/7quk8SRqJB8xnSlL1f/rsZ56z/71NJb09F6zt/gf5q2xz8OGq35jOk6TNx7aYz5Sk2oZfmM90/9FvOy+dnPJaHrEAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMVhbjRut8nVVFh22nOB2abzpOkRMMV85mSFFel+czwBfsfZSrhbTzCf3EKBJ3pTN9V+6480pA2nylJlcO+kpjpUuYjM8IX0woEbe/fkXjAdJ4kuQrbnI5JzLH/ecnZz0xVerDPG9T+a5UClWHTmZuPbTGdJ0lDS6+az5SkP728yHzm7HO2P7NUYlTaObW1PGIBAADMUCwAAIAZigUAADBDsQAAAGYoFgAAwExOxaKjo0MrVqxQJBLR/PnztWHDBp09e9arvaEIkQFI5ABkANnlVCyOHDmilpYWHTt2TIcOHdKVK1f0+OOPa2RkxKv9ociQAUjkAGQA2eX0RgUHDx4c9/WePXs0f/589fT0aO3ataYbQ3EiA5DIAcgAsrul11jE43FJUm1trclmUHrIACRyADKA66b91orpdFptbW1as2aNli1bNuGaRCKhRCKR+XpoaGi6N4ciNJUMSOSg3HEuAOcC3Gjaj1i0tLTo1KlT6uzszLqmo6ND0Wg0c2lsbJzuzaEITSUDEjkod5wLwLkAN5pWsWhtbdWBAwf029/+Vg0NDVnXbdu2TfF4PHOJxWLT3iiKy1QzIJGDcsa5AJwL8HM5PRXinNPLL7+sffv26fDhw1qyZMmk60OhkEKh0C1tEMUl1wxI5KAccS4A5wJkk1OxaGlp0ccff6z9+/crEomov79fkhSNRjVr1ixPNojiQgYgkQOQAWSX01Mh7733nuLxuB599FEtXLgwc9m7d69X+0ORIQOQyAHIALLL+akQzGxkABI5ABlAdnxWCAAAMEOxAAAAZigWAADADMUCAACYoVgAAAAz0/6skFvhAj65Cp/pzOifTcdJkuKqtB8qKVF31Xxmsta+I6Z/SpnPvFEqKCloO9OXtp0nSZXDtlkdk/bgt8+TmR7GIFXpkypt79/q/7T/eSXmeJOB0Tr7O9dVBMxnpkbNR44T+eSEKny259vahl+YzpOkP728yHymJD297ivzmUcH/tZ03tWRhLRzamt5xAIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFQU4kZTQb98lbadJnQpZTpPksIXvLl7krX2fe622y+az0xdTihmPvWG+SGfFPKZznQeVOXKYds9jkl7EK9kjTOfmR61nzlmZKFPAeMMVJ/zYL/Omwy4ioD5zNG6q+Yz0z/ZzxxnxT1SRdh0pPuPftN5kjT7nDc5ODrwt+Yzl9WeN52XDCb1b1NcyyMWAADADMUCAACYoVgAAAAzFAsAAGCGYgEAAMxQLAAAgJlbKhZvvPGGfD6f2trajLaDUkMGQAYgkQNcN+1iceLECe3cuVP33Xef5X5QQsgAyAAkcoDxplUshoeH9eyzz+qDDz7QbbfdZr0nlAAyADIAiRzgZtMqFi0tLVq/fr2am5snXZdIJDQ0NDTugvIw1QxI5KBckQFI5AA3y/lNhTs7O9Xb26sTJ0781bUdHR3avn37tDaG4pVLBiRyUI7IACRygInl9IhFLBbTli1b9NFHHykc/uvv675t2zbF4/HMJRbz8pMnkA+5ZkAiB+WGDEAiB8gup0csenp6NDg4qIceeihzXSqV0tGjR/Xuu+8qkUgoELj+oTqhUEihUMhutyi4XDMgkYNyQwYgkQNkl1OxeOyxx3Ty5Mlx123atElLly7VK6+8clOIUH7IAMgAJHKA7HIqFpFIRMuWLRt33ezZszV37tybrkd5IgMgA5DIAbLjnTcBAICZnP8q5OcOHz5ssA2UMjIAMgCJHOAaHrEAAABmKBYAAMAMxQIAAJi55ddY5MI5J0lKXRk1n+37/7MtpRLe3D3pn1LmM1OXE57NdMb3bSYHSfscOA+qsrP/cUmS0h7MTY/a/x6kR6/9nCxzMDYrnbDPQCrpwbmg0mc+U5JS9oev9E9X7Wd6kIEb511N2Z+/XDppPjPlQV4l6eqI/fEng7bHnxy5ImlqGfA566RM4ty5c2psbMzXzcFILBZTQ0OD2TxyUJosc0AGShPnAkwlA3ktFul0Wn19fYpEIvL5Jv8/gKGhITU2NioWi6mmpiZPO/ReKR2Xc06XLl1SfX29/H67hwLIQWkdlxc5IAOldVyFPheU0n2Vi1I6rlwykNenQvx+f85tt6ampujv8OkoleOKRqPmM8nBdaVyXNY5IAPXlcpxFcO5oFTuq1yVynFNNQO8eBMAAJihWAAAADNFWyxCoZDa29vL7pPwyvW4vFKu91e5HpcXyvW+Ktfj8kK53lflelx5ffEmAAAob0X7iAUAACg9FAsAAGCGYgEAAMxQLAAAgJmCFosdO3Zo8eLFCofDWrVqlY4fPz7p+k8//VRLly5VOBzWvffeq88//zxPO52ajo4OrVixQpFIRPPnz9eGDRt09uzZSb9nz5498vl84y7hcDhPOy4O5IAckAEyQAbKKAOuQDo7O10wGHS7du1yp0+fdps3b3Zz5sxxAwMDE67/+uuvXSAQcL/+9a/dmTNn3K9+9StXWVnpTp48meedZ7du3Tq3e/dud+rUKffNN9+4J554wi1atMgNDw9n/Z7du3e7mpoad/78+cylv78/j7suLHJwzUzOARm4hgyQgXLJQMGKxcqVK11LS0vm61Qq5err611HR8eE65966im3fv36cdetWrXKvfjii57u81YMDg46Se7IkSNZ1+zevdtFo9H8barIkINrZnIOyMA1ZIAMlEsGCvJUSDKZVE9Pj5qbmzPX+f1+NTc3q7u7e8Lv6e7uHrdektatW5d1fTGIx+OSpNra2knXDQ8P6/bbb1djY6OefPJJnT59Oh/bKzhyMN5MzAEZGI8MXEMGSjsDBSkWFy5cUCqVUl1d3bjr6+rq1N/fP+H39Pf357S+0NLptNra2rRmzRotW7Ys67o777xTu3bt0v79+/Xhhx8qnU5r9erVOnfuXB53Wxjk4LqZmgMycB0ZIAPlkoG8frrpTNLS0qJTp07pq6++mnRdU1OTmpqaMl+vXr1ad911l3bu3KnXX3/d623CY+QAZAAzLQMFKRbz5s1TIBDQwMDAuOsHBga0YMGCCb9nwYIFOa0vpNbWVh04cEBHjx7N+aOhKysr9eCDD+q7777zaHfFgxxkN1NyQAayIwNkoFQzUJCnQoLBoJYvX66urq7Mdel0Wl1dXePa2o2amprGrZekQ4cOZV1fCM45tba2at++ffryyy+1ZMmSnGekUimdPHlSCxcu9GCHxYUcZDdTckAGsiMDZKBkM1CoV412dna6UCjk9uzZ486cOeNeeOEFN2fOnMyf1jz33HPu1Vdfzaz/+uuvXUVFhXvzzTfdt99+69rb24vuz4teeuklF41G3eHDh8f9udDly5cza35+XNu3b3dffPGF+/77711PT4975plnXDgcdqdPny7EIeQdObhmJueADFxDBshAuWSgYMXCOefeeecdt2jRIhcMBt3KlSvdsWPHMv/2yCOPuI0bN45b/8knn7g77rjDBYNBd88997jPPvsszzuenKQJL7t3786s+flxtbW1Ze6Duro698QTT7je3t78b76AyAE5IANkgAyUTwb42HQAAGCGzwoBAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFTk88bS6bT6+voUiUTk8/nyedOYBuecLl26pPr6evn9dh2UHJQWL3JABkqLV+cClKe8Fou+vj41Njbm8yZhIBaLqaGhwWweOShNljkgA6XJ+lyA8pTXYhGJRCRJf7/3eVVUBU1nn++yP0lFv0+Zz5SkZMS+8Y/Otf+/vlRyVH/e+c+Zn5uVsXlH/22eqqtt74un/qXNdJ4kNXRdMp8pSf4rafOZ59dEzWemkqP6029sczA264fexaoxzsDaN543nSdJ8/592HymJF2dXWk+8+LfhcxnppKj+vb/vm5+LkB5ymuxGHvIs6IqqIrZtuEPhMKm8ySpotKbYpEK2heLQMi7h5OtH6oem1dd7Ve1ccnyJAcVV8xnSpI/bZ8vL45/jGUOxmbVVPtVY52BoAcZCFw1n3ltsH2xCATti8UYnrbCVPBkGQAAMEOxAAAAZigWAADADMUCAACYmVax2LFjhxYvXqxwOKxVq1bp+PHj1vtCkSMDkMgBgJvlXCz27t2rrVu3qr29Xb29vbr//vu1bt06DQ4OerE/FCEyAIkcAJhYzsXirbfe0ubNm7Vp0ybdfffdev/991VVVaVdu3Z5sT8UITIAiRwAmFhOxSKZTKqnp0fNzc3XB/j9am5uVnd3903rE4mEhoaGxl1Q2nLNgEQOyhHnAgDZ5FQsLly4oFQqpbq6unHX19XVqb+//6b1HR0dikajmQtv4Vv6cs2ARA7KEecCANl4+lch27ZtUzwez1xisZiXN4ciRQ5ABoCZI6e39J43b54CgYAGBgbGXT8wMKAFCxbctD4UCikU8u7tZZF/uWZAIgfliHMBgGxyesQiGAxq+fLl6urqylyXTqfV1dWlpqYm882h+JABSOQAQHY5fwjZ1q1btXHjRj388MNauXKl3n77bY2MjGjTpk1e7A9FiAxAIgcAJpZzsXj66af1448/6rXXXlN/f78eeOABHTx48KYXcaF8kQFI5ADAxKb1semtra1qbW213gtKCBmARA4A3IzPCgEAAGYoFgAAwAzFAgAAmKFYAAAAM9N68eatOt/VqEAobDqz8u//YjpPki7655rPlKTqc2nzmVUDznxmKmm/zxs99S9t5jnwPxQ3nSdJP1RGzWdK0tzT9vdvdV/KfObVK/Yzx6x943kFgrYZuLjM/n4d+ptq85mSFPovn/lML84vqaT9+QXli0csAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMBMRSFuNPp9ShWVKdOZF/1zTedJ0uVf2O5xTDJq3+dCF33mM1MJb3tnQ9clVVRcMZ35Q2XUdJ4kJW9PmM+UpPO19r9+c07b/8xSSe9yMO/fh1URuGo6c+hvqk3nSdLVebY5zUhXmo+sHC6tDKD8kBYAAGCGYgEAAMxQLAAAgBmKBQAAMEOxAAAAZigWAADATE7FoqOjQytWrFAkEtH8+fO1YcMGnT171qu9oQiRAUjkAEB2ORWLI0eOqKWlRceOHdOhQ4d05coVPf744xoZGfFqfygyZAASOQCQXU7v0HPw4MFxX+/Zs0fz589XT0+P1q5da7oxFCcyAIkcAMjult76Lx6PS5Jqa2sn/PdEIqFE4vq7Fg4NDd3KzaEI/bUMSORgJuBcAGDMtF+8mU6n1dbWpjVr1mjZsmUTruno6FA0Gs1cGhsbp71RFJ+pZEAiB+WOcwGAG027WLS0tOjUqVPq7OzMumbbtm2Kx+OZSywWm+7NoQhNJQMSOSh3nAsA3GhaT4W0trbqwIEDOnr0qBoaGrKuC4VCCoVC094citdUMyCRg3LGuQDAz+VULJxzevnll7Vv3z4dPnxYS5Ys8WpfKFJkABI5AJBdTsWipaVFH3/8sfbv369IJKL+/n5JUjQa1axZszzZIIoLGYBEDgBkl9NrLN577z3F43E9+uijWrhwYeayd+9er/aHIkMGIJEDANnl/FQIZjYyAIkcAMiOzwoBAABmKBYAAMAMxQIAAJihWAAAADO39Fkh05WM+JUK2naa6nNp03mSlIx607uuVnnwwjcPRqZHvX2Bnv9KWv50ynTm3NP2OThf682vyaz6YfOZly7XmM9Mj5qPzLg6u1KqqDSdGfovn+k8SVLado9jUlX2eb1cZ3/e8vpcgPLCIxYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYKaiEDc6OtenQMhnOrNqwJnOk6TQRds9ZthvValZ9kPTPg82eoPza6IKhMKmM6v7UqbzJGnOaW/696XLNeYz0/Wj9jMv288cc/HvQgoEQ6Yzq8+lTedJUuWwNxm4XGc/N/U/SisDKD88YgEAAMxQLAAAgBmKBQAAMEOxAAAAZigWAADADMUCAACYuaVi8cYbb8jn86mtrc1oOyg1ZABkAMCNpl0sTpw4oZ07d+q+++6z3A9KCBkAGQDwc9MqFsPDw3r22Wf1wQcf6LbbbrPeE0oAGQAZADCRaRWLlpYWrV+/Xs3NzZOuSyQSGhoaGndBeZhqBiRyUK7IAICJ5PyW3p2dnert7dWJEyf+6tqOjg5t3759WhtD8colAxI5KEdkAEA2OT1iEYvFtGXLFn300UcKh//6Zzxs27ZN8Xg8c4nFYtPeKIpDrhmQyEG5IQMAJpPTIxY9PT0aHBzUQw89lLkulUrp6NGjevfdd5VIJBQIBDL/FgqFFArZfsAQCivXDEjkoNyQAQCTyalYPPbYYzp58uS46zZt2qSlS5fqlVdeuelkgvJDBkAGAEwmp2IRiUS0bNmycdfNnj1bc+fOvel6lCcyADIAYDK88yYAADCT81+F/Nzhw4cNtoFSRgZABgCM4RELAABghmIBAADMUCwAAICZW36NRS6cc5KkVHLUfHYqmbafmfCmd6VHnf1MnwczR6/9nMZ+bla8zMHVKynzmamkVznwYOZl+6HpnxKSbHPg7bnA/nfBuwx48HtbIhlA+fK5PCbl3LlzamxszNfNwUgsFlNDQ4PZPHJQmixzQAZKk/W5AOUpr8UinU6rr69PkUhEPp9v0rVDQ0NqbGxULBZTTU1NnnbovVI6LuecLl26pPr6evn9dv/HRg5K67i8yAEZKK3j8upcgPKU16dC/H5/zm23pqam6H/ppqNUjisajZrPJAfXlcpxWeeADFxXKsflxbkA5YnqCQAAzFAsAACAmaItFqFQSO3t7WX3iYjlelxeKdf7q1yPywvlel+V63EBeX3xJgAAKG9F+4gFAAAoPRQLAABghmIBAADMUCwAAICZghaLHTt2aPHixQqHw1q1apWOHz8+6fpPP/1US5cuVTgc1r333qvPP/88Tzudmo6ODq1YsUKRSETz58/Xhg0bdPbs2Um/Z8+ePfL5fOMu4XA4TzsuDuSAHJABMoDyUbBisXfvXm3dulXt7e3q7e3V/fffr3Xr1mlwcHDC9b/73e/0y1/+Us8//7z+8Ic/aMOGDdqwYYNOnTqV551nd+TIEbW0tOjYsWM6dOiQrly5oscff1wjIyOTfl9NTY3Onz+fufzwww952nHhkYPrZmoOyMB1MzUDKDOuQFauXOlaWloyX6dSKVdfX+86OjomXP/UU0+59evXj7tu1apV7sUXX/R0n7dicHDQSXJHjhzJumb37t0uGo3mb1NFhhxcM5NzQAaumckZQHkpyCMWyWRSPT09am5uzlzn9/vV3Nys7u7uCb+nu7t73HpJWrduXdb1xSAej0uSamtrJ103PDys22+/XY2NjXryySd1+vTpfGyv4MjBeDMxB2RgvJmYAZSfghSLCxcuKJVKqa6ubtz1dXV16u/vn/B7+vv7c1pfaOl0Wm1tbVqzZo2WLVuWdd2dd96pXbt2af/+/frwww+VTqe1evVqnTt3Lo+7LQxycN1MzQEZuG6mZgDlJ6+fbjqTtLS06NSpU/rqq68mXdfU1KSmpqbM16tXr9Zdd92lnTt36vXXX/d6m/AYOQAZwExTkGIxb948BQIBDQwMjLt+YGBACxYsmPB7FixYkNP6QmptbdWBAwd09OjRnD8aurKyUg8++KC+++47j3ZXPMhBdjMlB2Qgu5mSAZSfgjwVEgwGtXz5cnV1dWWuS6fT6urqGtfYb9TU1DRuvSQdOnQo6/pCcM6ptbVV+/bt05dffqklS5bkPCOVSunkyZNauHChBzssLuQgu5mSAzKQ3UzJAMpQoV412tnZ6UKhkNuzZ487c+aMe+GFF9ycOXNcf3+/c8655557zr366quZ9V9//bWrqKhwb775pvv2229de3u7q6ysdCdPnizUIdzkpZdectFo1B0+fNidP38+c7l8+XJmzc+Pa/v27e6LL75w33//vevp6XHPPPOMC4fD7vTp04U4hLwjB9fM5ByQgWtmcgZQXgpWLJxz7p133nGLFi1ywWDQrVy50h07dizzb4888ojbuHHjuPWffPKJu+OOO1wwGHT33HOP++yzz/K848lJmvCye/fuzJqfH1dbW1vmPqirq3NPPPGE6+3tzf/mC4gckAMyQAZQPvjYdAAAYIbPCgEAAGYoFgAAwAzFAgAAmKFYAAAAMxQLAABghmIBAADMUCwAAIAZigUAADBDsQAAAGYoFgAAwAzFAgAAmKFYAAAAM/8NCp4ys6K4ma8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 7 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# show all second moment matrices\n",
    "vmin = G_hat.min()\n",
    "vmax = G_hat.max()\n",
    "for i in range(N):\n",
    "    plt.subplot(2,4,i+1)\n",
    "    plt.imshow(G_hat[i,:,:],vmin=vmin,vmax=vmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice - up to a scaling factor (Subject 4 has especially high signal-to-noise) all seven subjects have a very similar structure of the reresentation of fingers in M1.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are more used to looking in representational dissimilarity matrices (RDMs), you can also transform the second momement matrix into this (see Diedrichsen & Kriegeskorte,2017)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x336ab6f10>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGdCAYAAAAv9mXmAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAEfJJREFUeJzt3W9o3YW9x/Fv2i6napPM1rUut8ksONy6LnW2VoJjczZTel1R7pM9EBY6GDjS0dInI08sezDSR0OZpZb988G1tGxQBbnalW5tkGtnmpJLdeid4CCjazMZJG3E05qc+2Asd53a5cR88zu/9vWC8+AcfvH34afkzTm/JDbVarVaAMA8W1T0AACuTQIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKZYs9Amnp6fj7Nmz0dLSEk1NTQt9egA+hlqtFhcuXIj29vZYtOjq71EWPDBnz56Njo6OhT4tAPNodHQ0Vq9efdVjFjwwLS0tERHx5fj3WBKfWOjTl8rZZz9f9IRSeLrrP4ueUApfqjQXPaEUHh/7YtETGtqlycux/8H/mvlefjULHpi/fyy2JD4RS5oE5moW31gpekIpLGtxK3E2Wiuu02xU3vV9aTZmc4vDf3EApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIo5BWbv3r1x2223xdKlS+Oee+6JV199db53AVBydQfm0KFDsWvXrti9e3ecPn061q9fHw8++GCMjY1l7AOgpOoOzI9+9KP4zne+E9u2bYu1a9fG008/HTfeeGP8/Oc/z9gHQEnVFZhLly7F8PBw9PT0/P8/YNGi6OnpiVdeeWXexwFQXkvqOfidd96JqampWLVq1RWvr1q1Kt54440P/ZpqtRrVanXm+cTExBxmAlA26T9FNjAwEG1tbTOPjo6O7FMC0ADqCswtt9wSixcvjvPnz1/x+vnz5+PWW2/90K/p7++P8fHxmcfo6Ojc1wJQGnUFprm5OTZs2BDHjh2beW16ejqOHTsW3d3dH/o1lUolWltbr3gAcO2r6x5MRMSuXbuit7c3Nm7cGJs2bYonnngiJicnY9u2bRn7ACipugPzzW9+M/7yl7/E448/HufOnYs777wzXnrppQ/c+Afg+lZ3YCIitm/fHtu3b5/vLQBcQ/wtMgBSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkGJJUSc+++znY/GNlaJOXwr/9h+vFz2hFB7btqPoCaXw7jcmip5QCmtW/LXoCQ3t8uSlWR/rHQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUtQdmMHBwdi6dWu0t7dHU1NTPPfccwmzACi7ugMzOTkZ69evj71792bsAeAasaTeL9iyZUts2bIlYwsA1xD3YABIUfc7mHpVq9WoVqszzycmJrJPCUADSH8HMzAwEG1tbTOPjo6O7FMC0ADSA9Pf3x/j4+Mzj9HR0exTAtAA0j8iq1QqUalUsk8DQIOpOzAXL16Mt956a+b522+/HSMjI7F8+fLo7Oyc13EAlFfdgTl16lR87Wtfm3m+a9euiIjo7e2NZ555Zt6GAVBudQfmvvvui1qtlrEFgGuI34MBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAAplhR14qe7/jOWtejb1Ty2bUfRE0ph+S9eKXpCKXzyD3cWPaEU/vjlNUVPaGhT1fdmfazv8ACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIUVdgBgYG4u67746WlpZYuXJlPPLII/Hmm29mbQOgxOoKzIkTJ6Kvry9OnjwZR48ejcuXL8cDDzwQk5OTWfsAKKkl9Rz80ksvXfH8mWeeiZUrV8bw8HB85StfmddhAJRbXYH5Z+Pj4xERsXz58o88plqtRrVanXk+MTHxcU4JQEnM+Sb/9PR07Ny5M+69995Yt27dRx43MDAQbW1tM4+Ojo65nhKAEplzYPr6+uK1116LgwcPXvW4/v7+GB8fn3mMjo7O9ZQAlMicPiLbvn17vPDCCzE4OBirV6++6rGVSiUqlcqcxgFQXnUFplarxfe+9704fPhwHD9+PNasWZO1C4CSqyswfX19ceDAgXj++eejpaUlzp07FxERbW1tccMNN6QMBKCc6roHs2/fvhgfH4/77rsvPv3pT888Dh06lLUPgJKq+yMyAJgNf4sMgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkWFLUib9UaY7Wir5dzbvfmCh6Qil88g93Fj2hFBa9PFL0hFJon15f9ISG9v7778X/zvJY3+EBSCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkKKuwOzbty+6urqitbU1Wltbo7u7O1588cWsbQCUWF2BWb16dezZsyeGh4fj1KlTcf/998fDDz8cr7/+etY+AEpqST0Hb9269YrnP/zhD2Pfvn1x8uTJ+MIXvjCvwwAot7oC84+mpqbil7/8ZUxOTkZ3d/dHHletVqNarc48n5iYmOspASiRum/ynzlzJpYtWxaVSiUee+yxOHz4cKxdu/Yjjx8YGIi2traZR0dHx8caDEA51B2YO+64I0ZGRuJ3v/tdfPe7343e3t74/e9//5HH9/f3x/j4+MxjdHT0Yw0GoBzq/oisubk5br/99oiI2LBhQwwNDcWTTz4Z+/fv/9DjK5VKVCqVj7cSgNL52L8HMz09fcU9FgCIqPMdTH9/f2zZsiU6OzvjwoULceDAgTh+/HgcOXIkax8AJVVXYMbGxuJb3/pW/PnPf462trbo6uqKI0eOxNe//vWsfQCUVF2B+dnPfpa1A4BrjL9FBkAKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUiwp6sSPj30xKu9+oqjTl8KaFX8tekIp/PHLa4qeUArt0+uLnlAKTf/9P0VPaGhNtcuzPtY7GABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACk+FiB2bNnTzQ1NcXOnTvnaQ4A14o5B2ZoaCj2798fXV1d87kHgGvEnAJz8eLFePTRR+MnP/lJ3HzzzfO9CYBrwJwC09fXFw899FD09PT8y2Or1WpMTExc8QDg2rek3i84ePBgnD59OoaGhmZ1/MDAQPzgBz+oexgA5VbXO5jR0dHYsWNHPPvss7F06dJZfU1/f3+Mj4/PPEZHR+c0FIByqesdzPDwcIyNjcVdd90189rU1FQMDg7GU089FdVqNRYvXnzF11QqlahUKvOzFoDSqCswmzdvjjNnzlzx2rZt2+Jzn/tcfP/73/9AXAC4ftUVmJaWlli3bt0Vr910002xYsWKD7wOwPXNb/IDkKLunyL7Z8ePH5+HGQBca7yDASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSLFnoE9ZqtYiIuDR5eaFPXTqXJy8VPaEUpqrvFT2hFN5/33Wajaaa701X83787fr8/Xv51TTVZnPUPPrTn/4UHR0dC3lKAObZ6OhorF69+qrHLHhgpqen4+zZs9HS0hJNTU0LeeqPNDExER0dHTE6Ohqtra1Fz2lIrtHsuE6z4zrNTiNep1qtFhcuXIj29vZYtOjqd1kW/COyRYsW/cvqFaW1tbVh/iU2Ktdodlyn2XGdZqfRrlNbW9usjnOTH4AUAgNACoGJiEqlErt3745KpVL0lIblGs2O6zQ7rtPslP06LfhNfgCuD97BAJBCYABIITAApBAYAFJc94HZu3dv3HbbbbF06dK455574tVXXy16UsMZHByMrVu3Rnt7ezQ1NcVzzz1X9KSGMzAwEHfffXe0tLTEypUr45FHHok333yz6FkNZ9++fdHV1TXzi4Pd3d3x4osvFj2r4e3Zsyeamppi586dRU+py3UdmEOHDsWuXbti9+7dcfr06Vi/fn08+OCDMTY2VvS0hjI5ORnr16+PvXv3Fj2lYZ04cSL6+vri5MmTcfTo0bh8+XI88MADMTk5WfS0hrJ69erYs2dPDA8Px6lTp+L++++Phx9+OF5//fWipzWsoaGh2L9/f3R1dRU9pX6169imTZtqfX19M8+npqZq7e3ttYGBgQJXNbaIqB0+fLjoGQ1vbGysFhG1EydOFD2l4d188821n/70p0XPaEgXLlyoffazn60dPXq09tWvfrW2Y8eOoifV5bp9B3Pp0qUYHh6Onp6emdcWLVoUPT098corrxS4jGvB+Ph4REQsX7684CWNa2pqKg4ePBiTk5PR3d1d9JyG1NfXFw899NAV36fKZMH/2GWjeOedd2JqaipWrVp1xeurVq2KN954o6BVXAump6dj586dce+998a6deuKntNwzpw5E93d3fHee+/FsmXL4vDhw7F27dqiZzWcgwcPxunTp2NoaKjoKXN23QYGsvT19cVrr70WL7/8ctFTGtIdd9wRIyMjMT4+Hr/61a+it7c3Tpw4ITL/YHR0NHbs2BFHjx6NpUuXFj1nzq7bwNxyyy2xePHiOH/+/BWvnz9/Pm699daCVlF227dvjxdeeCEGBwcb9n9LUbTm5ua4/fbbIyJiw4YNMTQ0FE8++WTs37+/4GWNY3h4OMbGxuKuu+6aeW1qaioGBwfjqaeeimq1GosXLy5w4exct/dgmpubY8OGDXHs2LGZ16anp+PYsWM+D6ZutVottm/fHocPH47f/OY3sWbNmqInlcb09HRUq9WiZzSUzZs3x5kzZ2JkZGTmsXHjxnj00UdjZGSkFHGJuI7fwURE7Nq1K3p7e2Pjxo2xadOmeOKJJ2JycjK2bdtW9LSGcvHixXjrrbdmnr/99tsxMjISy5cvj87OzgKXNY6+vr44cOBAPP/889HS0hLnzp2LiL/9j5luuOGGgtc1jv7+/tiyZUt0dnbGhQsX4sCBA3H8+PE4cuRI0dMaSktLywfu3910002xYsWKct3XK/rH2Ir24x//uNbZ2Vlrbm6ubdq0qXby5MmiJzWc3/72t7WI+MCjt7e36GkN48OuT0TUfvGLXxQ9raF8+9vfrn3mM5+pNTc31z71qU/VNm/eXPv1r39d9KxSKOOPKftz/QCkuG7vwQCQS2AASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUvwfKenZbi3TT7kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "RDM = pcm.util.G_to_dist(np.mean(G_hat,axis=0))\n",
    "plt.imshow(RDM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W,glam = pcm.util.classical_mds(np.mean(G_hat,axis=0))\n",
    "W\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the models\n",
    "Now we are building a list of models, using a list of second moment matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make an empty list\n",
    "M = []\n",
    "# Null model: All fingers are represented independently.\n",
    "# For RSA model that would mean that all distances are equivalent\n",
    "M.append(pcm.FixedModel('null',np.eye(5)))\n",
    "# Muscle model: Structure is given by covariance structure of EMG signals\n",
    "M.append(pcm.FixedModel('muscle',modelM[0]))\n",
    "# Usage model: Structure is given by covariance structure of EMG signals\n",
    "M.append(pcm.FixedModel('usage',modelM[1]))\n",
    "# Component model: Linear combination of the muscle and usage model\n",
    "M.append(pcm.ComponentModel('muscle+usage',[modelM[0],modelM[1]]))\n",
    "# Free noise ceiling model\n",
    "M.append(pcm.FreeModel('ceil',5)) # Noise ceiling model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at two underlying second moment matrices - these are pretty similar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAEjCAYAAABuGEhQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHTZJREFUeJzt3XtwVPXdx/HPbiAbLkkgAsGYgAIKcgkqERujgiQlwyAFqvVSqwErjk5QGWpro1WUWoN/eJsagXqjozKAKOBg5VIIpCqXAMUijnKR+kQUAoJJCLKQ3d/zRx/2cRuQnPBLTs7m/Zo5M56Ts+d8smy+fnJ2N+szxhgBAABY4Hc7AAAAiB0UCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAs0u+HDh2v48OFuxwAANAGKBQAAsIZiAQAArKFYAAAAaygWMeCxxx6Tz+fTjh079Ktf/UrJycnq2rWrHnnkERljVFFRobFjxyopKUndu3fX008/HbntnDlz5PP59O9//zvqmGvWrJHP59OaNWsi23bu3Knrr79e3bt3V0JCgtLT03XzzTerqqoq6rZvvPGGhg4dqvbt26tz58665pprtGLFih/9HoLBoKZNm6Y+ffooEAgoIyNDv/vd7xQMBs/6/gHQMBMmTND5559fb/vJGXPSypUrddVVV6lTp07q2LGj+vbtq4ceeijy9ePHj+vRRx/VkCFDlJycrA4dOujqq69WaWlpvWN/++23uu2225SUlKROnTqpoKBAH3/8sXw+n+bMmRO172effaYbbrhBKSkpSkhIUFZWlt59911r3z/saON2ANhz00036eKLL9aMGTP03nvv6YknnlBKSopmz56tESNG6KmnntKbb76pBx54QJdffrmuueaaBh/7+PHjys/PVzAY1L333qvu3btr7969Wrp0qb777jslJydLkh5//HE99thjuvLKKzV9+nTFx8drw4YNWr16tUaOHHnKY4fDYf3sZz/TBx98oLvuuksXX3yxtm3bpmeffVY7duzQ4sWLbdw9ACzYvn27rrvuOmVmZmr69OkKBALatWuXPvzww8g+1dXVevnll3XLLbdo0qRJqqmp0SuvvKL8/Hxt3LhRl1xyiaT//OyPGTNGGzdu1D333KN+/fppyZIlKigoOOV5c3JydN555+n3v/+9OnTooAULFmjcuHF6++23NX78+Oa6C3AmBp43bdo0I8ncddddkW11dXUmPT3d+Hw+M2PGjMj2w4cPm3bt2pmCggJjjDGvvfaakWT27NkTdczS0lIjyZSWlhpjjPnnP/9pJJm33nrrtDl27txp/H6/GT9+vAmFQlFfC4fDkf8eNmyYGTZsWGT99ddfN36/3/zjH/+Ius2sWbOMJPPhhx825G4AcJYKCgpMz549620/OWOMMebZZ581ksyBAwdOe5y6ujoTDAajth0+fNikpqaaO+64I7Lt7bffNpLMc889F9kWCoXMiBEjjCTz2muvRbbn5uaaQYMGmWPHjkW2hcNhc+WVV5oLL7zQ6beKJsRTITHkzjvvjPx3XFycsrKyZIzRr3/968j2Tp06qW/fvvriiy8cHfvkFYnly5fr6NGjp9xn8eLFCofDevTRR+X3Rz+0fngZ9b+99dZbuvjii9WvXz8dPHgwsowYMUKSTnn5FIA7OnXqJElasmSJwuHwKfeJi4tTfHy8pP9clTh06JDq6uqUlZWlLVu2RPZbtmyZ2rZtq0mTJkW2+f1+FRYWRh3v0KFDWr16tW688UbV1NREZsS3336r/Px87dy5U3v37rX8naKxKBYxpEePHlHrycnJSkhIUJcuXeptP3z4sKNjX3DBBZo6dapefvlldenSRfn5+SopKYl6fcXu3bvl9/vVv39/R8feuXOntm/frq5du0YtF110kSSpsrLS0fEANJ2bbrpJOTk5uvPOO5Wamqqbb75ZCxYsqFcy/vrXvyozM1MJCQk655xz1LVrV7333ntRM+PLL7/Uueeeq/bt20fdtk+fPlHru3btkjFGjzzySL05MW3aNEnMiZaE11jEkLi4uAZtkyRjjKTTX0kIhUL1tj399NOaMGGClixZohUrVui+++5TcXGx1q9fr/T09EbnDofDGjRokJ555plTfj0jI6PRxwbQcA2ZB+3atVNZWZlKS0v13nvvadmyZZo/f75GjBihFStWKC4uTm+88YYmTJigcePG6be//a26deumuLg4FRcXa/fu3Y5znSwtDzzwgPLz80+5z3+XEbiHYtHKde7cWZL03XffRW3/8ssvT7n/oEGDNGjQIP3hD3/QRx99pJycHM2aNUtPPPGEevfurXA4rE8//TTy4qyG6N27tz7++GPl5ub+6FMmAJpW586d680Cqf488Pv9ys3NVW5urp555hk9+eSTevjhh1VaWqq8vDwtXLhQvXr10jvvvBP1M33y6sJJPXv2VGlpqY4ePRp11WLXrl1R+/Xq1UuS1LZtW+Xl5Z3tt4kmxlMhrVzv3r0lSWVlZZFtoVBIf/nLX6L2q66uVl1dXdS2QYMGye/3R94SOm7cOPn9fk2fPr3eZdGTV0hO5cYbb9TevXv10ksv1fva999/r9raWmffFIBG6d27t6qqqvSvf/0rsu2bb77RokWLIuuHDh2qd7uTv0icnAUnr5T+8Od+w4YNWrduXdTt8vPzdeLEiaif/XA4rJKSkqj9unXrpuHDh2v27Nn65ptv6p3/wIEDDf0W0Qy4YtHKDRgwQD/5yU9UVFSkQ4cOKSUlRfPmzatXIlavXq3JkyfrF7/4hS666CLV1dXp9ddfV1xcnK6//npJ/7kU+fDDD+uPf/yjrr76av385z9XIBBQeXm50tLSVFxcfMoMt912mxYsWKC7775bpaWlysnJUSgU0meffaYFCxZo+fLlysrKavL7Amjtbr75Zj344IMaP3687rvvPh09elQzZ87URRddFHnR5fTp01VWVqbRo0erZ8+eqqys1Isvvqj09HRdddVVkqTrrrtO77zzjsaPH6/Ro0drz549mjVrlvr3768jR45Ezjdu3DgNHTpUv/nNb7Rr1y7169dP7777bqS8/PBqR0lJia666ioNGjRIkyZNUq9evbR//36tW7dOX331lT7++ONmvKfwo1x9TwqsOPlWsP9++1dBQYHp0KFDvf2HDRtmBgwYEFnfvXu3ycvLM4FAwKSmppqHHnrIrFy5Murtpl988YW54447TO/evU1CQoJJSUkx1157rfn73/9e7/ivvvqqufTSS00gEDCdO3c2w4YNMytXrow6/w/fbmqMMcePHzdPPfWUGTBgQOR2Q4YMMY8//ripqqo6i3sHgBMrVqwwAwcONPHx8aZv377mjTfeiHq76apVq8zYsWNNWlqaiY+PN2lpaeaWW24xO3bsiBwjHA6bJ5980vTs2dMEAgFz6aWXmqVLl57y7awHDhwwv/zlL01iYqJJTk42EyZMMB9++KGRZObNmxe17+7du83tt99uunfvbtq2bWvOO+88c91115mFCxc2+f2ChvMZ8yPXqAEAaGaLFy/W+PHj9cEHHygnJ8ftOHCIYgEAcM3333+vdu3aRdZDoZBGjhypTZs2ad++fVFfgzfwGgsAgGvuvfdeff/998rOzlYwGNQ777yjjz76SE8++SSlwqO4YgEAcM3cuXP19NNPa9euXTp27Jj69Omje+65R5MnT3Y7GhqJYgEAAKzh71gAAABrKBYAAMCaZn/xZjgc1tdff63ExET+fDPgAmOMampqlJaWVu9TaFsq5gbgvobOjmYvFl9//TUfKgW0ABUVFWf14XHNibkBtBxnmh3NXiwSExMlSV9uOV9JHb3x25Ik3ZA9zO0IjoW+rf83/Vs6f4f2Z96pBfIlJLgdocHqwse19vCbkZ9FL/Ds3Bg91u0IjoX/Z6/bERwz//cZJV7jTwi4HcGROnNCZcFFZ5wdzV4sTl7GTOroV1KidwZEG3+82xEc8/nauh3BMb/Pe/ezJPk8+fjwzlMKnp0bcd76H4ckhT04N4wvfOadWiDPzrszzA7v/IQCAIAWj2IBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrGlUsSkpKdP755yshIUFXXHGFNm7caDsXgBjE7ABin+NiMX/+fE2dOlXTpk3Tli1bNHjwYOXn56uysrIp8gGIEcwOoHVwXCyeeeYZTZo0SRMnTlT//v01a9YstW/fXq+++mpT5AMQI5gdQOvgqFgcP35cmzdvVl5e3v8fwO9XXl6e1q1bd8rbBINBVVdXRy0AWhens4O5AXiXo2Jx8OBBhUIhpaamRm1PTU3Vvn37Tnmb4uJiJScnR5aMjIzGpwXgSU5nB3MD8K4mf1dIUVGRqqqqIktFRUVTnxKAxzE3AO9q42TnLl26KC4uTvv374/avn//fnXv3v2UtwkEAgoEAo1PCMDznM4O5gbgXY6uWMTHx2vIkCFatWpVZFs4HNaqVauUnZ1tPRyA2MDsAFoPR1csJGnq1KkqKChQVlaWhg4dqueee061tbWaOHFiU+QDECOYHUDr4LhY3HTTTTpw4IAeffRR7du3T5dccomWLVtW70VZAPBDzA6gdXBcLCRp8uTJmjx5su0sAGIcswOIfXxWCAAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAa9q4deIbsoepjT/erdM7tuP5Hm5HcKzXixluR3DuSNDtBI1S1zHgdoQGq6s7Jn3rdorGuWH0WLWJ8859vWNaR7cjONblb5e6HcGxlA2VbkdonIOH3U7giM/EScfOvB9XLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgjeNiUVZWpjFjxigtLU0+n0+LFy9uglgAYglzA2g9HBeL2tpaDR48WCUlJU2RB0AMYm4ArUcbpzcYNWqURo0a1RRZAMQo5gbQejguFk4Fg0EFg8HIenV1dVOfEoDHMTcA72ryF28WFxcrOTk5smRkZDT1KQF4HHMD8K4mLxZFRUWqqqqKLBUVFU19SgAex9wAvKvJnwoJBAIKBAJNfRoAMYS5AXgXf8cCAABY4/iKxZEjR7Rr167I+p49e7R161alpKSoR48eVsMBiA3MDaD1cFwsNm3apGuvvTayPnXqVElSQUGB5syZYy0YgNjB3ABaD8fFYvjw4TLGNEUWADGKuQG0HrzGAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANW3cOnHo20Py+dq6dXrHer2Y4XYEx3be4do/b6P1XOidx8QPxVefcDtCqxD+n70Ke2hudPnbpW5HcOxQf5/bERxrW9vV7QiN0nF3vNsRHDGhoHT4zPtxxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANY6KRXFxsS6//HIlJiaqW7duGjdunD7//POmygYgRjA7gNbDUbFYu3atCgsLtX79eq1cuVInTpzQyJEjVVtb21T5AMQAZgfQerRxsvOyZcui1ufMmaNu3bpp8+bNuuaaa6wGAxA7mB1A6+GoWPy3qqoqSVJKSspp9wkGgwoGg5H16urqszklgBhwptnB3AC8q9Ev3gyHw5oyZYpycnI0cODA0+5XXFys5OTkyJKRkdHYUwKIAQ2ZHcwNwLsaXSwKCwv1ySefaN68eT+6X1FRkaqqqiJLRUVFY08JIAY0ZHYwNwDvatRTIZMnT9bSpUtVVlam9PT0H903EAgoEAg0KhyA2NLQ2cHcALzLUbEwxujee+/VokWLtGbNGl1wwQVNlQtADGF2AK2Ho2JRWFiouXPnasmSJUpMTNS+ffskScnJyWrXrl2TBATgfcwOoPVw9BqLmTNnqqqqSsOHD9e5554bWebPn99U+QDEAGYH0Ho4fioEAJxidgCtB58VAgAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAmjZundjfob38vni3Tu/ckaDbCRzrubCt2xEcq8iPcztCo6Sv8rkdocHqToTdjtBoJhiU8Xknf8qGSrcjONa2tqvbERw7cKk3f0euSU9xO4IjoeAx6ZMz7+fNfw0AANAiUSwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDWOisXMmTOVmZmppKQkJSUlKTs7W++//35TZQMQI5gdQOvhqFikp6drxowZ2rx5szZt2qQRI0Zo7Nix2r59e1PlAxADmB1A69HGyc5jxoyJWv/Tn/6kmTNnav369RowYIDVYABiB7MDaD0cFYsfCoVCeuutt1RbW6vs7OzT7hcMBhUMBiPr1dXVjT0lgBjQkNnB3AC8y/GLN7dt26aOHTsqEAjo7rvv1qJFi9S/f//T7l9cXKzk5OTIkpGRcVaBAXiTk9nB3AC8y3Gx6Nu3r7Zu3aoNGzbonnvuUUFBgT799NPT7l9UVKSqqqrIUlFRcVaBAXiTk9nB3AC8y/FTIfHx8erTp48kaciQISovL9fzzz+v2bNnn3L/QCCgQCBwdikBeJ6T2cHcALzrrP+ORTgcjnouFAAagtkBxCZHVyyKioo0atQo9ejRQzU1NZo7d67WrFmj5cuXN1U+ADGA2QG0Ho6KRWVlpW6//XZ98803Sk5OVmZmppYvX66f/vSnTZUPQAxgdgCth6Ni8corrzRVDgAxjNkBtB58VggAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGvauHViX0KCfP54t07vWF3HgNsRHIuvPuF2BMfSV/ncjtAoe4d5p6OHj/ml991O0Tj+hID8Pu/MDR087HYCxzru9tD9+39q0lPcjtAoR7KPuh3BkfDRY9LsM+/nnWkIAABaPIoFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALDmrIrFjBkz5PP5NGXKFEtxAMQ65gYQ2xpdLMrLyzV79mxlZmbazAMghjE3gNjXqGJx5MgR3XrrrXrppZfUuXNn25kAxCDmBtA6NKpYFBYWavTo0crLyzvjvsFgUNXV1VELgNaHuQG0Dm2c3mDevHnasmWLysvLG7R/cXGxHn/8ccfBAMQO5gbQeji6YlFRUaH7779fb775phISEhp0m6KiIlVVVUWWioqKRgUF4E3MDaB1cXTFYvPmzaqsrNRll10W2RYKhVRWVqYXXnhBwWBQcXFxUbcJBAIKBAJ20gLwHOYG0Lo4Kha5ubnatm1b1LaJEyeqX79+evDBB+sNBwBgbgCti6NikZiYqIEDB0Zt69Chg84555x62wFAYm4ArQ1/eRMAAFjj+F0h/23NmjUWYgBoTZgbQOziigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKw5649Nd8oYI0mqCx9v7lOflbq6Y25HaBXqToTdjtAo4WPe6ejhY/95LJ/8WfSCyNwwJ1xO4ozPxLkdwTETCrodwbFQ0JvzOXzUW7nD3//nsXGm2eEzzTxdvvrqK2VkZDTnKQGcQkVFhdLT092O0SDMDaDlONPsaPZiEQ6H9fXXXysxMVE+n8/acaurq5WRkaGKigolJSVZO25TInPz8GJmqelyG2NUU1OjtLQ0+f3euNLSVHND8ubjg8zNg8zRGjo7mv2pEL/f36S/JSUlJXnmAXASmZuHFzNLTZM7OTnZ6vGaWlPPDcmbjw8yNw8y/7+GzA5v/LoCAAA8gWIBAACsiZliEQgENG3aNAUCAbejNBiZm4cXM0veze01Xryfydw8yNw4zf7iTQAAELti5ooFAABwH8UCAABYQ7EAAADWUCwAAIA1FAsAAGBNzBSLkpISnX/++UpISNAVV1yhjRs3uh3ptMrKyjRmzBilpaXJ5/Np8eLFbkc6o+LiYl1++eVKTExUt27dNG7cOH3++edux/pRM2fOVGZmZuQv0GVnZ+v99993O5YjM2bMkM/n05QpU9yOEpO8NDck780OL84Nyfuzw+25ERPFYv78+Zo6daqmTZumLVu2aPDgwcrPz1dlZaXb0U6ptrZWgwcPVklJidtRGmzt2rUqLCzU+vXrtXLlSp04cUIjR45UbW2t29FOKz09XTNmzNDmzZu1adMmjRgxQmPHjtX27dvdjtYg5eXlmj17tjIzM92OEpO8Njck780OL84Nyduzo0XMDRMDhg4dagoLCyProVDIpKWlmeLiYhdTNYwks2jRIrdjOFZZWWkkmbVr17odxZHOnTubl19+2e0YZ1RTU2MuvPBCs3LlSjNs2DBz//33ux0p5nh5bhjjzdnh1blhjDdmR0uZG56/YnH8+HFt3rxZeXl5kW1+v195eXlat26di8liW1VVlSQpJSXF5SQNEwqFNG/ePNXW1io7O9vtOGdUWFio0aNHRz2uYQ9zwx1emxuSt2ZHS5kbzf7pprYdPHhQoVBIqampUdtTU1P12WefuZQqtoXDYU2ZMkU5OTkaOHCg23F+1LZt25Sdna1jx46pY8eOWrRokfr37+92rB81b948bdmyReXl5W5HiVnMjebnpbkheW92tKS54fligeZXWFioTz75RB988IHbUc6ob9++2rp1q6qqqrRw4UIVFBRo7dq1LXZAVFRU6P7779fKlSuVkJDgdhzAGi/NDclbs6OlzQ3PF4suXbooLi5O+/fvj9q+f/9+de/e3aVUsWvy5MlaunSpysrKlJ6e7nacM4qPj1efPn0kSUOGDFF5ebmef/55zZ492+Vkp7Z582ZVVlbqsssui2wLhUIqKyvTCy+8oGAwqLi4OBcTxgbmRvPy2tyQvDU7Wtrc8PxrLOLj4zVkyBCtWrUqsi0cDmvVqlUt/vkwLzHGaPLkyVq0aJFWr16tCy64wO1IjRIOhxUMBt2OcVq5ubnatm2btm7dGlmysrJ06623auvWrZQKS5gbzSNW5obUsmdHS5sbnr9iIUlTp05VQUGBsrKyNHToUD333HOqra3VxIkT3Y52SkeOHNGuXbsi63v27NHWrVuVkpKiHj16uJjs9AoLCzV37lwtWbJEiYmJ2rdvnyQpOTlZ7dq1czndqRUVFWnUqFHq0aOHampqNHfuXK1Zs0bLly93O9ppJSYm1nv+uUOHDjrnnHM88by0l3htbkjemx1enBuS92ZHi5sbrrwXpQn8+c9/Nj169DDx8fFm6NChZv369W5HOq3S0lIjqd5SUFDgdrTTOlVeSea1115zO9pp3XHHHaZnz54mPj7edO3a1eTm5poVK1a4Hcsx3m7adLw0N4zx3uzw4twwJjZmh5tzw2eMMc1ZZAAAQOzy/GssAABAy0GxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDX/C/xeRF0C8g89AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    plt.subplot(1,2,i+1)\n",
    "    plt.imshow(M[i+1].G)\n",
    "    plt.title(M[i+1].name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model fitting\n",
    "Now let's fit the models to individual data set. There are three ways to do this. We can fit the models\n",
    "\n",
    "* to each individual participant with it's own parameters $\\theta$\n",
    "* to each all participants together with shared parameters, but with an individual parameter for the signal strength and for group. \n",
    "* in a cross-participant crossvalidated fashion. The models are fit to N-1 subjects and evaluated on the Nth subject. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do the individual fits - suppress verbose printout\n",
    "T_in, theta_in = pcm.fit_model_individ(Y,M,fit_scale = True, verbose = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting group model 0\n",
      "Fitting group model 1\n",
      "Fitting group model 2\n",
      "Fitting group model 3\n",
      "Fitting group model 4\n"
     ]
    }
   ],
   "source": [
    "# Fit the model in to the full group, using a individual scaling parameter for each\n",
    "T_gr, theta_gr = pcm.fit_model_group(Y, M, fit_scale=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting group cross model 0\n",
      "Fitting group cross model 1\n",
      "Fitting group cross model 2\n",
      "Fitting group cross model 3\n",
      "Fitting group cross model 4\n"
     ]
    }
   ],
   "source": [
    "# crossvalidated likelihood\n",
    "T_cv, theta_cv = pcm.fit_model_group_crossval(Y, M, fit_scale=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting and interpreting the results\n",
    "The results are returned as a nested data frame with the likelihood, noise and scale parameter for each individuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th>variable</th>\n",
       "      <th colspan=\"5\" halign=\"left\">likelihood</th>\n",
       "      <th colspan=\"5\" halign=\"left\">noise</th>\n",
       "      <th colspan=\"5\" halign=\"left\">iterations</th>\n",
       "      <th colspan=\"5\" halign=\"left\">scale</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model</th>\n",
       "      <th>null</th>\n",
       "      <th>muscle</th>\n",
       "      <th>usage</th>\n",
       "      <th>muscle+usage</th>\n",
       "      <th>ceil</th>\n",
       "      <th>null</th>\n",
       "      <th>muscle</th>\n",
       "      <th>usage</th>\n",
       "      <th>muscle+usage</th>\n",
       "      <th>ceil</th>\n",
       "      <th>null</th>\n",
       "      <th>muscle</th>\n",
       "      <th>usage</th>\n",
       "      <th>muscle+usage</th>\n",
       "      <th>ceil</th>\n",
       "      <th>null</th>\n",
       "      <th>muscle</th>\n",
       "      <th>usage</th>\n",
       "      <th>muscle+usage</th>\n",
       "      <th>ceil</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-42231.412711</td>\n",
       "      <td>-41966.470799</td>\n",
       "      <td>-41786.672956</td>\n",
       "      <td>-41786.672927</td>\n",
       "      <td>-41689.860468</td>\n",
       "      <td>0.875853</td>\n",
       "      <td>0.871286</td>\n",
       "      <td>0.868482</td>\n",
       "      <td>0.868483</td>\n",
       "      <td>0.872297</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.109319</td>\n",
       "      <td>0.750145</td>\n",
       "      <td>0.786771</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.996688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-34965.171104</td>\n",
       "      <td>-34923.791342</td>\n",
       "      <td>-34915.406608</td>\n",
       "      <td>-34914.959612</td>\n",
       "      <td>-34889.042772</td>\n",
       "      <td>1.070401</td>\n",
       "      <td>1.067480</td>\n",
       "      <td>1.069075</td>\n",
       "      <td>1.068119</td>\n",
       "      <td>1.066987</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>0.045008</td>\n",
       "      <td>0.324407</td>\n",
       "      <td>0.322917</td>\n",
       "      <td>0.961247</td>\n",
       "      <td>0.997908</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-34767.538097</td>\n",
       "      <td>-34679.107626</td>\n",
       "      <td>-34632.643241</td>\n",
       "      <td>-34632.642946</td>\n",
       "      <td>-34571.750328</td>\n",
       "      <td>1.026408</td>\n",
       "      <td>1.021219</td>\n",
       "      <td>1.019122</td>\n",
       "      <td>1.019123</td>\n",
       "      <td>1.023312</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>0.059863</td>\n",
       "      <td>0.435484</td>\n",
       "      <td>0.463987</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.992282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-45697.970627</td>\n",
       "      <td>-45609.052395</td>\n",
       "      <td>-45448.518276</td>\n",
       "      <td>-45448.518254</td>\n",
       "      <td>-45225.784824</td>\n",
       "      <td>1.480699</td>\n",
       "      <td>1.479592</td>\n",
       "      <td>1.474026</td>\n",
       "      <td>1.474025</td>\n",
       "      <td>1.478701</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.173031</td>\n",
       "      <td>1.193770</td>\n",
       "      <td>1.235628</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.998173</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-31993.363827</td>\n",
       "      <td>-31866.288313</td>\n",
       "      <td>-31806.982719</td>\n",
       "      <td>-31806.982521</td>\n",
       "      <td>-31707.184232</td>\n",
       "      <td>0.808482</td>\n",
       "      <td>0.805621</td>\n",
       "      <td>0.805774</td>\n",
       "      <td>0.805774</td>\n",
       "      <td>0.807319</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.073935</td>\n",
       "      <td>0.516338</td>\n",
       "      <td>0.532421</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.999362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>-41817.234010</td>\n",
       "      <td>-41632.061473</td>\n",
       "      <td>-41543.438786</td>\n",
       "      <td>-41543.438769</td>\n",
       "      <td>-41439.111952</td>\n",
       "      <td>1.035696</td>\n",
       "      <td>1.031827</td>\n",
       "      <td>1.031649</td>\n",
       "      <td>1.031648</td>\n",
       "      <td>1.034879</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>0.116696</td>\n",
       "      <td>0.801114</td>\n",
       "      <td>0.828773</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.999324</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-50336.142592</td>\n",
       "      <td>-50201.799362</td>\n",
       "      <td>-50173.300358</td>\n",
       "      <td>-50173.300306</td>\n",
       "      <td>-50099.140711</td>\n",
       "      <td>1.479001</td>\n",
       "      <td>1.472401</td>\n",
       "      <td>1.474430</td>\n",
       "      <td>1.474428</td>\n",
       "      <td>1.476146</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>0.101477</td>\n",
       "      <td>0.714043</td>\n",
       "      <td>0.723969</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.979400</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "variable    likelihood                                            \\\n",
       "model             null        muscle         usage  muscle+usage   \n",
       "0        -42231.412711 -41966.470799 -41786.672956 -41786.672927   \n",
       "1        -34965.171104 -34923.791342 -34915.406608 -34914.959612   \n",
       "2        -34767.538097 -34679.107626 -34632.643241 -34632.642946   \n",
       "3        -45697.970627 -45609.052395 -45448.518276 -45448.518254   \n",
       "4        -31993.363827 -31866.288313 -31806.982719 -31806.982521   \n",
       "5        -41817.234010 -41632.061473 -41543.438786 -41543.438769   \n",
       "6        -50336.142592 -50201.799362 -50173.300358 -50173.300306   \n",
       "\n",
       "variable                   noise                                             \\\n",
       "model             ceil      null    muscle     usage muscle+usage      ceil   \n",
       "0        -41689.860468  0.875853  0.871286  0.868482     0.868483  0.872297   \n",
       "1        -34889.042772  1.070401  1.067480  1.069075     1.068119  1.066987   \n",
       "2        -34571.750328  1.026408  1.021219  1.019122     1.019123  1.023312   \n",
       "3        -45225.784824  1.480699  1.479592  1.474026     1.474025  1.478701   \n",
       "4        -31707.184232  0.808482  0.805621  0.805774     0.805774  0.807319   \n",
       "5        -41439.111952  1.035696  1.031827  1.031649     1.031648  1.034879   \n",
       "6        -50099.140711  1.479001  1.472401  1.474430     1.474428  1.476146   \n",
       "\n",
       "variable iterations                                     scale            \\\n",
       "model          null muscle usage muscle+usage  ceil      null    muscle   \n",
       "0               4.0    4.0   4.0          7.0  30.0  0.109319  0.750145   \n",
       "1               4.0    4.0   4.0         16.0  21.0  0.045008  0.324407   \n",
       "2               4.0    4.0   4.0          7.0  36.0  0.059863  0.435484   \n",
       "3               4.0    4.0   4.0          6.0  18.0  0.173031  1.193770   \n",
       "4               4.0    4.0   4.0          6.0  32.0  0.073935  0.516338   \n",
       "5               4.0    4.0   4.0          5.0  19.0  0.116696  0.801114   \n",
       "6               4.0    4.0   4.0         12.0  31.0  0.101477  0.714043   \n",
       "\n",
       "variable                                   \n",
       "model        usage muscle+usage      ceil  \n",
       "0         0.786771     1.000000  0.996688  \n",
       "1         0.322917     0.961247  0.997908  \n",
       "2         0.463987     1.000000  0.992282  \n",
       "3         1.235628     1.000000  0.998173  \n",
       "4         0.532421     1.000000  0.999362  \n",
       "5         0.828773     1.000000  0.999324  \n",
       "6         0.723969     1.000000  0.979400  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T_in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihoods are very negative and quite different across participants, which is expected (see documentation). What we need to interpret are the difference is the likelihood relative to a null model. \n",
    "We can visualized these using the model_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAANM1JREFUeJzt3XtUVPX+//HXAIIgt0AFSUA9ZOZJtNQIsyQl8ZL3+mZ51MrSFDQlzYN11Lwc+uq37FSmWR6xzKN1zCxLSy0ti9QwU7uoKYWmiGmAoILA/v3Rcv+avDEwOOP2+Vhrr8X+fD6z571lz/hi78/MthmGYQgAAMCiPFxdAAAAQE0i7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEvzcnUB7qCiokIHDx5UQECAbDabq8sBAACVYBiGjh8/roiICHl4nP/8DWFH0sGDBxUZGenqMgAAQBXs379fDRs2PG8/YUdSQECApN//sQIDA11cDQAAqIzCwkJFRkaa/4+fD2FHMi9dBQYGEnYAALjMXGwKChOUAQCApXFm5w+Ki4vtJimXlpbq9OnT8vLyko+Pj904SfL19TUnRJ0+fVqlpaXy9PRU7dq1qzT2xIkTMgxDtWvXlqenpySprKxMJSUl8vDwkK+vb5XGnjx5UhUVFfLx8ZGX1++/8vLycp06dcqhsTabTX5+fubYU6dOqby8XN7e3qpVq5bDYysqKnTy5ElJUp06dcyxJSUlKisrU61ateTt7e3wWMMwdOLECUmSn5/fWb9PR8ZW5nfvjOPkXL9PZxwnZ36f1T1O/vz7rO5xcr7fZ3WPkz/+Pqt7nJzv91nV44T3CN4jeI+omfeISjFgFBQUGJIMSUZeXp7ZPm3aNEOS8dBDD9mN9/PzMyQZ2dnZZtusWbMMScZ9991nN7Zu3bqGJGPnzp1m27x58wxJRq9evezGRkdHG5KMzZs3m22LFi0yJBmJiYl2Y5s3b25IMj755BOzbfny5YYko127dnZj27RpY0gyVq5cabZ99NFHhiSjZcuWdmM7dOhgSDLefPNNs23jxo2GJCMmJsZubLdu3QxJxoIFC8y2r7/+2pBkRERE2I296667DEnGiy++aLbt3r3bkGQEBQXZjR08eLAhyZgxY4bZduDAAUOS4eXlZTd2xIgRhiRj0qRJZttvv/1m/j5LS0vN9rFjxxqSjLFjx5ptpaWl5tjffvvNbJ80aZIhyRgxYoTd83l5eRmSjAMHDphtM2bMMCQZgwcPthsbFBRkSDJ2795ttr344ouGJOOuu+6yGxsREWFIMr7++muzbcGCBYYko1u3bnZjY2JiDEnGxo0bzbY333zTkGR06NDBbmzLli0NScZHH31ktq1cudKQZLRp08ZubLt27QxJxvLly822Tz75xJBkNG/e3G5sYmKiIclYtGiR2bZ582ZDkhEdHW03tlevXoYkY968eWbbzp07DUlG3bp17cbed999hiRj1qxZZlt2drYhyfDz87Mb+9BDDxmSjGnTpplteXl55u/zjx599FFDkjFhwgSzraioyBxbVFRktk+YMMGQZDz66KN22+A94ne8R/yO94jfufo94sz/3wUFBcaFcBkLAABYms0wDMPVRbhaYWGhgoKCdPDgQYWHh3OKmlPUnKJ2o1PUXMbiPaIyv0/eI67M94gz/38XFBRc8ANGhB2p0v9YAADAfVT2/28uYwEAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEsj7AAAAEvj3lgAAFwihmGYXyYo/f6FeRe7Yzeqz6VndubMmaPY2FgFBgYqMDBQ8fHxWrVqldl/6tQpJScnKzQ0VP7+/urXr58OHz5st42cnBx1795dfn5+ql+/vsaNG6eysrJLvSsAAFxUcXGxevXqZS5/DD6oOS4NOw0bNtTTTz+trKwsffXVV+rYsaN69eqlb7/9VpI0ZswYvffee3rrrbe0YcMGHTx4UH379jUfX15eru7du6u0tFRffPGFFi5cqIyMDE2cONFVuwQAANyM290uIiQkRDNnztRdd92levXqafHixbrrrrskST/88IOuu+46ZWZm6uabb9aqVat055136uDBgwoLC5MkzZ07V+PHj9eRI0fMe5tcDLeLAKyJSwZwN0VFRerVq5e5vmLFCvn7+7uwosvbZXe7iPLyci1ZskTFxcWKj49XVlaWTp8+rcTERHNMs2bNFBUVpczMTElSZmamWrRoYQYdSUpKSlJhYaF5duhcSkpKVFhYaLcAsB4uGQCQ3CDs7NixQ/7+/vLx8dEjjzyi5cuXq3nz5srNzZW3t7eCg4PtxoeFhSk3N1eSlJubaxd0zvSf6Tuf9PR0BQUFmUtkZKRzdwoAALgNl4eda6+9Vtu2bdOmTZs0fPhwDR48WN99912NPmdaWpoKCgrMZf/+/TX6fAAAwHVc/tFzb29vxcTESJJat26tLVu26F//+pfuuecelZaWKj8/3+7szuHDhxUeHi5JCg8P1+bNm+22d+bTWmfGnIuPj498fHycvCcAAMAdufzMzp9VVFSopKRErVu3Vq1atbRu3Tqzb9euXcrJyVF8fLwkKT4+Xjt27FBeXp45Zs2aNQoMDFTz5s0vee0AAMD9uPTMTlpamrp27aqoqCgdP35cixcv1vr16/Xhhx8qKChIQ4YMUWpqqkJCQhQYGKiRI0cqPj5eN998sySpc+fOat68uQYOHKgZM2YoNzdXTz75pJKTkzlzAwAAJLk47OTl5WnQoEE6dOiQgoKCFBsbqw8//FB33HGHJGnWrFny8PBQv379VFJSoqSkJL300kvm4z09PbVy5UoNHz5c8fHxqlOnjgYPHqwpU6a4apcAAICbcWnYmT9//gX7a9eurdmzZ2v27NnnHRMdHa0PPvjA2aUBAACLcLs5OwAAAM5E2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJZG2AEAAJbm0rCTnp6utm3bKiAgQPXr11fv3r21a9cuuzEJCQmy2Wx2yyOPPGI3JicnR927d5efn5/q16+vcePGqays7FLuCgAAcFNernzyDRs2KDk5WW3btlVZWZkmTJigzp0767vvvlOdOnXMcQ8//LCmTJlirvv5+Zk/l5eXq3v37goPD9cXX3yhQ4cOadCgQapVq5b++c9/XtL9AQAA7selYWf16tV26xkZGapfv76ysrJ02223me1+fn4KDw8/5zY++ugjfffdd1q7dq3CwsLUqlUrTZ06VePHj9fkyZPl7e1do/sAAADcm1vN2SkoKJAkhYSE2LW/8cYbqlu3rq6//nqlpaXpxIkTZl9mZqZatGihsLAwsy0pKUmFhYX69ttvz/k8JSUlKiwstFsAAIA1ufTMzh9VVFRo9OjRuuWWW3T99deb7ffdd5+io6MVERGh7du3a/z48dq1a5fefvttSVJubq5d0JFkrufm5p7zudLT0/XUU0/V0J4AAAB34jZhJzk5WTt37tTGjRvt2ocOHWr+3KJFCzVo0ECdOnXS3r179Ze//KVKz5WWlqbU1FRzvbCwUJGRkVUrHAAAuDW3uIyVkpKilStX6pNPPlHDhg0vODYuLk6S9OOPP0qSwsPDdfjwYbsxZ9bPN8/Hx8dHgYGBdgsAALAml4YdwzCUkpKi5cuX6+OPP1bjxo0v+pht27ZJkho0aCBJio+P144dO5SXl2eOWbNmjQIDA9W8efMaqRsAAFw+XHoZKzk5WYsXL9aKFSsUEBBgzrEJCgqSr6+v9u7dq8WLF6tbt24KDQ3V9u3bNWbMGN12222KjY2VJHXu3FnNmzfXwIEDNWPGDOXm5urJJ59UcnKyfHx8XLl7AADADbj0zM6cOXNUUFCghIQENWjQwFyWLl0qSfL29tbatWvVuXNnNWvWTI899pj69eun9957z9yGp6enVq5cKU9PT8XHx+tvf/ubBg0aZPe9PAAA4Mrl0jM7hmFcsD8yMlIbNmy46Haio6P1wQcfOKssAABgIW4xQRkAAKCmEHYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClEXYAAIClebm6AAA1p/W411xdgkvZykoV9If1hH8skeHl7bJ6XC1r5iBXl8AxyTFp51Idk5zZAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAluZQ2CkrK9OUKVN04MCBmqoHAADAqRwKO15eXpo5c6bKyspqqh4AAACncvgyVseOHbVhw4aaqAUAAMDpHP4G5a5du+rvf/+7duzYodatW6tOnTp2/T179nRacQAAANXlcNgZMWKEJOnZZ589q89ms6m8vLz6VQEAADiJw2GnoqKiJuoAAACoES796Hl6erratm2rgIAA1a9fX71799auXbvsxpw6dUrJyckKDQ2Vv7+/+vXrp8OHD9uNycnJUffu3eXn56f69etr3LhxTKIGAACSqhh2NmzYoB49eigmJkYxMTHq2bOnPvvssyptJzk5WV9++aXWrFmj06dPq3PnziouLjbHjBkzRu+9957eeustbdiwQQcPHlTfvn3N/vLycnXv3l2lpaX64osvtHDhQmVkZGjixIlV2TUAAGAxDoedRYsWKTExUX5+fho1apRGjRolX19fderUSYsXL3ZoW6tXr9b999+vv/71r2rZsqUyMjKUk5OjrKwsSVJBQYHmz5+vZ599Vh07dlTr1q21YMECffHFF/ryyy8lSR999JG+++47LVq0SK1atVLXrl01depUzZ49W6WlpY7uHgAAsBiHw8706dM1Y8YMLV261Aw7S5cu1dNPP62pU6dWq5iCggJJUkhIiCQpKytLp0+fVmJiojmmWbNmioqKUmZmpiQpMzNTLVq0UFhYmDkmKSlJhYWF+vbbb8/5PCUlJSosLLRbAACANTkcdvbt26cePXqc1d6zZ09lZ2dXuZCKigqNHj1at9xyi66//npJUm5urry9vRUcHGw3NiwsTLm5ueaYPwadM/1n+s4lPT1dQUFB5hIZGVnlugEAgHtzOOxERkZq3bp1Z7WvXbu2WqEhOTlZO3fu1JIlS6q8jcpKS0tTQUGBuezfv7/GnxMAALiGwx89f+yxxzRq1Cht27ZN7dq1kyR9/vnnysjI0L/+9a8qFZGSkqKVK1fq008/VcOGDc328PBwlZaWKj8/3+7szuHDhxUeHm6O2bx5s932znxa68yYP/Px8ZGPj0+VagUAAJcXh8/sDB8+XEuWLNGOHTs0evRojR49Wjt37tTSpUs1bNgwh7ZlGIZSUlK0fPlyffzxx2rcuLFdf+vWrVWrVi27M0m7du1STk6O4uPjJUnx8fHasWOH8vLyzDFr1qxRYGCgmjdv7ujuAQAAi3H4zI4k9enTR3369Kn2kycnJ2vx4sVasWKFAgICzDk2QUFB8vX1VVBQkIYMGaLU1FSFhIQoMDBQI0eOVHx8vG6++WZJUufOndW8eXMNHDhQM2bMUG5urp588kklJydz9gYAADh+ZqdJkyY6evToWe35+flq0qSJQ9uaM2eOCgoKlJCQoAYNGpjL0qVLzTGzZs3SnXfeqX79+um2225TeHi43n77bbPf09NTK1eulKenp+Lj4/W3v/1NgwYN0pQpUxzdNQAAYEEOn9n56aefznn/q5KSEv3yyy8ObcswjIuOqV27tmbPnq3Zs2efd0x0dLQ++OADh54bAABcGSoddt59913z5w8//FBBQUHmenl5udatW6dGjRo5tTgAAIDqqnTY6d27t6Tf72w+ePBgu75atWqpUaNGeuaZZ5xaHAAAQHVVOuycudt548aNtWXLFtWtW7fGigIAAHAWh+fsVOdbkgEAAC41hz+NNWrUKD3//PNntb/44osaPXq0M2oCAABwGofDzrJly3TLLbec1d6uXTv997//dUpRAAAAzuJw2Dl69KjdJ7HOCAwM1K+//uqUogAAAJzF4bATExOj1atXn9W+atUqh79UENZhGIaKiorMpTLfoQQAwKXg8ATl1NRUpaSk6MiRI+rYsaMkad26dXrmmWf03HPPObs+XCaKi4vVq1cvc33FihXy9/d3YUUAAPzO4bDz4IMPqqSkRNOnT9fUqVMlSY0aNdKcOXM0aNAgpxcIAABQHVW6Eejw4cM1fPhwHTlyRL6+vvwFDwAA3FaVws4Z9erVc1YdAAAANaJKYee///2v3nzzTeXk5Ki0tNSub+vWrU4pDAAAwBkc/jTW888/rwceeEBhYWH6+uuvddNNNyk0NFT79u1T165da6JGAACAKnM47Lz00kuaN2+eXnjhBXl7e+vxxx/XmjVrNGrUKBUUFNREjQAAAFXmcNjJyclRu3btJEm+vr46fvy4JGngwIH6z3/+49zqAAAAqsnhsBMeHq5jx45JkqKiovTll19K+v0GoXyRHAAAcDcOh52OHTvq3XfflSQ98MADGjNmjO644w7dc8896tOnj9MLBAAAqA6HP401b948VVRUSJKSk5MVGhqqL774Qj179tSwYcOcXiAAAEB1VDrsREVF6euvv1ZoaKg8PDz04osvatCgQerfv7/69+9fkzUCAABUWaUvYx04cEDl5eXm+oQJE7jLOQAAcHsOz9k5g8nIAADgclDlsAMAAHA5cGiC8quvvmre9LOsrEwZGRmqW7eu3ZhRo0Y5rzoAAIBqcmiC8iuvvGKuh4eH6/XXX7cbY7PZCDsAAMCtVDrs/PTTTzVYBgAAQM1gzg4AALA0wg4AALA0wg4AALA0wg4AALA0wg4AALA0h8PO1q1btWPHDnN9xYoV6t27tyZMmKDS0lKnFgcAAFBdDoedYcOGaffu3ZKkffv2qX///vLz89Nbb72lxx9/3OkFAgAAVIfDYWf37t1q1aqVJOmtt97SbbfdpsWLFysjI0PLli1zdn0AAADV4nDYMQxDFRUVkqS1a9eqW7dukqTIyEjugg4AANyOw2GnTZs2mjZtml5//XVt2LBB3bt3lyRlZ2crLCzM6QUCAABUh8Nh57nnntPWrVuVkpKiJ554QjExMZKk//73v2rXrp3TCwQAAKgOh+56LkmxsbF2n8Y6Y+bMmfL09HRKUQAAAM5Spe/Zyc/P16uvvqq0tDQdO3ZMkvTdd98pLy/PqcUBAABUl8NndrZv365OnTopODhYP/30kx5++GGFhITo7bffVk5Ojl577bWaqBMAAKBKHD6zk5qaqgceeEB79uxR7dq1zfZu3brp008/dWpxAAAA1eVw2NmyZYuGDRt2VvvVV1+t3NxcpxQFAADgLA6HHR8fHxUWFp7Vvnv3btWrV88pRQEAADiLw2GnZ8+emjJlik6fPi1JstlsysnJ0fjx49WvXz+nFwgAAFAdDoedZ555RkVFRapfv75OnjypDh06KCYmRgEBAZo+fXpN1AgAAFBlDn8aKygoSGvWrNHGjRu1fft2FRUV6cYbb1RiYmJN1AcAAFAtDoedffv2qUmTJmrfvr3at29fEzUBAAA4jcOXsWJiYnT77bdr0aJFOnXqVE3UBAAA4DQOh52tW7cqNjZWqampCg8P17Bhw7Rp06YqPfmnn36qHj16KCIiQjabTe+8845d//333y+bzWa3dOnSxW7MsWPHNGDAAAUGBio4OFhDhgxRUVFRleoBAADW43DYadWqlf71r3/p4MGD+ve//61Dhw7p1ltv1fXXX69nn31WR44cqfS2iouL1bJlS82ePfu8Y7p06aJDhw6Zy3/+8x+7/gEDBujbb7/VmjVrtHLlSn366acaOnSoo7sFAAAsqkr3xpIkLy8v9e3bV2+99Zb+93//Vz/++KPGjh2ryMhIDRo0SIcOHbroNrp27app06apT58+5x3j4+Oj8PBwc7nqqqvMvu+//16rV6/Wq6++qri4OLVv314vvPCClixZooMHD553myUlJSosLLRbAACANVU57Hz11VcaMWKEGjRooGeffVZjx47V3r17tWbNGh08eFC9evVySoHr169X/fr1de2112r48OE6evSo2ZeZmang4GC1adPGbEtMTJSHh8cFL62lp6crKCjIXCIjI51SKwAAcD8Ofxrr2Wef1YIFC7Rr1y5169ZNr732mrp16yYPj99zU+PGjZWRkaFGjRpVu7guXbqob9++aty4sfbu3asJEyaoa9euyszMlKenp3Jzc1W/fn37HfLyUkhIyAVvXZGWlqbU1FRzvbCwkMADAIBFORx25syZowcffFD333+/GjRocM4x9evX1/z586tdXP/+/c2fW7RoodjYWP3lL3/R+vXr1alTpypv18fHRz4+PtWuDwAAuD+Hw86ePXsuOsbb21uDBw+uUkEX0qRJE9WtW1c//vijOnXqpPDwcOXl5dmNKSsr07FjxxQeHu705wcAAJcfh8POGSdOnFBOTo5KS0vt2mNjY6td1PkcOHBAR48eNc8oxcfHKz8/X1lZWWrdurUk6eOPP1ZFRYXi4uJqrA4AAHD5cDjsHDlyRPfff79Wr159zv7y8vJKb6uoqEg//vijuZ6dna1t27YpJCREISEheuqpp9SvXz+Fh4dr7969evzxxxUTE6OkpCRJ0nXXXacuXbro4Ycf1ty5c3X69GmlpKSof//+ioiIcHTXAACABTn8aazRo0eroKBAmzZtkq+vr1avXq2FCxfqmmuu0bvvvuvQtr766ivdcMMNuuGGGyRJqampuuGGGzRx4kR5enpq+/bt6tmzp5o2baohQ4aodevW+uyzz+zm27zxxhtq1qyZOnXqpG7duql9+/aaN2+eo7sFAAAsyuEzOx9//LFWrFihNm3ayMPDQ9HR0brjjjsUGBio9PR0de/evdLbSkhIkGEY5+3/8MMPL7qNkJAQLV68uNLPCQAAriwOn9kpLi42P+591VVXmd+Y3KJFC23dutW51QEAAFSTw2Hn2muv1a5duyRJLVu21Msvv6xffvlFc+fOPe9H0QEAAFzF4ctYjz76qHkriEmTJqlLly5644035O3trYyMDGfXBwAAUC0Oh52//e1v5s+tW7fWzz//rB9++EFRUVGqW7euU4sDAACorip/z84Zfn5+uvHGG51RCwAAgNM5NGdnz549WrZsmbKzsyVJ77//vm677Ta1bdtW06dPv+AnqwAAAFyh0md2li9frv/5n/+Rh4eHbDab5s2bp2HDhikhIUGBgYGaPHmyvLy8NH78+JqsFwAAwCGVPrMzffp0Pf744zp16pTmzJmjRx55ROnp6Vq1apVWrlyp2bNnM0EZAAC4nUqHnV27dunBBx+UzWbT4MGDVVpaqsTERLO/c+fO+vnnn2ukSAAAgKqqdNgpLi5WQEDA7w/y8JCvr6/8/PzMfl9fX5WUlDi/QgAAgGqodNix2Wyy2WznXQcAAHBHlZ6gbBiGmjZtagacoqIi3XDDDfLw8DD7AQAA3E2lw86CBQtqso7LXutxr7m6BJeylZUq6A/rCf9YIsPL22X1uFrWzEGuLgGAGzI8a6kg9l67ddS8SoedwYMH12QdAABYn812Rf8h6CrV/gZlAHBX/BUNQCLsALAy/ooGIAdvFwEAAHC5IewAAABLI+wAAABLc3jOTmpq6jnbbTabateurZiYGPXq1UshISHVLg4AAKC6HA47X3/9tbZu3ary8nJde+21kqTdu3fL09NTzZo100svvaTHHntMGzduVPPmzZ1eMAAAgCMcvozVq1cvJSYm6uDBg8rKylJWVpYOHDigO+64Q/fee69++eUX3XbbbRozZkxN1AsAAOAQh8POzJkzNXXqVAUGBpptQUFBmjx5smbMmCE/Pz9NnDhRWVlZTi0UAACgKhwOOwUFBcrLyzur/ciRIyosLJQkBQcHq7S0tPrVAQAAVFOVLmM9+OCDWr58uQ4cOKADBw5o+fLlGjJkiHr37i1J2rx5s5o2bersWgEAABzm8ATll19+WWPGjFH//v1VVlb2+0a8vDR48GDNmjVLktSsWTO9+uqrzq0UAACgChwOO/7+/nrllVc0a9Ys7du3T5LUpEkT+fv7m2NatWrltAIBAACqo8r3xvL39ze/S+ePQQcAAMCdODxnp6KiQlOmTFFQUJCio6MVHR2t4OBgTZ06VRUVFTVRIwAAQJU5fGbniSee0Pz58/X000/rlltukSRt3LhRkydP1qlTpzR9+nSnFwkAAFBVDoedhQsX6tVXX1XPnj3NttjYWF199dUaMWIEYQcAALgVhy9jHTt2TM2aNTurvVmzZjp27JhTigIAAHAWh8NOy5Yt9eKLL57V/uKLL6ply5ZOKQoAAMBZHL6MNWPGDHXv3l1r165VfHy8JCkzM1P79+/XBx984PQCAQAAqsPhMzsdOnTQ7t271adPH+Xn5ys/P199+/bVrl27dOutt9ZEjQAAAFVWpe/ZiYiIOGsi8oEDBzR06FDNmzfPKYUBAAA4g8Nnds7n6NGjmj9/vrM2BwAA4BROCzsAAADuiLADAAAsjbADAAAsrdITlPv27XvB/vz8/OrWAgAA4HSVDjtBQUEX7R80aFC1CwIAAHCmSoedBQsW1GQdAAAANYI5OwAAwNIIOwAAwNIIOwAAwNIIOwAAwNJcGnY+/fRT9ejRQxEREbLZbHrnnXfs+g3D0MSJE9WgQQP5+voqMTFRe/bssRtz7NgxDRgwQIGBgQoODtaQIUNUVFR0CfcCAAC4M5eGneLiYrVs2VKzZ88+Z/+MGTP0/PPPa+7cudq0aZPq1KmjpKQknTp1yhwzYMAAffvtt1qzZo1WrlypTz/9VEOHDr1UuwAAANxcle567ixdu3ZV165dz9lnGIaee+45Pfnkk+rVq5ck6bXXXlNYWJjeeecd9e/fX99//71Wr16tLVu2qE2bNpKkF154Qd26ddP//d//KSIi4pzbLikpUUlJibleWFjo5D0DAADuwm3n7GRnZys3N1eJiYlmW1BQkOLi4pSZmSlJyszMVHBwsBl0JCkxMVEeHh7atGnTebednp6uoKAgc4mMjKy5HQEAAC7ltmEnNzdXkhQWFmbXHhYWZvbl5uaqfv36dv1eXl4KCQkxx5xLWlqaCgoKzGX//v1Orh4AALgLl17GchUfHx/5+Pi4ugwAAHAJuO2ZnfDwcEnS4cOH7doPHz5s9oWHhysvL8+uv6ysTMeOHTPHAACAK5vbhp3GjRsrPDxc69atM9sKCwu1adMmxcfHS5Li4+OVn5+vrKwsc8zHH3+siooKxcXFXfKaAQCA+3HpZayioiL9+OOP5np2dra2bdumkJAQRUVFafTo0Zo2bZquueYaNW7cWP/4xz8UERGh3r17S5Kuu+46denSRQ8//LDmzp2r06dPKyUlRf379z/vJ7EAAMCVxaVh56uvvtLtt99urqempkqSBg8erIyMDD3++OMqLi7W0KFDlZ+fr/bt22v16tWqXbu2+Zg33nhDKSkp6tSpkzw8PNSvXz89//zzl3xfAACAe3Jp2ElISJBhGOftt9lsmjJliqZMmXLeMSEhIVq8eHFNlAcAACzAbefsAAAAOANhBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWBphBwAAWJpLv0EZ1mF41lJB7L126wAAuAPCDpzDZpPh5e3qKgAAOAuXsQAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKURdgAAgKW5ddiZPHmybDab3dKsWTOz/9SpU0pOTlZoaKj8/f3Vr18/HT582IUVAwAAd+PWYUeS/vrXv+rQoUPmsnHjRrNvzJgxeu+99/TWW29pw4YNOnjwoPr27evCagEAgLvxcnUBF+Pl5aXw8PCz2gsKCjR//nwtXrxYHTt2lCQtWLBA1113nb788kvdfPPNl7pUAADghtz+zM6ePXsUERGhJk2aaMCAAcrJyZEkZWVl6fTp00pMTDTHNmvWTFFRUcrMzLzgNktKSlRYWGi3AAAAa3LrsBMXF6eMjAytXr1ac+bMUXZ2tm699VYdP35cubm58vb2VnBwsN1jwsLClJube8HtpqenKygoyFwiIyNrcC8AAIArufVlrK5du5o/x8bGKi4uTtHR0XrzzTfl6+tb5e2mpaUpNTXVXC8sLCTwAABgUW59ZufPgoOD1bRpU/34448KDw9XaWmp8vPz7cYcPnz4nHN8/sjHx0eBgYF2CwAAsKbLKuwUFRVp7969atCggVq3bq1atWpp3bp1Zv+uXbuUk5Oj+Ph4F1YJAADciVtfxho7dqx69Oih6OhoHTx4UJMmTZKnp6fuvfdeBQUFaciQIUpNTVVISIgCAwM1cuRIxcfH80ksAABgcuuwc+DAAd177706evSo6tWrp/bt2+vLL79UvXr1JEmzZs2Sh4eH+vXrp5KSEiUlJemll15ycdUAAMCduHXYWbJkyQX7a9eurdmzZ2v27NmXqCIAAHC5uazm7AAAADiKsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACyNsAMAACzNMmFn9uzZatSokWrXrq24uDht3rzZ1SUBAAA3YImws3TpUqWmpmrSpEnaunWrWrZsqaSkJOXl5bm6NAAA4GKWCDvPPvusHn74YT3wwANq3ry55s6dKz8/P/373/92dWkAAMDFvFxdQHWVlpYqKytLaWlpZpuHh4cSExOVmZl5zseUlJSopKTEXC8oKJAkFRYWVrmO8pKTVX4srKc6x5IzcVzij9zhuOSYxB9V95g883jDMC447rIPO7/++qvKy8sVFhZm1x4WFqYffvjhnI9JT0/XU089dVZ7ZGRkjdSIK0/QC4+4ugTgLByXcDfOOiaPHz+uoKCg8/Zf9mGnKtLS0pSammquV1RU6NixYwoNDZXNZnNhZZe3wsJCRUZGav/+/QoMDHR1OYAkjku4H45J5zEMQ8ePH1dERMQFx132Yadu3bry9PTU4cOH7doPHz6s8PDwcz7Gx8dHPj4+dm3BwcE1VeIVJzAwkBcw3A7HJdwNx6RzXOiMzhmX/QRlb29vtW7dWuvWrTPbKioqtG7dOsXHx7uwMgAA4A4u+zM7kpSamqrBgwerTZs2uummm/Tcc8+puLhYDzzwgKtLAwAALmaJsHPPPffoyJEjmjhxonJzc9WqVSutXr36rEnLqFk+Pj6aNGnSWZcIAVfiuIS74Zi89GzGxT6vBQAAcBm77OfsAAAAXAhhBwAAWBphBwAAWBphBy6XkJCg0aNHu7oMAHAY71+XB8IOAACwNMIOAAAukpGRoYSEBFeXYXmEnStQQkKCRo4cqdGjR+uqq65SWFiYXnnlFfOLGAMCAhQTE6NVq1ZJ+v3F+Ofbabzzzjt29xH75ptvdPvttysgIECBgYFq3bq1vvrqK7P/888/V0JCgvz8/HTVVVcpKSlJv/322znrKykp0dixY3X11VerTp06iouL0/r1653+7wD31ahRIz333HN2ba1atdLkyZNlGIYmT56sqKgo+fj4KCIiQqNGjTLHvf7662rTpo0CAgIUHh6u++67T3l5eXbbevfdd3XNNdeodu3auv3227Vw4ULZbDbl5+ebYzZu3Khbb71Vvr6+ioyM1KhRo1RcXFyTu41KuNLevy70WpDE66GSCDtXqIULF6pu3bravHmzRo4cqeHDh+vuu+9Wu3bttHXrVnXu3FkDBw7UiRMnKrW9AQMGqGHDhtqyZYuysrL097//XbVq1ZIkbdu2TZ06dVLz5s2VmZmpjRs3qkePHiovLz/ntlJSUpSZmaklS5Zo+/btuvvuu9WlSxft2bPHafuPy9eyZcs0a9Ysvfzyy9qzZ4/eeecdtWjRwuw/ffq0pk6dqm+++UbvvPOOfvrpJ91///1mf3Z2tu666y717t1b33zzjYYNG6YnnnjC7jn27t2rLl26qF+/ftq+fbuWLl2qjRs3KiUl5VLtJi6A96//j9dDJRm44nTo0MFo3769uV5WVmbUqVPHGDhwoNl26NAhQ5KRmZlpLFiwwAgKCrLbxvLly40/Hj4BAQFGRkbGOZ/v3nvvNW655ZYL1vPoo48ahmEYP//8s+Hp6Wn88ssvdmM6depkpKWlVXYXcZmLjo42Zs2aZdfWsmVLY9KkScYzzzxjNG3a1CgtLa3UtrZs2WJIMo4fP24YhmGMHz/euP766+3GPPHEE4Yk47fffjMMwzCGDBliDB061G7MZ599Znh4eBgnT56s2k7BKaz2/rVgwQKjQ4cO593+hV4LhmHweqgkzuxcoWJjY82fPT09FRoaavfXwJlbbfz5dOf5pKam6qGHHlJiYqKefvpp7d271+w785dRZezYsUPl5eVq2rSp/P39zWXDhg1228SV6+6779bJkyfVpEkTPfzww1q+fLnKysrM/qysLPXo0UNRUVEKCAhQhw4dJEk5OTmSpF27dqlt27Z227zpppvs1r/55htlZGTYHYNJSUmqqKhQdnZ2De8hLuZyfv/Kycmx63vkkUf02Wef2bX985//rNTzSbweKssS98aC486coj3DZrPZtZ25nl1RUSEPDw8Zf7qryOnTp+3WJ0+erPvuu0/vv/++Vq1apUmTJmnJkiXq06ePfH19K11XUVGRPD09lZWVJU9PT7s+f3//Sm8Hl7cLHXORkZHatWuX1q5dqzVr1mjEiBGaOXOmNmzYoNLSUiUlJSkpKUlvvPGG6tWrp5ycHCUlJam0tLTSz19UVKRhw4bZzX04Iyoqqno7h2q7nN+/IiIitG3bNrP97bff1rJly/TGG2+YbSEhIebPF6uf10PlEHZwUfXq1dPx48dVXFysOnXqSJLdi/WMpk2bqmnTphozZozuvfdeLViwQH369FFsbKzWrVunp5566qLPdcMNN6i8vFx5eXm69dZbnb0ruEzUq1dPhw4dMtcLCwvt/oL09fVVjx491KNHDyUnJ6tZs2basWOHDMPQ0aNH9fTTTysyMlKS7CaaStK1116rDz74wK5ty5Ytdus33nijvvvuO8XExDh713CJudv7l5eXl91xVb9+ffn6+p73WLvYa0Hi9VAZXMbCRcXFxcnPz08TJkzQ3r17tXjxYmVkZJj9J0+eVEpKitavX6+ff/5Zn3/+ubZs2aLrrrtOkpSWlqYtW7ZoxIgR2r59u3744QfNmTNHv/7661nP1bRpUw0YMECDBg3S22+/rezsbG3evFnp6el6//33L9Uuw8U6duyo119/XZ999pl27NihwYMHm38pZ2RkaP78+dq5c6f27dunRYsWydfXV9HR0YqKipK3t7deeOEF7du3T++++66mTp1qt+1hw4bphx9+0Pjx47V79269+eab5vF85ozA+PHj9cUXXyglJUXbtm3Tnj17tGLFistrQiYkXf7vXxd6LUi8HirNlROG4Bp/nFB3xrkmwUkyli9fbhjG7xP6YmJiDF9fX+POO+805s2bZ07wKykpMfr3729ERkYa3t7eRkREhJGSkmI3cW39+vVGu3btDB8fHyM4ONhISkoyJ7/9uZ7S0lJj4sSJRqNGjYxatWoZDRo0MPr06WNs377d2f8UcFMFBQXGPffcYwQGBhqRkZFGRkaGOSlz+fLlRlxcnBEYGGjUqVPHuPnmm421a9eaj128eLHRqFEjw8fHx4iPjzfeffddQ5Lx9ddfm2NWrFhhxMTEGD4+PkZCQoIxZ84cQ5LdMbt582bjjjvuMPz9/Y06deoYsbGxxvTp0y/lPwPOwWrvXxeboHyh18KZfeP1cHE2w/jTxUAAuMJMnz5dc+fO1f79+11dCuByVnw9MGcHwBXnpZdeUtu2bRUaGqrPP/9cM2fOvLxOyQNOdCW8Hgg7AK44e/bs0bRp03Ts2DFFRUXpscceU1pamqvLAlziSng9cBkLAABYGp/GAgAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAQAAlkbYAXBFSkhI0OjRoys9PiMjQ8HBwTVWD4CaQ9gBAACWRtgBAACWRtgB4FYSEhI0cuRIjR49WldddZXCwsL0yiuvqLi4WA888IACAgIUExOjVatWmY/ZsGGDbrrpJvn4+KhBgwb6+9//rrKyMrO/uLhYgwYNkr+/vxo0aKBnnnnmrOctKSnR2LFjdfXVV6tOnTqKi4vT+vXrL8UuA6hhhB0AbmfhwoWqW7euNm/erJEjR2r48OG6++671a5dO23dulWdO3fWwIEDdeLECf3yyy/q1q2b2rZtq2+++UZz5szR/PnzNW3aNHN748aN04YNG7RixQp99NFHWr9+vbZu3Wr3nCkpKcrMzNSSJUu0fft23X333erSpYv27NlzqXcfgJNxuwgAbiUhIUHl5eX67LPPJEnl5eUKCgpS37599dprr0mScnNz1aBBA2VmZuq9997TsmXL9P3338tms0n6/caG48ePV0FBgU6cOKHQ0FAtWrRId999tyTp2LFjatiwoYYOHarnnntOOTk5atKkiXJychQREWHWkpiYqJtuukn//Oc/lZGRodGjRys/P//S/oMAqDZuBArA7cTGxpo/e3p6KjQ0VC1atDDbwsLCJEl5eXn6/vvvFR8fbwYdSbrllltUVFSkAwcO6LffflNpaani4uLM/pCQEF177bXm+o4dO1ReXq6mTZva1VFSUqLQ0FCn7x+AS4uwA8Dt1KpVy27dZrPZtZ0JNhUVFU55vqKiInl6eiorK0uenp52ff7+/k55DgCuQ9gBcFm77rrrtGzZMhmGYYagzz//XAEBAWrYsKFCQkJUq1Ytbdq0SVFRUZKk3377Tbt371aHDh0kSTfccIPKy8uVl5enW2+91WX7AqBmMEEZwGVtxIgR2r9/v0aOHKkffvhBK1as0KRJk5SamioPDw/5+/tryJAhGjdunD7++GPt3LlT999/vzw8/v/bX9OmTTVgwAANGjRIb7/9trKzs7V582alp6fr/fffd+HeAXAGzuwAuKxdffXV+uCDDzRu3Di1bNlSISEhGjJkiJ588klzzMyZM1VUVKQePXooICBAjz32mAoKCuy2s2DBAk2bNk2PPfaYfvnlF9WtW1c333yz7rzzzku9SwCcjE9jAQAAS+MyFgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsDTCDgAAsLT/B8/88p62YOxtAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = pcm.model_plot(T_in.likelihood,\n",
    "                        null_model = 'null',\n",
    "                        noise_ceiling= 'ceil')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem with the noise ceiling is that it is individually fit to each subject. It has much more parameters than the models it is competing against, so it is overfitting. To compare models with different numbers of parameters directly, we need to look at our cross-validated group fits. The group fits can be used as an upper noise ceiling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMdtJREFUeJzt3XlcVfW+//E3oCAIGwJlSnCInG6i5oA4FCmJQ863m2Zq5k1T0KOUebCOmmZ09KQ2mGZ5xEqujWZZWWpXysIhzLQsNbXQFDANEFIQWL8/+rlvO4fYsHFvl6/n47EeD9b3+91rfRaPtfHtWt+9l5thGIYAAABMyt3ZBQAAANQkwg4AADA1wg4AADA1wg4AADA1wg4AADA1wg4AADA1wg4AADC1Ws4uwBVUVFTo2LFj8vPzk5ubm7PLAQAAlWAYhk6fPq3w8HC5u1/6+g1hR9KxY8cUERHh7DIAAEAVHDlyRA0aNLhkP2FHkp+fn6Tff1kWi8XJ1QAAgMooLCxURESE9d/xSyHsSNZbVxaLhbADAMBV5q+moDBBGQAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBpPPQcA4AoxDEPFxcXW9bp16/7lE7tRfYQdAACukOLiYg0YMMC6vnbtWvn6+jqxomsDt7EAAICpEXYAAICpEXYAAICpEXYAAICpEXYAAICpEXYAAICpEXYAAICpEXYAAICpEXYAAICpEXYAAICp8bgIAKbFc4gASIQdACbGc4gASE6+jbVkyRJFR0fLYrHIYrEoNjZWH374obX/7NmzSkxMVFBQkHx9fTVkyBDl5ubabCM7O1t9+/aVj4+PgoODNXXqVJWVlV3pQwEAAC7KqWGnQYMGevLJJ5WVlaUvv/xS3bt314ABA/Ttt99KkqZMmaL33ntPb7zxhjIyMnTs2DENHjzY+vry8nL17dtXpaWl+uKLL7Ry5UqlpaVpxowZzjokAADgYpx6G6tfv34263PnztWSJUu0detWNWjQQMuXL1d6erq6d+8uSVqxYoVatGihrVu3qlOnTvr444+1d+9ebdy4USEhIWrTpo3mzJmjadOmadasWfL09HTGYQEAABfiMp/GKi8v1+rVq1VcXKzY2FhlZWXp3Llzio+Pt45p3ry5IiMjlZmZKUnKzMxUq1atFBISYh2TkJCgwsJC69WhiykpKVFhYaHNAgAAzMnpYWfPnj3y9fWVl5eXHnjgAa1Zs0YtW7ZUTk6OPD09FRAQYDM+JCREOTk5kqScnByboHO+/3zfpaSmpsrf39+6REREOPagAACAy3B62GnWrJl27dqlbdu2afz48Ro1apT27t1bo/tMSUlRQUGBdTly5EiN7g8AADiP0z967unpqaioKElSu3bttGPHDj399NO66667VFpaqvz8fJurO7m5uQoNDZUkhYaGavv27TbbO/9prfNjLsbLy0teXl4OPhIAAOCKnH5l588qKipUUlKidu3aqXbt2tq0aZO1b9++fcrOzlZsbKwkKTY2Vnv27FFeXp51zIYNG2SxWNSyZcsrXjsAAHA9Tr2yk5KSot69eysyMlKnT59Wenq6Nm/erI8++kj+/v4aM2aMkpOTFRgYKIvFookTJyo2NladOnWSJPXs2VMtW7bUiBEjNG/ePOXk5OjRRx9VYmIiV24AAIAkJ4edvLw8jRw5UsePH5e/v7+io6P10Ucf6fbbb5ckLVy4UO7u7hoyZIhKSkqUkJCg559/3vp6Dw8PrVu3TuPHj1dsbKzq1q2rUaNGafbs2c46JAAA4GKcGnaWL19+2f46depo8eLFWrx48SXHNGzYUB988IGjSwMAACbhcnN2AAAAHImwAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATK2WswtwFWFhYcrNzVVRUZGzSwHgIGfOnJGvr691PTc3V4WFhU6sCNc6zknHKi4urtQ4ruwAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTc2rYSU1NVYcOHeTn56fg4GANHDhQ+/btsxkTFxcnNzc3m+WBBx6wGZOdna2+ffvKx8dHwcHBmjp1qsrKyq7koQAAABdVy5k7z8jIUGJiojp06KCysjJNnz5dPXv21N69e1W3bl3ruPvvv1+zZ8+2rvv4+Fh/Li8vV9++fRUaGqovvvhCx48f18iRI1W7dm098cQTV/R4AACA63Fq2Fm/fr3NelpamoKDg5WVlaVbbrnF2u7j46PQ0NCLbuPjjz/W3r17tXHjRoWEhKhNmzaaM2eOpk2bplmzZsnT07NGjwEAALg2l5qzU1BQIEkKDAy0aV+1apXq1aunm266SSkpKfrtt9+sfZmZmWrVqpVCQkKsbQkJCSosLNS333570f2UlJSosLDQZgEAAObk1Cs7f1RRUaHJkyerS5cuuummm6ztd999txo2bKjw8HDt3r1b06ZN0759+/T2229LknJycmyCjiTrek5OzkX3lZqaqscee8ymLSwszJGHAwAAXITLhJ3ExER988032rJli0372LFjrT+3atVKYWFh6tGjhw4ePKgbbrihSvtKSUlRcnKydb2wsFAdO3asWuEAAMCluUTYSUpK0rp16/Tpp5+qQYMGlx0bExMjSfrhhx90ww03KDQ0VNu3b7cZk5ubK0mXnOfj5eUlLy8vm7bjx48rJCREFoulqocBwMUUFRWpqKjIuh4SEiJfX18nVoRrHeekY1V2GopT5+wYhqGkpCStWbNGn3zyiRo3bvyXr9m1a5ek/7vtFBsbqz179igvL886ZsOGDbJYLGrZsmWN1A0AAK4eTr2yk5iYqPT0dK1du1Z+fn7WOTb+/v7y9vbWwYMHlZ6erj59+igoKEi7d+/WlClTdMsttyg6OlqS1LNnT7Vs2VIjRozQvHnzlJOTo0cffVSJiYkXXL0BAADXHqde2VmyZIkKCgoUFxensLAw6/Laa69Jkjw9PbVx40b17NlTzZs314MPPqghQ4bovffes27Dw8ND69atk4eHh2JjY3XPPfdo5MiRNt/LAwAArl1OvbJjGMZl+yMiIpSRkfGX22nYsKE++OADR5UFAABMxKW+ZwcAAMDRCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDU7Ao7ZWVlmj17to4ePVpT9QAAADiUXWGnVq1amj9/vsrKymqqHgAAAIey+zZW9+7dlZGRURO1AAAAOFwte1/Qu3dv/f3vf9eePXvUrl071a1b16a/f//+DisOQPW0m/qys0twKreyUvn/YT3uH6tl1PJ0Wj3OljV/pLNL4JzknLRxpc5Ju8POhAkTJEkLFiy4oM/NzU3l5eXVrwoAAMBB7A47FRUVNVEHAABAjeCj5wAAwNSqFHYyMjLUr18/RUVFKSoqSv3799dnn33m6NoAAACqze6w8+qrryo+Pl4+Pj6aNGmSJk2aJG9vb/Xo0UPp6ek1USMAAECV2T1nZ+7cuZo3b56mTJlibZs0aZIWLFigOXPm6O6773ZogQAAANVh95WdQ4cOqV+/fhe09+/fX4cPH3ZIUQAAAI5id9iJiIjQpk2bLmjfuHGjIiIiHFIUAACAo9h9G+vBBx/UpEmTtGvXLnXu3FmS9PnnnystLU1PP/20wwsEAACoDrvDzvjx4xUaGqqnnnpKr7/+uiSpRYsWeu211zRgwACHFwgAAFAddocdSRo0aJAGDRrk6FoAAAAczu45O02aNNHJkycvaM/Pz1eTJk3s2lZqaqo6dOggPz8/BQcHa+DAgdq3b5/NmLNnzyoxMVFBQUHy9fXVkCFDlJubazMmOztbffv2lY+Pj4KDgzV16lSezA4AACRVIez8+OOPF33+VUlJiX7++We7tpWRkaHExERt3bpVGzZs0Llz59SzZ08VFxdbx0yZMkXvvfee3njjDWVkZOjYsWMaPHiwtb+8vFx9+/ZVaWmpvvjiC61cuVJpaWmaMWOGvYcGAABMqNK3sd59913rzx999JH8/f/vua3l5eXatGmTGjVqZNfO169fb7Oelpam4OBgZWVl6ZZbblFBQYGWL1+u9PR0de/eXZK0YsUKtWjRQlu3blWnTp308ccfa+/evdq4caNCQkLUpk0bzZkzR9OmTdOsWbPk6XntPk0WAADYEXYGDhwo6fcnm48aNcqmr3bt2mrUqJGeeuqpahVTUFAgSQoMDJQkZWVl6dy5c4qPj7eOad68uSIjI5WZmalOnTopMzNTrVq1UkhIiHVMQkKCxo8fr2+//VZt27a9YD8lJSUqKSmxrhcWFlarbgAA4LoqHXbOP+28cePG2rFjh+rVq+fQQioqKjR58mR16dJFN910kyQpJydHnp6eCggIsBkbEhKinJwc65g/Bp3z/ef7LiY1NVWPPfaYQ+sHAACuye45O4cPH3Z40JGkxMREffPNN1q9erXDt/1nKSkpKigosC5Hjhyp8X0CAADnsDvsTJo0Sc8888wF7c8995wmT55cpSKSkpK0bt06/e///q8aNGhgbQ8NDVVpaany8/Ntxufm5io0NNQ65s+fzjq/fn7Mn3l5eclisdgsAADAnOwOO2+99Za6dOlyQXvnzp315ptv2rUtwzCUlJSkNWvW6JNPPlHjxo1t+tu1a6fatWvbPJ5i3759ys7OVmxsrCQpNjZWe/bsUV5ennXMhg0bZLFY1LJlS7vqAQAA5mP3lwqePHnS5pNY51ksFv3yyy92bSsxMVHp6elau3at/Pz8rHNs/P395e3tLX9/f40ZM0bJyckKDAyUxWLRxIkTFRsbq06dOkmSevbsqZYtW2rEiBGaN2+ecnJy9OijjyoxMVFeXl72Hh4AADAZu6/sREVFXfCRcUn68MMP7f5SwSVLlqigoEBxcXEKCwuzLq+99pp1zMKFC3XHHXdoyJAhuuWWWxQaGqq3337b2u/h4aF169bJw8NDsbGxuueeezRy5EjNnj3b3kMDAAAmZPeVneTkZCUlJenEiRPW777ZtGmTnnrqKS1atMiubRmG8Zdj6tSpo8WLF2vx4sWXHNOwYUN98MEHdu0bAABcG+wOO/fdd59KSko0d+5czZkzR5LUqFEjLVmyRCNHjnR4gQAAANVRpQeBjh8/XuPHj9eJEyfk7e0tX19fR9cFAADgEFUKO+fVr1/fUXUAAADUiCqFnTfffFOvv/66srOzVVpaatO3c+dOhxQGAADgCHZ/GuuZZ57R6NGjFRISoq+++kodO3ZUUFCQDh06pN69e9dEjbgKGIahoqIi61KZyecAAFwJdl/Zef7557Vs2TINGzZMaWlpevjhh9WkSRPNmDFDp06dqokacRUoLi7WgAEDrOtr165lLhcAwCXYfWUnOztbnTt3liR5e3vr9OnTkqQRI0bof/7nfxxbHQAAQDXZHXZCQ0OtV3AiIyO1detWSb8/IJRbFwAAwNXYHXa6d++ud999V5I0evRoTZkyRbfffrvuuusuDRo0yOEFAgAAVIfdc3aWLVumiooKSb8/2yooKEhffPGF+vfvr3Hjxjm8QAAAgOqodNiJjIzUV199paCgILm7u+u5557TyJEjNXToUA0dOrQmawQAAKiySt/GOnr0qMrLy63r06dPt/sp5wAAAFea3XN2zmMyMgAAuBpUOewAAABcDeyaoPzSSy9ZvyiurKxMaWlpqlevns2YSZMmOa46AACAarJrgvKLL75oXQ8NDdUrr7xiM8bNzY2wAwAAXEqlw86PP/5Yg2UAAADUDObsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAU7M77OzcuVN79uyxrq9du1YDBw7U9OnTVVpa6tDiAAAAqsvusDNu3Djt379fknTo0CENHTpUPj4+euONN/Twww87vEAAAIDqsDvs7N+/X23atJEkvfHGG7rllluUnp6utLQ0vfXWW46uDwAAoFrsDjuGYaiiokKStHHjRvXp00eSFBERwVPQAQCAy7E77LRv316PP/64XnnlFWVkZKhv376SpMOHDyskJMThBQIAAFSH3WFn0aJF2rlzp5KSkvTII48oKipKkvTmm2+qc+fODi8QAACgOux66rkkRUdH23wa67z58+fLw8PDIUUBAAA4SpW+Zyc/P18vvfSSUlJSdOrUKUnS3r17lZeX59DiAAAAqsvuKzu7d+9Wjx49FBAQoB9//FH333+/AgMD9fbbbys7O1svv/xyTdQJAABQJXZf2UlOTtbo0aN14MAB1alTx9rep08fffrppw4tDgAAoLrsDjs7duzQuHHjLmi//vrrlZOT45CiAAAAHMXusOPl5aXCwsIL2vfv36/69es7pCgAAABHsTvs9O/fX7Nnz9a5c+ckSW5ubsrOzta0adM0ZMgQhxcIAABQHXaHnaeeekpFRUUKDg7WmTNndOuttyoqKkp+fn6aO3duTdQIAABQZXZ/Gsvf318bNmzQli1btHv3bhUVFenmm29WfHx8TdQHAABQLXaHnUOHDqlJkybq2rWrunbtWhM1AQAAOIzdt7GioqJ022236dVXX9XZs2droiYAAACHsTvs7Ny5U9HR0UpOTlZoaKjGjRunbdu21URtAAAA1WZ32GnTpo2efvppHTt2TP/+9791/PhxdevWTTfddJMWLFigEydO1ESdAAAAVVKlZ2NJUq1atTR48GC98cYb+uc//6kffvhBDz30kCIiIjRy5EgdP37ckXUCAABUSZXDzpdffqkJEyYoLCxMCxYs0EMPPaSDBw9qw4YNOnbsmAYMGODIOgEAAKrE7k9jLViwQCtWrNC+ffvUp08fvfzyy+rTp4/c3X/PTY0bN1ZaWpoaNWrk6FoBAADsZnfYWbJkie677z7de++9CgsLu+iY4OBgLV++vNrFAQAAVJfdt7EOHDiglJSUSwYdSfL09NSoUaP+cluffvqp+vXrp/DwcLm5uemdd96x6b/33nvl5uZms/Tq1ctmzKlTpzR8+HBZLBYFBARozJgxKioqsvewAACASdl9Zee83377TdnZ2SotLbVpj46OrvQ2iouL1bp1a913330aPHjwRcf06tVLK1assK57eXnZ9A8fPlzHjx/Xhg0bdO7cOY0ePVpjx45Venq6HUcDAADMyu6wc+LECd17771av379RfvLy8srva3evXurd+/elx3j5eWl0NDQi/Z99913Wr9+vXbs2KH27dtLkp599ln16dNH//rXvxQeHl7pWgAAgDnZfRtr8uTJKigo0LZt2+Tt7a3169dr5cqVuvHGG/Xuu+86vMDNmzcrODhYzZo10/jx43Xy5ElrX2ZmpgICAqxBR5Li4+Pl7u5+2S86LCkpUWFhoc0CAADMye4rO5988onWrl2r9u3by93dXQ0bNtTtt98ui8Wi1NRU9e3b12HF9erVS4MHD1bjxo118OBBTZ8+Xb1791ZmZqY8PDyUk5Oj4OBg2wOqVUuBgYHKycm55HZTU1P12GOPOaxOAADguuwOO8XFxdaAcd111+nEiRNq2rSpWrVqpZ07dzq0uKFDh1p/btWqlaKjo3XDDTdo8+bN6tGjR5W3m5KSouTkZOt6YWGhIiIiqlUrAABwTXbfxmrWrJn27dsnSWrdurVeeOEF/fzzz1q6dOllP6HlCE2aNFG9evX0ww8/SJJCQ0OVl5dnM6asrEynTp265Dwf6fd5QBaLxWYBAADmZPeVnb/97W/WR0HMnDlTvXr10qpVq+Tp6am0tDRH12fj6NGjOnnypDVUxcbGKj8/X1lZWWrXrp2k32+zVVRUKCYmpkZrAQAAVwe7w84999xj/bldu3b66aef9P333ysyMlL16tWza1tFRUXWqzSSdPjwYe3atUuBgYEKDAzUY489piFDhig0NFQHDx7Uww8/rKioKCUkJEiSWrRooV69eun+++/X0qVLde7cOSUlJWno0KF8EgsAAEiqxrOxzvPx8dHNN99sd9CRfn++Vtu2bdW2bVtJUnJystq2basZM2bIw8NDu3fvVv/+/dW0aVONGTNG7dq102effWbzXTurVq1S8+bN1aNHD/Xp00ddu3bVsmXLqntYAADAJOy6snPgwAHt3r1bN998sxo3bqz3339f//znP3XmzBkNHDhQ06dPl5ubW6W3FxcXJ8MwLtn/0Ucf/eU2AgMD+QJBAABwSZUOO2vWrNF//dd/yd3dXW5ublq2bJnGjRunuLg4WSwWzZo1S7Vq1dK0adNqsl4AAAC7VPo21ty5c/Xwww/r7NmzWrJkiR544AGlpqbqww8/1Lp167R48eIan6AMAABgr0qHnX379um+++6Tm5ubRo0apdLSUsXHx1v7e/bsqZ9++qlGigQAAKiqSoed4uJi+fn5/f4id3d5e3vLx8fH2u/t7a2SkhLHVwgAAFANlQ47bm5uNpOP/7wOAADgiio9QdkwDDVt2tQacIqKitS2bVu5u7tb+wEAAFxNpcPOihUrarIOAACAGlHpsDNq1KiarAMAAKBGVPsblAEAAFwZYQcAAJgaYQcAAJgaYQcAAJiaXQ8CxaW1m/qys0twKreyUvn/YT3uH6tl1PJ0Wj3OljV/pLNLAAD8f3aHneTk5Iu2u7m5qU6dOoqKitKAAQMUGBhY7eIAAACqy+6w89VXX2nnzp0qLy9Xs2bNJEn79++Xh4eHmjdvrueff14PPvigtmzZopYtWzq8YAAArlaGR20VRA+zWUfNs3vOzoABAxQfH69jx44pKytLWVlZOnr0qG6//XYNGzZMP//8s2655RZNmTKlJuoFAODq5eYmo5andRGPXboi7A478+fP15w5c2SxWKxt/v7+mjVrlubNmycfHx/NmDFDWVlZDi0UAACgKuy+jVVQUKC8vLwLblGdOHFChYWFkqSAgACVlpY6pkIAqCJuGQCQqngb67777tOaNWt09OhRHT16VGvWrNGYMWM0cOBASdL27dvVtGlTR9cKAPbhlgEAVeHKzgsvvKApU6Zo6NChKisr+30jtWpp1KhRWrhwoSSpefPmeumllxxbKQAAQBXYHXZ8fX314osvauHChTp06JAkqUmTJvL19bWOadOmjcMKBAAAqI4qf6mgr6+v9bt0/hh0AAAAXIndc3YqKio0e/Zs+fv7q2HDhmrYsKECAgI0Z84cVVRU1ESNAAAAVWb3lZ1HHnlEy5cv15NPPqkuXbpIkrZs2aJZs2bp7Nmzmjt3rsOLBAAAqCq7w87KlSv10ksvqX///ta26OhoXX/99ZowYQJhBwAAuBS7b2OdOnVKzZs3v6C9efPmOnXqlEOKAgAAcBS7w07r1q313HPPXdD+3HPPqXXr1g4pCgAAwFHsvo01b9489e3bVxs3blRsbKwkKTMzU0eOHNEHH3zg8AIBAACqw+4rO7feeqv279+vQYMGKT8/X/n5+Ro8eLD27dunbt261USNAAAAVVal79kJDw+/YCLy0aNHNXbsWC1btswhhQEAADiC3Vd2LuXkyZNavny5ozYHAADgEA4LOwAAAK6IsAMAAEyNsAMAAEyt0hOUBw8efNn+/Pz86tYCAADgcJUOO/7+/n/ZP3LkyGoXBAAA4EiVDjsrVqyoyToAAABqBHN2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqRF2AACAqTk17Hz66afq16+fwsPD5ebmpnfeecem3zAMzZgxQ2FhYfL29lZ8fLwOHDhgM+bUqVMaPny4LBaLAgICNGbMGBUVFV3BowAAAK7MqWGnuLhYrVu31uLFiy/aP2/ePD3zzDNaunSptm3bprp16yohIUFnz561jhk+fLi+/fZbbdiwQevWrdOnn36qsWPHXqlDAAAALq7Sj4uoCb1791bv3r0v2mcYhhYtWqRHH31UAwYMkCS9/PLLCgkJ0TvvvKOhQ4fqu+++0/r167Vjxw61b99ekvTss8+qT58++te//qXw8PArdiwAAMA1ueycncOHDysnJ0fx8fHWNn9/f8XExCgzM1OSlJmZqYCAAGvQkaT4+Hi5u7tr27Ztl9x2SUmJCgsLbRYAAGBOLht2cnJyJEkhISE27SEhIda+nJwcBQcH2/TXqlVLgYGB1jEXk5qaKn9/f+sSERHh4OoBAICrcNmwU5NSUlJUUFBgXY4cOeLskgAAQA1x2bATGhoqScrNzbVpz83NtfaFhoYqLy/Ppr+srEynTp2yjrkYLy8vWSwWmwUAAJiTy4adxo0bKzQ0VJs2bbK2FRYWatu2bYqNjZUkxcbGKj8/X1lZWdYxn3zyiSoqKhQTE3PFawYAAK7HqZ/GKioq0g8//GBdP3z4sHbt2qXAwEBFRkZq8uTJevzxx3XjjTeqcePG+sc//qHw8HANHDhQktSiRQv16tVL999/v5YuXapz584pKSlJQ4cO5ZNYAABAkpPDzpdffqnbbrvNup6cnCxJGjVqlNLS0vTwww+ruLhYY8eOVX5+vrp27ar169erTp061tesWrVKSUlJ6tGjh9zd3TVkyBA988wzV/xYAACAa3Jq2ImLi5NhGJfsd3Nz0+zZszV79uxLjgkMDFR6enpNlAcAAEzAZefsAAAAOAJhBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmJpTv0EZ5mF41FZB9DCbdQAAXAFhB47h5iajlqezqwAA4ALcxgIAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKbm0mFn1qxZcnNzs1maN29u7T979qwSExMVFBQkX19fDRkyRLm5uU6sGAAAuBqXDjuS9B//8R86fvy4ddmyZYu1b8qUKXrvvff0xhtvKCMjQ8eOHdPgwYOdWC0AAHA1tZxdwF+pVauWQkNDL2gvKCjQ8uXLlZ6eru7du0uSVqxYoRYtWmjr1q3q1KnTlS4VAAC4IJe/snPgwAGFh4erSZMmGj58uLKzsyVJWVlZOnfunOLj461jmzdvrsjISGVmZl52myUlJSosLLRZAACAObl02ImJiVFaWprWr1+vJUuW6PDhw+rWrZtOnz6tnJwceXp6KiAgwOY1ISEhysnJuex2U1NT5e/vb10iIiJq8CgAAIAzufRtrN69e1t/jo6OVkxMjBo2bKjXX39d3t7eVd5uSkqKkpOTreuFhYUEHgAATMqlr+z8WUBAgJo2baoffvhBoaGhKi0tVX5+vs2Y3Nzci87x+SMvLy9ZLBabBQAAmNNVFXaKiop08OBBhYWFqV27dqpdu7Y2bdpk7d+3b5+ys7MVGxvrxCoBAIArcenbWA899JD69eunhg0b6tixY5o5c6Y8PDw0bNgw+fv7a8yYMUpOTlZgYKAsFosmTpyo2NhYPokFAACsXDrsHD16VMOGDdPJkydVv359de3aVVu3blX9+vUlSQsXLpS7u7uGDBmikpISJSQk6Pnnn3dy1QAAwJW4dNhZvXr1Zfvr1KmjxYsXa/HixVeoIgAAcLW5qubsAAAA2IuwAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATM00YWfx4sVq1KiR6tSpo5iYGG3fvt3ZJQEAABdgirDz2muvKTk5WTNnztTOnTvVunVrJSQkKC8vz9mlAQAAJzNF2FmwYIHuv/9+jR49Wi1bttTSpUvl4+Ojf//7384uDQAAOFktZxdQXaWlpcrKylJKSoq1zd3dXfHx8crMzLzoa0pKSlRSUmJdLygokCQVFhZWuY7ykjNVfi3MpzrnkiNxXuKPXOG85JzEH1X3nDz/esMwLjvuqg87v/zyi8rLyxUSEmLTHhISou+///6ir0lNTdVjjz12QXtERESN1Ihrj/+zDzi7BOACnJdwNY46J0+fPi1/f/9L9l/1YacqUlJSlJycbF2vqKjQqVOnFBQUJDc3NydWdnUrLCxURESEjhw5IovF4uxyAEmcl3A9nJOOYxiGTp8+rfDw8MuOu+rDTr169eTh4aHc3Fyb9tzcXIWGhl70NV5eXvLy8rJpCwgIqKkSrzkWi4U3MFwO5yVcDeekY1zuis55V/0EZU9PT7Vr106bNm2ytlVUVGjTpk2KjY11YmUAAMAVXPVXdiQpOTlZo0aNUvv27dWxY0ctWrRIxcXFGj16tLNLAwAATmaKsHPXXXfpxIkTmjFjhnJyctSmTRutX7/+gknLqFleXl6aOXPmBbcIAWfivISr4Zy88tyMv/q8FgAAwFXsqp+zAwAAcDmEHQAAYGqEHQAAYGqEHThdXFycJk+e7OwyAMBu/P26OhB2AACAqRF2AABwkrS0NMXFxTm7DNMj7FyD4uLiNHHiRE2ePFnXXXedQkJC9OKLL1q/iNHPz09RUVH68MMPJf3+Zvzz4zTeeecdm+eIff3117rtttvk5+cni8Widu3a6csvv7T2f/7554qLi5OPj4+uu+46JSQk6Ndff71ofSUlJXrooYd0/fXXq27duoqJidHmzZsd/nuA62rUqJEWLVpk09amTRvNmjVLhmFo1qxZioyMlJeXl8LDwzVp0iTruFdeeUXt27eXn5+fQkNDdffddysvL89mW++++65uvPFG1alTR7fddptWrlwpNzc35efnW8ds2bJF3bp1k7e3tyIiIjRp0iQVFxfX5GGjEq61v1+Xey9I4v1QSYSda9TKlStVr149bd++XRMnTtT48eN15513qnPnztq5c6d69uypESNG6LfffqvU9oYPH64GDRpox44dysrK0t///nfVrl1bkrRr1y716NFDLVu2VGZmprZs2aJ+/fqpvLz8ottKSkpSZmamVq9erd27d+vOO+9Ur169dODAAYcdP65eb731lhYuXKgXXnhBBw4c0DvvvKNWrVpZ+8+dO6c5c+bo66+/1jvvvKMff/xR9957r7X/8OHD+s///E8NHDhQX3/9tcaNG6dHHnnEZh8HDx5Ur169NGTIEO3evVuvvfaatmzZoqSkpCt1mLgM/n79H94PlWTgmnPrrbcaXbt2ta6XlZUZdevWNUaMGGFtO378uCHJyMzMNFasWGH4+/vbbGPNmjXGH08fPz8/Iy0t7aL7GzZsmNGlS5fL1vO3v/3NMAzD+OmnnwwPDw/j559/thnTo0cPIyUlpbKHiKtcw4YNjYULF9q0tW7d2pg5c6bx1FNPGU2bNjVKS0srta0dO3YYkozTp08bhmEY06ZNM2666SabMY888oghyfj1118NwzCMMWPGGGPHjrUZ89lnnxnu7u7GmTNnqnZQcAiz/f1asWKFceutt15y+5d7LxiGwfuhkriyc42Kjo62/uzh4aGgoCCb/w2cf9TGny93XkpycrL++7//W/Hx8XryySd18OBBa9/5/xlVxp49e1ReXq6mTZvK19fXumRkZNhsE9euO++8U2fOnFGTJk10//33a82aNSorK7P2Z2VlqV+/foqMjJSfn59uvfVWSVJ2drYkad++ferQoYPNNjt27Giz/vXXXystLc3mHExISFBFRYUOHz5cw0eIv3I1//3Kzs626XvggQf02Wef2bQ98cQTldqfxPuhskzxbCzY7/wl2vPc3Nxs2s7fz66oqJC7u7uMPz1V5Ny5czbrs2bN0t133633339fH374oWbOnKnVq1dr0KBB8vb2rnRdRUVF8vDwUFZWljw8PGz6fH19K70dXN0ud85FRERo37592rhxozZs2KAJEyZo/vz5ysjIUGlpqRISEpSQkKBVq1apfv36ys7OVkJCgkpLSyu9/6KiIo0bN85m7sN5kZGR1Ts4VNvV/PcrPDxcu3btsra//fbbeuutt7Rq1SprW2BgoPXnv6qf90PlEHbwl+rXr6/Tp0+ruLhYdevWlSSbN+t5TZs2VdOmTTVlyhQNGzZMK1as0KBBgxQdHa1Nmzbpscce+8t9tW3bVuXl5crLy1O3bt0cfSi4StSvX1/Hjx+3rhcWFtr8D9Lb21v9+vVTv379lJiYqObNm2vPnj0yDEMnT57Uk08+qYiICEmymWgqSc2aNdMHH3xg07Zjxw6b9Ztvvll79+5VVFSUow8NV5ir/f2qVauWzXkVHBwsb2/vS55rf/VekHg/VAa3sfCXYmJi5OPjo+nTp+vgwYNKT09XWlqatf/MmTNKSkrS5s2b9dNPP+nzzz/Xjh071KJFC0lSSkqKduzYoQkTJmj37t36/vvvtWTJEv3yyy8X7Ktp06YaPny4Ro4cqbfffluHDx/W9u3blZqaqvfff/9KHTKcrHv37nrllVf02Wefac+ePRo1apT1f8ppaWlavny5vvnmGx06dEivvvqqvL291bBhQ0VGRsrT01PPPvusDh06pHfffVdz5syx2fa4ceP0/fffa9q0adq/f79ef/116/l8/orAtGnT9MUXXygpKUm7du3SgQMHtHbt2qtrQiYkXf1/vy73XpB4P1SaMycMwTn+OKHuvItNgpNkrFmzxjCM3yf0RUVFGd7e3sYdd9xhLFu2zDrBr6SkxBg6dKgRERFheHp6GuHh4UZSUpLNxLXNmzcbnTt3Nry8vIyAgAAjISHBOvntz/WUlpYaM2bMMBo1amTUrl3bCAsLMwYNGmTs3r3b0b8KuKiCggLjrrvuMiwWixEREWGkpaVZJ2WuWbPGiImJMSwWi1G3bl2jU6dOxsaNG62vTU9PNxo1amR4eXkZsbGxxrvvvmtIMr766ivrmLVr1xpRUVGGl5eXERcXZyxZssSQZHPObt++3bj99tsNX19fo27dukZ0dLQxd+7cK/lrwEWY7e/XX01Qvtx74fyx8X74a26G8aebgQBwjZk7d66WLl2qI0eOOLsUwOnM+H5gzg6Aa87zzz+vDh06KCgoSJ9//rnmz59/dV2SBxzoWng/EHYAXHMOHDigxx9/XKdOnVJkZKQefPBBpaSkOLsswCmuhfcDt7EAAICp8WksAABgaoQdAABgaoQdAABgaoQdAABgaoQdAABgaoQdANekuLg4TZ48udLj09LSFBAQUGP1AKg5hB0AAGBqhB0AAGBqhB0ALiUuLk4TJ07U5MmTdd111ykkJEQvvviiiouLNXr0aPn5+SkqKkoffvih9TUZGRnq2LGjvLy8FBYWpr///e8qKyuz9hcXF2vkyJHy9fVVWFiYnnrqqQv2W1JSooceekjXX3+96tatq5iYGG3evPlKHDKAGkbYAeByVq5cqXr16mn79u2aOHGixo8frzvvvFOdO3fWzp071bNnT40YMUK//fabfv75Z/Xp00cdOnTQ119/rSVLlmj58uV6/PHHrdubOnWqMjIytHbtWn388cfavHmzdu7cabPPpKQkZWZmavXq1dq9e7fuvPNO9erVSwcOHLjShw/AwXhcBACXEhcXp/Lycn322WeSpPLycvn7+2vw4MF6+eWXJUk5OTkKCwtTZmam3nvvPb311lv67rvv5ObmJun3BxtOmzZNBQUF+u233xQUFKRXX31Vd955pyTp1KlTatCggcaOHatFixYpOztbTZo0UXZ2tsLDw621xMfHq2PHjnriiSeUlpamyZMnKz8//8r+QgBUGw8CBeByoqOjrT97eHgoKChIrVq1sraFhIRIkvLy8vTdd98pNjbWGnQkqUuXLioqKtLRo0f166+/qrS0VDExMdb+wMBANWvWzLq+Z88elZeXq2nTpjZ1lJSUKCgoyOHHB+DKIuwAcDm1a9e2WXdzc7NpOx9sKioqHLK/oqIieXh4KCsrSx4eHjZ9vr6+DtkHAOch7AC4qrVo0UJvvfWWDMOwhqDPP/9cfn5+atCggQIDA1W7dm1t27ZNkZGRkqRff/1V+/fv16233ipJatu2rcrLy5WXl6du3bo57VgA1AwmKAO4qk2YMEFHjhzRxIkT9f3332vt2rWaOXOmkpOT5e7uLl9fX40ZM0ZTp07VJ598om+++Ub33nuv3N3/789f06ZNNXz4cI0cOVJvv/22Dh8+rO3btys1NVXvv/++E48OgCNwZQfAVe3666/XBx98oKlTp6p169YKDAzUmDFj9Oijj1rHzJ8/X0VFRerXr5/8/Pz04IMPqqCgwGY7K1as0OOPP64HH3xQP//8s+rVq6dOnTrpjjvuuNKHBMDB+DQWAAAwNW5jAQAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAU/t/VEwvvCHjABEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = pcm.model_plot(T_cv.likelihood,\n",
    "                        null_model = 'null',\n",
    "                        noise_ceiling= 'ceil',\n",
    "                        upper_ceiling = T_gr.likelihood['ceil'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the likelihood for individual, group, and crossvalidated group fits for the fixed models (null, muscle + usage) are all identical, because these models do not have common group parameters - in all cases we are fitting an individual scale and noise parameter. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the model fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, it is very useful to visualize the model prediction in comparision to the fitted data. The model parameters are stored in the return argument `theta`. We can pass these to the `Model.predict()` function to get the predicted second moment matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x336ea9c10>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGdCAYAAAAv9mXmAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAEeFJREFUeJzt3V9o1Yfdx/Fv1OWkrUlo7LT1SVwLHR3OxVGtXVZYu5q1SJGW3eyisODGYCMOJTcjN5PBRoQHRssqTvavDCbKBraj0Dpnp2EPdcZIwHZU6Cgjw2rW5yLR0B5tznkuxvLMtXU5ab7+zomvF5yLc/id/j782ubtOSeJTdVqtRoAsMCWFD0AgMVJYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFsut9wkqlEufOnYvW1tZoamq63qcH4COoVqtx8eLFWL16dSxZcu3XKNc9MOfOnYuurq7rfVoAFtD4+Hh0dnZe85jrHpjW1taIiPjr6Tujbbl36K7lwe9/vegJDaHjlyeLnsAiUv3cZ4qeUNfee68c/3Pqv2e/ll/LdQ/MP98Wa1u+JNpaBeZalja3FD2hISxr+ljRE1hEqsv8fzcXc/mIw1d4AFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFvAKzZ8+euPPOO6OlpSXuv//+OHny5ELvAqDB1RyYgwcPxsDAQOzatStOnz4d69evj0cffTQmJiYy9gHQoGoOzA9/+MP4xje+Edu2bYu1a9fGj3/847j55pvj5z//ecY+ABpUTYG5fPlyjI6ORm9v7///A5Ysid7e3njllVcWfBwAjWtZLQe//fbbMTMzE6tWrbrq8VWrVsXrr7/+gc8pl8tRLpdn709NTc1jJgCNJv27yIaGhqK9vX321tXVlX1KAOpATYG57bbbYunSpXHhwoWrHr9w4ULcfvvtH/icwcHBmJycnL2Nj4/Pfy0ADaOmwDQ3N8eGDRvi6NGjs49VKpU4evRo9PT0fOBzSqVStLW1XXUDYPGr6TOYiIiBgYHo6+uLjRs3xqZNm+Kpp56K6enp2LZtW8Y+ABpUzYH5yle+En//+9/ju9/9bpw/fz4++9nPxksvvfS+D/4BuLHVHJiIiO3bt8f27dsXegsAi4jfRQZACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIsK+rED37/67G0uaWo0zeE6hP/W/SEhjDR8vmiJzSEO34/UfSEhtD0l7eKnlDXmiqX53ysVzAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASFFzYIaHh2Pr1q2xevXqaGpqiueeey5hFgCNrubATE9Px/r162PPnj0ZewBYJJbV+oQtW7bEli1bMrYAsIj4DAaAFDW/gqlVuVyOcrk8e39qair7lADUgfRXMENDQ9He3j576+rqyj4lAHUgPTCDg4MxOTk5exsfH88+JQB1IP0tslKpFKVSKfs0ANSZmgNz6dKleOONN2bvv/nmmzE2NhYdHR2xZs2aBR0HQOOqOTCnTp2KL37xi7P3BwYGIiKir68vnn322QUbBkBjqzkwDz30UFSr1YwtACwifg4GgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkWFbUiTt+eTKWNX2sqNM3hImWzxc9oSFc6qoWPaEh/PXLK4ue0BA6zq4oekJde+/KuxG/nduxXsEAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIEVNgRkaGor77rsvWltbY+XKlfHEE0/E2bNns7YB0MBqCszx48ejv78/Tpw4EUeOHIkrV67EI488EtPT01n7AGhQy2o5+KWXXrrq/rPPPhsrV66M0dHR+MIXvrCgwwBobDUF5t9NTk5GRERHR8eHHlMul6NcLs/en5qa+iinBKBBzPtD/kqlEjt37owHHngg1q1b96HHDQ0NRXt7++ytq6trvqcEoIHMOzD9/f3x6quvxoEDB6553ODgYExOTs7exsfH53tKABrIvN4i2759e7zwwgsxPDwcnZ2d1zy2VCpFqVSa1zgAGldNgalWq/Htb387Dh06FMeOHYu77roraxcADa6mwPT398f+/fvj+eefj9bW1jh//nxERLS3t8dNN92UMhCAxlTTZzB79+6NycnJeOihh+KOO+6YvR08eDBrHwANqua3yABgLvwuMgBSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkGJZ0QP4cHf8fqLoCQ3hr19eWfSEhvDOf80UPaEhTNzsz93XUnl3ScRv53asKwlACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFDUFZu/evdHd3R1tbW3R1tYWPT098eKLL2ZtA6CB1RSYzs7O2L17d4yOjsapU6fi4Ycfjscffzxee+21rH0ANKhltRy8devWq+7/4Ac/iL1798aJEyfi05/+9IIOA6Cx1RSYfzUzMxO//vWvY3p6Onp6ej70uHK5HOVyefb+1NTUfE8JQAOp+UP+M2fOxPLly6NUKsU3v/nNOHToUKxdu/ZDjx8aGor29vbZW1dX10caDEBjqDkw99xzT4yNjcWf/vSn+Na3vhV9fX3x5z//+UOPHxwcjMnJydnb+Pj4RxoMQGOo+S2y5ubmuPvuuyMiYsOGDTEyMhJPP/107Nu37wOPL5VKUSqVPtpKABrOR/45mEqlctVnLAAQUeMrmMHBwdiyZUusWbMmLl68GPv3749jx47F4cOHs/YB0KBqCszExER89atfjbfeeiva29uju7s7Dh8+HF/60pey9gHQoGoKzM9+9rOsHQAsMn4XGQApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASLGsqBNXP/eZqC5rKer0DaHpL28VPaEhdJxdUfSEhjBxsz9PzsXlzstFT6hrlXfmfn38FwdACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFB8pMLt3746mpqbYuXPnAs0BYLGYd2BGRkZi37590d3dvZB7AFgk5hWYS5cuxZNPPhk/+clP4tZbb13oTQAsAvMKTH9/fzz22GPR29v7H48tl8sxNTV11Q2AxW9ZrU84cOBAnD59OkZGRuZ0/NDQUHzve9+reRgAja2mVzDj4+OxY8eO+NWvfhUtLS1zes7g4GBMTk7O3sbHx+c1FIDGUtMrmNHR0ZiYmIh777139rGZmZkYHh6OZ555JsrlcixduvSq55RKpSiVSguzFoCGUVNgNm/eHGfOnLnqsW3btsWnPvWp+M53vvO+uABw46opMK2trbFu3bqrHrvllltixYoV73scgBubn+QHIEXN30X2744dO7YAMwBYbLyCASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSLLveJ6xWqxER8d575et96obTVLlc9ISG8N6Vd4ue0BAq7/rz5FxU3vH/3bVU3vnH1+5/fi2/lqbqXI5aQH/729+iq6vrep4SgAU2Pj4enZ2d1zzmugemUqnEuXPnorW1NZqamq7nqT/U1NRUdHV1xfj4eLS1tRU9py65RnPjOs2N6zQ39XidqtVqXLx4MVavXh1Lllz7VfF1f4tsyZIl/7F6RWlra6ubf4n1yjWaG9dpblynuam369Te3j6n47wpC0AKgQEghcBERKlUil27dkWpVCp6St1yjebGdZob12luGv06XfcP+QG4MXgFA0AKgQEghcAAkEJgAEhxwwdmz549ceedd0ZLS0vcf//9cfLkyaIn1Z3h4eHYunVrrF69OpqamuK5554relLdGRoaivvuuy9aW1tj5cqV8cQTT8TZs2eLnlV39u7dG93d3bM/ONjT0xMvvvhi0bPq3u7du6OpqSl27txZ9JSa3NCBOXjwYAwMDMSuXbvi9OnTsX79+nj00UdjYmKi6Gl1ZXp6OtavXx979uwpekrdOn78ePT398eJEyfiyJEjceXKlXjkkUdienq66Gl1pbOzM3bv3h2jo6Nx6tSpePjhh+Pxxx+P1157rehpdWtkZCT27dsX3d3dRU+pXfUGtmnTpmp/f//s/ZmZmerq1aurQ0NDBa6qbxFRPXToUNEz6t7ExEQ1IqrHjx8vekrdu/XWW6s//elPi55Rly5evFj95Cc/WT1y5Ej1wQcfrO7YsaPoSTW5YV/BXL58OUZHR6O3t3f2sSVLlkRvb2+88sorBS5jMZicnIyIiI6OjoKX1K+ZmZk4cOBATE9PR09PT9Fz6lJ/f3889thjV32daiTX/Zdd1ou33347ZmZmYtWqVVc9vmrVqnj99dcLWsViUKlUYufOnfHAAw/EunXrip5Td86cORM9PT3x7rvvxvLly+PQoUOxdu3aomfVnQMHDsTp06djZGSk6CnzdsMGBrL09/fHq6++Gn/84x+LnlKX7rnnnhgbG4vJycn4zW9+E319fXH8+HGR+Rfj4+OxY8eOOHLkSLS0tBQ9Z95u2MDcdtttsXTp0rhw4cJVj1+4cCFuv/32glbR6LZv3x4vvPBCDA8P1+1fS1G05ubmuPvuuyMiYsOGDTEyMhJPP/107Nu3r+Bl9WN0dDQmJibi3nvvnX1sZmYmhoeH45lnnolyuRxLly4tcOHc3LCfwTQ3N8eGDRvi6NGjs49VKpU4evSo94OpWbVaje3bt8ehQ4fi5ZdfjrvuuqvoSQ2jUqlEueyvUP9XmzdvjjNnzsTY2NjsbePGjfHkk0/G2NhYQ8Ql4gZ+BRMRMTAwEH19fbFx48bYtGlTPPXUUzE9PR3btm0relpduXTpUrzxxhuz9998880YGxuLjo6OWLNmTYHL6kd/f3/s378/nn/++WhtbY3z589HxD/+Yqabbrqp4HX1Y3BwMLZs2RJr1qyJixcvxv79++PYsWNx+PDhoqfVldbW1vd9fnfLLbfEihUrGutzvaK/ja1oP/rRj6pr1qypNjc3Vzdt2lQ9ceJE0ZPqzh/+8IdqRLzv1tfXV/S0uvFB1yciqr/4xS+KnlZXvva1r1U/8YlPVJubm6sf//jHq5s3b67+7ne/K3pWQ2jEb1P26/oBSHHDfgYDQC6BASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEjxf4j62wovVe3cAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "G,_ = M[4].predict(theta_gr[4][:M[4].n_param])\n",
    "plt.imshow(G)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}