Estimating Parameters1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# Estimating Parameters in Chaotic Systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Suppose we obseve a dynamical system and we wish to estimate its\n",
    "parameters. One way of doing this is to use a Monte Carlo\n",
    "method. But if the system is chaotic then this approach is highly\n",
    "unlikely to work as the small changes in the proposals for the\n",
    "parameters will result in large changes in the path.\n",
    "\n",
    "Here's a [well-known chaotic dynamical system](https://en.wikipedia.org/wiki/Lorenz_system):\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\mathrm{d}x}{\\mathrm{d}t} &= \\alpha (y - x), \\\\\n",
    "\\frac{\\mathrm{d}y}{\\mathrm{d}t} &= x (\\rho - z) - y, \\\\\n",
    "\\frac{\\mathrm{d}z}{\\mathrm{d}t} &= x y - \\beta z.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "And here is its representation in [libbi](http://libbi.org) together with some\n",
    "instructions to generate noisy values of the $x$ variable (we use\n",
    "the noisy values later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model Lorenz {\n",
      "  const rho = 45.92\n",
      "  const beta = 4.0\n",
      "  const alpha = 16.0\n",
      "\n",
      "  state X, Y, Z\n",
      "  obs X_obs\n",
      "\n",
      "  sub initial {\n",
      "    X ~ log_normal(log(1.0), 0.00002)\n",
      "    Y ~ log_normal(log(1.0), 0.00002)\n",
      "    Z ~ log_normal(log(1.0), 0.00002)\n",
      "  }\n",
      "\n",
      "  sub transition(delta = 0.0001) {\n",
      "    ode {\n",
      "      dX/dt = alpha * (Y - X)\n",
      "      dY/dt = X * (rho - Z) - Y\n",
      "      dZ/dt = X * Y - beta * Z\n",
      "    }\n",
      "  }\n",
      "\n",
      "  sub observation {\n",
      "    X_obs ~ normal(X, 0.2)\n",
      "  }\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(readr)\n",
    "model_file_name <- \"LorenzGenerate.bi\"\n",
    "writeLines(read_file(model_file_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "library('rbi')\n",
    "library(ggplot2)\n",
    "\n",
    "Lorenz <- bi_model(model_file_name)\n",
    "\n",
    "T <- 10.0\n",
    "nObs <- 100\n",
    "init_parameters <- list(X = 1, Y = 1, Z = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "synthetic_dataset <- bi_generate_dataset(end_time=T, model=Lorenz,\n",
    "                                         init=init_parameters,\n",
    "                                         noutputs = nObs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>X.time</th><th scope=col>X.value</th><th scope=col>Y.time</th><th scope=col>Y.value</th><th scope=col>Z.time</th><th scope=col>Z.value</th><th scope=col>X_obs.time</th><th scope=col>X_obs.value</th><th scope=col>clock</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>96</th><td> 9.5      </td><td>-14.764658</td><td> 9.5      </td><td>-19.536008</td><td> 9.5      </td><td>39.76409  </td><td> 9.5      </td><td>-15.044416</td><td>96309     </td></tr>\n",
       "\t<tr><th scope=row>97</th><td> 9.6      </td><td>-17.611354</td><td> 9.6      </td><td>-14.413804</td><td> 9.6      </td><td>53.75645  </td><td> 9.6      </td><td>-17.397395</td><td>96309     </td></tr>\n",
       "\t<tr><th scope=row>98</th><td> 9.7      </td><td> -9.657817</td><td> 9.7      </td><td> -5.782971</td><td> 9.7      </td><td>45.82000  </td><td> 9.7      </td><td>-10.099852</td><td>96309     </td></tr>\n",
       "\t<tr><th scope=row>99</th><td> 9.8      </td><td> -8.019435</td><td> 9.8      </td><td> -9.629767</td><td> 9.8      </td><td>35.46636  </td><td> 9.8      </td><td> -8.105424</td><td>96309     </td></tr>\n",
       "\t<tr><th scope=row>100</th><td> 9.9      </td><td>-14.220221</td><td> 9.9      </td><td>-19.872234</td><td> 9.9      </td><td>37.34805  </td><td> 9.9      </td><td>-14.511350</td><td>96309     </td></tr>\n",
       "\t<tr><th scope=row>101</th><td>10.0      </td><td>-18.793750</td><td>10.0      </td><td>-15.760650</td><td>10.0      </td><td>55.08816  </td><td>10.0      </td><td>-18.686090</td><td>96309     </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllllll}\n",
       "  & X.time & X.value & Y.time & Y.value & Z.time & Z.value & X\\_obs.time & X\\_obs.value & clock\\\\\n",
       "\\hline\n",
       "\t96 &  9.5       & -14.764658 &  9.5       & -19.536008 &  9.5       & 39.76409   &  9.5       & -15.044416 & 96309     \\\\\n",
       "\t97 &  9.6       & -17.611354 &  9.6       & -14.413804 &  9.6       & 53.75645   &  9.6       & -17.397395 & 96309     \\\\\n",
       "\t98 &  9.7       &  -9.657817 &  9.7       &  -5.782971 &  9.7       & 45.82000   &  9.7       & -10.099852 & 96309     \\\\\n",
       "\t99 &  9.8       &  -8.019435 &  9.8       &  -9.629767 &  9.8       & 35.46636   &  9.8       &  -8.105424 & 96309     \\\\\n",
       "\t100 &  9.9       & -14.220221 &  9.9       & -19.872234 &  9.9       & 37.34805   &  9.9       & -14.511350 & 96309     \\\\\n",
       "\t101 & 10.0       & -18.793750 & 10.0       & -15.760650 & 10.0       & 55.08816   & 10.0       & -18.686090 & 96309     \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | X.time | X.value | Y.time | Y.value | Z.time | Z.value | X_obs.time | X_obs.value | clock | \n",
       "|---|---|---|---|---|---|\n",
       "| 96 |  9.5       | -14.764658 |  9.5       | -19.536008 |  9.5       | 39.76409   |  9.5       | -15.044416 | 96309      | \n",
       "| 97 |  9.6       | -17.611354 |  9.6       | -14.413804 |  9.6       | 53.75645   |  9.6       | -17.397395 | 96309      | \n",
       "| 98 |  9.7       |  -9.657817 |  9.7       |  -5.782971 |  9.7       | 45.82000   |  9.7       | -10.099852 | 96309      | \n",
       "| 99 |  9.8       |  -8.019435 |  9.8       |  -9.629767 |  9.8       | 35.46636   |  9.8       |  -8.105424 | 96309      | \n",
       "| 100 |  9.9       | -14.220221 |  9.9       | -19.872234 |  9.9       | 37.34805   |  9.9       | -14.511350 | 96309      | \n",
       "| 101 | 10.0       | -18.793750 | 10.0       | -15.760650 | 10.0       | 55.08816   | 10.0       | -18.686090 | 96309      | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "    X.time X.value    Y.time Y.value    Z.time Z.value  X_obs.time X_obs.value\n",
       "96   9.5   -14.764658  9.5   -19.536008  9.5   39.76409  9.5       -15.044416 \n",
       "97   9.6   -17.611354  9.6   -14.413804  9.6   53.75645  9.6       -17.397395 \n",
       "98   9.7    -9.657817  9.7    -5.782971  9.7   45.82000  9.7       -10.099852 \n",
       "99   9.8    -8.019435  9.8    -9.629767  9.8   35.46636  9.8        -8.105424 \n",
       "100  9.9   -14.220221  9.9   -19.872234  9.9   37.34805  9.9       -14.511350 \n",
       "101 10.0   -18.793750 10.0   -15.760650 10.0   55.08816 10.0       -18.686090 \n",
       "    clock\n",
       "96  96309\n",
       "97  96309\n",
       "98  96309\n",
       "99  96309\n",
       "100 96309\n",
       "101 96309"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "synthetic_data <- bi_read(synthetic_dataset)\n",
    "synthetic_df <- as.data.frame(synthetic_data)\n",
    "tail(synthetic_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Saving 10 x 5 in image\n"
     ]
    }
   ],
   "source": [
    "p <- ggplot(synthetic_df, aes(X.time)) +\n",
    "    geom_path(aes(y = X.value, colour=\"alpha 16.0\")) +\n",
    "    theme(legend.position=\"bottom\") +\n",
    "    ggtitle(\"Lorenz\") +\n",
    "    theme(plot.title = element_text(hjust = 0.5)) +\n",
    "    xlab(\"Time\") +\n",
    "    ylab(\"X Value\")\n",
    "ggsave(filename = \"diagrams/xpath.svg\", plot = p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "![title](diagrams/xpath.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We can check that the system becomes chaotic, in the sense that small\n",
    "changes in initial conditions lead to qualitatively different\n",
    "behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "path0 <- ggplot() +\n",
    "    theme(legend.position=\"bottom\") +\n",
    "    ggtitle(\"Lorenz\") +\n",
    "    theme(plot.title = element_text(hjust = 0.5)) +\n",
    "    xlab(\"Time\") +\n",
    "    ylab(\"Value\")\n",
    "\n",
    "\n",
    "set.seed(42)\n",
    "\n",
    "T <- 20.0\n",
    "\n",
    "for (i in c(\"red\", \"blue\", \"green\")) {\n",
    "    init_parameters <- list(X = 1 + rnorm(1,0.0,0.01),\n",
    "                            Y = 1 + rnorm(1,0.0,0.01),\n",
    "                            Z = 1 + rnorm(1,0.0,0.01))\n",
    "\n",
    "    synthetic_dataset <- bi_generate_dataset(end_time=T, model=Lorenz,\n",
    "                                             init=init_parameters,\n",
    "                                             noutputs = nObs)\n",
    "\n",
    "    synthetic_data <- bi_read(synthetic_dataset)\n",
    "    synthetic_df <- as.data.frame(synthetic_data)\n",
    "\n",
    "    path0 <- path0 +\n",
    "        geom_line(data = synthetic_df, aes(x = X.time, y = X.value), color = i)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Saving 10 x 5 in image\n"
     ]
    }
   ],
   "source": [
    "ggsave(filename = \"diagrams/xpath4.svg\", plot = path0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "![title](diagrams/xpath4.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Parameters as State"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Alternatively we can model the parameter as part of the state and\n",
    "assume that it undergoes Brownian Motion. This seems reasonable: the\n",
    "further we go into the future, the less certain we are about its\n",
    "value. An improvement might be to model it as an Ornstein-Uhlenbeck\n",
    "process which is mean-reverting - after all we don't expect the\n",
    "parameter to take on arbitrarily large or small values but\n",
    "let's learn to walk before we learn to run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Since we expect our parameter to positive let's model it as geometric\n",
    "Brownian Motion.\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{d}\\alpha = \\alpha\\sigma\\mathrm{d}W_t\n",
    "\\end{equation}\n",
    "\n",
    "By [Itô's lemma](https://en.wikipedia.org/wiki/Itô%27s_lemma#Geometric_Brownian_motion),\n",
    "we have\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{d}(\\log \\alpha) = -\\frac{\\sigma^2}{2}\\mathrm{d}t + \\sigma\\mathrm{d}W_t\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We can model this in libbi as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model LorenzState {\n",
      "  const rho = 45.92\n",
      "  const beta = 4\n",
      "\n",
      "  const h         = 0.1;    // time step\n",
      "  const delta_abs = 1.0e-3; // absolute error tolerance\n",
      "  const delta_rel = 1.0e-6; // relative error tolerance\n",
      "\n",
      "  state X, Y, Z\n",
      "  state ln_alpha\n",
      "\n",
      "  param mu, sigma\n",
      "\n",
      "  noise w\n",
      "\n",
      "  obs X_obs\n",
      "\n",
      "  sub parameter {\n",
      "    mu ~ uniform(12.0, 20.0)\n",
      "    sigma ~ uniform(0.0, 0.5)\n",
      "  }\n",
      "\n",
      "  sub proposal_parameter {\n",
      "     mu ~ truncated_gaussian(mu, 0.02, 12.0, 20.0);\n",
      "     sigma ~ truncated_gaussian(sigma, 0.01, 0.0, 0.5);\n",
      "   }\n",
      "\n",
      "  sub initial {\n",
      "    X ~ log_normal(log(1.0), 0.2)\n",
      "    Y ~ log_normal(log(1.0), 0.2)\n",
      "    Z ~ log_normal(log(1.0), 0.2)\n",
      "    ln_alpha ~ gaussian(log(mu), sigma)\n",
      "  }\n",
      "\n",
      "  sub transition(delta = h) {\n",
      "    w ~ normal (0.0, sqrt(h))\n",
      "    ode(h = h, atoler = delta_abs, rtoler = delta_rel, alg =  'RK4(3)') {\n",
      "      dX/dt = exp(ln_alpha) * (Y - X)\n",
      "      dY/dt = X * (rho - Z) - Y\n",
      "      dZ/dt = X * Y - beta * Z\n",
      "      dln_alpha/dt = -sigma * sigma / 2 - sigma * w / h\n",
      "    }\n",
      "  }\n",
      "\n",
      "  sub observation {\n",
      "    X_obs ~ normal(X, 0.2)\n",
      "  }\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model_file_name <- \"LorenzState.bi\"\n",
    "writeLines(read_file(model_file_name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We can then run the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File /private/var/folders/h1/bkwn2ct12hvb6__dzrk5gwx80000gn/T/RtmpQ7zb3k/LorenzState1176a1bafc6c5/LorenzState_output1176a1d4a68c2.nc (NC_FORMAT_NETCDF4):\n",
      "\n",
      "     13 variables (excluding dimension variables):\n",
      "        double time[nr]   (Contiguous storage)  \n",
      "        double w[np,nr]   (Contiguous storage)  \n",
      "        double X[np,nr]   (Contiguous storage)  \n",
      "        double Y[np,nr]   (Contiguous storage)  \n",
      "        double Z[np,nr]   (Contiguous storage)  \n",
      "        double ln_alpha[np,nr]   (Contiguous storage)  \n",
      "        double mu[]   (Contiguous storage)  \n",
      "        double sigma[]   (Contiguous storage)  \n",
      "        8 byte int clock[]   (Contiguous storage)  \n",
      "        int ancestor[np,nr]   (Contiguous storage)  \n",
      "        double logweight[np,nr]   (Contiguous storage)  \n",
      "        double loglikelihood[]   (Contiguous storage)  \n",
      "        double X_obs[np,nr]   (Contiguous storage)  \n",
      "\n",
      "     2 dimensions:\n",
      "        np  Size:8192\n",
      "        nr  Size:101\n",
      "\n",
      "    3 global attributes:\n",
      "        libbi_schema: ParticleFilter\n",
      "        libbi_schema_version: 1\n",
      "        libbi_version: 1.4.1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>Min.</th><th scope=col>1st Qu.</th><th scope=col>Median</th><th scope=col>Mean</th><th scope=col>3rd Qu.</th><th scope=col>Max.</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>mu</th><td>14.1060017</td><td>14.1060017</td><td>14.1060017</td><td>14.1060017</td><td>14.1060017</td><td>14.1060017</td></tr>\n",
       "\t<tr><th scope=row>sigma</th><td> 0.2652258</td><td> 0.2652258</td><td> 0.2652258</td><td> 0.2652258</td><td> 0.2652258</td><td> 0.2652258</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llllll}\n",
       "  & Min. & 1st Qu. & Median & Mean & 3rd Qu. & Max.\\\\\n",
       "\\hline\n",
       "\tmu & 14.1060017 & 14.1060017 & 14.1060017 & 14.1060017 & 14.1060017 & 14.1060017\\\\\n",
       "\tsigma &  0.2652258 &  0.2652258 &  0.2652258 &  0.2652258 &  0.2652258 &  0.2652258\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. | \n",
       "|---|---|\n",
       "| mu | 14.1060017 | 14.1060017 | 14.1060017 | 14.1060017 | 14.1060017 | 14.1060017 | \n",
       "| sigma |  0.2652258 |  0.2652258 |  0.2652258 |  0.2652258 |  0.2652258 |  0.2652258 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "      Min.       1st Qu.    Median     Mean       3rd Qu.    Max.      \n",
       "mu    14.1060017 14.1060017 14.1060017 14.1060017 14.1060017 14.1060017\n",
       "sigma  0.2652258  0.2652258  0.2652258  0.2652258  0.2652258  0.2652258"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LorenzState <- bi_model(model_file_name)\n",
    "\n",
    "bi_state_model <- libbi(model=LorenzState)\n",
    "bi_state <- filter(bi_state_model,\n",
    "                   nparticles = 8192,\n",
    "                   nthreads = 1,\n",
    "                   end_time = T,\n",
    "                   obs = synthetic_dataset,\n",
    "                   init = init_parameters,\n",
    "                   ess_rel = 1,\n",
    "                   sample_obs = TRUE)\n",
    "\n",
    "bi_file_summary(bi_state$output_file_name)\n",
    "bi_state\n",
    "summary(bi_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "And then after some manipulation, draw the path of the \"state\" parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "output <- bi_read(bi_state)\n",
    "logw <- xtabs(value ~ time + np, data = output$logweight, addNA = TRUE)\n",
    "X <- output$X$value\n",
    "Y <- output$Y$value\n",
    "Z <- output$Z$value\n",
    "A <- output$ln_alpha$value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "log2normw <- function(lw){\n",
    "  w <- exp(lw - max(lw))\n",
    "  return(w / sum(w))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "w = t(apply(X=logw, MARGIN=1, FUN=log2normw))\n",
    "Xmeans = apply(X = X*w, MARGIN=1, FUN=sum)\n",
    "Ymeans = apply(X = X*w, MARGIN=1, FUN=sum)\n",
    "Zmeans = apply(X = Z*w, MARGIN=1, FUN=sum)\n",
    "Ameans = apply(X = A*w, MARGIN=1, FUN=sum)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "synthetic_data <- bi_read(synthetic_dataset)\n",
    "X_original <- synthetic_data$X$value\n",
    "Y_original <- synthetic_data$Y$value\n",
    "Z_original <- synthetic_data$Z$value\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "synthetic_df <- as.data.frame(synthetic_data)\n",
    "synthetic_df$Xmeans <- Xmeans\n",
    "synthetic_df$Ymeans <- Ymeans\n",
    "synthetic_df$Zmeans <- Zmeans\n",
    "synthetic_df$Ameans <- Ameans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "simpleWarning in if (class(res$value) == \"help_files_with_topic\") {: the condition has length > 1 and only the first element will be used\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pAmeans <- ggplot(synthetic_df, aes(X.time)) +\n",
    "           geom_path(aes(y = exp(Ameans), colour=\"Ameans\")) +\n",
    "           theme(legend.position=\"bottom\") +\n",
    "           ggtitle(\"Lorenz\") +\n",
    "           theme(plot.title = element_text(hjust = 0.5)) +\n",
    "           ylim(0.0, max(exp(Ameans))) +\n",
    "           xlab(\"Time\") +\n",
    "           ylab(\"Value\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Saving 10 x 5 in image\n"
     ]
    }
   ],
   "source": [
    "ggsave(filename = \"diagrams/xpath5.svg\", plot = pAmeans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "![title](diagrams/xpath5.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We can try inferring the parameter for the different chaotic solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dataset_list <- list()\n",
    "parameters_list <- list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for (i in c(1,2,3)) {\n",
    "    init_parameters <- list(X = 1 + rnorm(1,0.0,0.01),\n",
    "                            Y = 1 + rnorm(1,0.0,0.01),\n",
    "                            Z = 1 + rnorm(1,0.0,0.01))\n",
    "\n",
    "    parameters_list[[i]] <- init_parameters\n",
    "    synthetic_dataset <- bi_generate_dataset(end_time=T, model=Lorenz,\n",
    "                                             init=init_parameters,\n",
    "                                             noutputs = nObs)\n",
    "\n",
    "    dataset_list[[i]] <- synthetic_dataset\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "X_list <- list()\n",
    "Y_list <- list()\n",
    "Z_list <- list()\n",
    "A_list <- list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    " for (i in c(1,2,3)) {\n",
    "    bi_state <- filter(bi_state_model, nparticles = 8192, nthreads = 1, end_time = T, obs = dataset_list[[i]], init = parameters_list[[i]], ess_rel = 1, sample_obs = TRUE)\n",
    "    output <- bi_read(bi_state)\n",
    "    logw <- xtabs(value ~ time + np, data = output$logweight, addNA = TRUE)\n",
    "    w = t(apply(X=logw, MARGIN=1, FUN=log2normw))\n",
    "    X <- output$X$value\n",
    "    Y <- output$Y$value\n",
    "    Z <- output$Z$value\n",
    "    A <- output$ln_alpha$value\n",
    "    X_list[[i]] = apply(X = X*w, MARGIN=1, FUN=sum)\n",
    "    Y_list[[i]] = apply(X = X*w, MARGIN=1, FUN=sum)\n",
    "    Z_list[[i]] = apply(X = Z*w, MARGIN=1, FUN=sum)\n",
    "    A_list[[i]] = apply(X = A*w, MARGIN=1, FUN=sum)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "simpleWarning in if (class(res$value) == \"help_files_with_topic\") {: the condition has length > 1 and only the first element will be used\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "path2 <- ggplot() +\n",
    "    theme(legend.position=\"bottom\") +\n",
    "    ggtitle(\"Lorenz\") +\n",
    "    theme(plot.title = element_text(hjust = 0.5)) +\n",
    "    xlab(\"Time\") +\n",
    "    ylab(\"Value\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for (i in c(1,2,3)) {\n",
    "        synthetic_data <- bi_read(dataset_list[[i]])\n",
    "        synthetic_df <- as.data.frame(synthetic_data)\n",
    "        synthetic_df$Ameans <- exp(A_list[[i]])\n",
    "        path2 <- path2 + geom_line(data = synthetic_df,\n",
    "                                   aes(x = X.time, y = Ameans), color = \"blue\")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Saving 10 x 5 in image\n"
     ]
    }
   ],
   "source": [
    "ggsave(filename = \"diagrams/xpath7.svg\", plot = path2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "![title](diagrams/xpath7.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "And if we take the mean of the tails then we get pretty decent\n",
    "estimates of the parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x <- list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for (i in c(1:3)) {\n",
    "        x[[i]] <- tail(exp(A_list[[i]]), n = 50)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 15.92264\n",
      "[1] 16.00402\n",
      "[1] 16.19887\n"
     ]
    }
   ],
   "source": [
    "for (i in 1:3) print(mean(x[[i]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 1.170018\n",
      "[1] 0.6070526\n",
      "[1] 0.763862\n"
     ]
    }
   ],
   "source": [
    "for (i in 1:3) print(sd(x[[i]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook can be downloaded from"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 3.5.1 (Juniper)",
   "language": "R",
   "name": "juniper_r3.5.1"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".R",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.1"
  },
  "name": "Estimating Parameters1.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}