TD and Markov.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A wrong view of TD equation\n",
    "As from the TD equation $V(s_t) = V(s_t) +\\alpha (R_t - V(s_t))$, intuitively we see the term $R_t - V(s_t)$ as somewhat of supervised learning fashion, as we treat V(s_t) as the expectation of the reward given the current state is $s_t$ and we try to minimize the difference between estimated reward and the true reward collected at time t. This view along with the form of the equation remind us of the gradient based method, seeing $(R_t - V(s_t))$ as the gradient of some function E, and we are trying to optimize this funciton. However, later on we will prove that this gradient based view is wrong, and a more plausible way of viewing TD is through Markov process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov process view of TD equation\n",
    "Now consider that we are in the setting of Markov process. Let us define $w_s$ as the probability of state s giving rise to the terminal state z, and then $w_s$ can be estimated through emperical frequency, i.e.: $w_s = n_zs/n_s$, where n_zs is the frequency of seeing z after s, and n_s is just the frequency seeing s. Then we can derive an iterative update for $w_s$. After we encounter the state s, two things could happen, either we go directly into the terminal state z or we go to some other state r. Denote the transition probability from state s to state r as $P_sr$, and the probability of encountering terminal state after seeing s as $h_s$, we will have $w_s = h_s + \\sum_r P_sr w_r$ as our interative update rule. All these probability can be estimated from sample, i.e. $P_sr = m_sr/n_s$, $h_s = m_sz/n_s$ where $m_sr$ is the frequency of seeing r after seeing s. So now, plug these estimate back to the update euqation, we will have: $\\sum_r{m_sr w_s} - n_s w_s + m_sz = 0$, and then we can write it in matrix form, $(M - N)w + m = 0$. This is our familiar linear system again, denote M - N as A while m as b, and the iterative solution to this system should be $w = w + \\alpha (Aw +n)$. However if we view this form as a gradient update, then the function we try to optimize would be $(Aw + b - 0)^2$ and the true gradient update would be $w = w + \\alpha A^T(Aw +b)$ thus it's different from the iterative rule we had. \n",
    "<br />\n",
    "<br />\n",
    "Another justification of this would be looking into the second order of the derivative. Now let us define $x_t$ as a vector of all zeros with the only one appearing on the corresponding state that it is current in. For example, suppose we have 2 state, {0, 1}, in total and at time 1, the process enter state 0, then $x_1 = [1, 0]$. Also, define $\\sigma_{s_t, z}$ being one only if $s_t$ is actually $z$. Then our updating rule can be rewritten as $w = w + \\alpha ((x_t \\otimes x_{t+1} - x_t \\otimes x_{t})w + x_t \\sigma_{s_{t+1}, z})$. Now for this update to be a gradient based method, the gradient of the corresponding function E has to be $(x_t \\otimes x_{t+1} - x_t \\otimes x_{t})w + x_t \\sigma_{s_{t+1}, z}$. If we take the second order of derivative w.r.t $w_j$ and $w_i$, the result will be $\\frac{\\partial^2 E}{\\partial w_i \\partial w_j} = x_{ti}(x_{(t+1)j} - x_{tj})$ and  $\\frac{\\partial^2 E}{\\partial w_j \\partial w_i} = x_{tj}(x_{(t+1)i} - x_{ti})$. However, these two values in general will not equal to each other as  $x_{ti}x_{(t+1)j}$ will generally not equal to $x_{tj}x_{(t+1)i}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiment for the first order Markov Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we want to conduct experiment to see if this TD like update equation will actually recover the probability of termination given different states. The Markov model that we use is as shown in the following graph. \n",
    "<img src=\"files/Markov_model.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have in total seven states, let us denote them as 0 to 6, from left side to right. State 0 and 6 arethe termination state, and the process always start with state 3 and each transfer probability is 1/3. We generate different number of samples ranging from 5 to 1000 and we evaluate the performance based on the mean squared error between the predicted probability of termination given each state, which is $w_s$ and the true probability calculated through math. From the result we can observe that the error shrinks as we increase the sample size, also it reaches to a relatively small value, indicating this interative method works for the first order Markov Model. However, the result I got is not the same as in the paper, perhaps he used a differnent normalization method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def gen_first_order_samples(num_samples):\n",
    "    A = np.asarray([[0.,0.,0.,0.,0.,0.,0.],\n",
    "         [1/3,1/3,1/3,0.,0.,0.,0.],\n",
    "         [0.,1/3,1/3,1/3,0.,0.,0.],\n",
    "         [0.,0.,1/3,1/3,1/3,0.,0.],\n",
    "         [0.,0.,0.,1/3,1/3,1/3,0.],\n",
    "         [0.,0.,0.,0.,1/3,1/3,1/3],\n",
    "         [0.,0.,0.,0.,0.,0.,0.]\n",
    "        ])\n",
    "    seq_samples = []\n",
    "    for j in range(0, num_samples):\n",
    "        obs = [3]\n",
    "        while True:\n",
    "            thr = np.random.randn(1)\n",
    "            if thr <= 1/3:\n",
    "                obs.append(obs[-1] - 1)\n",
    "            elif thr <= 2/3:\n",
    "                obs.append(obs[-1])\n",
    "            else:\n",
    "                obs.append(obs[-1] + 1)\n",
    "            if obs[-1] == 0 or obs[-1] == 6:\n",
    "                break\n",
    "        seq_samples.append(np.asarray(obs))\n",
    "    seq_samples = np.asarray(seq_samples)\n",
    "    return seq_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def sum_to_one(x):\n",
    "    return x/np.sum(x)\n",
    "def cal_first_order_term_prob(index):\n",
    "    temp = (1./3)**index\n",
    "    temp2 = (1./3)**(6-index)\n",
    "    return temp/(1.-temp) + temp2/(1.-temp2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gen_vecsamples_from_seqsamples(seq_samples):\n",
    "    vec_samples = []\n",
    "    for i in range(0, len(seq_samples)):\n",
    "        vec = []\n",
    "        for j in range(0, len(seq_samples[i])):\n",
    "            temp = np.zeros(7)\n",
    "            if seq_samples[i][j] != 0 and seq_samples[i][j] != 6:\n",
    "                temp[seq_samples[i][j]] = 1.\n",
    "            \n",
    "            vec.append(temp)\n",
    "        vec_samples.append(np.asarray(vec))\n",
    "    return np.asarray(vec_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cal_prob_each_state(index):\n",
    "    diff = abs(index - 3)\n",
    "    if diff == 0:\n",
    "        temp = (1./3)**1\n",
    "        temp2 = (1./3)**2\n",
    "        temp3 = (1./3)**4\n",
    "        return temp/(1.-temp) + temp2/(1.-temp2) + temp3/(1.-temp3)\n",
    "    if diff == 1:\n",
    "        temp = (1./3)**1\n",
    "        temp2 = (1./3)**3\n",
    "        temp3 = (1./3)**5\n",
    "        return temp/(1.-temp) + temp2/(1.-temp2) + temp3/(1.-temp3)\n",
    "    if diff == 2:\n",
    "        temp = (1./3)**2\n",
    "        temp2 = (1./3)**4\n",
    "        temp3 = (1./3)**6\n",
    "        return temp/(1.-temp) + temp2/(1.-temp2) + temp3/(1.-temp3)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.07001067114703478, 0.07460664335664335, 0.04903846153846152, 0.07460664335664335, 0.07001067114703478]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGHCAYAAABrpPKuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xm8W1W9/vHP0xFKS5lb5lLaBBxQWkVRBL0qFRzu9YrD\ncQRHFIVfvV7F61XUq+IIjiiKgjhUcbw4oqhXEQG1FRxA6EgZC2UolJa2tN/fH2uH7qbJOclpkr1P\n87xfr7zSrOxhZSU9ebLW2nsrIjAzMzMri1FFV8DMzMwsz+HEzMzMSsXhxMzMzErF4cTMzMxKxeHE\nzMzMSsXhxMzMzErF4cTMzMxKxeHEzMzMSsXhxMzMzErF4cSsiyTNkbRJ0hG5sm9Juq5uuZ0lXSDp\n9mz5D2Xl+0j6oaS7JG2U9PpevwbbmqTx2fv00aLr0gpJYyWdLemmrN7fLLpOnSSpmr2uFxVdF+sM\nh5M+JulV2X/oTZKe1GSZ2h+zi+vKd5L0Pkl/k7Ra0kpJf5H0SUlTc8udkdtH/W2jpL26/TqHS9K7\nJT27A5uqv0ZEAJvqyt4LvAT4FPBy4FtZ+WeBo4H3A68AftWB+hSug21rrXkjcBrwDdLn67PFVsds\ncGOKroCVwlrgpcAf8oWSjgH2BR6sKx8DXAZUgK8CnwYmAo8EBoDvA7fnVgngZOCBBvu+tyOvoDve\nA3wJ+EmHt/tyQHVlTwN+FxFnNii/KCI+1eE6FK1bbWuNPQ1YFBGnF10Rs1Y4nBjAT4EXSjo1IvK/\n6F8K/BnYo2755wOPBQYi4tv5JySNA8Y12Mf3IuLuDtZ5m2Uhi4h4qJf7jYiNDYr3AhbmC7K2nAys\n6tS+i3rNNjySdoiIB4deckh7Ue4fAmZb8LCOBTAP2B14Zq1Q0ljgBOCbbP0rf3q23h/qyomI9RGx\nuhMVk/Q9SfPryn6UDQk9J1d2RFY2Z5BtHZwtc5qkt0paTOoxqjRYdrykTaT/HyfnhqHOGaK+B2b1\nW53NHfkoMLbBct+SdG327znZvqYCJ+SGu95B6rEK4G1Z+ZrcNnaT9Nls2G2dpBskvbVuP7Vx+FMk\nvU3Skuw1T8/NhXmepPdKukXSGkmXSDqwQZ2fLOmXklZlr+9X+Xk02TIfzra5v6SvS7pX0t2Szs2C\n1rDatp26Zu2+1bYkXSnppw22+a+SPiDp1uy1fVPSBEk7SPqcpDsk3Ze9hoY/5iSdmLX/WklXSXpi\ng2X2l3ShpBWSHpT0V0kvb/I6ny/pI5JuAVbX2q7JvidJ+rSkm7PtXifp1Nzz1ay9nwg8Lvf5OmKQ\nbR4p6VKlodoHJC2W9Pm6Zd4p6Q9Kc6HWZK/7eXXLPDwvR9JAVrc1ki6TdEi2zFuy7a/NPl/71G3j\nSkl/lPQESVdk6y+S9Opm9a9b/1GSflBXz2fVLTMu+wwszOpxp6TfSjq6lX1Yd7jnxACWAVeShmQu\nycqOB3YmzX04rW75G0mB5ZXAB1vcx+6S6kPOQxExWK/AZcDzJE3MBZ4nARuBpwA/zsqOzsoub6Ee\nryMFhi8A62n8a3I9aejlQuDXwPlZ+Q3NNippIvAbYE/gbGAl8Crg2AaL5+egXJPt63PA9aQhslr5\n7dm+f0IKkBtz+/o9sFv2Om4htcHHJe0REf9Vt783AqOBc4CH2LIn5gxgHfBhUkB9O3ABaRig9tqe\nBfwvcAVpOAbgtcD/SXpiRPw197oC+CGprd4BHJEteyvwPobRtu3Ula3n9wxVfgZwP/AB4FDgTaRQ\nuDPp7+N7gKOy17AQ+Hjd+nPYPIfjIeAU4BJJsyNiEYCkfYE/AmuATwJ3A88BLpQ0ISK+WLfN/yEN\ngX4E2Insfa8naRTwM1LwOBf4O/Bs4JOSpkTEu0jt/nJS2wt4d3a/qMk298m2eXPWJvcDB2XbzTsN\n+DbpfRyf7eP7ko6NiF/XLftM4AXA50lt+l/AxVmIfBVpntVewH8CX8zapiaAKcDFpPky3yD9nTpP\n0pqI+BZNSHos8DtgCfAhUjAfAH4s6TkR8fNs0TOz1/N54C+k3sojSL3Dv2u2feuyiPCtT2+kPwwb\ngVmkP8r3AuOz574NXJr9eylwcW69HYDrSJM6lwJfAU4C9mywjzOy5Rrdrh2ifrOz5eZkjx+VPf4W\n8Ifccj8E/jzEtg7O1r0L2KXF9tkAnNPisu/I2vL4XNmErH02AkfkyufVv3bgNtLcknzZ+KzOH60r\n/wBwD7B/XflZpC/WPbPH1Wz9O4Gd65adkz23ABidK//PrL7Ts8ejstfw/br1JwDLgR/mys7Mtvnp\numV/AizfhrZtqa65dtxqu6Rg9dMG2/wTMCpX/r1sm9+tW39+/j3LvTcbgENz5dNJAerrubKvZ204\nqW6b3wfuAMbU1ekftbIh2uXF2fL/r678f0khcN+61//HFre5Mf+amiw3vu7xWOCfwI8atNH9wNRc\n+Vuy8mX57QCfIAW8KXX13gi8oW67fwduzJXVPusvypX9Hriq7v1V9p5fnSu7jrr/e74Vf/OwjtVc\nRPrCeU72y/w5pF8pW4k0Bn4E8FHSL5tXAV8Gbsu6mOuHMoI0T+UZdbeThqjTX4DVpF4BSL0lN5F+\nrc2WtENWfhSpl6UVF0VEN8bejwOWRcTDQwcRsYbULp12AqnXYY2k3Ws34FLSfJ+j6pb/VkTc12Rb\n58WWc2AuI/0Bn549PgI4EJhXt68JpJ6ip7GlIP2Kz7sM2KfB56JdjepKrq7DcX5sOc/qquz+K3XL\nXQVMa7D+/0XEw4eFR8QS0hyu4wAkjQb+lRQYxtW14SWkHqBH123zK9HanKDjSL0BX6grP4vUQ9F0\nmHMQ95Le/3/N6t5QRKyr/VvSLqTehstJP3Tq/Swi8hPka2387fx2snKRemry1rK5h6227y8B+0k6\nrFH9lI4YfBLpR9auuTbfHfgF8GhJu+Ze82GS6vdrBfKwjgEQESslXUqaBLsT6RfzdwdZ/n7gdOB0\nSfsDTwfeRurWvpfN3f81l0WbE2IjYpOkK0ihhOz+MtIfwdHAEyXdQRreaDWcLGunDm04kDQsU69R\n2baaAcwkBb56Qeoiz1s2yLZuqnt8T3Zf+8M9M7v/NlsLICSNr/uSWd5gmwJ2IfXiDFejuipX105s\nc9Ug5eMl7RgRa3PljYZHbiANR+5M+tLeidRbcGqDZdt9v/IOBG6KrSfMXpd7vl2/IA2hfBB4h6Rf\nk3omvxURG2oLSXo+8E5SsBqfW38NW2vWxjc3Ka9/P2+KiPV1ZbVhwGnAX9la7XP7MbYeioPN7X4P\n8C5Sj9liSX8lDWt9LSKubbCe9YjDieV9k/SLZG/Sr537W1kpIm4CLpD0Q9L47svYOpwM1++B/5I0\nnhRO/iciVkn6e/b4DtIfmlbDydqhFymv3Lydn5DmLzTyz7rHg73mhvMZ2DwJehSpfU9l85devfov\njqG2OVytbLfZ3JJmvQDNttmp11Drnf4KaTivkavrHhf2GY2IAP5N0pGk3tM5pNMFnCbpyRGxTtIz\nSV/mvwTeQJob9RDpdAHPabDZdtpYtN7Ggy1Xa/cPkXr4GlkOEBG/lnQwqYfrWNJr+g9JJ0VEw95j\n6z6HE8v7AalL/gmksee2RMS9SkfBPLKDdbqMNFQxAOzD5hDyO9JwzwrghojYll/kzTT7omvkRjb/\nWss7pEN1AdKXh6RlwITYeuJhNywmfQms6vD+2mnbdtxD6qGpdyBpknKnNXrPK8C9EXGfpLWksKEu\nvF83Akc06Lk6NPf8sETEFaT5Hu+SdBJpePIFpB8w/07q5TguPyQm6ZTh7m8I+0saV9d7UiV9hpY1\nWWdxdr+ulXbPenXPB87PhrWvIM2XczgpiOec2MMi4gHSr5/3Aj9qtpykw7Lx2/ryA4FHsPUv921x\nFelX2TuAu3Pj+5eRjlI4mrpeE0ljskMop7S6k2z5feuKH6DxF10jPwWmKXfW0+yPXEuHPLbpIuCp\njQ51lLRrg6OimmklIFxJ6pZ/u6QdG+yv/hw4rWqnbaH1MLMYeFJ2JAsAkk4gHUU13G0O5hhJD4fx\n7Bf48aShAbKhkP8FBiQ1Omy9vv3aqdNPgR1J/2fz5pL+z/x8qzWGkJuHkXdNdl8bvtlImnz6cG+U\npJmk172tGr3+HYHX5PY1jnTU3S3A3xpuJOJm0mf3lEaf0XyZpN3q1l1N6gEeX7+e9Y57TmyLL7KI\n+FoL6zwTeJ/SKe2vJE1aPZg0wXUcKdzU7+OFkhqd/+QXg/V6RMRapXOdPJE0Fl7zO9JY/gS2HtI5\nkDQEcR4w5LVosiGj60h/zPN/YOcDxymdN2IF6Qyb8xtsAtJhum8Evi3pU6ThphNJ828OGKoObfoQ\n6dDOX0r6CmlYYBJwGOlX7V40HvuvN2SIiYiHJL2O9AX7N0kXkg5P3Y80qfkWhtHLRntt21JdM+eR\n5kj8TNL3Sb0YLyEdLTPcbQ7mH6T34bOkL+w3kYa53p9b5m2kScp/lvQl0mdtD+BxwJGkszAPp07f\nJc2/+ngWfGqHEh8HnBkRtw7j9bxe0qtIbbiEFCBfTzr8uXaagR+TXufPJX2b1KP5JtKPkuow9pnX\n6PXfBLw3C0C1YeNDgFdkw1DNnAz8Fvi7pPNIn4G9gSeT5rXUzkezWNLPSEeD3UN6T55DmvBvBXE4\nsVZ+qdXOX1HzXdLp6o8lHa2xG+k/9VXAJyKi/twAQfrybuRpDD1J8jLSUNPDISQiVkhaRDpSo9F8\nk/o6Nysb7LlTSUdCnEk6fPpc0pfq1itHrJb0VNL5Lk4jHT751axuFzdapcW6bVWe7evJwH+TutpP\nJHWzX0+apLi2bv1mWjonSET8QunaS+8mTezciXTI7hVsfaRIq7/8W27bNut6saTTs+0fRTq/yLOy\n7Tdq83Y0Wv8S4FrSe7EvaXLmiyLi4fO2RMStkh5Pmod1Aum8HStJYaL+dPIt1ymbMH4c6dDyE0i9\nC0tJhxZ/poX6N/Ir0vk9XprV8x7S+/zeWtiJiJ9LegPpUO5PknqrTiOF4/pwMtjnuuHLalB2B/Bm\n0jmATiaF49dFRP3FC+s/C3+V9DjS8MxrSIFkBSmEfCC36NmkUDeH1FuyNPfarCAaPHiamZkVIzta\nb3REND2jrW2fSjHnRNJTJF2sdFrqTao7DXKTdZ4qab7SKZtvyLoizczMbIQrRTghdRNfTRq3HLIr\nR9I00rjnr4DHkE5/fF52iJuZmZmNYKWYcxLpGgc/hy3O4zCYNwJLIuLt2ePrJR1FmqX+y+7U0szM\nCuC5B32oFOFkGJ5IOlV33iWkiU1mZrYdiIgji66DFaMswzrtmkqadZ23Atg5OyzUzMzMRqiR2nPS\ntuykYXNIZxSsvxaFmZmZNbcD6VpGl0TEXd3e2UgNJ7eTjsHPmwLcV3ca57w5+FTEZmZm2+JlpMsY\ndNVIDSdXkF2SPOfYrLyZZQBf//rXOfTQQwdZzDpp7ty5nH22pwL1ktu899zmvec2763rrruOl7/8\n5dC9K7tvoRThRNJOpMvA147UmS7pMaRrqdwk6Uxgn4ioncvkC6RrJnyEdLXPp5POkDjYtR0eBDj0\n0EOZNWtWN16GNTB58mS3d4+5zXvPbd57bvPC9GRaRFkmxD4O+Avp9NUBfIJ0iuH3Zc9PBfavLRwR\ny0inG34G6fwoc4HXRET9ETxmZmY2wpSi5yQifssgQSkiTmpQ9jtgdjfrZWZmZr1Xlp4TMzMzM8Dh\nxLpsYGCg6Cr0Hbd577nNe89tvn3rm6sSS5oFzJ8/f74nUZmZmbVhwYIFzJ49G2B2RCzo9v7cc2Jm\nZmal4nBiZmZmpeJwYmZmZqXicGJmZmal4nBiZmZmpeJwYmZmZqXSd+Hk5pth4cKia2FmZmbN9F04\nee974T3vKboWZmZm1kzfhZNp0+CGG4quhZmZmTXTd+HkgANSOOmTE+OamZmNOH0XTqZNg9Wr4bbb\niq6JmZmZNdJ34eSAA9L99dcXWw8zMzNrrO/CyX77wZgxDidmZmZl1XfhZMwYmD7dk2LNzMzKqu/C\nCUCl4p4TMzOzsurLcFKtuufEzMysrPoynFQqsHQprF9fdE3MzMysXl+Gk2oVNm6ExYuLromZmZnV\n69twAh7aMTMzK6O+DCdTpsCkSZ4Ua2ZmVkZ9GU6k1HvicGJmZlY+fRlOIE2K9bCOmZlZ+fRtOHHP\niZmZWTn1dTi58064556ia2JmZmZ5fRtOKpV076EdMzOzcunbcDJzZrr30I6ZmVm59G04mTgR9t3X\nPSdmZmZl07fhBDwp1szMrIz6Ppy458TMzKxc+jqcVCqwcCFs2lR0TczMzKymr8NJtQpr18JNNxVd\nEzMzM6vp+3ACHtoxMzMrk74OJwceCOPGeVKsmZlZmfR1OBk9GmbMcM+JmZlZmfR1OIE0KdY9J2Zm\nZuXR9+HE5zoxMzMrF4eTKixfno7aMTMzs+L1fTipVCACFi0quiZmZmYGDic+nNjMzKxk+j6c7L47\n7Lqr552YmZmVRd+HE8mTYs3MzMqk78MJ+AKAZmZmZeJwwuZznUQUXRMzMzNzOCH1nNxzD6xcWXRN\nzMzMzOGE1HMCHtoxMzMrA4cT0vV1JE+KNTMzKwOHE2DHHdMVit1zYmZmVjyHk4wvAGhmZlYOpQkn\nkk6RtFTSWklXSnr8EMu/TNLVkh6QdKukL0vabbj797lOzMzMyqEU4UTSi4FPAGcAhwPXAJdI2qPJ\n8k8Gvgp8CXgEcAJwBPDF4dahUknX19m4cbhbMDMzs04oRTgB5gLnRsSFEfFP4GRgDfDqJss/EVga\nEZ+LiBsj4g/AuaSAMizVKmzYAMuWDXcLZmZm1gmFhxNJY4HZwK9qZRERwKXAkU1WuwLYX9Jx2Tam\nAC8EfjLcevgCgGZmZuVQeDgB9gBGAyvqylcAUxutkPWUvBz4tqT1wG3APcCbh1uJ/fZLR+143omZ\nmVmxyhBO2ibpEcCngPcCs4A5wEGkoZ1hGTUKZs50ODEzMyvamKIrAKwENgJT6sqnALc3Wed04PKI\nOCt7/HdJbwIuk/SuiKjvhXnY3LlzmTx58hZlAwMDDAwMUKl4WMfMzPrbvHnzmDdv3hZlq1at6mkd\nCg8nEbFB0nzg6cDFAJKUPf50k9UmAOvryjYBAWiw/Z199tnMmjWr4XPVKlxwQctVNzMz2+7UfrDn\nLViwgNmzZ/esDmUZ1jkLeJ2kV0o6BPgCKYBcACDpTElfzS3/I+AFkk6WdFB2aPGngKsiollvy5Cq\nVbjlFli9etivw8zMzLZR4T0nABFxUXZOk/eThnOuBuZExJ3ZIlOB/XPLf1XSROAU4OPAvaSjfU7f\nlnrULgC4cCEcfvi2bMnMzMyGqxThBCAizgHOafLcSQ3KPgd8rpN1qIWT6693ODEzMytKWYZ1SmHX\nXWHPPT0p1szMrEgOJ3V8jR0zM7NiOZzUcTgxMzMrlsNJndq5TiKKromZmVl/cjipU63C/ffD7cM+\nINnMzMy2hcNJHV8A0MzMrFgOJ3WmT4fRoz3vxMzMrCgOJ3XGjYODDnI4MTMzK4rDSQO+AKCZmVlx\nHE4a8OHEZmZmxXE4aaBahSVLYMOGomtiZmbWfxxOGqhUYOPGFFDMzMystxxOGqgdTuyhHTMzs95z\nOGlg771h4kRPijUzMyuCw0kDUhracc+JmZlZ7zmcNFGtuufEzMysCA4nTbjnxMzMrBgOJ01Uq7Bi\nBaxaVXRNzMzM+ovDSROVSrr30I6ZmVlvOZw0UQsnHtoxMzPrLYeTJiZNgn32cTgxMzPrNYeTQfgC\ngGZmZr3ncDIIXwDQzMys9xxOBlGpwMKFsGlT0TUxMzPrHw4ng6hWYc0auOWWomtiZmbWPxxOBuEL\nAJqZmfWew8kgpk2DsWM9KdbMzKyXHE4GMWYMHHywe07MzMx6yeFkCD6c2MzMrLccTobgw4nNzMx6\ny+FkCNUqLFsGDz5YdE3MzMz6g8PJECoViIDFi4uuiZmZWX9wOBmCDyc2MzPrLYeTIey5J0ye7Emx\nZmZmveJwMgTJk2LNzMx6yeGkBQ4nZmZmveNw0gKf68TMzKx3HE5aUK3CXXelm5mZmXWXw0kLfMSO\nmZlZ7zictGDGjHTvoR0zM7PuczhpwYQJcMAB7jkxMzPrBYeTFnlSrJmZWW84nLTIhxObmZn1hsNJ\ni6pVWLQINm4suiZmZmbbN4eTFlUqsG4dLF9edE3MzMy2bw4nLfLhxGZmZr3hcNKi/feH8eM9KdbM\nzKzbHE5aNHo0zJzpnhMzM7Nuczhpg4/YMTMz6z6Hkzb4XCdmZmbd53DShmoVbroJHnig6JqYmZlt\nv0oTTiSdImmppLWSrpT0+CGWHyfpg5KWSXpQ0hJJJ3azjpVKul+0qJt7MTMz62+lCCeSXgx8AjgD\nOBy4BrhE0h6DrPYd4GnASUAFGAC6OiPEhxObmZl135iiK5CZC5wbERcCSDoZeDbwauCj9QtLehbw\nFGB6RNybFXf99Gi77QZ77OFwYmZm1k2F95xIGgvMBn5VK4uIAC4Fjmyy2nOBPwPvkHSzpOslfUzS\nDt2uryfFmpmZdVcZek72AEYDK+rKVwDVJutMJ/WcPAj8W7aNzwO7Aa/pTjWTahX+/vdu7sHMzKy/\nlSGcDMcoYBPw0ohYDSDprcB3JL0pItY1W3Hu3LlMnjx5i7KBgQEGBgZa2nGlAt/7HkSANOz6m5mZ\nldK8efOYN2/eFmWrVq3qaR3KEE5WAhuBKXXlU4Dbm6xzG3BLLZhkrgME7Acsbrazs88+m1mzZg27\nstUq3Hcf3HEHTKmvsZmZ2QjX6Af7ggULmD17ds/qUPick4jYAMwHnl4rk6Ts8R+arHY5sI+kCbmy\nKqk35eYuVTXtxEfsmJmZdVXh4SRzFvA6Sa+UdAjwBWACcAGApDMlfTW3/DeBu4DzJR0q6WjSUT1f\nHmxIpxMOPhhGjfKkWDMzs24pw7AOEXFRdk6T95OGc64G5kTEndkiU4H9c8s/IOmZwGeAP5GCyreB\nd3e7ruPHw7Rp7jkxMzPrllKEE4CIOAc4p8lzJzUouwGY0+16NVKpOJyYmZl1S1mGdUaUatXDOmZm\nZt3icDIM1SosXgwbNhRdEzMzs+2Pw8kwVCrw0EOwbFnRNTEzM9v+OJwMgw8nNjMz6x6Hk2HYd1+Y\nMMHhxMzMrBscToZB8gUAzczMusXhZJiqVfecmJmZdYPDyTC558TMzKw7HE6GqVqF225LFwE0MzOz\nznE4GabaETvuPTEzM+ssh5NhqlTSvcOJmZlZZzmcDNPOO8PUqZ4Ua2Zm1mkOJ9vAk2LNzMw6z+Fk\nG/hwYjMzs85zONkGtasTRxRdEzMzs+1HW+FE0tsl7Zh7/GRJ43OPJ0k6p5MVLLNKBR54AG69teia\nmJmZbT/a7Tk5E5iUe/wzYN/c4wnAG7a1UiOFLwBoZmbWee2GEw3xuK8cdBCMGeNwYmZm1kmec7IN\nxo6F6dN9xI6ZmVknOZxsIx+xY2Zm1lljhrHOayWtzq1/oqSV2eNJTdbZblUq8MMfFl0LMzOz7Ue7\n4WQ58Lrc49uBVzRYpm9Uq7B0KaxbB+PHD728mZmZDa6tcBIR07pUjxGrUoFNm2DxYnjEI4qujZmZ\n2cjnOSfbyFcnNjMz66x2T8J2pKTn1JW9UtJSSXdI+mL+pGz9YMqUdBFAT4o1MzPrjHZ7Tt4DPLL2\nQNKjgS8DlwIfBp4LvLNjtRsBJF8A0MzMrJPaDSePBX6Ve/wS4KqIeF1EnAWcCryoU5UbKXw4sZmZ\nWee0G052BVbkHh9DOoV9zZ+A/be1UiNNpeJwYmZm1inthpMVwEEAksYBs4Arc89PAjZ0pmojR7UK\nK1fC3XcXXRMzM7ORr91w8lPgw5KeQroI4BrgstzzhwGLO1S3EcNH7JiZmXVOu+Hk3cBDwG9JJ2N7\nfUSszz3/auAXHarbiDFzZrp3ODEzM9t27Z6EbSVwtKTJwOqI2Fi3yAuB+ztVuZFip51gv/0878TM\nzKwT2gonkr5S97jZoq8eboVGKk+KNTMz64x2r61zInAj8BegaTLpR9Uq/P73RdfCzMxs5Gs3nHwe\nGCAdsXM+8PWI8DEqpHBy/vnpOjujfFEAMzOzYWvrazQiTgH2Bj5KOhvsTZIukjRHg4zx9INKBR58\nEJb31TWZzczMOq/t3/gRsS4i5kXEM4FHAP8AzgGWSZrY6QqOFD6c2MzMrDO2dQBiExCk+Sejt706\nI9eBB8K4cZ4Ua2Zmtq3aDieSxksakPRL4Abg0cCbgQMiYnWnKzhSjB4NM2a458TMzGxbtXso8Tmk\ni/3dBHwFGMjOfWL4AoBmZmad0O7ROicDy4ElpIv+HdNoHmxE/Pu2V23kqVRg3ryia2FmZjaytRtO\nLiTNMbEGqtV0tM7atbDjjkXXxszMbGRq9/T1J3apHtuF2hE7CxfCYYcVWxczM7ORyqcL66BKJd17\nUqyZmdnwOZx00B57wG67eVKsmZnZtnA46TBfANDMzGzbOJx0WLXqYR0zM7Nt4XDSYbVznYSPaTIz\nMxsWh5MOq1Tg3nthpU9NZ2ZmNiwOJx1WO5zY807MzMyGx+Gkww4+GCSHEzMzs+EqTTiRdIqkpZLW\nSrpS0uNbXO/JkjZIWtDtOrZixx3TFYo9KdbMzGx4ShFOJL0Y+ARwBnA4cA1wiaQ9hlhvMvBV4NKu\nV7INvgCgmZnZ8JUinABzgXMj4sKI+CfpAoNrgFcPsd4XgG8AV3a5fm3xuU7MzMyGr/BwImksMBv4\nVa0sIoLUG3LkIOudBBwEvK/bdWxXtQqLF8NDDxVdEzMzs5Gn8HAC7AGMBlbUla8ApjZaQdJM4EPA\nyyJiU3er175KBTZsgGXLiq6JmZnZyFOGcNIWSaNIQzlnRMTiWnGBVdpK7XBiT4o1MzNr35iiKwCs\nBDYCU+rKpwC3N1h+EvA44LGSPpeVjQIkaT1wbET8X7OdzZ07l8mTJ29RNjAwwMDAwPBq38B++6Wj\ndq6/Ho4/vmObNTMz67p58+Yxb968LcpWrVrV0zooSnCedUlXAldFxGnZYwHLgU9HxMfqlhVwaN0m\nTgGeBrwz1hmQAAAYdklEQVQAWBYRaxvsYxYwf/78+cyaNasLr2JLj3kMHHkkfOELXd+VmZlZVy1Y\nsIDZs2cDzI6Irp+6oww9JwBnARdImg/8kXT0zgTgAgBJZwL7RMSrssmy1+ZXlnQH8GBEXNfTWg/C\nFwA0MzMbnlKEk4i4KDunyftJwzlXA3Mi4s5skanA/kXVbzgqFbj88qJrYWZmNvKUIpwARMQ5wDlN\nnjtpiHXfR8kOKa5W4dZbYfVqmDix6NqYmZmNHCPuaJ2RwkfsmJmZDY/DSZdUKuneZ4o1MzNrj8NJ\nl+yyC+y1l3tOzMzM2uVw0kW+xo6ZmVn7HE66yIcTm5mZtc/hpIuq1dRzUoLz3JmZmY0YDiddVKmk\nQ4lvu63ompiZmY0cDidd5MOJzczM2udw0kXTp8Po0Z4Ua2Zm1g6Hky4aNw4OOsjhxMzMrB0OJ13m\nI3bMzMza43DSZT7XiZmZWXscTrqsWoWlS2H9+qJrYmZmNjI4nHRZtQobN8KSJUXXxMzMbGRwOOky\nXwDQzMysPQ4nXbb33jBxoifFmpmZtcrhpMskT4o1MzNrh8NJD/hwYjMzs9Y5nPRA7QKAZmZmNjSH\nkx6oVOCOO+Dee4uuiZmZWfk5nPSALwBoZmbWOoeTHpg5M917aMfMzGxoDic9MGkS7LOPe07MzMxa\n4XDSI54Ua2Zm1hqHkx7xuU7MzMxa43DSI9UqLFwImzYVXRMzM7NyczjpkUoF1q6Fm28uuiZmZmbl\n5nDSIz6c2MzMrDUOJz0ybRqMHet5J2ZmZkNxOOmRMWPg4IMdTszMzIbicNJDvgCgmZnZ0BxOesiH\nE5uZmQ3N4aSHqlW48cZ01I6ZmZk15nDSQ9UqRMDixUXXxMzMrLwcTnqoUkn3HtoxMzNrzuGkh/bc\nE3bZxZNizczMBuNw0kOSJ8WamZkNxeGkx3x1YjMzs8E5nPSYz3ViZmY2OIeTHqtU4O67YeXKomti\nZmZWTg4nPeYLAJqZmQ3O4aTHZs5M9553YmZm1pjDSY/tuCMccIDDiZmZWTMOJwXwpFgzM7PmHE4K\n4HOdmJmZNedwUoBqFRYtgo0bi66JmZlZ+TicFKBahfXr0xWKzczMbEsOJwXwBQDNzMyaczgpwAEH\nwPjxnhRrZmbWiMNJAUaNSuc7cc+JmZnZ1hxOCuILAJqZmTVWmnAi6RRJSyWtlXSlpMcPsuzzJf1C\n0h2SVkn6g6Rje1nfbeVznZiZmTVWinAi6cXAJ4AzgMOBa4BLJO3RZJWjgV8AxwGzgN8AP5L0mB5U\ntyMqFbj5ZnjggaJrYmZmVi6lCCfAXODciLgwIv4JnAysAV7daOGImBsRH4+I+RGxOCLeBSwEntu7\nKm+b2gUAFy4sth5mZmZlU3g4kTQWmA38qlYWEQFcChzZ4jYETALu7kYdu8GHE5uZmTVWeDgB9gBG\nAyvqylcAU1vcxn8COwEXdbBeXbXbbrDHHg4nZmZm9cYUXYFtJemlwLuB50XEyqLr0w5PijUzM9ta\nGcLJSmAjMKWufApw+2ArSnoJ8EXghIj4TSs7mzt3LpMnT96ibGBggIGBgZYr3CmVCvztbz3frZmZ\nWVPz5s1j3rx5W5StWrWqp3VQmt5RLElXAldFxGnZYwHLgU9HxMearDMAnAe8OCJ+3MI+ZgHz58+f\nz6xZszpX+W3wkY/Ahz4E994LUtG1MTMza2zBggXMnj0bYHZELOj2/sow5wTgLOB1kl4p6RDgC8AE\n4AIASWdK+mpt4Wwo56vAfwB/kjQlu+3c+6oPX6UC990HK+pn25iZmfWxUoSTiLgIeBvwfuAvwGHA\nnIi4M1tkKrB/bpXXkSbRfg64NXf7ZK/q3Am1w4k9KdbMzGyzMsw5ASAizgHOafLcSXWPn9aTSnXZ\nwQen6+zccAMcc0zRtTEzMyuHUvSc9Kvx42HaNPecmJmZ5TmcFMyHE5uZmW3J4aRglYp7TszMzPIc\nTgpWrcKSJbBhQ9E1MTMzKweHk4JVq/DQQ7B0adE1MTMzKweHk4L5AoBmZmZbcjgp2L77woQJDidm\nZmY1DicFk1LviY/YMTMzSxxOSqBadc+JmZlZjcNJCfhcJ2ZmZps5nJRApQK3354uAmhmZtbvHE5K\nwBcANDMz28zhpARqhxN7aMfMzMzhpBR23hmmTnXPiZmZGTiclIYnxZqZmSUOJyXhCwCamZklDicl\nUes52bSp6JqYmZkVy+GkJKpVWLMGbr216JqYmZkVy+GkJHwBQDMzs8ThpCQOOgjGjPGkWDMzM4eT\nkhg7FqZPd8+JmZmZw0mJ+AKAZmZmDiel4nOdmJmZOZyUSqUCy5bBunVF18TMzKw4DiclUq2m85ws\nWlR0TczMzIrjcFIivgCgmZmZw0mpTJmSLgLoSbFmZtbPHE5KRPKkWDMzM4eTkvEFAM3MrN85nJSM\nz3ViZmb9zuGkZCoVuOuudDMzM+tHDiclU62me887MTOzfuVwUjIzZ6Z7hxMzM+tXDicls9NOsN9+\nnndiZmb9y+GkhDwp1szM+pnDSQlVKh7WMTOz/uVwUkLVKixcCBs3Fl0TMzOz3nM4KaFqNV2Z+Kab\niq6JmZlZ7zmclFDtAoCed2JmZv3I4aSEDjwQxo93ODEzs/7kcFJCo0fDjBmeFGtmZv3J4aSkfAFA\nMzPrVw4nJeVznZiZWb9yOCmpSiUdrbNmTdE1MTMz6y2Hk5KqXQBw4cJi62FmZtZrDicl5asTm5lZ\nv3I4Kandd4fddvO8EzMz6z8OJyXmSbFmZtaPxhRdAWuuUoHrriu6FkOLgDvvhEWLtrwtXw577ZXO\n2TJz5ub7ffaBUY7FZmbWhMNJiVWr8MMfpi9/qdi6bNoEt922dQCp3Vav3rzs3nunIDJtGtxxB3z3\nu3DjjWkbADvuCAcfvGVgqf3bwcXMzBxOSqxahVWrUq/EXnt1f38bN6bDlxuFjyVLYO3atJwE+++f\nwsQRR8BLX5r+PWMGTJ8OO+209bbXrYOlS9O2Fi7cfD9YcMmHFwcXM7P+UZpwIukU4G3AVOAa4C0R\n8adBln8q8AngkcBy4IMR8dUeVLVn8hcA7FQ42bABli1rHECWLk3PQzqF/kEHpaDw1KfCa1+7OYAc\ndFC69k8r5s2bx8DAAOPHwyGHpFu9ZsHlO9/ZOrjU6uDg0lytza133Oa95zbfvpUinEh6MSlovB74\nIzAXuERSJSJWNlh+GvBj4BzgpcAzgPMk3RoRv+xVvbttxozUS3H99fCUp7S+3oMPpp6ORgFk+fLU\nQwIwblwKHzNmwPHHb/7inzEDDjgAxo7d9tfQyh+QbgaXmTPTMNNIDy7r18P99w9+u+8+eOAB+MEP\n5nHHHQPsvDNMnkzD+x13LH6ocHviL8rec5tv30oRTkhh5NyIuBBA0snAs4FXAx9tsPwbgSUR8fbs\n8fWSjsq2s92Ekx12SFcobnSuk9WrYfHixgHkllvSPBWACRM2f2m/8IWbw8iMGbDvvqmHpMy6EVzy\n4aVbwSUiBYVmIWKwgNGofN26wfe3ww4waVIaUrvtNnjnOzcPwzUyZkwKKYMFmPr7RmXjxnW23czM\noAThRNJYYDbwoVpZRISkS4Ejm6z2RODSurJLgLO7UskCVavw61/DBz6wOXwsXgy33755mZ133vyF\n++Qnp/taCJk6dfv9hdzJ4FI/TLT77kP3Ugz23OrVmwNiI6NGpTDR6LbXXs2fmzQpvd/5xxMnbtnL\n9bznwcUXpyG6++5Lt1WrGt/Xl91yC1x77ebHq1bBQw81fx077NB+oKm/nzSp/CHZzHqr8HAC7AGM\nBlbUla8Aqk3Wmdpk+Z0ljY+IIX5njhyPf3wKJjfeuPnX/zOfuWUA2X337TeADFcngku9Wu9EfUDY\nc880EXiwAFF/mzCh++/Z2LHps7H77sPfRkQaJhwq5NTfL1q0ddlgYW3ixM1DTb26Qee2dc018C//\nkl5jo1utLXtdPtQ6tTYYNaoz953c1lD3V18Nr3nN0O9lr59rZd38+1J/3+5zndzWYM/dcQc9VYZw\n0is7AFw3Ek4ckvNv/wbHHtv4CBhIc0iWL+9tndqxatUqFixYUHQ1Gtpnn3Q75pjNZevXp96D1atT\ngJgwIbX9hAlpKGRb1HpWuq2bbT5qFOyyS7q1Y9OmNMz0wAOpbWtDXvnH69a19sXarKydZYcqyz9u\nZflRo1Yxfnxq80ZfTo3KW1kmH2Abfcm1Wt5s241CSy2cb9rU+nPtLLthQ2f2d++9q/jTnxYM+v60\n81yz97zZcs223co+69+r+vcn/z6VZZn16x/+7tyBHlDkW7QA2bDOGuAFEXFxrvwCYHJEPL/BOr8F\n5kfEW3NlJwJnR8SuTfbzUuAbna29mZlZX3lZRHyz2zspvOckIjZImg88HbgYQJKyx59ustoVwHF1\nZcdm5c1cArwMWAY8uA1VNjMz6zc7ANNI36VdV3jPCYCkFwEXACez+VDiE4BDIuJOSWcC+0TEq7Ll\npwF/Ix1K/BVSkPkkcHxE1E+UNTMzsxGk8J4TgIi4SNIewPuBKcDVwJyIuDNbZCqwf275ZZKeTTo6\n51TgZuA1DiZmZmYjXyl6TszMzMxqRvh5M83MzGx743BiZmZmpdIX4UTSKZKWSlor6UpJjy+6TiOR\npHdK+qOk+yStkPQDSZUGy71f0q2S1kj6paQZdc+Pl/Q5SSsl3S/pu5J6cN3lkU/S6ZI2STqrrtxt\n3kGS9pH0tay91ki6RtKsumXc5h0iaZSk/5G0JGvPRZL+u8FybvNhkvQUSRdLuiX7G/K8Bstsc/tK\n2lXSNyStknSPpPMkNTlTV3PbfTjJXVTwDOBw0hWPL8km4Fp7ngJ8BngC6WKLY4FfSNqxtoCkdwBv\nJl3E8QjgAVJ756/C8knStZNeABwN7AN8rxcvYCTLQvXrSZ/hfLnbvIMk7QJcDqwD5gCHAv8B3JNb\nxm3eWacDbwDeBBwCvB14u6Q31xZwm2+znUgHm7wJ2GqyaQfb95uk/zNPz5Y9Gji37dpGxHZ9A64E\nPpV7LNLRPW8vum4j/Ua69MAm4Khc2a3A3NzjnYG1wItyj9cBz88tU822c0TRr6msN2AicD3wL8Bv\ngLPc5l1r6w8Dvx1iGbd5Z9v8R8CX6sq+C1zoNu9Ke28CnldXts3tSwolm4DDc8vMAR4CprZTx+26\n5yR3UcFf1coitdZgFxW01u1CSuB3A0g6iHTYd7697wOuYnN7P450CHt+meuB5fg9GczngB9FxK/z\nhW7zrngu8GdJF2XDlwskvbb2pNu8K/4APF3STABJjwGeDPw0e+w276IOtu8TgXsi4i+5zV9K+p54\nQjt1KsV5TrpoOBcVtBZkZ/H9JPD7iLg2K55K+hA2au+p2b+nAOuzD36zZSxH0kuAx5L+ONRzm3fe\ndOCNpOHgD5K6uD8taV1EfA23eTd8mPTL/J+SNpKmHLwrIr6VPe82765Ote9UYItLBEbERkl30+Z7\nsL2HE+uec4BHkH7dWJdI2o8UAp8RERuKrk+fGAX8MSLenT2+RtKjSGew/lpx1dquvRh4KfAS4FpS\nGP+UpFuzQGh9Zrse1gFWAhtJiS9vCnB776uzfZD0WeB44KkRcVvuqdtJc3oGa+/bgXGSdh5kGdts\nNrAnsEDSBkkbgGOA0yStJ/1qcZt31m1A/eXLrwMOyP7tz3nnfRT4cER8JyL+ERHfIJ0B/J3Z827z\n7upU+94O1B+9MxrYjTbfg+06nGS/NGsXFQS2uKjgH4qq10iWBZN/BZ4WEcvzz0XEUtIHMN/eO5PG\nGmvtPZ80OSq/TJX0h3+wCzf2q0uBR5N+ST4mu/0Z+DrwmIhYgtu80y5n62HfKnAj+HPeJRNIPyTz\nNpF9R7nNu6uD7XsFsIukw3Obfzop+FzVbqW26xvwImAN8ErSIWrnAncBexZdt5F2Iw3l3EM6pHhK\n7rZDbpm3Z+37XNKX6g+BhcC4uu0sBZ5K6hm4HLis6Nc3Um5sfbSO27yz7fs40lEJ7wQOJg033A+8\nxG3etTY/nzSx8njgQOD5pLkLH3Kbd6yNdyL9uHksKfj9v+zx/p1sX9Ik5j8DjycN+18PfK3t+hbd\nYD16U94ELCMdFnUF8Lii6zQSb9kHemOD2yvrlnsv6bC0NaTLa8+oe3486XwpK7M/+t8B9ir69Y2U\nG/DrfDhxm3eljY8H/pq15z+AVzdYxm3eufbeCTgr++J7IPtSfB8wxm3esTY+psnf8K90sn1JR3F+\nHVhF+jH7JWBCu/X1hf/MzMysVLbrOSdmZmY28jicmJmZWak4nJiZmVmpOJyYmZlZqTicmJmZWak4\nnJiZmVmpOJyYmZlZqTicmJmZWak4nJj1MUm/kXRW0fXIk/RFSXdJ2ijpsKLr0ypJx0ja1ODCaGbW\npjFFV8DMrEbSs0jXwTqGdCrzlcXWqG0+5bZZBzicmFlHSRoFRAzv2hgzgNsior0rmJrZdsXDOmYF\ny4ZWPiXpI9lwxm2Szsg9f2A2XHBYrmxyVnZ09rg2pHCspAWS1ki6VNKeko6TdK2kVZK+IWmHuiqM\nkfQZSfdKulPS++vqN07SxyXdLGm1pCskHZN7/lWS7pH0XEn/AB4E9m/yWo+RdJWkByXdKunMLMwg\n6Xzg08AB2WtZ0mQbB0i6WNLdWX3+lvW4IGmUpPMkLcna4J+STq1b/3xJP5D0Tkm3Z3X/b0mjJX00\new9uknRig/fgxZIul7Q22+/Rzd9ZkHSUpN9ldbkxe58n5J5/k6Qbsu3dLumiwbZn1i8cTszK4ZXA\nauAI0qXL3yPp6bnnW+2FOIN0Fe4jgQOAi4BTgZeQrrR7LPCWunVOBDaQLnF+KvBWSa/JPf854AnA\ni0iXUv8O8DNJB+eWmZDV+zXAI0mXu9+CpH2AnwBXAYcBJ2fL/3e2yKnAe4CbgSlZfRo5BxgHHAU8\nCngHqe0g/U27CXgBcCjpyrYflHRC3Tb+BdgbeAowF3g/8GPgbtJ78AXg3KzOeR8FPka67PwVwI8k\n7dqokln7/IzUXo8CXky6hPxnsucfB3wqe/0VYA7wuyav2ay/FH0ZZ9986/cb8Bvgt3VlVwEfyv59\nIOlS54flnp+clR2dPT6GdPnzp+aWeUdWdmCu7PPAT+v2/fe6fZ9ZKyMFnA3A1Lplfgl8IPv3q7L9\nPGqI1/lB4Nq6sjcCq3KPTwOWDLGda4B3t9G+nwEuyj0+v34fwHXA/+UejyJdEv5Fde/B23LLjAaW\n18py78HO2eMvAZ+v289RwEOkcPV80iXldyr6M+ibb2W7ec6JWTn8te7xbcBew9jO33L/XgGsiYgb\n68rqeySurHt8Ban3RKRf/KOBG7LHNePYcrLq+oj4+xB1OyTbdt7lwERJ+0XEzUOsX/Np4POS5gCX\nAt+LiIdft6RTgJNIwWrHrK5/qdvGP+oeryDXdhGxSdJdbP0eXJlbZqOkP5N6aBp5DPBoSS/PldXa\n8CBSwFsOLJX0c+DnwA8iYm2T7Zn1DQ/rmJXDhrrHweb/n5uy+3w4GNvCdmKI7bZiIumX/izSl23t\ndiipl6OmZ1+oEfFl0pf7haTw9OcskCDpJaRhly8Bz8zqej4poOQ1apdtbat6E4FzSUNYtXY7jDSE\nszgiVgOHk4bcbiUNQV3jQ5HNHE7MRoI7s/u9c2WH07nDVp9Q9/hIYGFEBKnHYTQwJSKW1N22mlcy\nhOuybecdBdzfRq8JABFxS0R8MSJOAD4BvC576knA5RFxbkRcExFLgIObbqh9T6z9Q9JoYDZwbZNl\nFwCPiIilDdruoex1bIqIX0fE6aTwMo00H8asr3lYx6zkIuJBSVcCp0taRpos+j8NFlWDslYcIOnj\nwBdJX7ZvJk0SJSIWSvomcKGkt5HCyl6kL9BrIuJnbeznHOA0SZ8BPksa5nkvKVy0TNLZpImmNwC7\nAU9jc0BYCLxC0rGk86S8gjSM1fDIn2E4RdIiUtB6K7ALqWfm4erl/v0R4Irs9Z4HPECaLPyMiHiL\npGcD00mTYO8Bnp2tf32H6mo2YjmcmBWvlR6QV5O+4P5M+vJ6O/CLYWyn0b4vJM3N+CNpCOfsiDgv\nt8yJpCNKPg7sS5prciXwo7Z2FHGrpONJwy5Xk46M+RJpomw7RpPCzX7AfaSg8tbsuXNJR9J8K3tt\n80hHGx03VPVaLDs9uz0GWAQ8NyLubrRORPwtO+T6g6QAImAx8O1skXuBfycdYbUDKVi9JCKuG6Ku\nZts9pZ5bMzNrRtKBpN6XwyOifvKymXWY55yYmbVmuMNmZtYmhxMzs9a4m9msRzysY2ZmZqXinhMz\nMzMrFYcTMzMzKxWHEzMzMysVhxMzMzMrFYcTMzMzKxWHEzMzMysVhxMzMzMrFYcTMzMzKxWHEzMz\nMyuV/w+oxFvLPqaLtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x98365c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def sum_to_one(x):\n",
    "    #x = x-min(x)\n",
    "    #x = abs(x)\n",
    "    return x/np.sum(x)\n",
    "sample_size = [5, 100, 200, 300, 400, 500, 600,700,800,900,1000]\n",
    "true = []\n",
    "\n",
    "for i in range(1, 6):\n",
    "    temp = cal_first_order_term_prob(i)\n",
    "    #print temp\n",
    "    true.append(temp* cal_prob_each_state(i))\n",
    "#true = cal_first_order_term_prob()\n",
    "print true\n",
    "total_err = []\n",
    "for t in range(0, 10000):\n",
    "    error = []\n",
    "    for i in range(0, len(sample_size)):\n",
    "        seq_samples = gen_first_order_samples(sample_size[i])\n",
    "        vec_samples = gen_vecsamples_from_seqsamples(seq_samples)\n",
    "        sigma = [1, 0, 0, 0, 0, 0, 1]\n",
    "        w = np.random.randn(7)\n",
    "        for j in range(0, len(seq_samples)):\n",
    "            for k in range(0, len(seq_samples[j])-1):\n",
    "                temp1 = np.outer(vec_samples[j][k], vec_samples[j][k+1])\n",
    "                temp2 = np.outer(vec_samples[j][k], vec_samples[j][k])\n",
    "                if seq_samples[j][k+1] == 0 or seq_samples[j][k+1] == 6:\n",
    "                    temp3 = vec_samples[j][k]\n",
    "                else: temp3 = 0.\n",
    "                w = w+0.01*(np.dot((temp1 - temp2), w)+temp3)\n",
    "        w = sum_to_one(w)\n",
    "        #print w\n",
    "        error.append(np.mean((true - w[1:6])**2))\n",
    "    #print '******'\n",
    "    total_err.append(np.asarray(error))\n",
    "total =  sum_to_one(np.mean(total_err, axis = 0))\n",
    "from matplotlib import pyplot as plt\n",
    "plt.title('MSE w.r.t differnent number of samples')\n",
    "plt.xlabel('number of samples')\n",
    "plt.ylabel('MSE')\n",
    "#plt.axis([1, 1000, 0, 1])\n",
    "plt.plot(sample_size, total)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Reference\n",
    "[1] Barnard E. Temporal-difference methods and Markov models[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1993, 23(2): 357-365."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}