Simple_Model1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import cPickle as pickle\n",
    "import plotly.plotly as py\n",
    "from plotly.graph_objs import *\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "#import scipy.optimize\n",
    "import statsmodels.api as sm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from statsmodels.regression.linear_model import OLSResults\n",
    "from statsmodels.regression.linear_model import OLS as ols\n",
    "import scipy.stats as stat\n",
    "from numpy import linalg as LA\n",
    "from scipy.misc import comb\n",
    "from pylab import *\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import basinhopping\n",
    "\n",
    "pd.set_option('display.width', 200)\n",
    "pd.set_option('display.max_rows', 30)\n",
    "pd.set_option('display.notebook_repr_html',True)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Size</th>\n",
       "      <th>Image</th>\n",
       "      <th>Ori</th>\n",
       "      <th>PosY</th>\n",
       "      <th>trial</th>\n",
       "      <th>time</th>\n",
       "      <th>terraincode</th>\n",
       "      <th>pausetime</th>\n",
       "      <th>choice</th>\n",
       "      <th>selectiontime</th>\n",
       "      <th>istarget</th>\n",
       "      <th>speed</th>\n",
       "      <th>lapdistance</th>\n",
       "      <th>timebeforezone</th>\n",
       "      <th>slowing_f</th>\n",
       "      <th>mouse</th>\n",
       "      <th>date</th>\n",
       "      <th>total_trial_ind</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td> 400</td>\n",
       "      <td> bird</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 1</td>\n",
       "      <td>  1.840244</td>\n",
       "      <td> 53</td>\n",
       "      <td> 0.267095</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2152.301946</td>\n",
       "      <td> 3937</td>\n",
       "      <td> 0.099917</td>\n",
       "      <td>  1.151120</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td> 550</td>\n",
       "      <td> bird</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 2</td>\n",
       "      <td>  3.024151</td>\n",
       "      <td> 36</td>\n",
       "      <td> 0.268031</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2155.681332</td>\n",
       "      <td> 2455</td>\n",
       "      <td> 0.100434</td>\n",
       "      <td>  0.995480</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td> 350</td>\n",
       "      <td> bird</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 3</td>\n",
       "      <td>  4.273621</td>\n",
       "      <td> 54</td>\n",
       "      <td> 0.399960</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 1231.647159</td>\n",
       "      <td> 2661</td>\n",
       "      <td> 0.115117</td>\n",
       "      <td>  1.148952</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td> 550</td>\n",
       "      <td> bird</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 4</td>\n",
       "      <td>  7.623390</td>\n",
       "      <td> 92</td>\n",
       "      <td> 0.667990</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 1416.589407</td>\n",
       "      <td> 4376</td>\n",
       "      <td> 0.200129</td>\n",
       "      <td> 11.667522</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td> 450</td>\n",
       "      <td> bird</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 5</td>\n",
       "      <td> 12.638539</td>\n",
       "      <td> 52</td>\n",
       "      <td> 0.232075</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2709.459253</td>\n",
       "      <td> 3349</td>\n",
       "      <td> 0.115528</td>\n",
       "      <td>  0.995757</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Size Image  Ori  PosY  trial       time  terraincode  pausetime choice  selectiontime istarget        speed  lapdistance  timebeforezone  slowing_f    mouse                        date  \\\n",
       "0   400  bird    0     0      1   1.840244           53   0.267095  False              1    False  2152.301946         3937        0.099917   1.151120  M130587  2014-10-03 13:29:50.625000   \n",
       "1   550  bird  -15     0      2   3.024151           36   0.268031  False              1    False  2155.681332         2455        0.100434   0.995480  M130587  2014-10-03 13:29:50.625000   \n",
       "2   350  bird    0     0      3   4.273621           54   0.399960  False              1    False  1231.647159         2661        0.115117   1.148952  M130587  2014-10-03 13:29:50.625000   \n",
       "3   550  bird   45     0      4   7.623390           92   0.667990  False              1    False  1416.589407         4376        0.200129  11.667522  M130587  2014-10-03 13:29:50.625000   \n",
       "4   450  bird    0     0      5  12.638539           52   0.232075  False              1    False  2709.459253         3349        0.115528   0.995757  M130587  2014-10-03 13:29:50.625000   \n",
       "\n",
       "   total_trial_ind  \n",
       "0                0  \n",
       "1                1  \n",
       "2                2  \n",
       "3                3  \n",
       "4                4  "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#open file with data -  actually there are two important files!\n",
    "fpath = r\"C:\\Users\\dorisv\\Desktop\\AllenInstWork\\Behavior Data for Doris\\avg_response_df.pkl\"\n",
    "f = open(fpath,'rb')\n",
    "\n",
    "dd = pickle.load(f)\n",
    "    \n",
    "gpath = r\"C:\\Users\\dorisv\\Desktop\\AllenInstWork\\Behavior Data for Doris\\trial_df.pkl\"\n",
    "g = open(gpath,'rb')\n",
    "df = pickle.load(g)\n",
    "\n",
    "df[\"total_trial_ind\"] = np.arange(len(df))\n",
    "dd.Image[dd.Image==\"CoxImage1.png\"] = \"mushroom\"\n",
    "dd.Image[dd.Image==\"CoxImage2.png\"] = \"bird\"\n",
    "\n",
    "df.Image[df.Image==\"CoxImage1.png\"] = \"mushroom\"\n",
    "df.Image[df.Image==\"CoxImage2.png\"] = \"bird\"\n",
    "\n",
    "sizes = np.sort(df.Size.unique())\n",
    "orientations = np.sort(df.Ori.unique())\n",
    "\n",
    "#for the purposes of the experiment, should be m=n=7\n",
    "m = len(sizes)\n",
    "n = len(orientations)\n",
    "\n",
    "evenTest1 = isinstance((m/float(2)), (int, long))\n",
    "evenTest2 = isinstance((n/float(2)), (int, long))\n",
    "    \n",
    "if evenTest1:\n",
    "    s_T = sizes[m/2-1]\n",
    "else:\n",
    "    s_T = sizes[m/2]\n",
    "        \n",
    "if evenTest2:\n",
    "    o_T = orientations[n/2-1]    \n",
    "else:\n",
    "    o_T = orientations[n/2]\n",
    "                \n",
    "s_max = sizes[m-1]\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Size</th>\n",
       "      <th>Image</th>\n",
       "      <th>Ori</th>\n",
       "      <th>PosY</th>\n",
       "      <th>trial</th>\n",
       "      <th>time</th>\n",
       "      <th>terraincode</th>\n",
       "      <th>pausetime</th>\n",
       "      <th>choice</th>\n",
       "      <th>selectiontime</th>\n",
       "      <th>istarget</th>\n",
       "      <th>speed</th>\n",
       "      <th>lapdistance</th>\n",
       "      <th>timebeforezone</th>\n",
       "      <th>slowing_f</th>\n",
       "      <th>mouse</th>\n",
       "      <th>date</th>\n",
       "      <th>total_trial_ind</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td> 500</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 101</td>\n",
       "      <td> 272.104879</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0.382876</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2000.578984</td>\n",
       "      <td> 2819</td>\n",
       "      <td> 0.115008</td>\n",
       "      <td>  1.155547</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>101</th>\n",
       "      <td> 350</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 102</td>\n",
       "      <td> 274.454671</td>\n",
       "      <td> 96</td>\n",
       "      <td> 0.265104</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2378.636351</td>\n",
       "      <td> 4341</td>\n",
       "      <td> 0.100000</td>\n",
       "      <td>  0.999937</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>102</th>\n",
       "      <td> 350</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 103</td>\n",
       "      <td> 276.754522</td>\n",
       "      <td> 96</td>\n",
       "      <td> 0.249983</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2269.571773</td>\n",
       "      <td> 5221</td>\n",
       "      <td> 0.099999</td>\n",
       "      <td>  1.000006</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>103</th>\n",
       "      <td> 600</td>\n",
       "      <td>     bird</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 104</td>\n",
       "      <td> 278.989269</td>\n",
       "      <td> 49</td>\n",
       "      <td> 0.335187</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2122.946592</td>\n",
       "      <td> 4922</td>\n",
       "      <td> 0.100393</td>\n",
       "      <td>  1.166350</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 103</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>104</th>\n",
       "      <td> 600</td>\n",
       "      <td> mushroom</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 105</td>\n",
       "      <td> 280.271210</td>\n",
       "      <td> 84</td>\n",
       "      <td> 2.735020</td>\n",
       "      <td>  True</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td>  573.372541</td>\n",
       "      <td> 2432</td>\n",
       "      <td> 0.083366</td>\n",
       "      <td>  1.199498</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 104</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>105</th>\n",
       "      <td> 600</td>\n",
       "      <td> mushroom</td>\n",
       "      <td> 15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 106</td>\n",
       "      <td> 288.404010</td>\n",
       "      <td> 56</td>\n",
       "      <td> 0.231864</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2440.298151</td>\n",
       "      <td> 4071</td>\n",
       "      <td> 0.083307</td>\n",
       "      <td>  0.982307</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>106</th>\n",
       "      <td> 500</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 107</td>\n",
       "      <td> 289.888518</td>\n",
       "      <td> 79</td>\n",
       "      <td> 5.134661</td>\n",
       "      <td>  True</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td>  996.488420</td>\n",
       "      <td> 3569</td>\n",
       "      <td> 0.100432</td>\n",
       "      <td> 80.824983</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 106</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107</th>\n",
       "      <td> 300</td>\n",
       "      <td>     bird</td>\n",
       "      <td>-45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 108</td>\n",
       "      <td> 299.587250</td>\n",
       "      <td> 13</td>\n",
       "      <td> 0.232888</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2018.135034</td>\n",
       "      <td> 3332</td>\n",
       "      <td> 0.082828</td>\n",
       "      <td>  1.206668</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 107</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108</th>\n",
       "      <td> 400</td>\n",
       "      <td>     bird</td>\n",
       "      <td>-30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 109</td>\n",
       "      <td> 300.853070</td>\n",
       "      <td> 25</td>\n",
       "      <td> 0.249982</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2102.311724</td>\n",
       "      <td> 2692</td>\n",
       "      <td> 0.100364</td>\n",
       "      <td>  0.996328</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 108</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>109</th>\n",
       "      <td> 400</td>\n",
       "      <td>     bird</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 110</td>\n",
       "      <td> 302.337151</td>\n",
       "      <td> 39</td>\n",
       "      <td> 0.249979</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2240.089092</td>\n",
       "      <td> 3107</td>\n",
       "      <td> 0.100420</td>\n",
       "      <td>  1.173070</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td> 600</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 111</td>\n",
       "      <td> 304.652576</td>\n",
       "      <td> 42</td>\n",
       "      <td> 1.917979</td>\n",
       "      <td>  True</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td>  785.294610</td>\n",
       "      <td> 5067</td>\n",
       "      <td> 0.099992</td>\n",
       "      <td>  1.180810</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>111</th>\n",
       "      <td> 600</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 112</td>\n",
       "      <td> 310.336608</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0.235153</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2539.956953</td>\n",
       "      <td> 3665</td>\n",
       "      <td> 0.100481</td>\n",
       "      <td>  0.848274</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>112</th>\n",
       "      <td> 400</td>\n",
       "      <td>     bird</td>\n",
       "      <td>-30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 113</td>\n",
       "      <td> 311.468634</td>\n",
       "      <td> 25</td>\n",
       "      <td> 0.215024</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2323.958837</td>\n",
       "      <td> 2615</td>\n",
       "      <td> 0.083233</td>\n",
       "      <td>  0.997026</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>113</th>\n",
       "      <td> 600</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 114</td>\n",
       "      <td> 313.785088</td>\n",
       "      <td> 28</td>\n",
       "      <td> 2.282738</td>\n",
       "      <td>  True</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td>  996.217790</td>\n",
       "      <td> 5392</td>\n",
       "      <td> 0.099999</td>\n",
       "      <td>  1.180852</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>114</th>\n",
       "      <td> 500</td>\n",
       "      <td> mushroom</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 115</td>\n",
       "      <td> 320.568518</td>\n",
       "      <td> 72</td>\n",
       "      <td> 0.234974</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2215.131643</td>\n",
       "      <td> 5009</td>\n",
       "      <td> 0.083345</td>\n",
       "      <td>  1.020778</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>115</th>\n",
       "      <td> 600</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 116</td>\n",
       "      <td> 321.701552</td>\n",
       "      <td> 28</td>\n",
       "      <td> 2.517901</td>\n",
       "      <td>  True</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td>  365.109852</td>\n",
       "      <td> 2642</td>\n",
       "      <td> 0.100466</td>\n",
       "      <td>  0.848074</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>116</th>\n",
       "      <td> 350</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 117</td>\n",
       "      <td> 335.100693</td>\n",
       "      <td> 33</td>\n",
       "      <td> 0.200007</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2486.574838</td>\n",
       "      <td> 5037</td>\n",
       "      <td> 0.082829</td>\n",
       "      <td>  0.992744</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117</th>\n",
       "      <td> 450</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 118</td>\n",
       "      <td> 336.967432</td>\n",
       "      <td> 80</td>\n",
       "      <td> 0.249982</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2031.691896</td>\n",
       "      <td> 4601</td>\n",
       "      <td> 0.100414</td>\n",
       "      <td>  1.174133</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 117</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118</th>\n",
       "      <td> 350</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 119</td>\n",
       "      <td> 338.366489</td>\n",
       "      <td> 82</td>\n",
       "      <td> 0.282900</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 1906.267066</td>\n",
       "      <td> 2840</td>\n",
       "      <td> 0.099966</td>\n",
       "      <td>  1.000347</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>119</th>\n",
       "      <td> 300</td>\n",
       "      <td>     bird</td>\n",
       "      <td>-30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 120</td>\n",
       "      <td> 340.716972</td>\n",
       "      <td> 27</td>\n",
       "      <td> 0.284779</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2509.558577</td>\n",
       "      <td> 4523</td>\n",
       "      <td> 0.115157</td>\n",
       "      <td>  0.868235</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>120</th>\n",
       "      <td> 300</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 121</td>\n",
       "      <td> 342.716413</td>\n",
       "      <td> 83</td>\n",
       "      <td> 0.267103</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 2098.499372</td>\n",
       "      <td> 4602</td>\n",
       "      <td> 0.100003</td>\n",
       "      <td>  0.999264</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 120</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>121</th>\n",
       "      <td> 500</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 122</td>\n",
       "      <td> 344.485596</td>\n",
       "      <td> 30</td>\n",
       "      <td> 0.385039</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td>  286.002002</td>\n",
       "      <td> 3682</td>\n",
       "      <td> 0.100371</td>\n",
       "      <td>  1.167157</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 121</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>122</th>\n",
       "      <td> 550</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>-30</td>\n",
       "      <td> 0</td>\n",
       "      <td> 123</td>\n",
       "      <td> 349.485320</td>\n",
       "      <td> 15</td>\n",
       "      <td> 0.199973</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2845.386312</td>\n",
       "      <td> 2501</td>\n",
       "      <td> 0.100149</td>\n",
       "      <td>  0.647218</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>123</th>\n",
       "      <td> 300</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 124</td>\n",
       "      <td> 351.266169</td>\n",
       "      <td> 48</td>\n",
       "      <td> 0.249960</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2022.695246</td>\n",
       "      <td> 4878</td>\n",
       "      <td> 0.083213</td>\n",
       "      <td>  0.996017</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>124</th>\n",
       "      <td> 300</td>\n",
       "      <td> mushroom</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 125</td>\n",
       "      <td> 353.636024</td>\n",
       "      <td> 90</td>\n",
       "      <td> 0.249992</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2298.511523</td>\n",
       "      <td> 4878</td>\n",
       "      <td> 0.099612</td>\n",
       "      <td>  0.825685</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>125</th>\n",
       "      <td> 500</td>\n",
       "      <td>     bird</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 126</td>\n",
       "      <td> 355.798941</td>\n",
       "      <td> 93</td>\n",
       "      <td> 0.517877</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td> False</td>\n",
       "      <td> 1460.909244</td>\n",
       "      <td> 4861</td>\n",
       "      <td> 0.100361</td>\n",
       "      <td>  1.494564</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>126</th>\n",
       "      <td> 500</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 127</td>\n",
       "      <td> 359.348928</td>\n",
       "      <td> 44</td>\n",
       "      <td> 0.317918</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 1732.767801</td>\n",
       "      <td> 5043</td>\n",
       "      <td> 0.082942</td>\n",
       "      <td>  1.205699</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 126</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>127</th>\n",
       "      <td> 400</td>\n",
       "      <td> mushroom</td>\n",
       "      <td> 45</td>\n",
       "      <td> 0</td>\n",
       "      <td> 128</td>\n",
       "      <td> 362.048781</td>\n",
       "      <td> 88</td>\n",
       "      <td> 0.299984</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2028.999396</td>\n",
       "      <td> 4812</td>\n",
       "      <td> 0.100033</td>\n",
       "      <td>  0.999692</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>128</th>\n",
       "      <td> 300</td>\n",
       "      <td> mushroom</td>\n",
       "      <td> 15</td>\n",
       "      <td> 0</td>\n",
       "      <td> 129</td>\n",
       "      <td> 364.065155</td>\n",
       "      <td> 62</td>\n",
       "      <td> 0.284890</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2290.439556</td>\n",
       "      <td> 4050</td>\n",
       "      <td> 0.100390</td>\n",
       "      <td>  0.844853</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 128</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>129</th>\n",
       "      <td> 300</td>\n",
       "      <td> mushroom</td>\n",
       "      <td>  0</td>\n",
       "      <td> 0</td>\n",
       "      <td> 130</td>\n",
       "      <td> 365.485172</td>\n",
       "      <td> 48</td>\n",
       "      <td> 0.199987</td>\n",
       "      <td> False</td>\n",
       "      <td> 1</td>\n",
       "      <td>  True</td>\n",
       "      <td> 2443.832035</td>\n",
       "      <td> 3075</td>\n",
       "      <td> 0.082521</td>\n",
       "      <td>  1.004444</td>\n",
       "      <td> M130587</td>\n",
       "      <td> 2014-10-03 13:29:50.625000</td>\n",
       "      <td> 129</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Size     Image  Ori  PosY  trial        time  terraincode  pausetime choice  selectiontime istarget        speed  lapdistance  timebeforezone  slowing_f    mouse                        date  \\\n",
       "100   500  mushroom  -15     0    101  272.104879           30   0.382876  False              1     True  2000.578984         2819        0.115008   1.155547  M130587  2014-10-03 13:29:50.625000   \n",
       "101   350      bird   45     0    102  274.454671           96   0.265104  False              1    False  2378.636351         4341        0.100000   0.999937  M130587  2014-10-03 13:29:50.625000   \n",
       "102   350      bird   45     0    103  276.754522           96   0.249983  False              1    False  2269.571773         5221        0.099999   1.000006  M130587  2014-10-03 13:29:50.625000   \n",
       "103   600      bird    0     0    104  278.989269           49   0.335187  False              1    False  2122.946592         4922        0.100393   1.166350  M130587  2014-10-03 13:29:50.625000   \n",
       "104   600  mushroom   45     0    105  280.271210           84   2.735020   True              1     True   573.372541         2432        0.083366   1.199498  M130587  2014-10-03 13:29:50.625000   \n",
       "105   600  mushroom   15     0    106  288.404010           56   0.231864  False              1     True  2440.298151         4071        0.083307   0.982307  M130587  2014-10-03 13:29:50.625000   \n",
       "106   500      bird   30     0    107  289.888518           79   5.134661   True              1    False   996.488420         3569        0.100432  80.824983  M130587  2014-10-03 13:29:50.625000   \n",
       "107   300      bird  -45     0    108  299.587250           13   0.232888  False              1    False  2018.135034         3332        0.082828   1.206668  M130587  2014-10-03 13:29:50.625000   \n",
       "108   400      bird  -30     0    109  300.853070           25   0.249982  False              1    False  2102.311724         2692        0.100364   0.996328  M130587  2014-10-03 13:29:50.625000   \n",
       "109   400      bird  -15     0    110  302.337151           39   0.249979  False              1    False  2240.089092         3107        0.100420   1.173070  M130587  2014-10-03 13:29:50.625000   \n",
       "110   600  mushroom    0     0    111  304.652576           42   1.917979   True              1     True   785.294610         5067        0.099992   1.180810  M130587  2014-10-03 13:29:50.625000   \n",
       "111   600  mushroom  -45     0    112  310.336608            0   0.235153  False              1     True  2539.956953         3665        0.100481   0.848274  M130587  2014-10-03 13:29:50.625000   \n",
       "112   400      bird  -30     0    113  311.468634           25   0.215024  False              1    False  2323.958837         2615        0.083233   0.997026  M130587  2014-10-03 13:29:50.625000   \n",
       "113   600  mushroom  -15     0    114  313.785088           28   2.282738   True              1     True   996.217790         5392        0.099999   1.180852  M130587  2014-10-03 13:29:50.625000   \n",
       "114   500  mushroom   30     0    115  320.568518           72   0.234974  False              1     True  2215.131643         5009        0.083345   1.020778  M130587  2014-10-03 13:29:50.625000   \n",
       "115   600  mushroom  -15     0    116  321.701552           28   2.517901   True              1     True   365.109852         2642        0.100466   0.848074  M130587  2014-10-03 13:29:50.625000   \n",
       "116   350  mushroom  -15     0    117  335.100693           33   0.200007  False              1     True  2486.574838         5037        0.082829   0.992744  M130587  2014-10-03 13:29:50.625000   \n",
       "117   450      bird   30     0    118  336.967432           80   0.249982  False              1    False  2031.691896         4601        0.100414   1.174133  M130587  2014-10-03 13:29:50.625000   \n",
       "118   350      bird   30     0    119  338.366489           82   0.282900  False              1    False  1906.267066         2840        0.099966   1.000347  M130587  2014-10-03 13:29:50.625000   \n",
       "119   300      bird  -30     0    120  340.716972           27   0.284779  False              1    False  2509.558577         4523        0.115157   0.868235  M130587  2014-10-03 13:29:50.625000   \n",
       "120   300      bird   30     0    121  342.716413           83   0.267103  False              1    False  2098.499372         4602        0.100003   0.999264  M130587  2014-10-03 13:29:50.625000   \n",
       "121   500  mushroom  -15     0    122  344.485596           30   0.385039  False              1     True   286.002002         3682        0.100371   1.167157  M130587  2014-10-03 13:29:50.625000   \n",
       "122   550  mushroom  -30     0    123  349.485320           15   0.199973  False              1     True  2845.386312         2501        0.100149   0.647218  M130587  2014-10-03 13:29:50.625000   \n",
       "123   300  mushroom    0     0    124  351.266169           48   0.249960  False              1     True  2022.695246         4878        0.083213   0.996017  M130587  2014-10-03 13:29:50.625000   \n",
       "124   300  mushroom   45     0    125  353.636024           90   0.249992  False              1     True  2298.511523         4878        0.099612   0.825685  M130587  2014-10-03 13:29:50.625000   \n",
       "125   500      bird   45     0    126  355.798941           93   0.517877  False              1    False  1460.909244         4861        0.100361   1.494564  M130587  2014-10-03 13:29:50.625000   \n",
       "126   500  mushroom    0     0    127  359.348928           44   0.317918  False              1     True  1732.767801         5043        0.082942   1.205699  M130587  2014-10-03 13:29:50.625000   \n",
       "127   400  mushroom   45     0    128  362.048781           88   0.299984  False              1     True  2028.999396         4812        0.100033   0.999692  M130587  2014-10-03 13:29:50.625000   \n",
       "128   300  mushroom   15     0    129  364.065155           62   0.284890  False              1     True  2290.439556         4050        0.100390   0.844853  M130587  2014-10-03 13:29:50.625000   \n",
       "129   300  mushroom    0     0    130  365.485172           48   0.199987  False              1     True  2443.832035         3075        0.082521   1.004444  M130587  2014-10-03 13:29:50.625000   \n",
       "\n",
       "     total_trial_ind  \n",
       "100              100  \n",
       "101              101  \n",
       "102              102  \n",
       "103              103  \n",
       "104              104  \n",
       "105              105  \n",
       "106              106  \n",
       "107              107  \n",
       "108              108  \n",
       "109              109  \n",
       "110              110  \n",
       "111              111  \n",
       "112              112  \n",
       "113              113  \n",
       "114              114  \n",
       "115              115  \n",
       "116              116  \n",
       "117              117  \n",
       "118              118  \n",
       "119              119  \n",
       "120              120  \n",
       "121              121  \n",
       "122              122  \n",
       "123              123  \n",
       "124              124  \n",
       "125              125  \n",
       "126              126  \n",
       "127              127  \n",
       "128              128  \n",
       "129              129  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[100:130]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def transform_intoPercept(prob_misID):\n",
    "    \n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "    percept = prob_misID*(p_b_real-p_m_real)+p_m_real  #notation: p(i_p=m) = percept !!\n",
    "    \n",
    "    percept_givenM = 1-prob_misID\n",
    "    percept_givenB = prob_misID\n",
    "    \n",
    "    return percept, percept_givenM, percept_givenB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def transform_into_prob_misID(percept):\n",
    "\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "    prob_misID = (percept-p_m_real)/(p_b_real-p_m_real)\n",
    "    \n",
    "    return prob_misID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def probability_misidentify(x):\n",
    "    \n",
    "    #parameters of the model\n",
    "    #k_s = x[0]\n",
    "    #k_o = x[1]\n",
    "    #k_d = x[2]    \n",
    "    #alpha = x[3]\n",
    "    alpha = x[0]\n",
    "    k_o = x[1]\n",
    "    k_s = x[2]\n",
    "    k_d = x[3]\n",
    "    x_0 = 1\n",
    "    \n",
    "    print alpha, k_o, k_s, k_d\n",
    "    \n",
    "    #probabilities of misidentification, from which we infer p-stop. The prob_misID depend on s,o, s_T,o_T\n",
    "    prob_misID = np.zeros((len(orientations), len(sizes)))\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):               \n",
    "                    #our model for prob_msID is a hill function\n",
    "            \n",
    "            if s!=0:\n",
    "                prob_misID[ind2,ind1] = 0.5*1/(1+x_0/(k_s*math.fabs(s-s_T)+k_o*math.fabs(o-o_T)+k_d*(s_max-s)+alpha))\n",
    "                #prob_misID[ind2,ind1] = k_s*math.fabs(s-s_T)+k_o*math.fabs(o-o_T)+k_d*(s_max-s)+alpha\n",
    "            else:\n",
    "                prob_misID[ind2,ind1] = 0.5\n",
    "                 #another version of the model   \n",
    "                #prob_misID[i,ind1,ind2] = 0.5*1/(1+x_0/(k_s*math.fabs(s-s_T)+k_o*math.fabs(o-o_T)+k_d*math.fabs*(s_max-s)+alpha))\n",
    "    \n",
    "    return prob_misID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def p_stop_ProbabilityMatching(x):\n",
    "    \n",
    "    p_stop_model1 = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_givenM_model = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_givenB_model = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    prob_misID = probability_misidentify(x)\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "                \n",
    "            p_stop_model1[ind2,ind1] = (1-prob_misID[ind2,ind1])*p_m_real[ind2,ind1]+prob_misID[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "            p_stop_givenM_model[ind2,ind1] = 1-prob_misID[ind2,ind1]\n",
    "            p_stop_givenB_model[ind2,ind1] = prob_misID[ind2,ind1]\n",
    "            \n",
    "    return p_stop_model1, p_stop_givenM_model, p_stop_givenB_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.56338028  0.57352941  0.49253731  0.65217391  0.609375    0.76811594\n",
      "   0.82539683]\n",
      " [ 0.6         0.55172414  0.70769231  0.6         0.77419355  0.68656716\n",
      "   0.83333333]\n",
      " [ 0.65671642  0.59375     0.64179104  0.63333333  0.61538462  0.6875\n",
      "   0.71875   ]\n",
      " [ 0.54411765  0.53731343  0.65151515  0.59677419  0.59677419  0.68571429\n",
      "   0.81428571]\n",
      " [ 0.60273973  0.64705882  0.65217391  0.70769231  0.83823529  0.80645161\n",
      "   0.75384615]\n",
      " [ 0.6875      0.64705882  0.7704918   0.74193548  0.77941176  0.78125\n",
      "   0.79104478]\n",
      " [ 0.75714286  0.62318841  0.71212121  0.80327869  0.64285714  0.79310345\n",
      "   0.75      ]] [[ 0.33333333  0.32631579  0.27368421  0.19587629  0.26530612  0.31132075\n",
      "   0.30107527]\n",
      " [ 0.23762376  0.20408163  0.2244898   0.1875      0.21428571  0.19148936\n",
      "   0.34090909]\n",
      " [ 0.27619048  0.25        0.17        0.14583333  0.17647059  0.29347826\n",
      "   0.31632653]\n",
      " [ 0.10204082  0.08247423  0.10416667  0.12380952  0.10784314  0.12745098\n",
      "   0.20618557]\n",
      " [ 0.25        0.10784314  0.12        0.0990099   0.12121212  0.16666667\n",
      "   0.26262626]\n",
      " [ 0.18888889  0.16842105  0.21782178  0.20618557  0.18556701  0.17894737\n",
      "   0.24175824]\n",
      " [ 0.34693878  0.28712871  0.22680412  0.3         0.29        0.37373737\n",
      "   0.33695652]]\n",
      "49\n"
     ]
    }
   ],
   "source": [
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "\n",
    "print trials_stop_m/trials_m, trials_stop_b/trials_b\n",
    "print sum(trials_stop_m/trials_m > trials_stop_b/trials_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def p_real():\n",
    "    \n",
    "    p_m_real = np.zeros((len(orientations),len(sizes)))\n",
    "    p_b_real = np.zeros((len(orientations),len(sizes)))\n",
    "    trials_tot = np.zeros((len(orientations),len(sizes)))\n",
    "    trials_stop = np.zeros((len(orientations),len(sizes)))\n",
    "    trials_m = np.zeros((len(orientations),len(sizes)))\n",
    "    trials_b = np.zeros((len(orientations),len(sizes)))\n",
    "    trials_stop_m = np.zeros((len(orientations),len(sizes)))\n",
    "    trials_stop_b = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "                \n",
    "            trials_tot[ind2,ind1] = len(df[(df.Size==s)&(df.Ori==o)].total_trial_ind)\n",
    "            trials_stop[ind2,ind1] = len(df[(df.Size==s)&(df.Ori==o)&(df.choice==True)].total_trial_ind)\n",
    "            \n",
    "            trials_m[ind2,ind1] = len(df[(df.Size==s)&(df.Ori==o)&(df.Image=='mushroom')].total_trial_ind)\n",
    "            trials_b[ind2,ind1] = len(df[(df.Size==s)&(df.Ori==o)&(df.Image=='bird')].total_trial_ind)   \n",
    "            \n",
    "            trials_stop_m[ind2,ind1] = len(df[(df.Size==s)&(df.Ori==o)&(df.Image=='mushroom')&(df.choice==True)].total_trial_ind)\n",
    "            trials_stop_b[ind2,ind1] = len(df[(df.Size==s)&(df.Ori==o)&(df.Image=='bird')&(df.choice==True)].total_trial_ind)   \n",
    "            \n",
    "            p_m_real[ind2, ind1] = trials_m[ind2,ind1]/trials_tot[ind2,ind1]\n",
    "            p_b_real[ind2, ind1] = trials_b[ind2,ind1]/trials_tot[ind2,ind1]   \n",
    "        \n",
    "    return p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m 0.683680029514 0.089285528193\n",
      "b 0.780284188962 0.0777394874638\n"
     ]
    }
   ],
   "source": [
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b\n",
    "\n",
    "correct_m = np.zeros((len(orientations), len(sizes)))\n",
    "correct_b = np.zeros((len(orientations), len(sizes)))\n",
    "\n",
    "for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):        \n",
    "            correct_m[ind2,ind1] = float(len(df.total_trial_ind[(df.Image == 'mushroom') & (df.Size==s) & (df.Ori==o) & (df.choice == 1)]))/trials_m[ind2,ind1]\n",
    "            correct_b[ind2,ind1] = float(len(df.total_trial_ind[(df.Image == 'bird') & (df.Size==s) & (df.Ori==o) & (df.choice == 0)]))/trials_b[ind2,ind1]\n",
    "            \n",
    "print \"m\", np.mean(correct_m), np.std(correct_m)    \n",
    "print \"b\", np.mean(correct_b), np.std(correct_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def p_stop_real_function():\n",
    "    \n",
    "    p_stop_real = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_toFit = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_toFit2 = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenM2 = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenB2 = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenM_toFit = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenB_toFit = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenM_toFit2 = np.zeros((len(orientations),len(sizes)))\n",
    "    p_stop_real_givenB_toFit2 = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b  = p_real()\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            \n",
    "            i1 = dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"mushroom\")].yp.index[0]\n",
    "            i2 = dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"bird\")].yp.index[0]\n",
    "            \n",
    "            p_stop_real[ind2,ind1] = dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"mushroom\")].yp[i1]*p_m_real[ind2,ind1]+dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"bird\")].yp[i2]*p_b_real[ind2,ind1]\n",
    "            p_stop_real_givenM[ind2,ind1] = dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"mushroom\")].yp[i1]\n",
    "            p_stop_real_givenB[ind2,ind1] = dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"bird\")].yp[i2]\n",
    "        \n",
    "            p_stop_real_toFit[ind2,ind1]  = (p_stop_real[ind2,ind1]-p_m_real[ind2,ind1])/(p_b_real[ind2,ind1]-p_m_real[ind2,ind1])\n",
    "            p_stop_real_toFit2[ind2,ind1] = 2*p_stop_real_toFit[ind2,ind1]/(1-2*p_stop_real_toFit[ind2,ind1])\n",
    "            \n",
    "            p_stop_real_givenM_toFit[ind2,ind1]  = 1-p_stop_real_givenM[ind2,ind1]\n",
    "            p_stop_real_givenM_toFit2[ind2,ind1] = 2*p_stop_real_givenM_toFit[ind2,ind1]/(1-2*p_stop_real_givenM_toFit[ind2,ind1])\n",
    "            p_stop_real_givenB_toFit[ind2,ind1]  = p_stop_real_givenB[ind2,ind1]\n",
    "            p_stop_real_givenB_toFit2[ind2,ind1] = 2*p_stop_real_givenB_toFit[ind2,ind1]/(1-2*p_stop_real_givenB_toFit[ind2,ind1])\n",
    "            \n",
    "    return p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_real_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def p_stop_AND_correct_function():\n",
    "    \n",
    "    p_stop_AND_correct = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            i1 = dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"mushroom\")].yp.index[0]\n",
    "            p_stop_AND_correct[ind2,ind1] =  dd[(dd.Size==s)&(dd.Ori==o)&(dd.Image==\"mushroom\")].yp[i1]\n",
    "            \n",
    "    return p_stop_AND_correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "global name 'orientations' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-18-362fafba16ad>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m#plot: p(i_p != i_r) or hill^{-1}(2*p(i_p != i_r)). One is better fit to the data, the other has outliers and is more \"flat\" (non-statistically significant coefficients found when regressing)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mp_stop_real\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_toFit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_toFit2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenM\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_givenB\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenM_toFit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenM_toFit2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenB_toFit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenB_toFit2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_stop_real_function\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;31m#print p_stop_real_toFit\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-16-343018d8fa23>\u001b[0m in \u001b[0;36mp_stop_real_function\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mp_stop_real_function\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m     \u001b[0mp_stop_real\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0morientations\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msizes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m     \u001b[0mp_stop_real_toFit\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0morientations\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msizes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0mp_stop_real_toFit2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0morientations\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msizes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: global name 'orientations' is not defined"
     ]
    }
   ],
   "source": [
    "#plot: p(i_p != i_r) or hill^{-1}(2*p(i_p != i_r)). One is better fit to the data, the other has outliers and is more \"flat\" (non-statistically significant coefficients found when regressing)\n",
    "\n",
    "p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "\n",
    "#print p_stop_real_toFit\n",
    "#print p_stop_real_toFit2\n",
    "\n",
    "or1 = np.zeros((len(orientations)*len(sizes),1))\n",
    "s1 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp1 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp2 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp3 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp4 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp5 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp6 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp7 = np.zeros((len(orientations)*len(sizes),1))\n",
    "pp8 = np.zeros((len(orientations)*len(sizes),1))\n",
    "or5 = np.zeros((len(orientations)*len(sizes),1))\n",
    "s5 = np.zeros((len(orientations)*len(sizes),1))\n",
    "\n",
    "evenTest1 = isinstance((m/float(2)), (int, long))\n",
    "evenTest2 = isinstance((n/float(2)), (int, long))\n",
    "    \n",
    "if evenTest1:\n",
    "    s_T = sizes[m/2-1]\n",
    "else:\n",
    "    s_T = sizes[m/2]\n",
    "        \n",
    "if evenTest2:\n",
    "    o_T = orientations[n/2-1]    \n",
    "else:\n",
    "    o_T = orientations[n/2]\n",
    "    \n",
    "s_max = sizes[m-1]\n",
    "\n",
    "for ind2,o in enumerate(orientations):\n",
    "    for ind1,s in enumerate(sizes):\n",
    "        or1[ind2*len(orientations)+ind1] = math.fabs(orientations[ind2]-o_T)\n",
    "        s1[ind2*len(orientations)+ind1] = math.fabs(sizes[ind1]-s_T)\n",
    "        pp1[ind2*len(orientations)+ind1] = p_stop_real_toFit[ind2,ind1]\n",
    "        pp2[ind2*len(orientations)+ind1] = p_stop_real_toFit2[ind2,ind1]\n",
    "        pp3[ind2*len(orientations)+ind1] = p_stop_real_givenM_toFit2[ind2,ind1]\n",
    "        pp4[ind2*len(orientations)+ind1] = p_stop_real_givenB_toFit2[ind2,ind1]\n",
    "        pp6[ind2*len(orientations)+ind1] = p_stop_real_givenM_toFit[ind2,ind1]\n",
    "        pp7[ind2*len(orientations)+ind1] = p_stop_real_givenB_toFit[ind2,ind1]\n",
    "        \n",
    "pp5 = np.concatenate((pp3,pp4)) \n",
    "pp8 = np.concatenate((pp6,pp7))\n",
    "\n",
    "or5 = np.concatenate((or1,or1))\n",
    "s5 = np.concatenate((s1,s1))\n",
    "\n",
    "# create matplotlib 3d axes\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = Axes3D(fig, azim=-115, elev=15)\n",
    "\n",
    "#plot the actual data\n",
    "ax.scatter(or1, s1, pp1, color='black', alpha=1.0)\n",
    "\n",
    "ax.set_xlabel('|size-s_T|')\n",
    "ax.set_ylabel('|orientation-o_T|')\n",
    "ax.set_zlabel('p(i_p != i_r)')\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = Axes3D(fig, azim=-115, elev=15)\n",
    "\n",
    "ax.scatter(or1, s1, pp2, color='black', alpha=1.0)\n",
    "   \n",
    "# set axis labels\n",
    "ax.set_xlabel('|size-s_T|')\n",
    "ax.set_ylabel('|orientation-o_T|')\n",
    "ax.set_zlabel('hill^{-1}(p(i_p != i_r))')\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = Axes3D(fig, azim=-115, elev=15)\n",
    "\n",
    "ax.scatter(or5, s5, pp5, color='black', alpha=1.0)\n",
    "# set axis labels\n",
    "ax.set_xlabel('|size-s_T|')\n",
    "ax.set_ylabel('|orientation-o_T|')\n",
    "ax.set_zlabel('hill^{-1}(p(i_p != i_r))')\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "ax = Axes3D(fig, azim=-115, elev=15)\n",
    "\n",
    "ax.scatter(or5, s5, pp8, color='black', alpha=1.0)\n",
    "# set axis labels\n",
    "ax.set_xlabel('|size-s_T|')\n",
    "ax.set_ylabel('|orientation-o_T|')\n",
    "ax.set_zlabel('hill^{-1}(p(i_p != i_r))')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#determine outliers in the data, in the limits limit_low and limit_high\n",
    "\n",
    "def determine_outliers():\n",
    "\n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "\n",
    "    outliers_orientationsM = []\n",
    "    outliers_sizeM = []\n",
    "    \n",
    "    outliers_orientationsB = []\n",
    "    outliers_sizeB = []\n",
    "    \n",
    "    limit_low = -1\n",
    "    limit_high = 1\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            \n",
    "            if p_stop_real_givenM_toFit2[ind2,ind1]<limit_low or p_stop_real_givenM_toFit2[ind2,ind1]>limit_high:\n",
    "                outliers_orientationsM.append(ind2)\n",
    "                outliers_sizeM.append(ind1)\n",
    "                \n",
    "            if p_stop_real_givenB_toFit2[ind2,ind1]<limit_low or p_stop_real_givenB_toFit2[ind2,ind1]>limit_high:\n",
    "                outliers_orientationsB.append(ind2)\n",
    "                outliers_sizeB.append(ind1)    \n",
    "                    \n",
    "    return outliers_orientationsM, outliers_sizeM, outliers_orientationsB, outliers_sizeB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fitting_function_withoutOutliers(x0):\n",
    "   \n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_given, Bp_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    \n",
    "    outliers_orientationsM, outliers_sizeM, outliers_orientationsB, outliers_sizeB = determine_outliers()\n",
    "    \n",
    "    evenTest1 = isinstance((m/float(2)), (int, long))\n",
    "    evenTest2 = isinstance((n/float(2)), (int, long))\n",
    "    \n",
    "    if evenTest1:\n",
    "        s_T = sizes[m/2-1]\n",
    "    else:\n",
    "        s_T = sizes[m/2]\n",
    "        \n",
    "    if evenTest2:\n",
    "        o_T = orientations[n/2-1]    \n",
    "    else:\n",
    "        o_T = orientations[n/2]\n",
    "    \n",
    "    s_max = sizes[m-1]\n",
    "    x1 = []\n",
    "    x2 = []\n",
    "    x3 = []\n",
    "    ytemp = []\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "        \n",
    "            ok=0\n",
    "            for i in range(0,len(outliers_orientationsM)):\n",
    "                if (ind2 == outliers_orientationsM[i]) and (ind1 == outliers_sizeM[i]):\n",
    "                    ok=1\n",
    "                    \n",
    "            if ok==0:        \n",
    "                x1.append(math.fabs(orientations[ind2]-o_T))\n",
    "                x2.append(math.fabs(sizes[ind1]-s_T))\n",
    "                x3.append(s_max-sizes[ind1])\n",
    "                #x2.append(1/sizes[ind1])\n",
    "                ytemp.append(x0*p_stop_real_givenM_toFit[ind2,ind1])\n",
    "\n",
    "            ok=0\n",
    "            for i in range(0,len(outliers_orientationsB)):\n",
    "                if (ind2 == outliers_orientationsB[i]) and (ind1 == outliers_sizeB[i]):\n",
    "                    ok=1\n",
    "                    \n",
    "            if ok==0:        \n",
    "                x1.append(math.fabs(orientations[ind2]-o_T))\n",
    "                x2.append(math.fabs(sizes[ind1]-s_T))\n",
    "                x3.append(s_max-sizes[ind1])\n",
    "                #x2.append(1/sizes[ind1])\n",
    "                ytemp.append(x0*p_stop_real_givenB_toFit[ind2,ind1])\n",
    "    \n",
    "    x_i = np.zeros((len(x1),3))\n",
    "    y_i = np.zeros((len(ytemp),1))\n",
    "    \n",
    "    x_i[:,0] = x1\n",
    "    x_i[:,1] = x2\n",
    "    x_i[:,2] = x3\n",
    "    y_i = ytemp\n",
    "    \n",
    "    print len(x_i[:,0])\n",
    "    print len(y_i)\n",
    "    \n",
    "    x_i = sm.add_constant(x_i)\n",
    "    est = sm.OLS(y_i, x_i).fit()\n",
    "    y_est = est.predict(x_i)\n",
    "    \n",
    "    print est.summary()\n",
    "    \n",
    "    #print sm.regression.linear_model.OLSResults.aic()\n",
    "    #print OLSResults.aic(est.params)\n",
    "    \n",
    "    ## Create the 3d plot -- skip reading this\n",
    "    xx1, xx2 = np.meshgrid(np.linspace(min(x_i[:,1]), max(x_i[:,1]), 100), np.linspace(min(x_i[:,2]), max(x_i[:,2]), 100))\n",
    "    \n",
    "    Z = est.params[0] + est.params[1]*xx1 + est.params[2]*xx2\n",
    "    \n",
    "    # create matplotlib 3d axes\n",
    "    fig = plt.figure(figsize=(12, 8))\n",
    "    ax = Axes3D(fig, azim=-115, elev=15)\n",
    "\n",
    "    # plot hyperplane\n",
    "    surf = ax.plot_surface(xx1, xx2, Z, cmap=plt.cm.RdBu_r, alpha=0.6, linewidth=0)\n",
    "    \n",
    "    #plot the actual data\n",
    "    ax.scatter(x_i[:,1], x_i[:,2], y_i, color='black', alpha=1.0)\n",
    "\n",
    "    # set axis labels\n",
    "    ax.set_xlabel('|size-s_T|')\n",
    "    ax.set_ylabel('|orientation-o_T|')\n",
    "    ax.set_zlabel('hill^{-1}(p(i_p != i_r))')\n",
    "    \n",
    "    #for ii in xrange(0,360,1):\n",
    "    #    ax.view_init(elev=10., azim=ii)\n",
    "    #    plt.savefig(\"movie_v3_%d\"%ii+\".png\")\n",
    "        \n",
    "            #f = f+stop_trials*math.log(p_stop_model1[ind2,ind1]) + (total_trials-stop_trials)*math.log(1-p_stop_model1[ind2,ind1])     \n",
    "            #err = err + np.square(p_stop_model1[ind2,ind1]-p_stop_real[ind2,ind1])\n",
    "                \n",
    "    #return err\n",
    "    #return -f\n",
    "    \n",
    "    p_stop_model, p_stop_givenM_model, p_stop_givenB_model = p_stop_ProbabilityMatching([est.params[0], est.params[1], est.params[2], est.params[3]])\n",
    "    mle=0\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "                 #mle = mle + trials_stop[ind2,ind1]*log(p_stop_model[ind2,ind1])+(trials_tot[ind2,ind1]-trials_stop[ind2,ind1])*log(1-p_stop_model[ind2,ind1])\n",
    "                 mle = mle + trials_stop_m[ind2,ind1]*log(p_stop_givenM_model[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_stop_givenM_model[ind2,ind1])+ trials_stop_b[ind2,ind1]*log(p_stop_givenB_model[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_stop_givenB_model[ind2,ind1])\n",
    "                    \n",
    "    print LA.norm(p_stop_real-p_stop_model, 'fro')\n",
    "    print mle\n",
    "    \n",
    "    return est.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    #def fitting_function(x0):\n",
    "    x0=1\n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    \n",
    "    x_i = np.zeros((len(orientations)*len(sizes),3))\n",
    "    y_i = np.zeros((len(orientations)*len(sizes),1))\n",
    "    x_i_givenM = np.zeros((len(orientations)*len(sizes),3))\n",
    "    y_i_givenM = np.zeros((len(orientations)*len(sizes),1))\n",
    "    x_i_givenB = np.zeros((len(orientations)*len(sizes),3))\n",
    "    y_i_givenB = np.zeros((len(orientations)*len(sizes),1))\n",
    "    x_i2_givenM = np.zeros((len(orientations)*len(sizes),3))\n",
    "    y_i2_givenM = np.zeros((len(orientations)*len(sizes),1))\n",
    "    x_i2_givenB = np.zeros((len(orientations)*len(sizes),3))\n",
    "    y_i2_givenB = np.zeros((len(orientations)*len(sizes),1))\n",
    "    x_i_withPrior = np.zeros((2*len(orientations)*len(sizes),3))\n",
    "    y_i_withPrior = np.zeros((2*len(orientations)*len(sizes),1))\n",
    "    modeled_curve = np.zeros((len(orientations)*len(sizes),1))\n",
    "    modeled_curve_givenM1 = np.zeros((len(orientations)*len(sizes),1))\n",
    "    modeled_curve_givenM2 = np.zeros((len(orientations)*len(sizes),1))\n",
    "    hill_modeled_curve = []\n",
    "    hill_modeled_curve_givenM1 = []\n",
    "    hill_modeled_curve_givenM2 = []\n",
    "    hill_modeled_curve_Z = []\n",
    "    hill_modeled_curve_noHill_Z = np.zeros((100,100))\n",
    "    \n",
    "    ideal_params = [0.3679, 0.0129, 0.0024, 0.001314]\n",
    "        \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            x_i[ind2*len(sizes)+ind1,0] = math.fabs(orientations[ind2]-o_T)\n",
    "            x_i[ind2*len(sizes)+ind1,1] = math.fabs(sizes[ind1]-s_T)\n",
    "            x_i[ind2*len(sizes)+ind1,2] = s_max-sizes[ind1]\n",
    "            #x_i[ind2*len(sizes)+ind1,2] = 1/sizes[ind1]\n",
    "            y_i[ind2*len(sizes)+ind1] = x0*p_stop_real_toFit[ind2,ind1]\n",
    "            modeled_curve[ind2*len(sizes)+ind1] = ideal_params[0] + ideal_params[1]*x_i[ind2*len(sizes)+ind1,0] + ideal_params[2]*x_i[ind2*len(sizes)+ind1,1]+ ideal_params[3]*x_i[ind2*len(sizes)+ind1,2]\n",
    "            \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            #x_i_givenM[ind2*len(sizes)+ind1,0] = math.fabs(orientations[ind2]-o_T)\n",
    "            x_i_givenM[ind2*len(sizes)+ind1,0] = orientations[ind2]\n",
    "            #x_i_givenM[ind2*len(sizes)+ind1,1] = math.fabs(sizes[ind1]-s_T)\n",
    "            x_i_givenM[ind2*len(sizes)+ind1,1] = sizes[ind1]\n",
    "            #x_i_givenM[ind2*len(sizes)+ind1,2] = s_max-sizes[ind1]\n",
    "            #x_i_givenM[ind2*len(sizes)+ind1,2] = 1/sizes[ind1]\n",
    "            y_i_givenM[ind2*len(sizes)+ind1] = x0*p_stop_real_givenM_toFit[ind2,ind1]\n",
    "            #modeled_curve_givenM1[ind2*len(sizes)+ind1] = ideal_params[0] + ideal_params[1]*x_i_givenM[ind2*len(sizes)+ind1,0] + ideal_params[2]*x_i_givenM[ind2*len(sizes)+ind1,1]+ ideal_params[3]*x_i_givenM[ind2*len(sizes)+ind1,2]\n",
    "            \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            #x_i2_givenM[ind2*len(sizes)+ind1,0] = math.fabs(orientations[ind2]-o_T)\n",
    "            x_i2_givenM[ind2*len(sizes)+ind1,0] = orientations[ind2]\n",
    "            #x_i2_givenM[ind2*len(sizes)+ind1,1] = math.fabs(sizes[ind1]-s_T)\n",
    "            x_i2_givenM[ind2*len(sizes)+ind1,1] = sizes[ind1]\n",
    "            #x_i2_givenM[ind2*len(sizes)+ind1,2] = s_max-sizes[ind1]\n",
    "            #x_i2_givenM[ind2*len(sizes)+ind1,2] = 1/sizes[ind1]\n",
    "            y_i2_givenM[ind2*len(sizes)+ind1] = x0*p_stop_real_givenB_toFit[ind2,ind1]        \n",
    "            #modeled_curve_givenM2[ind2*len(sizes)+ind1] = ideal_params[0] + ideal_params[1]*x_i2_givenM[ind2*len(sizes)+ind1,0] + ideal_params[2]*x_i2_givenM[ind2*len(sizes)+ind1,1]+ ideal_params[3]*x_i2_givenM[ind2*len(sizes)+ind1,2]\n",
    "            \n",
    "    #for element in modeled_curve:\n",
    "    #    hill_modeled_curve.append(1/(1+1/element))\n",
    "\n",
    "    #for element in modeled_curve_givenM1:\n",
    "    #    hill_modeled_curve_givenM1.append(1/(1+1/element))\n",
    "    #    hill_modeled_curve_givenM2.append(1/(1+1/element))\n",
    "    \n",
    "    #x_i_withPrior = np.concatenate((x_i_givenM, x_i2_givenM), axis=0)\n",
    "    #y_i_withPrior = np.concatenate((y_i_givenM, y_i2_givenM), axis=0)\n",
    "    \n",
    "    #x_i = sm.add_constant(x_i)\n",
    "    #est = sm.OLS(y_i, x_i).fit()\n",
    "    #y_est = est.predict(x_i)\n",
    "    \n",
    "    #x_i_withPrior = sm.add_constant(x_i_withPrior)\n",
    "    #est = sm.OLS(y_i_withPrior, x_i_withPrior).fit()\n",
    "    #y_est = est.predict(x_i_withPrior)\n",
    "    \n",
    "    #print est.summary()\n",
    "\n",
    "    #print sm.regression.linear_model.OLSResults.aic()\n",
    "    #print OLSResults.aic(est.params)\n",
    "    \n",
    "    ## Create the 3d plot -- skip reading this\n",
    "    #xx1, xx2 = np.meshgrid(np.linspace(min(x_i[:,1]), max(x_i[:,1]), 100), np.linspace(min(x_i[:,2]), max(x_i[:,2]), 100))\n",
    "    #xx1, xx2 = np.meshgrid(np.linspace(min(x_i_withPrior[:,1]), max(x_i_withPrior[:,1]), 100), np.linspace(min(x_i_withPrior[:,2]), max(x_i_withPrior[:,2]), 100))\n",
    "    xx1, xx2 = np.meshgrid(np.linspace(min(orientations),max(orientations),100), np.linspace(min(sizes),max(sizes),100))\n",
    "    print \"xx1\", xx1, xx1.shape\n",
    "    print \"xx2\", xx2, xx2.shape\n",
    "    \n",
    "    for i in range(0,100):\n",
    "        for j in range(0,100):\n",
    "            hill_modeled_curve_noHill_Z[i,j] = ideal_params[0]+ideal_params[1]*math.fabs(xx1[i,j]-o_T)+ideal_params[2]*math.fabs(xx2[i,j]-s_T) + ideal_params[3]*(s_max - xx2[i,j])\n",
    "    \n",
    "    for element in hill_modeled_curve_noHill_Z:\n",
    "        hill_modeled_curve_Z.append(0.5/(1+1/element))\n",
    "    \n",
    "    #Z = est.params[0] + est.params[1]*xx1 + est.params[2]*xx2 \n",
    "    #+ est.params[3]*(min(x_i[:,3])+max(x_i[:,3]))/2\n",
    "    \n",
    "    # create matplotlib 3d axes\n",
    "    fig = plt.figure(figsize=(12, 8))\n",
    "    ax = Axes3D(fig, azim=-115, elev=15)\n",
    "\n",
    "    # plot hyperplane\n",
    "    surf = ax.plot_surface(xx1, xx2, hill_modeled_curve_Z, cmap=plt.cm.RdBu_r, alpha=0.6, linewidth=0)\n",
    "    \n",
    "    #plot the actual data\n",
    "    #ax.scatter(x_i[:,1], x_i[:,2], y_i, color='black', alpha=1.0)\n",
    "    ax.scatter(x_i_givenM[:,0], x_i_givenM[:,1], y_i_givenM, color='black', alpha=1.0)\n",
    "    \n",
    "    # set axis labels\n",
    "    ax.set_xlabel('orientations')\n",
    "    ax.set_ylabel('sizes')\n",
    "    ax.set_zlabel('uncertainty')\n",
    "    \n",
    "    #p_stop_model, p_stop_givenM_model, p_stop_givenB_model = p_stop_ProbabilityMatching([est.params[0], est.params[1], est.params[2], est.params[3]])\n",
    "    mle=0\n",
    "    \n",
    "    #for ind2,o in enumerate(orientations):\n",
    "    #    for ind1,s in enumerate(sizes):\n",
    "                 #mle = mle + trials_stop[ind2,ind1]*log(p_stop_model[ind2,ind1])+(trials_tot[ind2,ind1]-trials_stop[ind2,ind1])*log(1-p_stop_model[ind2,ind1])\n",
    "    #             mle = mle + trials_stop_m[ind2,ind1]*log(p_stop_givenM_model[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_stop_givenM_model[ind2,ind1])+ trials_stop_b[ind2,ind1]*log(p_stop_givenB_model[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_stop_givenB_model[ind2,ind1])\n",
    "                    \n",
    "    #print LA.norm(p_stop_real-p_stop_model, 'fro')\n",
    "    #print mle\n",
    "    \n",
    "    #for ii in xrange(0,360,1):\n",
    "    #    ax.view_init(elev=10., azim=ii)\n",
    "    #    plt.savefig(\"movie_v110_%d\"%ii+\".png\")\n",
    "    \n",
    "    #   return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitting_function(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 0.0 0.0 0.0\n",
      "x [ 0.  0.  0.  0.]\n",
      "3.93219922085 inf\n",
      "0.00025 0.0 0.0 0.0\n",
      "x [ 0.00025  0.       0.       0.     ]\n",
      "3.93053268988 18667.6710921\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-30-b273f4ac2571>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mx_initial\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mminimize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfitting_probMatching\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx_initial\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmethod\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'Nelder-Mead'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[0mp_stop_model\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_givenM_model\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_givenB_model\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_stop_ProbabilityMatching\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mresults\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresults\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresults\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresults\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\scipy\\optimize\\_minimize.pyc\u001b[0m in \u001b[0;36mminimize\u001b[1;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[0;32m    411\u001b[0m                       callback=callback, **options)\n\u001b[0;32m    412\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'nelder-mead'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 413\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    414\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'powell'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    415\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0m_minimize_powell\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\scipy\\optimize\\optimize.pyc\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[1;34m(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)\u001b[0m\n\u001b[0;32m    447\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    448\u001b[0m         \u001b[0msim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 449\u001b[1;33m         \u001b[0mf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    450\u001b[0m         \u001b[0mfsim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    451\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\scipy\\optimize\\optimize.pyc\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[1;34m(*wrapper_args)\u001b[0m\n\u001b[0;32m    279\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    280\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 281\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    282\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    283\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-29-f61f47b8c89b>\u001b[0m in \u001b[0;36mfitting_probMatching\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m         \u001b[0mp_stop_real\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_toFit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_toFit2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenM\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenB\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenM_toFit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenM_toFit2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenB_toFit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_real_givenB_toFit2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_stop_real_function\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m         \u001b[0mp_m_real\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_b_real\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrials_tot\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrials_stop\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrials_m\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrials_b\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrials_stop_m\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrials_stop_b\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_real\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m         \u001b[0mp_base\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_givenM_model\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_stop_givenB_model\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mp_stop_ProbabilityMatching\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-27-19226db2aa28>\u001b[0m in \u001b[0;36mp_real\u001b[1;34m()\u001b[0m\n\u001b[0;32m     20\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     21\u001b[0m             \u001b[0mtrials_stop_m\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mind2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mind1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSize\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m&\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mOri\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0mo\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m&\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mImage\u001b[0m\u001b[1;33m==\u001b[0m\u001b[1;34m'mushroom'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m&\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mchoice\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtotal_trial_ind\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 22\u001b[1;33m             \u001b[0mtrials_stop_b\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mind2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mind1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSize\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m&\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mOri\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0mo\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m&\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mImage\u001b[0m\u001b[1;33m==\u001b[0m\u001b[1;34m'bird'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m&\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mchoice\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtotal_trial_ind\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     23\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     24\u001b[0m             \u001b[0mp_m_real\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mind2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mind1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtrials_m\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mind2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mind1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mtrials_tot\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mind2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mind1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\frame.pyc\u001b[0m in \u001b[0;36m__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m   1676\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mSeries\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1677\u001b[0m             \u001b[1;31m# either boolean or fancy integer index\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1678\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_getitem_array\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1679\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mDataFrame\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1680\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_getitem_frame\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\frame.pyc\u001b[0m in \u001b[0;36m_getitem_array\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m   1718\u001b[0m             \u001b[0mkey\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_check_bool_indexer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1719\u001b[0m             \u001b[0mindexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnonzero\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1720\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtake\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mconvert\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1721\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1722\u001b[0m             \u001b[0mindexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mix\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_convert_to_indexer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\generic.pyc\u001b[0m in \u001b[0;36mtake\u001b[1;34m(self, indices, axis, convert, is_copy)\u001b[0m\n\u001b[0;32m   1233\u001b[0m         new_data = self._data.take(indices,\n\u001b[0;32m   1234\u001b[0m                                    \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_block_manager_axis\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1235\u001b[1;33m                                    convert=True, verify=True)\n\u001b[0m\u001b[0;32m   1236\u001b[0m         \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_constructor\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnew_data\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__finalize__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1237\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\internals.pyc\u001b[0m in \u001b[0;36mtake\u001b[1;34m(self, indexer, axis, verify, convert)\u001b[0m\n\u001b[0;32m   2932\u001b[0m         \u001b[0mnew_labels\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0maxes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtake\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2933\u001b[0m         return self.reindex_indexer(new_axis=new_labels, indexer=indexer,\n\u001b[1;32m-> 2934\u001b[1;33m                                     axis=axis, allow_dups=True)\n\u001b[0m\u001b[0;32m   2935\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2936\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mmerge\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlsuffix\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m''\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrsuffix\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m''\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\internals.pyc\u001b[0m in \u001b[0;36mreindex_indexer\u001b[1;34m(self, new_axis, indexer, axis, fill_value, allow_dups, copy)\u001b[0m\n\u001b[0;32m   2821\u001b[0m                                       fill_tuple=(fill_value if fill_value is not None else\n\u001b[0;32m   2822\u001b[0m                                                   blk.fill_value,))\n\u001b[1;32m-> 2823\u001b[1;33m                           for blk in self.blocks]\n\u001b[0m\u001b[0;32m   2824\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2825\u001b[0m         \u001b[0mnew_axes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0maxes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\internals.pyc\u001b[0m in \u001b[0;36mtake_nd\u001b[1;34m(self, indexer, axis, new_mgr_locs, fill_tuple)\u001b[0m\n\u001b[0;32m    750\u001b[0m             \u001b[0mfill_value\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfill_tuple\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    751\u001b[0m             new_values = com.take_nd(self.get_values(), indexer, axis=axis,\n\u001b[1;32m--> 752\u001b[1;33m                                      allow_fill=True, fill_value=fill_value)\n\u001b[0m\u001b[0;32m    753\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    754\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mnew_mgr_locs\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda\\lib\\site-packages\\pandas\\core\\common.pyc\u001b[0m in \u001b[0;36mtake_nd\u001b[1;34m(arr, indexer, axis, out, fill_value, mask_info, allow_fill)\u001b[0m\n\u001b[0;32m    651\u001b[0m         \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfill_value\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    652\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 653\u001b[1;33m         \u001b[0mindexer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_ensure_int64\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    654\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mallow_fill\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    655\u001b[0m             \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfill_value\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0marr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "    # fitting in the case of probability matching\n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_real_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    x_initial = [0,0,0,0]\n",
    "\n",
    "    results = minimize(fitting_probMatching, x_initial, method = 'Nelder-Mead')\n",
    "    \n",
    "    p_stop_model, p_stop_givenM_model, p_stop_givenB_model = p_stop_ProbabilityMatching([results.x[0], results.x[1], results.x[2], results.x[3]])\n",
    "    mle=0\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "                 #mle = mle + trials_stop[ind2,ind1]*log(p_stop_model[ind2,ind1])+(trials_tot[ind2,ind1]-trials_stop[ind2,ind1])*log(1-p_stop_model[ind2,ind1])\n",
    "                 mle = mle + trials_stop_m[ind2,ind1]*log(p_stop_givenM_model[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_stop_givenM_model[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_stop_givenB_model[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_stop_givenB_model[ind2,ind1])\n",
    "                    \n",
    "    #print LA.norm(p_stop_real-p_stop_model, 'fro')\n",
    "    print LA.norm(p_stop_real - p_m_real*p_stop_givenM_model - p_b_real*p_stop_givenB_model, 'fro')\n",
    "    print mle\n",
    "    print \"real\", p_stop_real\n",
    "    print \"model\", p_stop_model\n",
    "    print \"1\", p_stop_givenM_model, p_stop_real_givenM\n",
    "    print \"2\", p_stop_givenB_model, p_stop_real_givenB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fitting_probMatching(x):\n",
    "        \n",
    "        p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_real_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "        p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "        \n",
    "        p_base, p_stop_givenM_model, p_stop_givenB_model = p_stop_ProbabilityMatching(x)\n",
    "        mle = 0\n",
    "        \n",
    "        print \"x\", x\n",
    "        \n",
    "        for ind2,o in enumerate(orientations):\n",
    "            for ind1,s in enumerate(sizes):    \n",
    "                mle = mle + trials_stop_m[ind2,ind1]*log(p_stop_givenM_model[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_stop_givenM_model[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_stop_givenB_model[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_stop_givenB_model[ind2,ind1])\n",
    "                \n",
    "        #norm = LA.norm(p_stop_real - p_base, 'fro')\n",
    "        norm = LA.norm(p_stop_real_givenM-p_stop_givenM_model, 'fro') + LA.norm(p_stop_real_givenB-p_stop_givenB_model, 'fro')\n",
    "        \n",
    "        print norm, -mle\n",
    "        return -mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# compute average T1, T0 for each combination (s,o)\n",
    "def pause_and_laptimes():\n",
    "\n",
    "    df['laptime'] = np.diff(np.concatenate([[0], df.time]))\n",
    "    \n",
    "    matrix = df.values\n",
    "\n",
    "    a,b = matrix.shape\n",
    "    matrix2=matrix\n",
    "\n",
    "    for i in range(0, a):\n",
    "        if matrix2[i,18]<0:\n",
    "            matrix2[i,18]=matrix2[i,5]\n",
    "    \n",
    "    dh = pd.DataFrame(matrix2, columns = list(df.columns.values))       \n",
    "        \n",
    "    print \"dh\", dh[dh.laptime<0].laptime\n",
    "    \n",
    "    pausetimes = np.zeros((len(orientations),len(sizes)))\n",
    "    pausetimes_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "    pausetimes_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    laptimes = np.zeros((len(orientations),len(sizes)))\n",
    "    laptimes_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "    laptimes_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    p_Correct_GivenStop = np.zeros((len(orientations),len(sizes)))\n",
    "\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            \n",
    "                pausetimes[ind2,ind1] = dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='mushroom')].pausetime.mean()*p_m_real[ind2,ind1]+dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='bird')].pausetime.mean()*p_b_real[ind2, ind1]\n",
    "                laptimes[ind2,ind1] = dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='mushroom')].laptime.mean()*p_m_real[ind2, ind1]+dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='bird')].laptime.mean()*p_b_real[ind2, ind1]\n",
    "\n",
    "                pausetimes_givenM[ind2,ind1] = dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='mushroom')].pausetime.mean()\n",
    "                pausetimes_givenB[ind2,ind1] = dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='bird')].pausetime.mean()\n",
    "                \n",
    "                laptimes_givenM[ind2,ind1] = dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='mushroom')].laptime.mean()\n",
    "                laptimes_givenB[ind2,ind1] = dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='bird')].laptime.mean()\n",
    "                \n",
    "                p_Correct_GivenStop[ind2,ind1] = len(dh[(dh.Size==s)&(dh.Ori==o)&(dh.Image=='mushroom')&(dh.choice==True)].index)/(float)(len(dh[(dh.Size==s)&(dh.Ori==o)&(dh.choice==True)].index))\n",
    "   \n",
    "    return pausetimes, laptimes, pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#show heatmaps of p(stop) = p(m|m)*p(m) + p(m|b)*p(b) and p(m|m) + p(m|b)\n",
    "\n",
    "sum_term_avg = 0\n",
    "sum_term = np.zeros((len(orientations),len(sizes)))\n",
    "sum_term2 = np.zeros((len(orientations),len(sizes)))\n",
    "sum_term3 = np.zeros((len(orientations),len(sizes)))\n",
    "\n",
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "\n",
    "for ind2,o in enumerate(orientations):\n",
    "    for ind1,s in enumerate(sizes):\n",
    "    \n",
    "        #testing how much p(m|m)+p(b/m)<>1\n",
    "    \n",
    "        i1 = dd[(dd.Size==s)&(dd.Ori==o)].yp.index[0]\n",
    "        i2 = dd[(dd.Size==s)&(dd.Ori==o)].yp.index[1]\n",
    "        sum_term_avg = sum_term_avg + dd[(dd.Size==s)&(dd.Ori==o)].yp[i1]+dd[(dd.Size==s)&(dd.Ori==o)].yp[i2]\n",
    "        sum_term[ind2,ind1] = dd[(dd.Size==s)&(dd.Ori==o)].yp[i1]+dd[(dd.Size==s)&(dd.Ori==o)].yp[i2]\n",
    "        sum_term2[ind2,ind1] = dd[(dd.Size==s)&(dd.Ori==o)].yp[i1]*p_m_real[ind2,ind1]+dd[(dd.Size==s)&(dd.Ori==o)].yp[i2]*p_b_real[ind2,ind1]\n",
    "        sum_term3[ind2,ind1] = (1-dd[(dd.Size==s)&(dd.Ori==o)].yp[i1])*p_m_real[ind2,ind1] + dd[(dd.Size==s)&(dd.Ori==o)].yp[i2]*p_b_real[ind2,ind1]      \n",
    "            \n",
    "sum_term_avg = sum_term_avg/(len(sizes)*len(orientations))       \n",
    "\n",
    "a,b = sum_term.shape\n",
    "\n",
    "plt.imshow(sum_term, interpolation='none')\n",
    "\n",
    "plt.xticks(range(b), sizes)\n",
    "plt.yticks(range(a), orientations)\n",
    "plt.suptitle('Heat map showing p(m|m)+p(m|b)')\n",
    "\n",
    "plt.xlabel('Sizes')\n",
    "plt.ylabel('Orientations')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(sum_term2, interpolation='none')\n",
    "\n",
    "plt.xticks(range(b), sizes)\n",
    "plt.yticks(range(a), orientations)\n",
    "plt.suptitle('Heat map showing p(m|m)*p(m)+p(m|b)*p(b) = p(stop)')\n",
    "\n",
    "plt.xlabel('Sizes')\n",
    "plt.ylabel('Orientations')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(sum_term3, interpolation='none')\n",
    "\n",
    "plt.xticks(range(b), sizes)\n",
    "plt.yticks(range(a), orientations)\n",
    "plt.suptitle('Uncertainty = p(m|b)*p(b)+p(b|m)*p(m)')\n",
    "\n",
    "plt.xlabel('Sizes')\n",
    "plt.ylabel('Orientations')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_real_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "\n",
    "plt.imshow(1-p_stop_real_givenM, interpolation='none')\n",
    "\n",
    "plt.xticks(range(b), sizes)\n",
    "plt.yticks(range(a), orientations)\n",
    "plt.suptitle('Heat map showing p(b|m), if probability matching')\n",
    "\n",
    "plt.xlabel('Sizes')\n",
    "plt.ylabel('Orientations')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(p_stop_real_givenB, interpolation='none')\n",
    "\n",
    "plt.xticks(range(b), sizes)\n",
    "plt.yticks(range(a), orientations)\n",
    "plt.suptitle('Heat map showing p(m|b), if probability matching')\n",
    "\n",
    "plt.xlabel('Sizes')\n",
    "plt.ylabel('Orientations')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def find_mu_sigma_delta():\n",
    "\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    \n",
    "    #parameters = fitting_function(1)  \n",
    "    #x = [parameters[0], parameters[1], parameters[2], parameters[3]]\n",
    "    x = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "    prob_misID = probability_misidentify(x)\n",
    "    percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "\n",
    "    sigma = 1\n",
    "    tau =0.1\n",
    "\n",
    "    delta = sigma*math.sqrt(tau)\n",
    "    mu = (2*percept-1)*sigma/(math.sqrt(tau))\n",
    "\n",
    "    mu_givenM = (2*percept_givenM-1)*sigma/(math.sqrt(tau))\n",
    "    mu_givenB = (2*percept_givenB-1)*sigma/(math.sqrt(tau))\n",
    "    \n",
    "    return percept, mu, mu_givenM, mu_givenB, sigma, delta, tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#B-2(cont)-1: contionuous case, end task at time M (get Gaussian) and then threshold\n",
    "\n",
    "percept, mu, mu_givenM, mu_givenB, sigma, delta, tau = find_mu_sigma_delta()\n",
    "p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "p_stop_AND_correct = p_stop_AND_correct_function()\n",
    "pausetimes, laptimes, pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop = pause_and_laptimes()\n",
    "\n",
    "p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "\n",
    "start_N = delta\n",
    "end_N = 31*delta\n",
    "no_N = (end_N/delta-start_N/delta)\n",
    "array_N = np.linspace(start_N,end_N,no_N,endpoint=False)\n",
    "\n",
    "print start_N, end_N, no_N, array_N\n",
    "\n",
    "start_M = tau\n",
    "end_M = 51*tau\n",
    "no_M = (end_M-start_M)/tau\n",
    "array_M = np.linspace(start_M,end_M,no_M, endpoint =False)\n",
    "\n",
    "p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenM2 = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenB2 = np.zeros((len(orientations),len(sizes)))\n",
    "\n",
    "norm = np.zeros((no_N, no_M))\n",
    "mle = np.zeros((no_N, no_M))\n",
    "total_len = len(df)\n",
    "\n",
    "avg_reward = np.zeros((no_N,no_M))\n",
    "avg_reward_numerator = np.zeros((no_N,no_M))\n",
    "avg_reward_denominator = np.zeros((no_N,no_M))\n",
    "avg_reward2 = np.zeros((no_N,no_M))\n",
    "\n",
    "avg_reward_modified = np.zeros((no_N,no_M))\n",
    "avg_reward_numerator_modified = np.zeros((no_N,no_M))\n",
    "avg_reward_denominator_modified = np.zeros((no_N,no_M))\n",
    "avg_reward2_modified = np.zeros((no_N,no_M))\n",
    "\n",
    "for i,N in enumerate(array_N):\n",
    "    for j,M in enumerate(array_M):\n",
    "        \n",
    "        for ind2,o in enumerate(orientations):\n",
    "            for ind1,s in enumerate(sizes):\n",
    "                \n",
    "                 mydist = stat.norm(mu[ind2,ind1]*M, sigma*M)\n",
    "                 p_modeled_stopping[ind2,ind1] = 1-mydist.cdf(N)\n",
    "                 \n",
    "                 mydist_givenM = stat.norm(mu_givenM[ind2,ind1]*M, sigma*M)\n",
    "                 mydist_givenB = stat.norm(mu_givenB[ind2,ind1]*M, sigma*M)\n",
    "                 print mydist_givenM\n",
    "                    \n",
    "                 p_modeled_stopping_givenM[ind2,ind1] = 1-mydist_givenM.cdf(N)   \n",
    "                 p_modeled_stopping_givenB[ind2,ind1] = 1-mydist_givenB.cdf(N)   \n",
    "                 \n",
    "                 p_modeled_stopping_givenM2[ind2,ind1] = 1-stat.norm.cdf(N, mu_givenM[ind2,ind1]*M, sigma*M)   \n",
    "                 p_modeled_stopping_givenB2[ind2,ind1] = 1-stat.norm.cdf(N, mu_givenB[ind2,ind1]*M, sigma*M)  \n",
    "                    \n",
    "                 p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                    \n",
    "                 if o<o_T or s<s_T:\n",
    "                        \n",
    "                     avg_reward[i,j] = avg_reward[i,j]+ (p_stop_AND_correct[ind2,ind1]/(laptimes[ind2,ind1]+p_modeled_stopping[ind2,ind1]*pausetimes[ind2,ind1]))*trials_tot[ind2,ind1]/total_len \n",
    "                     avg_reward_numerator[i,j] = avg_reward_numerator[i,j] + p_stop_AND_correct[ind2,ind1]*trials_tot[ind2,ind1]/total_len\n",
    "                     avg_reward_denominator[i,j] = avg_reward_denominator[i,j] + (laptimes[ind2,ind1]+p_modeled_stopping[ind2,ind1]*pausetimes[ind2,ind1])*trials_tot[ind2,ind1]/total_len   \n",
    "                     avg_reward_modified[i,j] = avg_reward_modified[i,j] + ((p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1])/(laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1]))*trials_tot[ind2,ind1]/total_len           \n",
    "                        \n",
    "                 #mle[i,j] = mle[i,j] + trials_stop[ind2,ind1]*log(p_modeled_stopping[ind2,ind1])+(trials_tot[ind2,ind1]-trials_stop[ind2,ind1])*log(1-p_modeled_stopping[ind2,ind1])\n",
    "                 mle[i,j] = mle[i,j] + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])\n",
    "                        \n",
    "        if (i==0) & (j==0):\n",
    "            print '1', p_modeled_stopping2, p_stop_real\n",
    "            print '2', p_modeled_stopping_givenM, p_stop_real_givenM\n",
    "            print '3', p_modeled_stopping_givenB, p_stop_real_givenB\n",
    "            \n",
    "        norm[i,j] = LA.norm(p_stop_real - p_modeled_stopping2, 'fro')\n",
    "        avg_reward2[i,j] = avg_reward_numerator[i,j]/avg_reward_denominator[i,j]\n",
    "   \n",
    "minimum = 10000\n",
    "maximum = -10000\n",
    "maximum_reward = -10000\n",
    "maximum_reward2 = -10000\n",
    "maximum_reward_modified = -10000\n",
    "\n",
    "for i,N in enumerate(array_N):\n",
    "    for j,M in enumerate(array_M):\n",
    "        \n",
    "        if norm[i,j]<minimum:\n",
    "            minimum = norm[i,j]\n",
    "            i_min = i\n",
    "            j_min = j\n",
    "        \n",
    "        if mle[i,j]>maximum:\n",
    "            maximum = mle[i,j]\n",
    "            i_max_mle = i\n",
    "            j_max_mle = j\n",
    "            \n",
    "        if avg_reward[i,j]>maximum_reward:\n",
    "            maximum_reward = avg_reward[i,j]\n",
    "            N_rew = i\n",
    "            M_rew = j\n",
    "            \n",
    "        if avg_reward2[i,j]>maximum_reward2:\n",
    "            maximum_reward2 = avg_reward2[i,j]\n",
    "            N_rew2 = i\n",
    "            M_rew2 = j\n",
    "            \n",
    "        if avg_reward_modified[i,j]>maximum_reward_modified:\n",
    "            maximum_reward_modified = avg_reward_modified[i,j]\n",
    "            N_rew_mod = i\n",
    "            M_rew_mod = j     \n",
    "                \n",
    "print minimum, i_min, j_min, array_N[i_min], array_M[j_min]\n",
    "print maximum, i_max_mle, j_max_mle, array_N[i_max_mle], array_M[j_max_mle]\n",
    "\n",
    "print \"norm\", norm[0,:]\n",
    "print \"mle\", mle\n",
    "\n",
    "a,b = norm.shape\n",
    "\n",
    "aic = 4 - mle[i_max_mle,j_max_mle]\n",
    "print \"aic =\", aic\n",
    "\n",
    "print \"avg reward rate for best fit is\", avg_reward[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward, \"at N,M\", N_rew, M_rew\n",
    "\n",
    "print \"avg reward rate for best fit is II\", avg_reward2[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward2[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward2, \"at N,M\", N_rew2, M_rew2\n",
    "\n",
    "print \"avg reward rate for best fit is II\", avg_reward_modified[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward_modified[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward_modified, \"at N,M\", N_rew_mod, M_rew_mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " #B - 1(discrete)1-+1,0 step, with complicated formula   \n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    percept, mu, mu_givenM, mu_givenB, sigma, delta, tau = find_mu_sigma_delta()\n",
    "    p_stop_AND_correct = p_stop_AND_correct_function()\n",
    "    pausetimes, laptimes, pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop = pause_and_laptimes()\n",
    "    \n",
    "    #x = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "    x = [1.62881729, 0.03056224, 0.00569523, 0.00306576]\n",
    "    \n",
    "    prob_misID = probability_misidentify(x)\n",
    "    percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "    \n",
    "    start_N = delta\n",
    "    end_N = delta*10\n",
    "    no_N = (end_N-start_N+delta)/delta\n",
    "    array_N = np.linspace(start_N,end_N,no_N,endpoint=False)\n",
    "    exp1N = (int)(start_N/delta)\n",
    "    exp2N = (int)(end_N/delta)\n",
    "    \n",
    "    start_M = tau\n",
    "    end_M = tau*15\n",
    "    no_M = (end_M-start_M+tau)/tau\n",
    "    array_M = np.linspace(start_M,end_M,no_M, endpoint =False)\n",
    "    exp1M = (int)(start_M/tau)\n",
    "    exp2M = (int)(end_M/tau)\n",
    "    \n",
    "    norm = np.zeros((no_N, no_M))\n",
    "    mle = np.zeros((no_N, no_M))\n",
    "    \n",
    "    total_len = len(df)\n",
    "    avg_reward = np.zeros((no_N,no_M))\n",
    "    avg_reward_numerator = np.zeros((no_N,no_M))\n",
    "    avg_reward_denominator = np.zeros((no_N,no_M))\n",
    "    \n",
    "    for exp1 in range(exp1N, exp2N):\n",
    "        for exp2 in range(exp1,exp2M):       \n",
    "            \n",
    "            p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "            p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "            p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "            p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))\n",
    "            \n",
    "            for ind2,o in enumerate(orientations):\n",
    "                for ind1,s in enumerate(sizes):\n",
    "\n",
    "                    for k in range(0,exp2-exp1+1):\n",
    "                        p_modeled_stopping[ind2,ind1] = p_modeled_stopping[ind2,ind1] + ((percept[ind2,ind1])**exp1)*((1-percept[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        p_modeled_stopping_givenM[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1] + ((percept_givenM[ind2,ind1])**exp1)*((1-percept_givenM[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        p_modeled_stopping_givenB[ind2,ind1] = p_modeled_stopping_givenB[ind2,ind1] + ((percept_givenB[ind2,ind1])**exp1)*((1-percept_givenB[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        \n",
    "                    p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                        \n",
    "                    #mle[exp1,exp2] = mle[exp1,exp2] + trials_stop[ind2,ind1]*log(p_modeled_stopping2[ind2,ind1])+(trials_tot[ind2,ind1]-trials_stop[ind2,ind1])*log(1-p_modeled_stopping2[ind2,ind1])            \n",
    "                    mle[exp1,exp2] = mle[exp1,exp2] + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])            \n",
    "                    \n",
    "                    if isnan(mle[exp1,exp2]):\n",
    "                        mle[exp1,exp2] = -np.inf\n",
    "                        \n",
    "                    if o<o_T or s<s_T:\n",
    "                        avg_reward_numerator[exp1,exp2] = avg_reward_numerator[exp1,exp2] + p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1]*trials_tot[ind2,ind1]/total_len\n",
    "                    avg_reward_denominator[exp1,exp2] = avg_reward_denominator[exp1,exp2] + (laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1])*trials_tot[ind2,ind1]/total_len\n",
    "            \n",
    "            norm[exp1,exp2] = math.sqrt(LA.norm(p_stop_real_givenM - p_modeled_stopping_givenM, 'fro')**2 + LA.norm(p_stop_real_givenM - p_modeled_stopping_givenM,'fro')**2)\n",
    "            avg_reward[exp1,exp2] = avg_reward_numerator[exp1,exp2]/avg_reward_denominator[exp1,exp2]\n",
    "            \n",
    "    i_min=0\n",
    "    j_min=0\n",
    "    \n",
    "    minimum = np.inf\n",
    "    maximum =-np.inf\n",
    "    maximum_reward = -np.inf\n",
    "    \n",
    "    for exp1 in range(exp1N, exp2N):\n",
    "        for exp2 in range(exp1,exp2M):  \n",
    "        \n",
    "            if norm[exp1,exp2]<minimum:\n",
    "                minimum = norm[exp1,exp2]\n",
    "                i_min = exp1\n",
    "                j_min = exp2     \n",
    "            \n",
    "            if mle[exp1,exp2]>maximum:\n",
    "                maximum = mle[exp1,exp2]\n",
    "                i_max_mle = exp1\n",
    "                j_max_mle = exp2\n",
    "                \n",
    "            if avg_reward[exp1,exp2]>maximum_reward:\n",
    "                maximum_reward = avg_reward[exp1,exp2]\n",
    "                N_rew = exp1\n",
    "                M_rew = exp2 \n",
    "                \n",
    "    print minimum, i_min, j_min\n",
    "    print minimum, i_min, j_min, array_N[i_min], array_M[j_min]\n",
    "    print maximum, i_max_mle, j_max_mle, array_N[i_max_mle], array_M[j_max_mle]    \n",
    "    \n",
    "    a,b = norm.shape\n",
    "    print \"a,b\", a,b\n",
    "    print \"norm\", norm[1,1]\n",
    "    print \"winning norm\", minimum\n",
    "    \n",
    "    aic = 4 - 2*mle[i_max_mle,j_max_mle]\n",
    "    print \"aic =\", aic\n",
    "    print \"log mle =\", mle[i_max_mle,j_max_mle]\n",
    "    \n",
    "    print \"avg reward rate for best fit is\", avg_reward[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward[i_min,j_min]\n",
    "    print \"maximum reward rate is\",  maximum_reward, \"at N,M\", N_rew, M_rew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#B - 1(discrete)-2 - with start point & threshold\n",
    "percept, mu, mu_givenM, mu_givenB, sigma, delta, tau = find_mu_sigma_delta()\n",
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "      \n",
    "x = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "prob_misID = probability_misidentify(x)\n",
    "\n",
    "p_stop_AND_correct = p_stop_AND_correct_function()\n",
    "pausetimes, laptimes,  pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop = pause_and_laptimes()\n",
    "\n",
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "    \n",
    "p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "\n",
    "start_N = delta\n",
    "end_N = delta*30\n",
    "no_N = (end_N-start_N+delta)/delta\n",
    "array_N = np.linspace(start_N,end_N,no_N,endpoint=False)\n",
    "\n",
    "norm = np.zeros((no_N, 2*no_N-1))\n",
    "mle = np.zeros((no_N, 2*no_N-1))\n",
    "mle2 = np.zeros((no_N, 2*no_N-1))\n",
    "total_len = len(df)\n",
    "\n",
    "avg_reward = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward_numerator = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward_denominator = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward2 = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward_modified = np.zeros((no_N,no_M))\n",
    "avg_reward_numerator_modified = np.zeros((no_N,no_M))\n",
    "avg_reward_denominator_modified = np.zeros((no_N,no_M))\n",
    "avg_reward2_modified = np.zeros((no_N,no_M))\n",
    "\n",
    "for i,N in enumerate(array_N):           \n",
    "    \n",
    "    exp1 = (int)(math.floor(N/delta))\n",
    "    for j,z in enumerate(range(-exp1+1, exp1)):\n",
    "        \n",
    "        if i==0:\n",
    "            break\n",
    "        \n",
    "        exp2 = z\n",
    "        \n",
    "        for ind2,o in enumerate(orientations):\n",
    "            for ind1,s in enumerate(sizes):\n",
    "\n",
    "                if percept[ind2,ind1]!=0.5:\n",
    "                    \n",
    "                    factor = percept[ind2, ind1]/(1-percept[ind2,ind1])\n",
    "                    factor_givenM = percept_givenM[ind2,ind1]/(1-percept_givenM[ind2,ind1])\n",
    "                    factor_givenB = percept_givenB[ind2,ind1]/(1-percept_givenB[ind2,ind1])\n",
    "                    \n",
    "                    p_modeled_stopping[ind2,ind1] = (factor**(2*exp1)- factor**(exp1-exp2))/(factor**(2*exp1)-1)\n",
    "                    p_modeled_stopping_givenM[ind2,ind1] = (factor_givenM**(2*exp1)- factor_givenM**(exp1-exp2))/(factor_givenM**(2*exp1)-1)\n",
    "                    p_modeled_stopping_givenB[ind2,ind1] = (factor_givenB**(2*exp1)- factor_givenB**(exp1-exp2))/(factor_givenB**(2*exp1)-1)\n",
    "                    p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                    \n",
    "                else:\n",
    "                    \n",
    "                    p_modeled_stopping[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    p_modeled_stopping_givenM[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    p_modeled_stopping_givenB[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                \n",
    "                mle[i,j] = mle[i,j] + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])                        \n",
    "                mle2[i,j] = mle2[i,j] + trials_stop_m[ind2,ind1]*log(percept_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-percept_givenM[ind2,ind1]) + trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])            \n",
    "                \n",
    "                if o<o_T or s<s_T:\n",
    "                     avg_reward_modified[i,j] = avg_reward_modified[i,j]+ (p_modeled_stopping_givenM/(laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1]))*trials_tot[ind2,ind1]/total_len\n",
    "                     avg_reward_numerator[i,j] = avg_reward_numerator[i,j] + p_modeled_stopping_givenM[ind2,ind1]*trials_tot[ind2,ind1]/total_len\n",
    "                     avg_reward_denominator[i,j] = avg_reward_denominator[i,j] + (laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1])*trials_tot[ind2,ind1]/total_len    \n",
    "                        \n",
    "        norm[i,j] = LA.norm(p_stop_real - p_modeled_stopping2, 'fro')             \n",
    "        avg_reward2[i,j] = avg_reward_numerator[i,j]/avg_reward_denominator[i,j]    \n",
    "#print norm\n",
    "minimum = 100000\n",
    "maximum = -100000\n",
    "maximum_reward=-10000\n",
    "maximum_reward2 = -10000\n",
    "\n",
    "for i in range(1,(int)(no_N)):\n",
    "    for j in range(0,2*i-1):\n",
    "        \n",
    "        if norm[i,j]<minimum:\n",
    "            minimum = norm[i,j]\n",
    "            i_min = i\n",
    "            j_min = j-i+1  \n",
    "            \n",
    "        if mle[i,j]>maximum:\n",
    "            maximum = mle[i,j]\n",
    "            i_max_mle = i\n",
    "            j_max_mle = j-i+1\n",
    "            \n",
    "        if avg_reward_modified[i,j]>maximum_reward:\n",
    "            maximum_reward = avg_reward[i,j]\n",
    "            N_rew = i\n",
    "            M_rew = j  \n",
    "            \n",
    "        if avg_reward2[i,j]>maximum_reward2:\n",
    "            maximum_reward2 = avg_reward2[i,j]\n",
    "            N_rew2 = i\n",
    "            M_rew2 = j    \n",
    "                \n",
    "print minimum, i_min, j_min, array_N[i_min]\n",
    "print maximum, i_max_mle, j_max_mle, array_N[i_max_mle]  \n",
    "\n",
    "#print \"mle2\", mle2\n",
    "aic = 4 - mle[i_max_mle,j_max_mle+i_max_mle-1]\n",
    "print \"aic =\", aic\n",
    "\n",
    "print \"avg reward rate for best fit is\", avg_reward_modified[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward_modified[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward, \"at N,M\", N_rew, M_rew\n",
    "\n",
    "print \"avg reward rate for best fit is II\", avg_reward2[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward2[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward2, \"at N,M\", N_rew2, M_rew2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#B-2(cont. case)2- looking for threshold and start point\n",
    "\n",
    "percept, mu, mu_givenM, mu_givenB, sigma, delta, tau = find_mu_sigma_delta()\n",
    "p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "p_stop_AND_correct = p_stop_AND_correct_function()\n",
    "pausetimes, laptimes,  pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop = pause_and_laptimes()\n",
    "\n",
    "total_len = len(df)\n",
    "avg_reward = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward_numerator = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward_denominator = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward2 = np.zeros((no_N,2*no_N-1))\n",
    "avg_reward_modified = np.zeros((no_N,no_M))\n",
    "avg_reward_numerator_modified = np.zeros((no_N,no_M))\n",
    "avg_reward_denominator_modified = np.zeros((no_N,no_M))\n",
    "avg_reward2_modified = np.zeros((no_N,no_M))\n",
    "\n",
    "p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))\n",
    "\n",
    "start_N = delta\n",
    "end_N = 30*delta\n",
    "no_N = (end_N/delta-start_N/delta)+1\n",
    "array_N = np.linspace(start_N,end_N,no_N,endpoint=False)\n",
    "exp1N = (int)(start_N/delta)\n",
    "exp2N = (int)(end_N/delta)\n",
    "    \n",
    "norm = np.zeros((no_N, 2*no_N-1))\n",
    "mle = np.zeros((no_N, 2*no_N-1))\n",
    "    \n",
    "for exp1 in range(exp1N, exp2N):\n",
    "    \n",
    "    for j,exp2 in enumerate(range(-exp1+1,exp1)):\n",
    "\n",
    "        for ind2,o in enumerate(orientations):\n",
    "            for ind1,s in enumerate(sizes):\n",
    "                 \n",
    "                 p_modeled_stopping[ind2,ind1] = (math.exp(-2*mu_givenM[ind2,ind1]*(exp1+exp2)/(sigma**2))-1)/(math.exp(-2*mu_givenM[ind2,ind1]*2*exp1/(sigma**2))-1)\n",
    "                 p_modeled_stopping_givenM[ind2,ind1] = (math.exp(-2*mu_givenM[ind2,ind1]*(exp1+exp2)/(sigma**2))-1)/(math.exp(-2*mu_givenM[ind2,ind1]*2*exp1/(sigma**2))-1)\n",
    "                 p_modeled_stopping_givenB[ind2,ind1] = (math.exp(-2*mu_givenB[ind2,ind1]*(exp1+exp2)/(sigma**2))-1)/(math.exp(-2*mu_givenB[ind2,ind1]*2*exp1/(sigma**2))-1)\n",
    "                 p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                    \n",
    "                 mle[exp1,j+1] = mle[exp1,j+1] + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])                        \n",
    "                    \n",
    "                 if o<o_T or s<s_T:   \n",
    "                     avg_reward_modified[exp1,j+1] = avg_reward[exp1,j+1]+ ((p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1])/(laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1]))*trials_tot[ind2,ind1]/total_len\n",
    "                     avg_reward_numerator[exp1,j+1] = avg_reward_numerator[exp1,j+1] + (p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1])*trials_tot[ind2,ind1]/total_len\n",
    "                     avg_reward_denominator[exp1,j+1] = avg_reward_denominator[exp1,j+1] + (laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1])*trials_tot[ind2,ind1]/total_len           \n",
    "                        \n",
    "        norm[exp1,j+1] = LA.norm(p_stop_real - p_modeled_stopping2, 'fro')           \n",
    "        avg_reward2[exp1,j+1] = avg_reward_numerator[exp1,j+1]/avg_reward_denominator[exp1,j+1]       \n",
    "        \n",
    "minimum = 10000\n",
    "maximum = -10000\n",
    "maximum_reward = -10000\n",
    "maximum_reward2 = -10000\n",
    "\n",
    "for i in range(exp1N, exp2N):\n",
    "    for j in range(1,2*i):       \n",
    "\n",
    "        if norm[i,j]<minimum:\n",
    "            minimum = norm[i,j]\n",
    "            i_min = i\n",
    "            j_min = j\n",
    "        \n",
    "        if mle[i,j]>maximum:\n",
    "            maximum = mle[i,j]\n",
    "            i_max_mle = i\n",
    "            j_max_mle = j\n",
    "            \n",
    "        if avg_reward_modified[i,j]>maximum_reward:\n",
    "            maximum_reward = avg_reward_modified[i,j]\n",
    "            N_rew = i\n",
    "            M_rew = j     \n",
    "            \n",
    "        if avg_reward2[i,j]>maximum_reward2:\n",
    "            maximum_reward2 = avg_reward2[i,j]\n",
    "            N_rew2 = i\n",
    "            M_rew2 = j\n",
    "            \n",
    "print minimum, i_min, j_min\n",
    "print maximum, i_max_mle, j_max_mle\n",
    "\n",
    "print \"norm\", norm\n",
    "print \"mle\", mle\n",
    "\n",
    "a,b = norm.shape\n",
    "\n",
    "aic = 4 - 2*mle[i_max_mle,j_max_mle]\n",
    "print \"aic =\", aic\n",
    "\n",
    "print \"avg reward rate for best fit is\", avg_reward_modified[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward_modified[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward, \"at N,M\", N_rew, M_rew\n",
    "\n",
    "print \"avg reward rate for best fit is II\", avg_reward2[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward2[i_min,j_min]\n",
    "print \"maximum reward rate is\",  maximum_reward2, \"at N,M\", N_rew2, M_rew2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fitting_accumulator2(x):    \n",
    "    \n",
    "    mle = 0\n",
    "    p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    prob_misID = probability_misidentify(x)\n",
    "    percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "    \n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "            for ind1,s in enumerate(sizes):\n",
    "\n",
    "                if percept[ind2,ind1]!=0.5:\n",
    "                     factor = percept[ind2, ind1]/(1-percept[ind2,ind1])\n",
    "                     p_modeled_stopping[ind2,ind1] = (factor**(2*exp1)- factor**(exp1+exp2))/(factor**(2*exp1)-1)\n",
    "                else:\n",
    "                     p_modeled_stopping[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                        \n",
    "                if percept_givenM[ind2,ind1]!=0.5:        \n",
    "                     factor_givenM = percept_givenM[ind2,ind1]/(1-percept_givenM[ind2,ind1])\n",
    "                     p_modeled_stopping_givenM[ind2,ind1] = (factor_givenM**(2*exp1)- factor_givenM**(exp1+exp2))/(factor_givenM**(2*exp1)-1)    \n",
    "                else:        \n",
    "                     p_modeled_stopping_givenM[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    \n",
    "                if percept_givenB[ind2,ind1]!=0.5:     \n",
    "                    factor_givenB = percept_givenB[ind2,ind1]/(1-percept_givenB[ind2,ind1])\n",
    "                    p_modeled_stopping_givenB[ind2,ind1] = (factor_givenB**(2*exp1)- factor_givenB**(exp1+exp2))/(factor_givenB**(2*exp1)-1)\n",
    "                else:\n",
    "                    p_modeled_stopping_givenB[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    \n",
    "                mle = mle + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])                        \n",
    "                        \n",
    "    if isnan(mle):\n",
    "        mle = -np.inf\n",
    "        \n",
    "    print mle\n",
    "    return -mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " #Ca -1(discrete)1 - long formula for start point and threshold -- looking for M,N by brute force\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    percept, mu, mu_givenM, mu_givenB, sigma, delta, tau = find_mu_sigma_delta()  \n",
    "    pausetimes, laptimes,  pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop = pause_and_laptimes()\n",
    "    total_len = len(df)\n",
    "    \n",
    "    x_initial = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "    print \"starting points\", x_initial\n",
    "    \n",
    "    p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    start_N = delta\n",
    "    end_N = delta*10\n",
    "    no_N = (end_N-start_N+delta)/delta\n",
    "    array_N = np.linspace(start_N,end_N,no_N,endpoint=False)\n",
    "    exp1N = (int)(start_N/delta)\n",
    "    exp2N = (int)(end_N/delta)\n",
    "    \n",
    "    norm2 = np.zeros((no_N, 2*no_N-1))\n",
    "    mle2 = np.zeros((no_N, 2*no_N-1))\n",
    "    xopt2 = np.zeros((no_N*(2*no_N-1), 4))\n",
    "    avg_reward3 = np.zeros((no_N, 2*no_N-1))\n",
    "    avg_reward2 = np.zeros((no_N, 2*no_N-1))\n",
    "    avg_reward_modified = np.zeros((no_N, 2*no_N-1))\n",
    "    avg_reward_numerator = np.zeros((no_N, 2*no_N-1))\n",
    "    avg_reward_denominator = np.zeros((no_N, 2*no_N-1)) \n",
    "        \n",
    "    increment = 0\n",
    "    \n",
    "    for exp1 in range(exp1N, exp2N):\n",
    "        for exp2 in range(-exp1+1, exp1):            \n",
    "            \n",
    "            i=exp1\n",
    "            j=exp2+exp1-1\n",
    "            \n",
    "            print \"exp1,exp2 v1\", exp1, exp2,i,j\n",
    "            #results = minimize(fitting_accumulator2, x_initial)\n",
    "            #xopt[increment,:] = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "            #print \"actual results\", xopt[increment, :] \n",
    "            \n",
    "            var = xopt[2,:]\n",
    "            prob_misID = probability_misidentify(var)\n",
    "            percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "            \n",
    "            for ind2,o in enumerate(orientations):\n",
    "                for ind1,s in enumerate(sizes):\n",
    "\n",
    "                    if percept[ind2,ind1]!=0.5:\n",
    "                        factor = percept[ind2, ind1]/(1-percept[ind2,ind1])\n",
    "                        p_modeled_stopping[ind2,ind1] = (factor**(2*exp1)- factor**(exp1-exp2))/(factor**(2*exp1)-1)\n",
    "                    else:\n",
    "                        p_modeled_stopping[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                        \n",
    "                    if percept_givenM[ind2,ind1]!=0.5:        \n",
    "                        factor_givenM = percept_givenM[ind2,ind1]/(1-percept_givenM[ind2,ind1])\n",
    "                        p_modeled_stopping_givenM[ind2,ind1] = (factor_givenM**(2*exp1)- factor_givenM**(exp1-exp2))/(factor_givenM**(2*exp1)-1)    \n",
    "                    else:        \n",
    "                        p_modeled_stopping_givenM[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    \n",
    "                    if percept_givenB[ind2,ind1]!=0.5:     \n",
    "                        factor_givenB = percept_givenB[ind2,ind1]/(1-percept_givenB[ind2,ind1])\n",
    "                        p_modeled_stopping_givenB[ind2,ind1] = (factor_givenB**(2*exp1)- factor_givenB**(exp1-exp2))/(factor_givenB**(2*exp1)-1)\n",
    "                    else:\n",
    "                        p_modeled_stopping_givenB[ind2,ind1] = (exp1+exp2)/(2*exp1)\n",
    "                    \n",
    "                    p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                    mle2[i,j] = mle2[i,j] + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])                        \n",
    "               \n",
    "                    if o<o_T or s<s_T:\n",
    "                         avg_reward_numerator[i,j] = avg_reward_numerator[i,j] + p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1]*trials_tot[ind2,ind1]/total_len\n",
    "                    avg_reward_denominator[i,j] = avg_reward_denominator[i,j] + (laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1])*trials_tot[ind2,ind1]/total_len    \n",
    "                  \n",
    "            norm2[i,j] = math.sqrt(LA.norm(p_modeled_stopping_givenM - p_stop_real_givenM,'fro')**2 + LA.norm(p_modeled_stopping_givenB - p_stop_real_givenB,'fro')**2)       \n",
    "            avg_reward3[i,j] = avg_reward_numerator[i,j]/avg_reward_denominator[i,j] \n",
    "            increment = increment+1\n",
    "        \n",
    "    minimum = np.inf\n",
    "    maximum =-np.inf\n",
    "    maximum_reward = -np.inf\n",
    "    \n",
    "    for exp1 in range(exp1N, exp2N):\n",
    "        for exp2 in range(-exp1+1, exp1):            \n",
    "            \n",
    "            i=exp1\n",
    "            j=exp2+exp1-1    \n",
    "            \n",
    "            if norm2[i,j]<minimum:\n",
    "                minimum = norm2[i,j]\n",
    "                i_min = [i,exp1]\n",
    "                j_min = [j,exp2]        \n",
    "            \n",
    "            if mle2[i,j]>maximum:\n",
    "                maximum = mle2[i,j]\n",
    "                i_max_mle = [i,exp1]\n",
    "                j_max_mle = [j,exp2]\n",
    "                \n",
    "            if avg_reward3[i,j]>maximum_reward:\n",
    "                maximum_reward = avg_reward3[i,j]\n",
    "                N_rew = [i,exp1]\n",
    "                M_rew = [j,exp2]     \n",
    "            \n",
    "    print \"norm\", minimum, i_min, j_min\n",
    "    print \"mle\", maximum, i_max_mle, j_max_mle\n",
    "        \n",
    "    a,b = norm2.shape\n",
    "    print \"a,b\", a,b\n",
    "    print \"norm\", norm\n",
    "    \n",
    "    aic = 4 - 2*mle2[i_max_mle,j_max_mle]\n",
    "    print \"aic =\", aic\n",
    "    \n",
    "    print \"avg reward rate for best fit is II\", avg_reward3[i_max_mle[0], j_max_mle[0]], \"or (if minimizing norm)\", avg_reward3[i_min[0],j_min[0]]\n",
    "    print \"maximum reward rate is\",  maximum_reward, \"at N,M\", N_rew[1], M_rew[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fitting_accumulator1(x):\n",
    "    \n",
    "    mle = 0\n",
    "    p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "    p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))\n",
    "    \n",
    "    prob_misID = probability_misidentify(x)\n",
    "    percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "    \n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real()\n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    \n",
    "    for ind2,o in enumerate(orientations):\n",
    "        for ind1,s in enumerate(sizes):\n",
    "            \n",
    "            for k in range(0,exp2-exp1+1):\n",
    "                p_modeled_stopping[ind2,ind1] = p_modeled_stopping[ind2,ind1] + ((percept[ind2,ind1])**exp1)*((1-percept[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                p_modeled_stopping_givenM[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1] + ((percept_givenM[ind2,ind1])**exp1)*((1-percept_givenM[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                p_modeled_stopping_givenB[ind2,ind1] = p_modeled_stopping_givenB[ind2,ind1] + ((percept_givenB[ind2,ind1])**exp1)*((1-percept_givenB[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        \n",
    "            p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "            mle = mle + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])            \n",
    "            \n",
    "    norm = math.sqrt(LA.norm(p_stop_real_givenM - p_modeled_stopping_givenM, 'fro')**2 + LA.norm(p_stop_real_givenB - p_modeled_stopping_givenB, 'fro')**2) \n",
    "    \n",
    "    if isnan(mle):\n",
    "        mle = -np.inf\n",
    "    \n",
    "    print norm, mle\n",
    "    return -mle\n",
    "    #return norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "  #Ca -1(discrete)1 - long formula for steps +1,0 -- looking for M,N by brute force\n",
    "    p_m_real, p_b_real, trials_tot, trials_stop, trials_m, trials_b, trials_stop_m, trials_stop_b = p_real() \n",
    "    p_stop_real, p_stop_real_toFit, p_stop_real_toFit2, p_stop_real_givenM, p_stop_givenB, p_stop_real_givenM_toFit, p_stop_real_givenM_toFit2, p_stop_real_givenB_toFit, p_stop_real_givenB_toFit2 = p_stop_real_function()\n",
    "    percept, mu, mu_givenM, mu_givenB, sigma, delta, tau = find_mu_sigma_delta()  \n",
    "    total_len = len(df)\n",
    "    pausetimes, laptimes,  pausetimes_givenM, laptimes_givenM, pausetimes_givenB, laptimes_givenB, p_Correct_GivenStop = pause_and_laptimes()\n",
    "    \n",
    "    x_initial = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "    print \"starting point\", x_initial\n",
    "    \n",
    "    start_N = delta\n",
    "    end_N = delta*10\n",
    "    no_N = (end_N-start_N+delta)/delta\n",
    "    array_N = np.linspace(start_N,end_N,no_N,endpoint=False)\n",
    "    exp1N = (int)(start_N/delta)\n",
    "    exp2N = (int)(end_N/delta)\n",
    "    \n",
    "    start_M = tau\n",
    "    end_M = tau*15\n",
    "    no_M = (end_M-start_M+tau)/tau\n",
    "    array_M = np.linspace(start_M,end_M,no_M, endpoint =False)\n",
    "    exp1M = (int)(start_M/tau)\n",
    "    exp2M = (int)(end_M/tau)\n",
    "    \n",
    "    norm = np.zeros((no_N, no_M))\n",
    "    mle = np.zeros((no_N, no_M))\n",
    "    xopt = np.zeros((no_N*no_M, 4))\n",
    "    avg_reward = np.zeros((no_N, no_M))\n",
    "    avg_reward_numerator = np.zeros((no_N,no_M))\n",
    "    avg_reward_denominator = np.zeros((no_N,no_M))\n",
    "    \n",
    "    increment = 0\n",
    "    \n",
    "    for exp1 in range(exp1N, exp2N):\n",
    "        for exp2 in range(exp1,exp2M):            \n",
    "             \n",
    "            p_modeled_stopping = np.zeros((len(orientations),len(sizes)))\n",
    "            p_modeled_stopping_givenM = np.zeros((len(orientations),len(sizes)))\n",
    "            p_modeled_stopping_givenB = np.zeros((len(orientations),len(sizes)))\n",
    "            p_modeled_stopping2 = np.zeros((len(orientations),len(sizes)))    \n",
    "            \n",
    "            print \"exp1,exp2 v1\", exp1, exp2\n",
    "            results = minimize(fitting_accumulator1, x_initial)\n",
    "            xopt[increment,:] = [results.x[0], results.x[1], results.x[2], results.x[3]]\n",
    "            print \"actual results found\", xopt[increment,: ], \"for \", exp1, exp2\n",
    "            \n",
    "            prob_misID = probability_misidentify(xopt[increment,:])\n",
    "            percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "            \n",
    "            for ind2,o in enumerate(orientations):\n",
    "                for ind1,s in enumerate(sizes):\n",
    "                    \n",
    "                    if (ind2==2)&(ind1==4)&(exp1==5)&(exp2==6):\n",
    "                            print \"wait, seriously?!\", percept_givenM\n",
    "                            for k in range(0,exp2-exp1+1):\n",
    "                                print ((percept_givenM[ind2,ind1])**exp1)*((1-percept_givenM[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                                \n",
    "                    for k in range(0,exp2-exp1+1):\n",
    "                        p_modeled_stopping[ind2,ind1] = p_modeled_stopping[ind2,ind1] + ((percept[ind2,ind1])**exp1)*((1-percept[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        p_modeled_stopping_givenM[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1] + ((percept_givenM[ind2,ind1])**exp1)*((1-percept_givenM[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        p_modeled_stopping_givenB[ind2,ind1] = p_modeled_stopping_givenB[ind2,ind1] + ((percept_givenB[ind2,ind1])**exp1)*((1-percept_givenB[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        \n",
    "                    p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]\n",
    "                    mle[exp1,exp2] = mle[exp1,exp2] + trials_stop_m[ind2,ind1]*log(p_modeled_stopping_givenM[ind2,ind1])+(trials_m[ind2,ind1]-trials_stop_m[ind2,ind1])*log(1-p_modeled_stopping_givenM[ind2,ind1])+trials_stop_b[ind2,ind1]*log(p_modeled_stopping_givenB[ind2,ind1])+(trials_b[ind2,ind1]-trials_stop_b[ind2,ind1])*log(1-p_modeled_stopping_givenB[ind2,ind1])            \n",
    "            \n",
    "                    if o<o_T or s<s_T:\n",
    "                       avg_reward_numerator[exp1,exp2] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1]*trials_tot[ind2,ind1]/total_len\n",
    "                    avg_reward_denominator[exp1,exp2] = (laptimes[ind2,ind1]+p_modeled_stopping2[ind2,ind1]*pausetimes[ind2,ind1])*trials_tot[ind2,ind1]/total_len \n",
    "                    avg_reward[exp1,exp2] = avg_reward_numerator[exp1,exp2]/avg_reward_denominator[exp1,exp2]\n",
    "                        \n",
    "            norm[exp1,exp2] = math.sqrt(LA.norm(p_stop_real_givenM - p_modeled_stopping_givenM, 'fro')**2 + LA.norm(p_stop_real_givenB - p_modeled_stopping_givenB, 'fro')**2)\n",
    "            print \"norm\", norm[exp1, exp2], \"mle\", mle[exp1, exp2]\n",
    "            \n",
    "            if (exp1==5) & (exp2==6):\n",
    "                print \"here we go:\"\n",
    "                print p_modeled_stopping_givenM\n",
    "                print p_modeled_stopping_givenB\n",
    "                print math.sqrt(LA.norm(p_modeled_stopping_givenM - p_stop_real_givenM,'fro')**2 + LA.norm(p_modeled_stopping_givenB - p_stop_real_givenB,'fro')**2)       \n",
    "                print mle[5,6]\n",
    "                print \"mle,\", mle\n",
    "            \n",
    "            increment = increment+1\n",
    "        \n",
    "    minimum = np.inf\n",
    "    maximum =-np.inf\n",
    "    maximum_reward = -np.inf\n",
    "    \n",
    "    for i in range(exp1N,exp2N):\n",
    "        for j in range(i, exp2M):\n",
    "        \n",
    "            if norm[i,j]<minimum:\n",
    "                minimum = norm[i,j]\n",
    "                i_min = i\n",
    "                j_min = j        \n",
    "            \n",
    "            if mle[i,j]>maximum:\n",
    "                maximum = mle[i,j]\n",
    "                i_max_mle = i\n",
    "                j_max_mle = j\n",
    "                \n",
    "            if avg_reward[i,j]>maximum_reward:\n",
    "                maximum_reward = avg_reward[i,j]\n",
    "                N_rew = i\n",
    "                M_rew = j     \n",
    "            \n",
    "    print minimum, i_min, j_min\n",
    "    print maximum, i_max_mle, j_max_mle\n",
    "        \n",
    "    a,b = norm.shape\n",
    "    print \"a,b\", a,b\n",
    "    print \"norm\", norm\n",
    "    \n",
    "    aic = 4 - 2*mle[i_max_mle,j_max_mle]\n",
    "    print \"aic =\", aic\n",
    "    \n",
    "    print \"avg reward rate for best fit is II\", avg_reward[i_max_mle,j_max_mle], \"or (if minimizing norm)\", avg_reward[i_min,j_min]\n",
    "    print \"maximum reward rate is\",  maximum_reward, \"at N,M\", N_rew, M_rew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_modeled_stopping = np.zeros((len(orientations), len(sizes)))\n",
    "p_modeled_stopping_givenM = np.zeros((len(orientations), len(sizes)))\n",
    "p_modeled_stopping_givenB = np.zeros((len(orientations), len(sizes)))\n",
    "\n",
    "exp1 = 4\n",
    "exp2 = 7\n",
    "increment = 42\n",
    "\n",
    "prob_misID = probability_misidentify(xopt[increment,:])\n",
    "percept, percept_givenM, percept_givenB = transform_intoPercept(prob_misID)\n",
    "\n",
    "for ind2,o in enumerate(orientations):\n",
    "    for ind1,s in enumerate(sizes):\n",
    "                     \n",
    "        for k in range(0,exp2-exp1+1):\n",
    "            p_modeled_stopping[ind2,ind1] = p_modeled_stopping[ind2,ind1] + ((percept[ind2,ind1])**exp1)*((1-percept[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "            p_modeled_stopping_givenM[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1] + ((percept_givenM[ind2,ind1])**exp1)*((1-percept_givenM[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "            p_modeled_stopping_givenB[ind2,ind1] = p_modeled_stopping_givenB[ind2,ind1] + ((percept_givenB[ind2,ind1])**exp1)*((1-percept_givenB[ind2,ind1])**k)*comb(exp1+k-1, exp1-1)\n",
    "                        \n",
    "        p_modeled_stopping2[ind2,ind1] = p_modeled_stopping_givenM[ind2,ind1]*p_m_real[ind2,ind1] + p_modeled_stopping_givenB[ind2,ind1]*p_b_real[ind2,ind1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(p_modeled_stopping2, interpolation='none')\n",
    "\n",
    "plt.suptitle('Probability of stopping')\n",
    "\n",
    "plt.xlabel('M')\n",
    "plt.ylabel('N')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(p_modeled_stopping_givenM, interpolation='none')\n",
    "\n",
    "plt.suptitle('Probability of stopping given mushroom')\n",
    "\n",
    "plt.xlabel('M')\n",
    "plt.ylabel('N')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(p_modeled_stopping_givenB, interpolation='none')\n",
    "\n",
    "plt.suptitle('Probability of stopping given bird')\n",
    "\n",
    "plt.xlabel('M')\n",
    "plt.ylabel('N')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(avg_reward, interpolation='none')\n",
    "\n",
    "plt.suptitle('Average Reward Rate')\n",
    "\n",
    "plt.xlabel('M')\n",
    "plt.ylabel('N')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "print avg_reward.max()\n",
    "print avg_reward.argmax()\n",
    "print avg_reward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print mle\n",
    "plt.imshow(mle, interpolation='none')\n",
    "\n",
    "plt.suptitle('MLE')\n",
    "\n",
    "plt.xlabel('M')\n",
    "plt.ylabel('N')\n",
    "\n",
    "plt.jet()\n",
    "plt.colorbar()\n",
    "\n",
    "print mle.max()\n",
    "print mle.argmax()\n",
    "print mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "#save object:\n",
    "with open('MN_swipe_vars_model3.pickle', 'w') as f:\n",
    "    pickle.dump([xopt], f)\n",
    "    \n",
    "#get back object:\n",
    "with open('MN_swipe_vars_model3.pickle') as f:\n",
    "    new_var = pickle.load(f)\n",
    "    \n",
    "new_var2 = np.asarray(new_var)\n",
    "new_var2 = np.squeeze(new_var2)\n",
    "   \n",
    "with open('model_mle_norm_model3.pickle', 'w') as f:\n",
    "    pickle.dump([mle, norm], f)    \n",
    "    \n",
    "with open('model_mle_norm_model3.pickle') as f:\n",
    "    norm2 = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}