{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: distrax in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (0.0.1)\n", "Collecting distrax\n", " Downloading distrax-0.1.2-py3-none-any.whl (272 kB)\n", "\u001b[K |████████████████████████████████| 272 kB 6.9 MB/s eta 0:00:01\n", "\u001b[?25hRequirement already satisfied: jax>=0.1.55 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from distrax) (0.2.11)\n", "Requirement already satisfied: absl-py>=0.9.0 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from distrax) (0.12.0)\n", "Requirement already satisfied: chex>=0.0.7 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from distrax) (0.0.8)\n", "Requirement already satisfied: jaxlib>=0.1.67 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from distrax) (0.1.70)\n", "Requirement already satisfied: numpy>=1.18.0 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from distrax) (1.19.5)\n", "Collecting tensorflow-probability>=0.15.0\n", " Using cached tensorflow_probability-0.16.0-py2.py3-none-any.whl (6.3 MB)\n", "Requirement already satisfied: six in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from absl-py>=0.9.0->distrax) (1.15.0)\n", "Requirement already satisfied: dm-tree>=0.1.5 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from chex>=0.0.7->distrax) (0.1.6)\n", "Requirement already satisfied: toolz>=0.9.0 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from chex>=0.0.7->distrax) (0.11.1)\n", "Requirement already satisfied: opt-einsum in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from jax>=0.1.55->distrax) (3.3.0)\n", "Requirement already satisfied: flatbuffers<3.0,>=1.12 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from jaxlib>=0.1.67->distrax) (1.12)\n", "Requirement already satisfied: scipy in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from jaxlib>=0.1.67->distrax) (1.6.3)\n", "Requirement already satisfied: cloudpickle>=1.3 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from tensorflow-probability>=0.15.0->distrax) (1.6.0)\n", "Requirement already satisfied: decorator in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from tensorflow-probability>=0.15.0->distrax) (4.4.2)\n", "Requirement already satisfied: gast>=0.3.2 in /opt/anaconda3/envs/spyder-dev/lib/python3.9/site-packages (from tensorflow-probability>=0.15.0->distrax) (0.4.0)\n", "Installing collected packages: tensorflow-probability, distrax\n", " Attempting uninstall: tensorflow-probability\n", " Found existing installation: tensorflow-probability 0.13.0\n", " Uninstalling tensorflow-probability-0.13.0:\n", " Successfully uninstalled tensorflow-probability-0.13.0\n", " Attempting uninstall: distrax\n", " Found existing installation: distrax 0.0.1\n", " Uninstalling distrax-0.0.1:\n", " Successfully uninstalled distrax-0.0.1\n", "\u001b[31mERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.\n", "jsl 0.0.0 requires dataclasses, which is not installed.\u001b[0m\n", "Successfully installed distrax-0.1.2 tensorflow-probability-0.16.0\n", "\u001b[33mWARNING: You are using pip version 21.2.4; however, version 22.0.4 is available.\n", "You should consider upgrading via the '/opt/anaconda3/envs/spyder-dev/bin/python -m pip install --upgrade pip' command.\u001b[0m\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "# meta-data does not work yet in VScode\n", "# https://github.com/microsoft/vscode-jupyter/issues/1121\n", "\n", "{\n", " \"tags\": [\n", " \"hide-cell\"\n", " ]\n", "}\n", "\n", "\n", "### Install necessary libraries\n", "\n", "try:\n", " import jax\n", "except:\n", " # For cuda version, see https://github.com/google/jax#installation\n", " %pip install --upgrade \"jax[cpu]\" \n", " import jax\n", "\n", "try:\n", " import distrax\n", "except:\n", " %pip install --upgrade distrax\n", " import distrax\n", "\n", "try:\n", " import jsl\n", "except:\n", " %pip install git+https://github.com/probml/jsl\n", " import jsl\n", "\n", "try:\n", " import rich\n", "except:\n", " %pip install rich\n", " import rich\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "{\n", " \"tags\": [\n", " \"hide-cell\"\n", " ]\n", "}\n", "\n", "\n", "### Import standard libraries\n", "\n", "import abc\n", "from dataclasses import dataclass\n", "import functools\n", "import itertools\n", "\n", "from typing import Any, Callable, NamedTuple, Optional, Union, Tuple\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "\n", "import jax\n", "import jax.numpy as jnp\n", "from jax import lax, vmap, jit, grad\n", "from jax.scipy.special import logit\n", "from jax.nn import softmax\n", "from functools import partial\n", "from jax.random import PRNGKey, split\n", "\n", "import inspect\n", "import inspect as py_inspect\n", "import rich\n", "from rich import inspect as r_inspect\n", "from rich import print as r_print\n", "\n", "def print_source(fname):\n", " r_print(py_inspect.getsource(fname))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{math}\n", "\n", "\\newcommand\\floor[1]{\\lfloor#1\\rfloor}\n", "\n", "\\newcommand{\\real}{\\mathbb{R}}\n", "\n", "% Numbers\n", "\\newcommand{\\vzero}{\\boldsymbol{0}}\n", "\\newcommand{\\vone}{\\boldsymbol{1}}\n", "\n", "% Greek https://www.latex-tutorial.com/symbols/greek-alphabet/\n", "\\newcommand{\\valpha}{\\boldsymbol{\\alpha}}\n", "\\newcommand{\\vbeta}{\\boldsymbol{\\beta}}\n", "\\newcommand{\\vchi}{\\boldsymbol{\\chi}}\n", "\\newcommand{\\vdelta}{\\boldsymbol{\\delta}}\n", "\\newcommand{\\vDelta}{\\boldsymbol{\\Delta}}\n", "\\newcommand{\\vepsilon}{\\boldsymbol{\\epsilon}}\n", "\\newcommand{\\vzeta}{\\boldsymbol{\\zeta}}\n", "\\newcommand{\\vXi}{\\boldsymbol{\\Xi}}\n", "\\newcommand{\\vell}{\\boldsymbol{\\ell}}\n", "\\newcommand{\\veta}{\\boldsymbol{\\eta}}\n", "%\\newcommand{\\vEta}{\\boldsymbol{\\Eta}}\n", "\\newcommand{\\vgamma}{\\boldsymbol{\\gamma}}\n", "\\newcommand{\\vGamma}{\\boldsymbol{\\Gamma}}\n", "\\newcommand{\\vmu}{\\boldsymbol{\\mu}}\n", "\\newcommand{\\vmut}{\\boldsymbol{\\tilde{\\mu}}}\n", "\\newcommand{\\vnu}{\\boldsymbol{\\nu}}\n", "\\newcommand{\\vkappa}{\\boldsymbol{\\kappa}}\n", "\\newcommand{\\vlambda}{\\boldsymbol{\\lambda}}\n", "\\newcommand{\\vLambda}{\\boldsymbol{\\Lambda}}\n", "\\newcommand{\\vLambdaBar}{\\overline{\\vLambda}}\n", "%\\newcommand{\\vnu}{\\boldsymbol{\\nu}}\n", "\\newcommand{\\vomega}{\\boldsymbol{\\omega}}\n", "\\newcommand{\\vOmega}{\\boldsymbol{\\Omega}}\n", "\\newcommand{\\vphi}{\\boldsymbol{\\phi}}\n", "\\newcommand{\\vvarphi}{\\boldsymbol{\\varphi}}\n", "\\newcommand{\\vPhi}{\\boldsymbol{\\Phi}}\n", "\\newcommand{\\vpi}{\\boldsymbol{\\pi}}\n", "\\newcommand{\\vPi}{\\boldsymbol{\\Pi}}\n", "\\newcommand{\\vpsi}{\\boldsymbol{\\psi}}\n", "\\newcommand{\\vPsi}{\\boldsymbol{\\Psi}}\n", "\\newcommand{\\vrho}{\\boldsymbol{\\rho}}\n", "\\newcommand{\\vtheta}{\\boldsymbol{\\theta}}\n", "\\newcommand{\\vthetat}{\\boldsymbol{\\tilde{\\theta}}}\n", "\\newcommand{\\vTheta}{\\boldsymbol{\\Theta}}\n", "\\newcommand{\\vsigma}{\\boldsymbol{\\sigma}}\n", "\\newcommand{\\vSigma}{\\boldsymbol{\\Sigma}}\n", "\\newcommand{\\vSigmat}{\\boldsymbol{\\tilde{\\Sigma}}}\n", "\\newcommand{\\vsigmoid}{\\vsigma}\n", "\\newcommand{\\vtau}{\\boldsymbol{\\tau}}\n", "\\newcommand{\\vxi}{\\boldsymbol{\\xi}}\n", "\n", "\n", "% Lower Roman (Vectors)\n", "\\newcommand{\\va}{\\mathbf{a}}\n", "\\newcommand{\\vb}{\\mathbf{b}}\n", "\\newcommand{\\vBt}{\\mathbf{\\tilde{B}}}\n", "\\newcommand{\\vc}{\\mathbf{c}}\n", "\\newcommand{\\vct}{\\mathbf{\\tilde{c}}}\n", "\\newcommand{\\vd}{\\mathbf{d}}\n", "\\newcommand{\\ve}{\\mathbf{e}}\n", "\\newcommand{\\vf}{\\mathbf{f}}\n", "\\newcommand{\\vg}{\\mathbf{g}}\n", "\\newcommand{\\vh}{\\mathbf{h}}\n", "%\\newcommand{\\myvh}{\\mathbf{h}}\n", "\\newcommand{\\vi}{\\mathbf{i}}\n", "\\newcommand{\\vj}{\\mathbf{j}}\n", "\\newcommand{\\vk}{\\mathbf{k}}\n", "\\newcommand{\\vl}{\\mathbf{l}}\n", "\\newcommand{\\vm}{\\mathbf{m}}\n", "\\newcommand{\\vn}{\\mathbf{n}}\n", "\\newcommand{\\vo}{\\mathbf{o}}\n", "\\newcommand{\\vp}{\\mathbf{p}}\n", "\\newcommand{\\vq}{\\mathbf{q}}\n", "\\newcommand{\\vr}{\\mathbf{r}}\n", "\\newcommand{\\vs}{\\mathbf{s}}\n", "\\newcommand{\\vt}{\\mathbf{t}}\n", "\\newcommand{\\vu}{\\mathbf{u}}\n", "\\newcommand{\\vv}{\\mathbf{v}}\n", "\\newcommand{\\vw}{\\mathbf{w}}\n", "\\newcommand{\\vws}{\\vw_s}\n", "\\newcommand{\\vwt}{\\mathbf{\\tilde{w}}}\n", "\\newcommand{\\vWt}{\\mathbf{\\tilde{W}}}\n", "\\newcommand{\\vwh}{\\hat{\\vw}}\n", "\\newcommand{\\vx}{\\mathbf{x}}\n", "%\\newcommand{\\vx}{\\mathbf{x}}\n", "\\newcommand{\\vxt}{\\mathbf{\\tilde{x}}}\n", "\\newcommand{\\vy}{\\mathbf{y}}\n", "\\newcommand{\\vyt}{\\mathbf{\\tilde{y}}}\n", "\\newcommand{\\vz}{\\mathbf{z}}\n", "%\\newcommand{\\vzt}{\\mathbf{\\tilde{z}}}\n", "\n", "\n", "% Upper Roman (Matrices)\n", "\\newcommand{\\vA}{\\mathbf{A}}\n", "\\newcommand{\\vB}{\\mathbf{B}}\n", "\\newcommand{\\vC}{\\mathbf{C}}\n", "\\newcommand{\\vD}{\\mathbf{D}}\n", "\\newcommand{\\vE}{\\mathbf{E}}\n", "\\newcommand{\\vF}{\\mathbf{F}}\n", "\\newcommand{\\vG}{\\mathbf{G}}\n", "\\newcommand{\\vH}{\\mathbf{H}}\n", "\\newcommand{\\vI}{\\mathbf{I}}\n", "\\newcommand{\\vJ}{\\mathbf{J}}\n", "\\newcommand{\\vK}{\\mathbf{K}}\n", "\\newcommand{\\vL}{\\mathbf{L}}\n", "\\newcommand{\\vM}{\\mathbf{M}}\n", "\\newcommand{\\vMt}{\\mathbf{\\tilde{M}}}\n", "\\newcommand{\\vN}{\\mathbf{N}}\n", "\\newcommand{\\vO}{\\mathbf{O}}\n", "\\newcommand{\\vP}{\\mathbf{P}}\n", "\\newcommand{\\vQ}{\\mathbf{Q}}\n", "\\newcommand{\\vR}{\\mathbf{R}}\n", "\\newcommand{\\vS}{\\mathbf{S}}\n", "\\newcommand{\\vT}{\\mathbf{T}}\n", "\\newcommand{\\vU}{\\mathbf{U}}\n", "\\newcommand{\\vV}{\\mathbf{V}}\n", "\\newcommand{\\vW}{\\mathbf{W}}\n", "\\newcommand{\\vX}{\\mathbf{X}}\n", "%\\newcommand{\\vXs}{\\vX_{\\vs}}\n", "\\newcommand{\\vXs}{\\vX_{s}}\n", "\\newcommand{\\vXt}{\\mathbf{\\tilde{X}}}\n", "\\newcommand{\\vY}{\\mathbf{Y}}\n", "\\newcommand{\\vZ}{\\mathbf{Z}}\n", "\\newcommand{\\vZt}{\\mathbf{\\tilde{Z}}}\n", "\\newcommand{\\vzt}{\\mathbf{\\tilde{z}}}\n", "\n", "\n", "%%%%\n", "\\newcommand{\\hidden}{\\vz}\n", "\\newcommand{\\hid}{\\hidden}\n", "\\newcommand{\\observed}{\\vy}\n", "\\newcommand{\\obs}{\\observed}\n", "\\newcommand{\\inputs}{\\vu}\n", "\\newcommand{\\input}{\\inputs}\n", "\n", "\\newcommand{\\hmmTrans}{\\vA}\n", "\\newcommand{\\hmmObs}{\\vB}\n", "\\newcommand{\\hmmInit}{\\vpi}\n", "\\newcommand{\\hmmhid}{\\hidden}\n", "\\newcommand{\\hmmobs}{\\obs}\n", "\n", "\\newcommand{\\ldsDyn}{\\vA}\n", "\\newcommand{\\ldsObs}{\\vC}\n", "\\newcommand{\\ldsDynIn}{\\vB}\n", "\\newcommand{\\ldsObsIn}{\\vD}\n", "\\newcommand{\\ldsDynNoise}{\\vQ}\n", "\\newcommand{\\ldsObsNoise}{\\vR}\n", "\n", "\\newcommand{\\ssmDynFn}{f}\n", "\\newcommand{\\ssmObsFn}{h}\n", "\n", "\n", "%%%\n", "\\newcommand{\\gauss}{\\mathcal{N}}\n", "\n", "\\newcommand{\\diag}{\\mathrm{diag}}\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(sec:lds-intro)=\n", "# Linear Gaussian SSMs\n", "\n", "\n", "Consider the state space model in \n", "{eq}`eq:SSM-ar`\n", "where we assume the observations are conditionally iid given the\n", "hidden states and inputs (i.e. there are no auto-regressive dependencies\n", "between the observables).\n", "We can rewrite this model as \n", "a stochastic nonlinear dynamical system (NLDS)\n", "by defining the distribution of the next hidden state \n", "as a deterministic function of the past state\n", "plus random process noise $\\vepsilon_t$ \n", "\\begin{align}\n", "\\hmmhid_t &= \\ssmDynFn(\\hmmhid_{t-1}, \\inputs_t, \\vepsilon_t) \n", "\\end{align}\n", "where $\\vepsilon_t$ is drawn from the distribution such\n", "that the induced distribution\n", "on $\\hmmhid_t$ matches $p(\\hmmhid_t|\\hmmhid_{t-1}, \\inputs_t)$.\n", "Similarly we can rewrite the observation distributions\n", "as a deterministic function of the hidden state\n", "plus observation noise $\\veta_t$:\n", "\\begin{align}\n", "\\hmmobs_t &= \\ssmObsFn(\\hmmhid_{t}, \\inputs_t, \\veta_t)\n", "\\end{align}\n", "\n", "\n", "If we assume additive Gaussian noise,\n", "the model becomes\n", "\\begin{align}\n", "\\hmmhid_t &= \\ssmDynFn(\\hmmhid_{t-1}, \\inputs_t) + \\vepsilon_t \\\\\n", "\\hmmobs_t &= \\ssmObsFn(\\hmmhid_{t}, \\inputs_t) + \\veta_t\n", "\\end{align}\n", "where $\\vepsilon_t \\sim \\gauss(\\vzero,\\vQ_t)$\n", "and $\\veta_t \\sim \\gauss(\\vzero,\\vR_t)$.\n", "We will call these Gaussian SSMs.\n", "\n", "If we additionally assume\n", "the transition function $\\ssmDynFn$\n", "and the observation function $\\ssmObsFn$ are both linear,\n", "then we can rewrite the model as follows:\n", "\\begin{align}\n", "p(\\hmmhid_t|\\hmmhid_{t-1},\\inputs_t) &= \\gauss(\\hmmhid_t|\\ldsDyn_t \\hmmhid_{t-1}\n", "+ \\ldsDynIn_t \\inputs_t, \\vQ_t)\n", "\\\\\n", "p(\\hmmobs_t|\\hmmhid_t,\\inputs_t) &= \\gauss(\\hmmobs_t|\\ldsObs_t \\hmmhid_{t}\n", "+ \\ldsObsIn_t \\inputs_t, \\vR_t)\n", "\\end{align}\n", "This is called a \n", "linear-Gaussian state space model\n", "(LG-SSM),\n", "or a\n", "linear dynamical system (LDS).\n", "We usually assume the parameters are independent of time, in which case\n", "the model is said to be time-invariant or homogeneous.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(sec:tracking-lds)=\n", "(sec:kalman-tracking)=\n", "## Example: tracking a 2d point\n", "\n", "\n", "\n", "% Sarkkar p43\n", "Consider an object moving in $\\real^2$.\n", "Let the state be\n", "the position and velocity of the object,\n", "$\\vz_t =\\begin{pmatrix} u_t & \\dot{u}_t & v_t & \\dot{v}_t \\end{pmatrix}$.\n", "(We use $u$ and $v$ for the two coordinates,\n", "to avoid confusion with the state and observation variables.)\n", "If we use Euler discretization,\n", "the dynamics become\n", "\\begin{align}\n", "\\underbrace{\\begin{pmatrix} u_t\\\\ \\dot{u}_t \\\\ v_t \\\\ \\dot{v}_t \\end{pmatrix}}_{\\vz_t}\n", " = \n", "\\underbrace{\n", "\\begin{pmatrix}\n", "1 & 0 & \\Delta & 0 \\\\\n", "0 & 1 & 0 & \\Delta\\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{pmatrix}\n", "}_{\\ldsDyn}\n", "\\\n", "\\underbrace{\\begin{pmatrix} u_{t-1} \\\\ \\dot{u}_{t-1} \\\\ v_{t-1} \\\\ \\dot{v}_{t-1} \\end{pmatrix}}_{\\vz_{t-1}}\n", "+ \\vepsilon_t\n", "\\end{align}\n", "where $\\vepsilon_t \\sim \\gauss(\\vzero,\\vQ)$ is\n", "the process noise.\n", "\n", "Let us assume\n", "that the process noise is \n", "a white noise process added to the velocity components\n", "of the state, but not to the location.\n", "(This is known as a random accelerations model.)\n", "We can approximate the resulting process in discrete time by assuming\n", "$\\vQ = \\diag(0, q, 0, q)$.\n", "(See {cite}`Sarkka13` p60 for a more accurate way\n", "to convert the continuous time process to discrete time.)\n", "\n", "\n", "Now suppose that at each discrete time point we\n", "observe the location,\n", "corrupted by Gaussian noise.\n", "Thus the observation model becomes\n", "\\begin{align}\n", "\\underbrace{\\begin{pmatrix} y_{1,t} \\\\ y_{2,t} \\end{pmatrix}}_{\\vy_t}\n", " &=\n", " \\underbrace{\n", " \\begin{pmatrix}\n", "1 & 0 & 0 & 0 \\\\\n", "0 & 0 & 1 & 0\n", " \\end{pmatrix}\n", " }_{\\ldsObs}\n", " \\\n", "\\underbrace{\\begin{pmatrix} u_t\\\\ \\dot{u}_t \\\\ v_t \\\\ \\dot{v}_t \\end{pmatrix}}_{\\vz_t} \n", " + \\veta_t\n", "\\end{align}\n", "where $\\veta_t \\sim \\gauss(\\vzero,\\vR)$ is the \\keywordDef{observation noise}.\n", "We see that the observation matrix $\\ldsObs$ simply ``extracts'' the\n", "relevant parts of the state vector.\n", "\n", "Suppose we sample a trajectory and corresponding set\n", "of noisy observations from this model,\n", "$(\\vz_{1:T}, \\vy_{1:T}) \\sim p(\\vz,\\vy|\\vtheta)$.\n", "(We use diagonal observation noise,\n", "$\\vR = \\diag(\\sigma_1^2, \\sigma_2^2)$.)\n", "The results are shown below. \n" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LDS(A=DeviceArray([[1., 0., 1., 0.],\n", " [0., 1., 0., 1.],\n", " [0., 0., 1., 0.],\n", " [0., 0., 0., 1.]], dtype=float32), C=DeviceArray([[1, 0, 0, 0],\n", " [0, 1, 0, 0]], dtype=int32), Q=DeviceArray([[0.001, 0. , 0. , 0. ],\n", " [0. , 0.001, 0. , 0. ],\n", " [0. , 0. , 0.001, 0. ],\n", " [0. , 0. , 0. , 0.001]], dtype=float32), R=DeviceArray([[1., 0.],\n", " [0., 1.]], dtype=float32), mu=DeviceArray([ 8., 10., 1., 0.], dtype=float32), Sigma=DeviceArray([[1., 0., 0., 0.],\n", " [0., 1., 0., 0.],\n", " [0., 0., 1., 0.],\n", " [0., 0., 0., 1.]], dtype=float32), state_offset=None, obs_offset=None, nstates=4, nobs=2)\n" ] } ], "source": [ "key = jax.random.PRNGKey(314)\n", "timesteps = 15\n", "delta = 1.0\n", "A = jnp.array([\n", " [1, 0, delta, 0],\n", " [0, 1, 0, delta],\n", " [0, 0, 1, 0],\n", " [0, 0, 0, 1]\n", "])\n", "\n", "C = jnp.array([\n", " [1, 0, 0, 0],\n", " [0, 1, 0, 0]\n", "])\n", "\n", "state_size, _ = A.shape\n", "observation_size, _ = C.shape\n", "\n", "Q = jnp.eye(state_size) * 0.001\n", "R = jnp.eye(observation_size) * 1.0\n", "# Prior parameter distribution\n", "mu0 = jnp.array([8, 10, 1, 0]).astype(float)\n", "Sigma0 = jnp.eye(state_size) * 1.0\n", "\n", "from jsl.lds.kalman_filter import LDS, smooth, filter\n", "\n", "lds = LDS(A, C, Q, R, mu0, Sigma0)\n", "print(lds)\n", "\n" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "from jsl.demos.plot_utils import plot_ellipse\n", "\n", "def plot_tracking_values(observed, filtered, cov_hist, signal_label, ax):\n", " timesteps, _ = observed.shape\n", " ax.plot(observed[:, 0], observed[:, 1], marker=\"o\", linewidth=0,\n", " markerfacecolor=\"none\", markeredgewidth=2, markersize=8, label=\"observed\", c=\"tab:green\")\n", " ax.plot(*filtered[:, :2].T, label=signal_label, c=\"tab:red\", marker=\"x\", linewidth=2)\n", " for t in range(0, timesteps, 1):\n", " covn = cov_hist[t][:2, :2]\n", " plot_ellipse(covn, filtered[t, :2], ax, n_std=2.0, plot_center=False)\n", " ax.axis(\"equal\")\n", " ax.legend()" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(7.24486608505249, 23.857812213897706, 8.0420747756958, 11.636079216003418)" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAd3klEQVR4nO3df3xU9Z3v8dcnP0jCj0ZNKFdFAXuLAjGiRKqrtroR662uXpVaeahbFVBsuy3edltsK9bHY/deVqnFChZ4UJfetloFtd12t11+9LJqL0gj5UeorL11gaK4EqojgTD59bl/zCSGkDCTmZOZc5L38/HIY2bOmTnzySTnPd/5nu/5jrk7IiISPQX5LkBERDKjABcRiSgFuIhIRCnARUQiSgEuIhJRRbl8ssrKSh87dmwun1JEJPJeffXVBncf2X15TgN87Nix1NXV5fIpRUQiz8z29LRcXSgiIhGlABcRiSgFuIhIROW0D1xEBq6Wlhb27dvH0aNH811KZJWWljJ69GiKi4vTur8CXEQCsW/fPkaMGMHYsWMxs3yXEznuzsGDB9m3bx/jxo1L6zHqQhGRQBw9epSKigqFd4bMjIqKij59glGAi0hgFN7Z6evrpy4UEcmLWDzG2j1rOdh0kIqyCqaNmUZ5SXm+y4oUtcBFJKfcnaXbllK7qpaHNj7E4q2LeWjjQ9SuqmXptqUE+R0Fu3fvpqqqKrDtBeXyyy8P5KRGtcBFJKeWbV/Gkq1LALj41IuZVDmJnQ072bh/Y+fyOefNyWeJJ9Ta2kpRUTiiUy1wEcmZWDzGih0rAFh0xSKWX7WcL13wJZZftZxFly8CYMWOFcTisYy2/+ijj1JVVUVVVRWLFiW219rayq233sqECROYPn06R44cAWDevHlMnDiR6upqvvKVrwBw4MABbrrpJi688EIuvPBCfvOb3wDwrW99i9tvv51LLrmE22+/nYsuuoidO3d2Pm9Hi/rw4cPcddddTJ06lfPPP5+f/exnADQ1NXHLLbcwYcIEbrjhBpqamjL6/Y7j7jn7mTJliovIwPT73/8+5X1W/fsqr1pZ5bP/dXaP62f96yyvWlnlq/99dZ+fv66uzquqqryxsdEPHTrkEydO9C1btjjgL7/8sru733nnnf7II494Q0ODjx8/3tvb293d/d1333V39xkzZvhLL73k7u579uzxc845x93dH3zwQb/gggv8yJEj7u7+6KOP+vz5893d/a233vLx48e7u/v999/vP/zhDzu3+dGPftQbGxv929/+tt95553u7r5t2zYvLCz03/72tz3+Hj29jkCd95CpaoGLSM4cbDoIwKTKST2ur6pM9Fc3NDX0edsvv/wyN9xwA8OGDWP48OHceOONvPTSS5xxxhlccsklANx22228/PLLlJeXU1paysyZM3n++ecZOnQoAOvWreMLX/gCkydP5rrrruP999+nsbERgOuuu46ysjIAbr75ZlavXg3As88+y/Tp0wFYs2YNCxYsYPLkyVx++eUcPXqUvXv38uKLL3LbbbcBUF1dTXV1dZ9/v56EoyNHRAaFirIKAHY27OxxfX1DPQCVZZWBPWf3oXlmRlFREZs3b2b9+vWsXr2axYsX8+tf/5r29nY2bdpEaWnpcdsZNmxY5/XTTz+diooKtm/fzjPPPMPSpUuBRI/Gc889x9lnnx1Y/SeiFriI5My0MdMoKSxh4/6NrN+z/ph16/esZ9P+TZQWlnLlmCv7vO3LLruMn/70pxw5coTDhw/zwgsvcNlll7F37142btwIwFNPPcWll15KY2MjsViMT33qU3znO99h27ZtAFx11VU8/vjjndvcunVrr8/3mc98hocffphYLNbZov7kJz/J448/3jmS5ne/+x0AH//4x3nqqacAqK+vZ/v27X3+/XqiABeRnCkvKWfWubMAmLthLrPXzOaxLY8xe81s5m6YC8DMc2dmNB78ggsu4I477mDq1Kl87GMfY9asWZx88smcffbZLFmyhAkTJvDuu+9y7733cujQIa699lqqq6u59NJLefTRRwH47ne/S11dHdXV1UycOLGzZd2T6dOn85Of/ISbb765c9kDDzxAS0sL1dXVTJo0iQceeACAe++9l8bGRiZMmMD8+fOZMmVKn3+/nljHO0Uu1NTUuL7QQWRgeu2115gwYULK+7k7y7YvY8WOFcTb4p3LSwpLmHXuLO6pvmdQn9HZ0+toZq+6e033+6oPXERyysyYc94cZpwzg3V71tHQ1EBlWSVXjrlSZ2L2kQJcRPKivKScm8bflO8yIk194CIiEZUywM3sSTN7x8zqe1j3ZTNzMwtuzI+IiKQlnRb4SuDq7gvN7AzgKmBvwDWJiEgaUga4u78I/LmHVd8BvgrkbhiLiIh0yuggppldD7zp7ttSDfcxs7uBuwHOPPPMTJ5ORAaQmr9bS0Njc6/rK4cPoe6b0/q83ffee4+nnnqKz33uc3163MqVK7nqqqs47bTTABg7dix1dXVUVoa/Z7jPBzHNbCjwdWB+Ovd39+XuXuPuNSNHjuzr04nIAHOi8E5nfW/ee+89nnjiieOWt7a2nvBxK1eu5K233sroOfMtkxb4R4BxQEfrezSwxcymuvvbQRYnIpKuefPm8cc//pHJkydTXFxMaWkpJ598Mrt27WLNmjVce+211NcnxmIsXLiQxsZGqqqqqKur49Zbb6WsrKzzlPvHH3+cn//857S0tLBq1SrOOeecfP5qvepzgLv7DuDDHbfNbDdQ4+59nz5MRAaksfP+uV8ev3vBNb0+ZsGCBdTX17N161Y2bNjANddcQ319PePGjWP37t09Pmb69OksXryYhQsXUlPzwYmOlZWVbNmyhSeeeIKFCxeyYsWKrH6f/pLOMMKngY3A2Wa2z8xm9n9ZIiLZmTp1KuPGjcvosTfeeCMAU6ZM6TX8wyBlC9zdZ6RYPzawakRkQDhRSzmd1vmJHp+urtO/FhUV0d7e3nn76NGjJ3xsSUkJAIWFhSn70PNJZ2KKyIAwYsQIDh061OO6UaNG8c4773Dw4EHi8Ti/+MUv0npc2GkuFBEZECoqKrjkkkuoqqqirKyMUaNGda4rLi5m/vz5TJ06ldNPP/2Yg5J33HEHc+bMOeYgZlRoOlkRCUS608n21zjwgULTyYpIaA3mcA6a+sBFRCJKAS4igclll+xA1NfXTwEuIoEoLS3l4MGDCvEMuTsHDx6ktLQ07ceoD1xEAjF69Gj27dvHgQMH8l1KZJWWljJ69Oi0768AF5FAFBcXZ3zmo2RGXSgiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiEoZ4Gb2pJm9Y2b1XZY9Yma7zGy7mb1gZif1a5UiInKcdFrgK4Gruy1bC1S5ezXwOnB/wHWJiEgKKQPc3V8E/txt2Rp3b03e3ASM7ofaRETkBILoA78L+GUA2xERkT7IKsDN7BtAK/DjE9znbjOrM7O6AwcOZPN0IiLSRcYBbmZ3ANcCt7q793Y/d1/u7jXuXjNy5MhMn05ERLopyuRBZnY18FXgE+5+JNiSREQkHekMI3wa2AicbWb7zGwmsBgYAaw1s61mtrSf6xQRkW5StsDdfUYPi7/fD7WIiEgfZNSFIiLREYvHWLtnLQebDlJRVsG0MdMoLynPd1kSAAW4yADl7izbvowVO1YQb4t3Ll+weQGzzp3FPdX3YGZ5rFCypQAXGaCWbV/Gkq1LALj41IuZVDmJnQ072bh/Y+fyOefNyWeJkiUFuESeugiOF4vHWLFjBQCLrlhE7Zm1nevW71nP3A1zWbFjBTPOmTHoX6soU4BLZKmLoHdr96wl3hbn4lMvPia8AWrH1HLRqRexaf8m1u1Zx03jb8pTlZItTScrkdXRRdARVLPOncXFp15MvC3Okq1LWLZ9Wb5LzJuDTQcBmFQ5qcf1VZVVADQ0NeSsJgmeWuASSeoiOLGKsgoAdjbs7HF9fUNidujKssqc1STBUwtcIimdLoJ4W5x1e9blqcL8mjZmGiWFJWzcv5H1e9Yfs279nvVs2r+J0sJSrhxzZZ4qlCCoBS6RlE4Xwab9mwZtF0F5STmzzp3Fkq1LmLthLhedehFVlVXUN9Szaf8mAGaeO3NQfjoZSBTgEknqIkjtnup7AFixYwWb9m/qDO6SwpLOg7wSbQpwiaRpY6axYPOCzi6C2jHH9oGriwDMjDnnzWHGOTNYt2cdDU0NVJZVcuWYK9XyHiAU4BJJ6iJIX3lJuYYKDlAKcIksdRHIYGcn+C6GwNXU1HhdXV3Onk8Gh1g8pi4CGdDM7FV3r+m+XC1wiTx1EUSbpkLInAI8IPonFOkbTYWQPQV4N30N4sH4T6g3KwmCZkvMnvrAk3oL4q4HxHoK4qXblvb6Twjw+cmfHzD/hJm+RiLdxeIxalfVEm+L9zoVQklhCes/vV6NA9QHnlImrYHBNh+HWkwSFM2WGAwFOJkH8WD6Jxxsb1bSvzQVQjA0mRWZT4w0WKbsjMVjLNi8gHhbnLPKz6Jm1LGf5DR5lPSVpkIIhgKczIN4oP8TujtLty2ldlUtv3jjFwC8EXuD2lW1LN22lK7HTwbKm5XkhmZLDIYCnMyDeKD/E3b9woSzys8CoHxIeY9fmBD1NyvJrY6pEADmbpjL7DWzeWzLY8xeM5u5G+YCmgohHeoDJ/OJkQbyfBzd+7xrRtVQu6qWWHOMuybdxZM7n+zs8657uy7yb1aSe/mYCqHm79bS0Njc6/rK4UOo++a0vG8zXQpwsgvigTofR0/HBTpeoyd3Pkn5kHJizTH++pd/zRuxN4DovllJfqSaLbE/gvFE20tnfa62mS4FeFKmQTxQp+zs6bhA19co1hwDEn3iUX+zCpP+as2Febu9TYWQz2CMCgV4UrZBPNDm4+jpuEDX16ij5f1XZ/0VX5v6tUi/WYVJf4VW1Lab7gmGP9m8l3hrO/HWNuIt7R9cb21P3k5eTy5Px4V/vw53p93pvGx3x3u43Z7DEyF7ogDvZqAFcaZOdFyg7u063oi9QWlh6aAN73z2e4ZFuiH7wu/2cTjexuF4K4ebE5dHmltp7FgWb+VwcytH4m00xls50tzG4ebWtLY97/kd2fwKPTpwKJ76TiGhAJceDeQDtEEIquXZ3u6819TCwcY4BxrTC47/8exWSGZnR4S6e5frx69Lx03f+7+0trXT0ua0trfT2ua0dFx2XdbWTmu709ae3nbve2ZbWvfLxGdqzqCkuICSogJKigoTl8W9XC8q5Lbvv5Jym5u/XouZYQYFZhQYx90uSN42jPHf/GW//X6pKMClV7k8QBvmPtpMvbb/fRoa44mfQ800JEO6obGZhkOJ5X8+3ExrmkHY4fktb/ZLva/uebdftnv95NMYVlLEsCGFycuixGVJ4bHXk+uGDingR7ue5OFnzki57X+YXh14vR/+UGng2+wvCnDp1YmOC9Q+spkFT/9Lr4/tazCGsY/W3Tna0k6sqYX3mpqJHWnhvaYWYk0taT33f3vspbTu96HSIipHlFA5rITNu/+c8v4LP30eAB3ThnXMH9bRIuy6rMOXfrI15XZXzbmYogKjuLCAokKjqKCA4sIPbhcXJC8LCygqMAoLjHH39/4/0OGxW85PeZ+ulm5byve2LwEW9Olxg1FoA1x9jMHor1EC+R4h0HEwqbW9nbZ2p7XdaU9etnX5ScdDP99JrKmF2JGWZFi3dN5ubmvPuMbxo4ZTMawkEc7Dh1A5vISRw0uoHJG4Xjm8hIrhQygpKux8zNh5/5xyu9OnjO5zLekE+IVjT+nzdoPW9fyD/lA5fEjK/SEM20xXygA3syeBa4F33L0quewU4BlgLLAbuNndA/381R8BEYU3haBrzGfQvvP+UQ4lD1I1Hm2lMZ74ORxv5VByWcf1dEyc/6tjAjoo//ib3b2uG1JUwEllxZSXFXPS0MRledkQntuyL+V219z3icBqHCy6nn+w5U/BB2N/7N/5zIx0WuArgcXA/+6ybB6w3t0XmNm85O2vBV9esIIOs7CeaODuNLW00Xg0vWBc/Os/0Nruxxy0am1rp6U9cdnalmjZtrYnD3Cl2Sqd+j/Xp75THxxpPnYYWIFBUUEBBQWJy8Lkx/rCAuv8iL/v3aaU2/3mNROSwdwR1EM6A7u0uLDHx6QT4Jnor9ZcVLbb9fyD5T3sO49teYwVO1bwhclf4J7zdN5BygB39xfNbGy3xdcDlyev/wDYQB4C/JIFv2Z4SRHDSxMHQkaUFHXeHt79emnwvUVBvyGkO1rgH361i0NHW2g82sqhjp94K4eOtnAo2dLtSwt14ZrX+1RnuiqHlzCiNHGAqvPvUZL4Ww0vLWL4kA/+dl9dvT3l9nY+9MkPQtqMgoLUXx6RTpfErMvOSuv3yYX+as1FZbsDfYK4oGWaaqPcfX/y+tvAqN7uaGZ3A3cDnHnmmRk+Xc/efC9166qvrl70ImVDCikrTv4krw8dUkjpkEKGFhdRNqSAsiHpvXTL/u2Px4x/7bzMYuzr9zb8MeV9SooKGFFaTEMaQ9PuvfwjFBcYRd0OVhUVFnQuL04e1EocxDLuWpn6m5Xqvpn+vCjpBPiwkvAcsslnv+dAlum8RGHWn19BmPUe4e5uZr0299x9ObAcEl+plu3zdfXSV6/o7Fdt7NrPevSDPtbGeAuH420cirfy4usHUm5z19uHgiyR//XLXYFuD+BvP3k2w0uKGFFaxIjS4i7XP7g9pCgx0WQ6LdCvXX1O4DUOdPk+VjJQDaTzD3LxfbmZBvh/mtmp7r7fzE4F3smqigydccrQPt0/nTD7ly9eRlNLK03N7RxpbqWppY2m5rbEZcf15jaOtLTx1Ct7U25v9mXjehz7OjTZrTB0SKJbYWhy+Ue+nnpY1uev+K9p/b5REpU+Wul/A2WCuFx8BWGmAf5PwGdJDNT8LPCzrKoIkYmnfSjt+6YT4N+4ZmI25YRW0MEYlT5ayUxfuhEGwgRxufoKwnSGET5N4oBlpZntAx4kEdzPmtlMYA9wc8YV9EItp2D01+uoYJR0ZNONEOV5iXL1fbnpjEKZ0cuq2l6WB6I/AiIKbwpRadmKpCMX3QhhlKsvbQ7PYf0cCDrMUoXt8DInFo/16SOSAlcGilx1I4RRroZDDqoAD1rXsO3to2LtqpLAjjiLREmuuhHCKFfDIRXgARmsHxVFepOrboQwytVwSAV4AAbzR0WR3gz2sypzMRxSAR6AqHxU7M8zwkS6G4hnVfZFLoZDKsADEPaPirk4I0yku4F0VmU2+nM4pAI8AGH/qKj+ecmXgXJWZVhZujPgBaGmpsbr6lJPghQ1sXiM2lW1xNviLLp80XEfFedumEtpYSnrPr0u562NY2rrpX++pLCE9Z9eP6BaQuouCpdYPBbZsyrDwMxedfea7svVAg9AmD8qRqV/PijqLgqnKJ9VGWYK8ICE9aNi2Pvng6buIhlMFOABCesEPGHvnw+ShnPKYFOQ7wIGmo6Pivecdw83jb8p70Exbcw0SgpLOodydZXPoVyxeIzVr69m2bZlrH59NbF4LOttptNdFG+Ls27PuqyfSyQM1AIf4MLWP9+ffdSDrbtIRAE+CISpf74/+6gHU3eRCGgY4aCS76Fc/T2kMczDOUWyoWGEkvehXP09pDFs3UUi/U0BLjmTiz7qMHUXifQ3BbjkTC76qMM6nFOkPyjAJWdyOTtdvruLRHJBAS45oz5qkWApwCWn1EctEhwNI5S8yPeQRpEo0TBCCRX1UYtkT3OhiIhElAJcRCSi1IUiIqGgb1HqOwW4iOSVvkUpcwpwEckrfYtS5hTgIpI3+hal7Oggpojkjb5FKTsKcBHJm3RmqAT0LUq9yCrAzew+M9tpZvVm9rSZlQZVmIgMfPoWpexkHOBmdjrwRaDG3auAQuCWoAoTkYEvrF+6HRXZHsQsAsrMrAUYCryVfUkiMlhohsrsZBzg7v6mmS0E9gJNwBp3X9P9fmZ2N3A3wJlnnpnp04nIAKUZKjOX8WyEZnYy8BzwGeA9YBWw2t1/1NtjNBuhiPRGM1T2rj9mI7wS+A93P5B8gueBvwB6DXARkd5ohsq+y2YUyl7gIjMbaonzXGuB14IpS0REUsk4wN39FWA1sAXYkdzW8oDqEhGRFLIaheLuDwIPBlSLiIj0gc7EFBGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSisgpwMzvJzFab2S4ze83MLg6qMBERObGiLB//GPArd59uZkOAoQHUJCIiacg4wM2sHPg4cAeAuzcDzcGUJSIiqWTThTIOOAD8o5n9zsxWmNmw7ncys7vNrM7M6g4cOJDF04mISFfZBHgRcAHwPXc/HzgMzOt+J3df7u417l4zcuTILJ5ORES6yibA9wH73P2V5O3VJAJdRERyIOMAd/e3gT+Z2dnJRbXA7wOpSkREUsp2FMrfAD9OjkB5A7gz+5JERCQdWQW4u28FaoIpRURE+kJnYoqIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkoszdc/dkZgeAPRk+vBJoCLCcoIW9Pgh/jaove2GvUfVlZoy7j+y+MKcBng0zq3P3mnzX0Zuw1wfhr1H1ZS/sNaq+YKkLRUQkohTgIiIRFaUAX57vAlIIe30Q/hpVX/bCXqPqC1Bk+sBFRORYUWqBi4hIFwpwEZGIikSAm9l9ZrbTzOrN7GkzK81zPU+a2TtmVt9l2SlmttbM/pC8PDmENT5iZrvMbLuZvWBmJ4Wpvi7rvmxmbmaV+agtWUOP9ZnZ3yRfw51m9nC+6kvW0tPfeLKZbTKzrWZWZ2ZT81TbGWb2f8zs98nX6kvJ5aHZT05QY2j2k1RCH+BmdjrwRaDG3auAQuCW/FbFSuDqbsvmAevd/aPA+uTtfFrJ8TWuBarcvRp4Hbg/10V1sZLj68PMzgCuAvbmuqBuVtKtPjO7ArgeOM/dJwEL81BXVys5/jV8GHjI3ScD85O386EV+LK7TwQuAj5vZhMJ137SW41h2k9OKPQBnlQElJlZETAUeCufxbj7i8Cfuy2+HvhB8voPgP+ey5q666lGd1/j7q3Jm5uA0Tkv7INaenoNAb4DfBXI69H1Xuq7F1jg7vHkfd7JeWFd9FKjAx9KXi8nT/uKu+939y3J64eA14DTCdF+0luNYdpPUgl9gLv7myRaOnuB/UDM3dfkt6oejXL3/cnrbwOj8llMGu4CfpnvIroys+uBN919W75r6cV44DIze8XM/s3MLsx3QT2YCzxiZn8isd/kvfVoZmOB84FXCOl+0q3GrkK3n3QV+gBP9pFdD4wDTgOGmdlt+a3qxDwxNjO04zPN7BskPj7+ON+1dDCzocDXSXzsD6si4BQSH7f/FnjWzCy/JR3nXuA+dz8DuA/4fj6LMbPhwHPAXHd/v+u6sOwnvdUYxv2ku9AHOHAl8B/ufsDdW4Dngb/Ic009+U8zOxUgeZnXj9e9MbM7gGuBWz1cJwF8hMSb9DYz203iY+sWM/svea3qWPuA5z1hM9BOYvKjMPksiX0EYBWQl4OYAGZWTCIYf+zuHTWFaj/ppcYw7yfHiEKA7wUuMrOhydZOLYm+qrD5JxI7D8nLn+Wxlh6Z2dUk+pevc/cj+a6nK3ff4e4fdvex7j6WRFhe4O5v57m0rn4KXAFgZuOBIYRv5rq3gE8kr/8l8Id8FJHcV78PvObuj3ZZFZr9pLcaw7yfHMfdQ/8DPATsAuqBHwIlea7naRL98S0kgmYmUEHiqPofgHXAKSGs8f8BfwK2Jn+Whqm+but3A5Vhqo9EYP8o+X+4BfjLEP6NLwVeBbaR6M+dkqfaLiXRPbK9y//bp8K0n5ygxtDsJ6l+dCq9iEhERaELRUREeqAAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hE1P8HfYldV3+GHYsAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\n", "z_hist, x_hist = lds.sample(key, timesteps)\n", "\n", "fig_truth, axs = plt.subplots()\n", "axs.plot(x_hist[:, 0], x_hist[:, 1],\n", " marker=\"o\", linewidth=0, markerfacecolor=\"none\",\n", " markeredgewidth=2, markersize=8,\n", " label=\"observed\", c=\"tab:green\")\n", "\n", "axs.plot(z_hist[:, 0], z_hist[:, 1],\n", " linewidth=2, label=\"truth\",\n", " marker=\"s\", markersize=8)\n", "axs.legend()\n", "axs.axis(\"equal\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main task is to infer the hidden states given the noisy\n", "observations, i.e., $p(\\vz|\\vy,\\vtheta)$. We discuss the topic of inference in {ref}`sec:inference`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.2" } }, "nbformat": 4, "nbformat_minor": 4 }