lds.html 50 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985
  1. <!DOCTYPE html>
  2. <html>
  3. <head>
  4. <meta charset="utf-8" />
  5. <meta name="viewport" content="width=device-width, initial-scale=1.0" />
  6. <title>Linear Gaussian SSMs &#8212; State Space Models: A Modern Approach</title>
  7. <link href="../../_static/css/theme.css" rel="stylesheet">
  8. <link href="../../_static/css/index.ff1ffe594081f20da1ef19478df9384b.css" rel="stylesheet">
  9. <link rel="stylesheet"
  10. href="../../_static/vendor/fontawesome/5.13.0/css/all.min.css">
  11. <link rel="preload" as="font" type="font/woff2" crossorigin
  12. href="../../_static/vendor/fontawesome/5.13.0/webfonts/fa-solid-900.woff2">
  13. <link rel="preload" as="font" type="font/woff2" crossorigin
  14. href="../../_static/vendor/fontawesome/5.13.0/webfonts/fa-brands-400.woff2">
  15. <link rel="stylesheet" type="text/css" href="../../_static/pygments.css" />
  16. <link rel="stylesheet" type="text/css" href="../../_static/sphinx-book-theme.css?digest=c3fdc42140077d1ad13ad2f1588a4309" />
  17. <link rel="stylesheet" type="text/css" href="../../_static/togglebutton.css" />
  18. <link rel="stylesheet" type="text/css" href="../../_static/copybutton.css" />
  19. <link rel="stylesheet" type="text/css" href="../../_static/mystnb.css" />
  20. <link rel="stylesheet" type="text/css" href="../../_static/sphinx-thebe.css" />
  21. <link rel="stylesheet" type="text/css" href="../../_static/panels-main.c949a650a448cc0ae9fd3441c0e17fb0.css" />
  22. <link rel="stylesheet" type="text/css" href="../../_static/panels-variables.06eb56fa6e07937060861dad626602ad.css" />
  23. <link rel="preload" as="script" href="../../_static/js/index.be7d3bbb2ef33a8344ce.js">
  24. <script data-url_root="../../" id="documentation_options" src="../../_static/documentation_options.js"></script>
  25. <script src="../../_static/jquery.js"></script>
  26. <script src="../../_static/underscore.js"></script>
  27. <script src="../../_static/doctools.js"></script>
  28. <script src="../../_static/clipboard.min.js"></script>
  29. <script src="../../_static/copybutton.js"></script>
  30. <script>let toggleHintShow = 'Click to show';</script>
  31. <script>let toggleHintHide = 'Click to hide';</script>
  32. <script>let toggleOpenOnPrint = 'true';</script>
  33. <script src="../../_static/togglebutton.js"></script>
  34. <script>var togglebuttonSelector = '.toggle, .admonition.dropdown, .tag_hide_input div.cell_input, .tag_hide-input div.cell_input, .tag_hide_output div.cell_output, .tag_hide-output div.cell_output, .tag_hide_cell.cell, .tag_hide-cell.cell';</script>
  35. <script src="../../_static/sphinx-book-theme.d59cb220de22ca1c485ebbdc042f0030.js"></script>
  36. <script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"
  37. const thebe_selector = ".thebe,.cell"
  38. const thebe_selector_input = "pre"
  39. const thebe_selector_output = ".output, .cell_output"
  40. </script>
  41. <script async="async" src="../../_static/sphinx-thebe.js"></script>
  42. <script>window.MathJax = {"TeX": {"Macros": {"N": "\\mathbb{N}", "floor": ["\\lfloor#1\\rfloor", 1], "bmat": ["\\left[\\begin{array}"], "emat": ["\\end{array}\\right]"]}}, "options": {"processHtmlClass": "tex2jax_process|mathjax_process|math|output_area"}}</script>
  43. <script defer="defer" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
  44. <link rel="index" title="Index" href="../../genindex.html" />
  45. <link rel="search" title="Search" href="../../search.html" />
  46. <link rel="next" title="Nonlinear Gaussian SSMs" href="nlds.html" />
  47. <link rel="prev" title="Hidden Markov Models" href="hmm.html" />
  48. <meta name="viewport" content="width=device-width, initial-scale=1" />
  49. <meta name="docsearch:language" content="None">
  50. <!-- Google Analytics -->
  51. </head>
  52. <body data-spy="scroll" data-target="#bd-toc-nav" data-offset="80">
  53. <div class="container-fluid" id="banner"></div>
  54. <div class="container-xl">
  55. <div class="row">
  56. <div class="col-12 col-md-3 bd-sidebar site-navigation show" id="site-navigation">
  57. <div class="navbar-brand-box">
  58. <a class="navbar-brand text-wrap" href="../../index.html">
  59. <h1 class="site-logo" id="site-title">State Space Models: A Modern Approach</h1>
  60. </a>
  61. </div><form class="bd-search d-flex align-items-center" action="../../search.html" method="get">
  62. <i class="icon fas fa-search"></i>
  63. <input type="search" class="form-control" name="q" id="search-input" placeholder="Search this book..." aria-label="Search this book..." autocomplete="off" >
  64. </form><nav class="bd-links" id="bd-docs-nav" aria-label="Main">
  65. <div class="bd-toc-item active">
  66. <ul class="nav bd-sidenav">
  67. <li class="toctree-l1">
  68. <a class="reference internal" href="../../root.html">
  69. State Space Models: A Modern Approach
  70. </a>
  71. </li>
  72. </ul>
  73. <ul class="current nav bd-sidenav">
  74. <li class="toctree-l1">
  75. <a class="reference internal" href="../scratch.html">
  76. Scratchpad
  77. </a>
  78. </li>
  79. <li class="toctree-l1 current active has-children">
  80. <a class="reference internal" href="ssm_index.html">
  81. State Space Models
  82. </a>
  83. <input checked="" class="toctree-checkbox" id="toctree-checkbox-1" name="toctree-checkbox-1" type="checkbox"/>
  84. <label for="toctree-checkbox-1">
  85. <i class="fas fa-chevron-down">
  86. </i>
  87. </label>
  88. <ul class="current">
  89. <li class="toctree-l2">
  90. <a class="reference internal" href="ssm_intro.html">
  91. What are State Space Models?
  92. </a>
  93. </li>
  94. <li class="toctree-l2">
  95. <a class="reference internal" href="hmm.html">
  96. Hidden Markov Models
  97. </a>
  98. </li>
  99. <li class="toctree-l2 current active">
  100. <a class="current reference internal" href="#">
  101. Linear Gaussian SSMs
  102. </a>
  103. </li>
  104. <li class="toctree-l2">
  105. <a class="reference internal" href="nlds.html">
  106. Nonlinear Gaussian SSMs
  107. </a>
  108. </li>
  109. <li class="toctree-l2">
  110. <a class="reference internal" href="inference.html">
  111. Inferential goals
  112. </a>
  113. </li>
  114. </ul>
  115. </li>
  116. <li class="toctree-l1 has-children">
  117. <a class="reference internal" href="../hmm/hmm_index.html">
  118. Hidden Markov Models
  119. </a>
  120. <input class="toctree-checkbox" id="toctree-checkbox-2" name="toctree-checkbox-2" type="checkbox"/>
  121. <label for="toctree-checkbox-2">
  122. <i class="fas fa-chevron-down">
  123. </i>
  124. </label>
  125. <ul>
  126. <li class="toctree-l2">
  127. <a class="reference internal" href="../hmm/hmm_filter.html">
  128. HMM filtering (forwards algorithm)
  129. </a>
  130. </li>
  131. <li class="toctree-l2">
  132. <a class="reference internal" href="../hmm/hmm_smoother.html">
  133. HMM smoothing (forwards-backwards algorithm)
  134. </a>
  135. </li>
  136. <li class="toctree-l2">
  137. <a class="reference internal" href="../hmm/hmm_viterbi.html">
  138. Viterbi algorithm
  139. </a>
  140. </li>
  141. <li class="toctree-l2">
  142. <a class="reference internal" href="../hmm/hmm_parallel.html">
  143. Parallel HMM smoothing
  144. </a>
  145. </li>
  146. <li class="toctree-l2">
  147. <a class="reference internal" href="../hmm/hmm_sampling.html">
  148. Forwards-filtering backwards-sampling algorithm
  149. </a>
  150. </li>
  151. </ul>
  152. </li>
  153. <li class="toctree-l1 has-children">
  154. <a class="reference internal" href="../lgssm/lgssm_index.html">
  155. Linear-Gaussian SSMs
  156. </a>
  157. <input class="toctree-checkbox" id="toctree-checkbox-3" name="toctree-checkbox-3" type="checkbox"/>
  158. <label for="toctree-checkbox-3">
  159. <i class="fas fa-chevron-down">
  160. </i>
  161. </label>
  162. <ul>
  163. <li class="toctree-l2">
  164. <a class="reference internal" href="../lgssm/kalman_filter.html">
  165. Kalman filtering
  166. </a>
  167. </li>
  168. <li class="toctree-l2">
  169. <a class="reference internal" href="../lgssm/kalman_smoother.html">
  170. Kalman (RTS) smoother
  171. </a>
  172. </li>
  173. <li class="toctree-l2">
  174. <a class="reference internal" href="../lgssm/kalman_parallel.html">
  175. Parallel Kalman Smoother
  176. </a>
  177. </li>
  178. <li class="toctree-l2">
  179. <a class="reference internal" href="../lgssm/kalman_sampling.html">
  180. Forwards-filtering backwards sampling
  181. </a>
  182. </li>
  183. </ul>
  184. </li>
  185. <li class="toctree-l1 has-children">
  186. <a class="reference internal" href="../extended/extended_index.html">
  187. Extended (linearized) methods
  188. </a>
  189. <input class="toctree-checkbox" id="toctree-checkbox-4" name="toctree-checkbox-4" type="checkbox"/>
  190. <label for="toctree-checkbox-4">
  191. <i class="fas fa-chevron-down">
  192. </i>
  193. </label>
  194. <ul>
  195. <li class="toctree-l2">
  196. <a class="reference internal" href="../extended/extended_filter.html">
  197. Extended Kalman filtering
  198. </a>
  199. </li>
  200. <li class="toctree-l2">
  201. <a class="reference internal" href="../extended/extended_smoother.html">
  202. Extended Kalman smoother
  203. </a>
  204. </li>
  205. <li class="toctree-l2">
  206. <a class="reference internal" href="../extended/extended_parallel.html">
  207. Parallel extended Kalman smoothing
  208. </a>
  209. </li>
  210. </ul>
  211. </li>
  212. <li class="toctree-l1 has-children">
  213. <a class="reference internal" href="../unscented/unscented_index.html">
  214. Unscented methods
  215. </a>
  216. <input class="toctree-checkbox" id="toctree-checkbox-5" name="toctree-checkbox-5" type="checkbox"/>
  217. <label for="toctree-checkbox-5">
  218. <i class="fas fa-chevron-down">
  219. </i>
  220. </label>
  221. <ul>
  222. <li class="toctree-l2">
  223. <a class="reference internal" href="../unscented/unscented_filter.html">
  224. Unscented filtering
  225. </a>
  226. </li>
  227. <li class="toctree-l2">
  228. <a class="reference internal" href="../unscented/unscented_smoother.html">
  229. Unscented smoothing
  230. </a>
  231. </li>
  232. </ul>
  233. </li>
  234. <li class="toctree-l1">
  235. <a class="reference internal" href="../quadrature/quadrature_index.html">
  236. Quadrature and cubature methods
  237. </a>
  238. </li>
  239. <li class="toctree-l1">
  240. <a class="reference internal" href="../postlin/postlin_index.html">
  241. Posterior linearization
  242. </a>
  243. </li>
  244. <li class="toctree-l1">
  245. <a class="reference internal" href="../adf/adf_index.html">
  246. Assumed Density Filtering
  247. </a>
  248. </li>
  249. <li class="toctree-l1">
  250. <a class="reference internal" href="../vi/vi_index.html">
  251. Variational inference
  252. </a>
  253. </li>
  254. <li class="toctree-l1">
  255. <a class="reference internal" href="../pf/pf_index.html">
  256. Particle filtering
  257. </a>
  258. </li>
  259. <li class="toctree-l1">
  260. <a class="reference internal" href="../smc/smc_index.html">
  261. Sequential Monte Carlo
  262. </a>
  263. </li>
  264. <li class="toctree-l1">
  265. <a class="reference internal" href="../learning/learning_index.html">
  266. Offline parameter estimation (learning)
  267. </a>
  268. </li>
  269. <li class="toctree-l1">
  270. <a class="reference internal" href="../tracking/tracking_index.html">
  271. Multi-target tracking
  272. </a>
  273. </li>
  274. <li class="toctree-l1">
  275. <a class="reference internal" href="../ensemble/ensemble_index.html">
  276. Data assimilation using Ensemble Kalman filter
  277. </a>
  278. </li>
  279. <li class="toctree-l1">
  280. <a class="reference internal" href="../bnp/bnp_index.html">
  281. Bayesian non-parametric SSMs
  282. </a>
  283. </li>
  284. <li class="toctree-l1">
  285. <a class="reference internal" href="../changepoint/changepoint_index.html">
  286. Changepoint detection
  287. </a>
  288. </li>
  289. <li class="toctree-l1">
  290. <a class="reference internal" href="../timeseries/timeseries_index.html">
  291. Timeseries forecasting
  292. </a>
  293. </li>
  294. <li class="toctree-l1">
  295. <a class="reference internal" href="../gp/gp_index.html">
  296. Markovian Gaussian processes
  297. </a>
  298. </li>
  299. <li class="toctree-l1">
  300. <a class="reference internal" href="../ode/ode_index.html">
  301. Differential equations and SSMs
  302. </a>
  303. </li>
  304. <li class="toctree-l1">
  305. <a class="reference internal" href="../control/control_index.html">
  306. Optimal control
  307. </a>
  308. </li>
  309. <li class="toctree-l1">
  310. <a class="reference internal" href="../../bib.html">
  311. Bibliography
  312. </a>
  313. </li>
  314. </ul>
  315. </div>
  316. </nav> <!-- To handle the deprecated key -->
  317. <div class="navbar_extra_footer">
  318. Powered by <a href="https://jupyterbook.org">Jupyter Book</a>
  319. </div>
  320. </div>
  321. <main class="col py-md-3 pl-md-4 bd-content overflow-auto" role="main">
  322. <div class="topbar container-xl fixed-top">
  323. <div class="topbar-contents row">
  324. <div class="col-12 col-md-3 bd-topbar-whitespace site-navigation show"></div>
  325. <div class="col pl-md-4 topbar-main">
  326. <button id="navbar-toggler" class="navbar-toggler ml-0" type="button" data-toggle="collapse"
  327. data-toggle="tooltip" data-placement="bottom" data-target=".site-navigation" aria-controls="navbar-menu"
  328. aria-expanded="true" aria-label="Toggle navigation" aria-controls="site-navigation"
  329. title="Toggle navigation" data-toggle="tooltip" data-placement="left">
  330. <i class="fas fa-bars"></i>
  331. <i class="fas fa-arrow-left"></i>
  332. <i class="fas fa-arrow-up"></i>
  333. </button>
  334. <div class="dropdown-buttons-trigger">
  335. <button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn" aria-label="Download this page"><i
  336. class="fas fa-download"></i></button>
  337. <div class="dropdown-buttons">
  338. <!-- ipynb file if we had a myst markdown file -->
  339. <!-- Download raw file -->
  340. <a class="dropdown-buttons" href="../../_sources/chapters/ssm/lds.ipynb"><button type="button"
  341. class="btn btn-secondary topbarbtn" title="Download source file" data-toggle="tooltip"
  342. data-placement="left">.ipynb</button></a>
  343. <!-- Download PDF via print -->
  344. <button type="button" id="download-print" class="btn btn-secondary topbarbtn" title="Print to PDF"
  345. onclick="printPdf(this)" data-toggle="tooltip" data-placement="left">.pdf</button>
  346. </div>
  347. </div>
  348. <!-- Source interaction buttons -->
  349. <div class="dropdown-buttons-trigger">
  350. <button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn"
  351. aria-label="Connect with source repository"><i class="fab fa-github"></i></button>
  352. <div class="dropdown-buttons sourcebuttons">
  353. <a class="repository-button"
  354. href="https://github.com/probml/ssm-book"><button type="button" class="btn btn-secondary topbarbtn"
  355. data-toggle="tooltip" data-placement="left" title="Source repository"><i
  356. class="fab fa-github"></i>repository</button></a>
  357. <a class="issues-button"
  358. href="https://github.com/probml/ssm-book/issues/new?title=Issue%20on%20page%20%2Fchapters/ssm/lds.html&body=Your%20issue%20content%20here."><button
  359. type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip" data-placement="left"
  360. title="Open an issue"><i class="fas fa-lightbulb"></i>open issue</button></a>
  361. </div>
  362. </div>
  363. <!-- Full screen (wrap in <a> to have style consistency -->
  364. <a class="full-screen-button"><button type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip"
  365. data-placement="bottom" onclick="toggleFullScreen()" aria-label="Fullscreen mode"
  366. title="Fullscreen mode"><i
  367. class="fas fa-expand"></i></button></a>
  368. <!-- Launch buttons -->
  369. <div class="dropdown-buttons-trigger">
  370. <button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn"
  371. aria-label="Launch interactive content"><i class="fas fa-rocket"></i></button>
  372. <div class="dropdown-buttons">
  373. <a class="binder-button" href="https://mybinder.org/v2/gh/probml/ssm-book/main?urlpath=tree/chapters/ssm/lds.ipynb"><button type="button"
  374. class="btn btn-secondary topbarbtn" title="Launch Binder" data-toggle="tooltip"
  375. data-placement="left"><img class="binder-button-logo"
  376. src="../../_static/images/logo_binder.svg"
  377. alt="Interact on binder">Binder</button></a>
  378. <a class="colab-button" href="https://colab.research.google.com/github/probml/ssm-book/blob/main/chapters/ssm/lds.ipynb"><button type="button" class="btn btn-secondary topbarbtn"
  379. title="Launch Colab" data-toggle="tooltip" data-placement="left"><img class="colab-button-logo"
  380. src="../../_static/images/logo_colab.png"
  381. alt="Interact on Colab">Colab</button></a>
  382. </div>
  383. </div>
  384. </div>
  385. <!-- Table of contents -->
  386. <div class="d-none d-md-block col-md-2 bd-toc show noprint">
  387. <div class="tocsection onthispage pt-5 pb-3">
  388. <i class="fas fa-list"></i> Contents
  389. </div>
  390. <nav id="bd-toc-nav" aria-label="Page">
  391. <ul class="visible nav section-nav flex-column">
  392. <li class="toc-h2 nav-item toc-entry">
  393. <a class="reference internal nav-link" href="#example-tracking-a-2d-point">
  394. Example: tracking a 2d point
  395. </a>
  396. </li>
  397. </ul>
  398. </nav>
  399. </div>
  400. </div>
  401. </div>
  402. <div id="main-content" class="row">
  403. <div class="col-12 col-md-9 pl-md-3 pr-md-0">
  404. <!-- Table of contents that is only displayed when printing the page -->
  405. <div id="jb-print-docs-body" class="onlyprint">
  406. <h1>Linear Gaussian SSMs</h1>
  407. <!-- Table of contents -->
  408. <div id="print-main-content">
  409. <div id="jb-print-toc">
  410. <div>
  411. <h2> Contents </h2>
  412. </div>
  413. <nav aria-label="Page">
  414. <ul class="visible nav section-nav flex-column">
  415. <li class="toc-h2 nav-item toc-entry">
  416. <a class="reference internal nav-link" href="#example-tracking-a-2d-point">
  417. Example: tracking a 2d point
  418. </a>
  419. </li>
  420. </ul>
  421. </nav>
  422. </div>
  423. </div>
  424. </div>
  425. <div>
  426. <div class="cell docutils container">
  427. <div class="cell_input docutils container">
  428. <div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># meta-data does not work yet in VScode</span>
  429. <span class="c1"># https://github.com/microsoft/vscode-jupyter/issues/1121</span>
  430. <span class="p">{</span>
  431. <span class="s2">&quot;tags&quot;</span><span class="p">:</span> <span class="p">[</span>
  432. <span class="s2">&quot;hide-cell&quot;</span>
  433. <span class="p">]</span>
  434. <span class="p">}</span>
  435. <span class="c1">### Install necessary libraries</span>
  436. <span class="k">try</span><span class="p">:</span>
  437. <span class="kn">import</span> <span class="nn">jax</span>
  438. <span class="k">except</span><span class="p">:</span>
  439. <span class="c1"># For cuda version, see https://github.com/google/jax#installation</span>
  440. <span class="o">%</span><span class="k">pip</span> install --upgrade &quot;jax[cpu]&quot;
  441. <span class="kn">import</span> <span class="nn">jax</span>
  442. <span class="k">try</span><span class="p">:</span>
  443. <span class="kn">import</span> <span class="nn">distrax</span>
  444. <span class="k">except</span><span class="p">:</span>
  445. <span class="o">%</span><span class="k">pip</span> install --upgrade distrax
  446. <span class="kn">import</span> <span class="nn">distrax</span>
  447. <span class="k">try</span><span class="p">:</span>
  448. <span class="kn">import</span> <span class="nn">jsl</span>
  449. <span class="k">except</span><span class="p">:</span>
  450. <span class="o">%</span><span class="k">pip</span> install git+https://github.com/probml/jsl
  451. <span class="kn">import</span> <span class="nn">jsl</span>
  452. <span class="k">try</span><span class="p">:</span>
  453. <span class="kn">import</span> <span class="nn">rich</span>
  454. <span class="k">except</span><span class="p">:</span>
  455. <span class="o">%</span><span class="k">pip</span> install rich
  456. <span class="kn">import</span> <span class="nn">rich</span>
  457. </pre></div>
  458. </div>
  459. </div>
  460. </div>
  461. <div class="cell docutils container">
  462. <div class="cell_input docutils container">
  463. <div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="p">{</span>
  464. <span class="s2">&quot;tags&quot;</span><span class="p">:</span> <span class="p">[</span>
  465. <span class="s2">&quot;hide-cell&quot;</span>
  466. <span class="p">]</span>
  467. <span class="p">}</span>
  468. <span class="c1">### Import standard libraries</span>
  469. <span class="kn">import</span> <span class="nn">abc</span>
  470. <span class="kn">from</span> <span class="nn">dataclasses</span> <span class="kn">import</span> <span class="n">dataclass</span>
  471. <span class="kn">import</span> <span class="nn">functools</span>
  472. <span class="kn">import</span> <span class="nn">itertools</span>
  473. <span class="kn">from</span> <span class="nn">typing</span> <span class="kn">import</span> <span class="n">Any</span><span class="p">,</span> <span class="n">Callable</span><span class="p">,</span> <span class="n">NamedTuple</span><span class="p">,</span> <span class="n">Optional</span><span class="p">,</span> <span class="n">Union</span><span class="p">,</span> <span class="n">Tuple</span>
  474. <span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
  475. <span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
  476. <span class="kn">import</span> <span class="nn">jax</span>
  477. <span class="kn">import</span> <span class="nn">jax.numpy</span> <span class="k">as</span> <span class="nn">jnp</span>
  478. <span class="kn">from</span> <span class="nn">jax</span> <span class="kn">import</span> <span class="n">lax</span><span class="p">,</span> <span class="n">vmap</span><span class="p">,</span> <span class="n">jit</span><span class="p">,</span> <span class="n">grad</span>
  479. <span class="kn">from</span> <span class="nn">jax.scipy.special</span> <span class="kn">import</span> <span class="n">logit</span>
  480. <span class="kn">from</span> <span class="nn">jax.nn</span> <span class="kn">import</span> <span class="n">softmax</span>
  481. <span class="kn">from</span> <span class="nn">functools</span> <span class="kn">import</span> <span class="n">partial</span>
  482. <span class="kn">from</span> <span class="nn">jax.random</span> <span class="kn">import</span> <span class="n">PRNGKey</span><span class="p">,</span> <span class="n">split</span>
  483. <span class="kn">import</span> <span class="nn">inspect</span>
  484. <span class="kn">import</span> <span class="nn">inspect</span> <span class="k">as</span> <span class="nn">py_inspect</span>
  485. <span class="kn">import</span> <span class="nn">rich</span>
  486. <span class="kn">from</span> <span class="nn">rich</span> <span class="kn">import</span> <span class="n">inspect</span> <span class="k">as</span> <span class="n">r_inspect</span>
  487. <span class="kn">from</span> <span class="nn">rich</span> <span class="kn">import</span> <span class="nb">print</span> <span class="k">as</span> <span class="n">r_print</span>
  488. <span class="k">def</span> <span class="nf">print_source</span><span class="p">(</span><span class="n">fname</span><span class="p">):</span>
  489. <span class="n">r_print</span><span class="p">(</span><span class="n">py_inspect</span><span class="o">.</span><span class="n">getsource</span><span class="p">(</span><span class="n">fname</span><span class="p">))</span>
  490. </pre></div>
  491. </div>
  492. </div>
  493. </div>
  494. <div class="math notranslate nohighlight">
  495. \[ \begin{align}\begin{aligned}\newcommand\floor[1]{\lfloor#1\rfloor}\\\newcommand{\real}{\mathbb{R}}\\% Numbers
  496. \newcommand{\vzero}{\boldsymbol{0}}
  497. \newcommand{\vone}{\boldsymbol{1}}\\% Greek https://www.latex-tutorial.com/symbols/greek-alphabet/
  498. \newcommand{\valpha}{\boldsymbol{\alpha}}
  499. \newcommand{\vbeta}{\boldsymbol{\beta}}
  500. \newcommand{\vchi}{\boldsymbol{\chi}}
  501. \newcommand{\vdelta}{\boldsymbol{\delta}}
  502. \newcommand{\vDelta}{\boldsymbol{\Delta}}
  503. \newcommand{\vepsilon}{\boldsymbol{\epsilon}}
  504. \newcommand{\vzeta}{\boldsymbol{\zeta}}
  505. \newcommand{\vXi}{\boldsymbol{\Xi}}
  506. \newcommand{\vell}{\boldsymbol{\ell}}
  507. \newcommand{\veta}{\boldsymbol{\eta}}
  508. %\newcommand{\vEta}{\boldsymbol{\Eta}}
  509. \newcommand{\vgamma}{\boldsymbol{\gamma}}
  510. \newcommand{\vGamma}{\boldsymbol{\Gamma}}
  511. \newcommand{\vmu}{\boldsymbol{\mu}}
  512. \newcommand{\vmut}{\boldsymbol{\tilde{\mu}}}
  513. \newcommand{\vnu}{\boldsymbol{\nu}}
  514. \newcommand{\vkappa}{\boldsymbol{\kappa}}
  515. \newcommand{\vlambda}{\boldsymbol{\lambda}}
  516. \newcommand{\vLambda}{\boldsymbol{\Lambda}}
  517. \newcommand{\vLambdaBar}{\overline{\vLambda}}
  518. %\newcommand{\vnu}{\boldsymbol{\nu}}
  519. \newcommand{\vomega}{\boldsymbol{\omega}}
  520. \newcommand{\vOmega}{\boldsymbol{\Omega}}
  521. \newcommand{\vphi}{\boldsymbol{\phi}}
  522. \newcommand{\vvarphi}{\boldsymbol{\varphi}}
  523. \newcommand{\vPhi}{\boldsymbol{\Phi}}
  524. \newcommand{\vpi}{\boldsymbol{\pi}}
  525. \newcommand{\vPi}{\boldsymbol{\Pi}}
  526. \newcommand{\vpsi}{\boldsymbol{\psi}}
  527. \newcommand{\vPsi}{\boldsymbol{\Psi}}
  528. \newcommand{\vrho}{\boldsymbol{\rho}}
  529. \newcommand{\vtheta}{\boldsymbol{\theta}}
  530. \newcommand{\vthetat}{\boldsymbol{\tilde{\theta}}}
  531. \newcommand{\vTheta}{\boldsymbol{\Theta}}
  532. \newcommand{\vsigma}{\boldsymbol{\sigma}}
  533. \newcommand{\vSigma}{\boldsymbol{\Sigma}}
  534. \newcommand{\vSigmat}{\boldsymbol{\tilde{\Sigma}}}
  535. \newcommand{\vsigmoid}{\vsigma}
  536. \newcommand{\vtau}{\boldsymbol{\tau}}
  537. \newcommand{\vxi}{\boldsymbol{\xi}}\\
  538. % Lower Roman (Vectors)
  539. \newcommand{\va}{\mathbf{a}}
  540. \newcommand{\vb}{\mathbf{b}}
  541. \newcommand{\vBt}{\mathbf{\tilde{B}}}
  542. \newcommand{\vc}{\mathbf{c}}
  543. \newcommand{\vct}{\mathbf{\tilde{c}}}
  544. \newcommand{\vd}{\mathbf{d}}
  545. \newcommand{\ve}{\mathbf{e}}
  546. \newcommand{\vf}{\mathbf{f}}
  547. \newcommand{\vg}{\mathbf{g}}
  548. \newcommand{\vh}{\mathbf{h}}
  549. %\newcommand{\myvh}{\mathbf{h}}
  550. \newcommand{\vi}{\mathbf{i}}
  551. \newcommand{\vj}{\mathbf{j}}
  552. \newcommand{\vk}{\mathbf{k}}
  553. \newcommand{\vl}{\mathbf{l}}
  554. \newcommand{\vm}{\mathbf{m}}
  555. \newcommand{\vn}{\mathbf{n}}
  556. \newcommand{\vo}{\mathbf{o}}
  557. \newcommand{\vp}{\mathbf{p}}
  558. \newcommand{\vq}{\mathbf{q}}
  559. \newcommand{\vr}{\mathbf{r}}
  560. \newcommand{\vs}{\mathbf{s}}
  561. \newcommand{\vt}{\mathbf{t}}
  562. \newcommand{\vu}{\mathbf{u}}
  563. \newcommand{\vv}{\mathbf{v}}
  564. \newcommand{\vw}{\mathbf{w}}
  565. \newcommand{\vws}{\vw_s}
  566. \newcommand{\vwt}{\mathbf{\tilde{w}}}
  567. \newcommand{\vWt}{\mathbf{\tilde{W}}}
  568. \newcommand{\vwh}{\hat{\vw}}
  569. \newcommand{\vx}{\mathbf{x}}
  570. %\newcommand{\vx}{\mathbf{x}}
  571. \newcommand{\vxt}{\mathbf{\tilde{x}}}
  572. \newcommand{\vy}{\mathbf{y}}
  573. \newcommand{\vyt}{\mathbf{\tilde{y}}}
  574. \newcommand{\vz}{\mathbf{z}}
  575. %\newcommand{\vzt}{\mathbf{\tilde{z}}}\\
  576. % Upper Roman (Matrices)
  577. \newcommand{\vA}{\mathbf{A}}
  578. \newcommand{\vB}{\mathbf{B}}
  579. \newcommand{\vC}{\mathbf{C}}
  580. \newcommand{\vD}{\mathbf{D}}
  581. \newcommand{\vE}{\mathbf{E}}
  582. \newcommand{\vF}{\mathbf{F}}
  583. \newcommand{\vG}{\mathbf{G}}
  584. \newcommand{\vH}{\mathbf{H}}
  585. \newcommand{\vI}{\mathbf{I}}
  586. \newcommand{\vJ}{\mathbf{J}}
  587. \newcommand{\vK}{\mathbf{K}}
  588. \newcommand{\vL}{\mathbf{L}}
  589. \newcommand{\vM}{\mathbf{M}}
  590. \newcommand{\vMt}{\mathbf{\tilde{M}}}
  591. \newcommand{\vN}{\mathbf{N}}
  592. \newcommand{\vO}{\mathbf{O}}
  593. \newcommand{\vP}{\mathbf{P}}
  594. \newcommand{\vQ}{\mathbf{Q}}
  595. \newcommand{\vR}{\mathbf{R}}
  596. \newcommand{\vS}{\mathbf{S}}
  597. \newcommand{\vT}{\mathbf{T}}
  598. \newcommand{\vU}{\mathbf{U}}
  599. \newcommand{\vV}{\mathbf{V}}
  600. \newcommand{\vW}{\mathbf{W}}
  601. \newcommand{\vX}{\mathbf{X}}
  602. %\newcommand{\vXs}{\vX_{\vs}}
  603. \newcommand{\vXs}{\vX_{s}}
  604. \newcommand{\vXt}{\mathbf{\tilde{X}}}
  605. \newcommand{\vY}{\mathbf{Y}}
  606. \newcommand{\vZ}{\mathbf{Z}}
  607. \newcommand{\vZt}{\mathbf{\tilde{Z}}}
  608. \newcommand{\vzt}{\mathbf{\tilde{z}}}\\
  609. %%%%
  610. \newcommand{\hidden}{\vz}
  611. \newcommand{\hid}{\hidden}
  612. \newcommand{\observed}{\vy}
  613. \newcommand{\obs}{\observed}
  614. \newcommand{\inputs}{\vu}
  615. \newcommand{\input}{\inputs}\\\newcommand{\hmmTrans}{\vA}
  616. \newcommand{\hmmObs}{\vB}
  617. \newcommand{\hmmInit}{\vpi}
  618. \newcommand{\hmmhid}{\hidden}
  619. \newcommand{\hmmobs}{\obs}\\\newcommand{\ldsDyn}{\vA}
  620. \newcommand{\ldsObs}{\vC}
  621. \newcommand{\ldsDynIn}{\vB}
  622. \newcommand{\ldsObsIn}{\vD}
  623. \newcommand{\ldsDynNoise}{\vQ}
  624. \newcommand{\ldsObsNoise}{\vR}\\\newcommand{\ssmDynFn}{f}
  625. \newcommand{\ssmObsFn}{h}\\
  626. %%%
  627. \newcommand{\gauss}{\mathcal{N}}\\\newcommand{\diag}{\mathrm{diag}}\end{aligned}\end{align} \]</div>
  628. <div class="tex2jax_ignore mathjax_ignore section" id="linear-gaussian-ssms">
  629. <span id="sec-lds-intro"></span><h1>Linear Gaussian SSMs<a class="headerlink" href="#linear-gaussian-ssms" title="Permalink to this headline">¶</a></h1>
  630. <p>Consider the state space model in
  631. <a class="reference internal" href="ssm_intro.html#equation-eq-ssm-ar">(2)</a>
  632. where we assume the observations are conditionally iid given the
  633. hidden states and inputs (i.e. there are no auto-regressive dependencies
  634. between the observables).
  635. We can rewrite this model as
  636. a stochastic nonlinear dynamical system (NLDS)
  637. by defining the distribution of the next hidden state
  638. as a deterministic function of the past state
  639. plus random process noise <span class="math notranslate nohighlight">\(\vepsilon_t\)</span></p>
  640. <div class="amsmath math notranslate nohighlight" id="equation-e31147c5-366d-4025-ad02-0d0886a1dffa">
  641. <span class="eqno">(5)<a class="headerlink" href="#equation-e31147c5-366d-4025-ad02-0d0886a1dffa" title="Permalink to this equation">¶</a></span>\[\begin{align}
  642. \hmmhid_t &amp;= \ssmDynFn(\hmmhid_{t-1}, \inputs_t, \vepsilon_t)
  643. \end{align}\]</div>
  644. <p>where <span class="math notranslate nohighlight">\(\vepsilon_t\)</span> is drawn from the distribution such
  645. that the induced distribution
  646. on <span class="math notranslate nohighlight">\(\hmmhid_t\)</span> matches <span class="math notranslate nohighlight">\(p(\hmmhid_t|\hmmhid_{t-1}, \inputs_t)\)</span>.
  647. Similarly we can rewrite the observation distributions
  648. as a deterministic function of the hidden state
  649. plus observation noise <span class="math notranslate nohighlight">\(\veta_t\)</span>:</p>
  650. <div class="amsmath math notranslate nohighlight" id="equation-f582ebd3-bbc8-4436-8a13-41e99905dfbd">
  651. <span class="eqno">(6)<a class="headerlink" href="#equation-f582ebd3-bbc8-4436-8a13-41e99905dfbd" title="Permalink to this equation">¶</a></span>\[\begin{align}
  652. \hmmobs_t &amp;= \ssmObsFn(\hmmhid_{t}, \inputs_t, \veta_t)
  653. \end{align}\]</div>
  654. <p>If we assume additive Gaussian noise,
  655. the model becomes</p>
  656. <div class="amsmath math notranslate nohighlight" id="equation-09c96fd9-5478-4aa7-ac38-43504373febc">
  657. <span class="eqno">(7)<a class="headerlink" href="#equation-09c96fd9-5478-4aa7-ac38-43504373febc" title="Permalink to this equation">¶</a></span>\[\begin{align}
  658. \hmmhid_t &amp;= \ssmDynFn(\hmmhid_{t-1}, \inputs_t) + \vepsilon_t \\
  659. \hmmobs_t &amp;= \ssmObsFn(\hmmhid_{t}, \inputs_t) + \veta_t
  660. \end{align}\]</div>
  661. <p>where <span class="math notranslate nohighlight">\(\vepsilon_t \sim \gauss(\vzero,\vQ_t)\)</span>
  662. and <span class="math notranslate nohighlight">\(\veta_t \sim \gauss(\vzero,\vR_t)\)</span>.
  663. We will call these Gaussian SSMs.</p>
  664. <p>If we additionally assume
  665. the transition function <span class="math notranslate nohighlight">\(\ssmDynFn\)</span>
  666. and the observation function <span class="math notranslate nohighlight">\(\ssmObsFn\)</span> are both linear,
  667. then we can rewrite the model as follows:</p>
  668. <div class="amsmath math notranslate nohighlight" id="equation-44914eac-1d98-497d-87e7-1d6b5346a8f2">
  669. <span class="eqno">(8)<a class="headerlink" href="#equation-44914eac-1d98-497d-87e7-1d6b5346a8f2" title="Permalink to this equation">¶</a></span>\[\begin{align}
  670. p(\hmmhid_t|\hmmhid_{t-1},\inputs_t) &amp;= \gauss(\hmmhid_t|\ldsDyn_t \hmmhid_{t-1}
  671. + \ldsDynIn_t \inputs_t, \vQ_t)
  672. \\
  673. p(\hmmobs_t|\hmmhid_t,\inputs_t) &amp;= \gauss(\hmmobs_t|\ldsObs_t \hmmhid_{t}
  674. + \ldsObsIn_t \inputs_t, \vR_t)
  675. \end{align}\]</div>
  676. <p>This is called a
  677. linear-Gaussian state space model
  678. (LG-SSM),
  679. or a
  680. linear dynamical system (LDS).
  681. We usually assume the parameters are independent of time, in which case
  682. the model is said to be time-invariant or homogeneous.</p>
  683. <div class="section" id="example-tracking-a-2d-point">
  684. <span id="sec-kalman-tracking"></span><span id="sec-tracking-lds"></span><h2>Example: tracking a 2d point<a class="headerlink" href="#example-tracking-a-2d-point" title="Permalink to this headline">¶</a></h2>
  685. <p>Consider an object moving in <span class="math notranslate nohighlight">\(\real^2\)</span>.
  686. Let the state be
  687. the position and velocity of the object,
  688. <span class="math notranslate nohighlight">\(\vz_t =\begin{pmatrix} u_t &amp; \dot{u}_t &amp; v_t &amp; \dot{v}_t \end{pmatrix}\)</span>.
  689. (We use <span class="math notranslate nohighlight">\(u\)</span> and <span class="math notranslate nohighlight">\(v\)</span> for the two coordinates,
  690. to avoid confusion with the state and observation variables.)
  691. If we use Euler discretization,
  692. the dynamics become</p>
  693. <div class="amsmath math notranslate nohighlight" id="equation-56d91323-8e7f-4b13-8792-ca0ca145de95">
  694. <span class="eqno">(9)<a class="headerlink" href="#equation-56d91323-8e7f-4b13-8792-ca0ca145de95" title="Permalink to this equation">¶</a></span>\[\begin{align}
  695. \underbrace{\begin{pmatrix} u_t\\ \dot{u}_t \\ v_t \\ \dot{v}_t \end{pmatrix}}_{\vz_t}
  696. =
  697. \underbrace{
  698. \begin{pmatrix}
  699. 1 &amp; 0 &amp; \Delta &amp; 0 \\
  700. 0 &amp; 1 &amp; 0 &amp; \Delta\\
  701. 0 &amp; 0 &amp; 1 &amp; 0 \\
  702. 0 &amp; 0 &amp; 0 &amp; 1
  703. \end{pmatrix}
  704. }_{\ldsDyn}
  705. \
  706. \underbrace{\begin{pmatrix} u_{t-1} \\ \dot{u}_{t-1} \\ v_{t-1} \\ \dot{v}_{t-1} \end{pmatrix}}_{\vz_{t-1}}
  707. + \vepsilon_t
  708. \end{align}\]</div>
  709. <p>where <span class="math notranslate nohighlight">\(\vepsilon_t \sim \gauss(\vzero,\vQ)\)</span> is
  710. the process noise.</p>
  711. <p>Let us assume
  712. that the process noise is
  713. a white noise process added to the velocity components
  714. of the state, but not to the location.
  715. (This is known as a random accelerations model.)
  716. We can approximate the resulting process in discrete time by assuming
  717. <span class="math notranslate nohighlight">\(\vQ = \diag(0, q, 0, q)\)</span>.
  718. (See <span id="id1">[<a class="reference internal" href="../../bib.html#id18" title="Simo Sarkka. Bayesian Filtering and Smoothing. Cambridge University Press, 2013. URL: https://users.aalto.fi/~ssarkka/pub/cup_book_online_20131111.pdf.">Sar13</a>]</span> p60 for a more accurate way
  719. to convert the continuous time process to discrete time.)</p>
  720. <p>Now suppose that at each discrete time point we
  721. observe the location,
  722. corrupted by Gaussian noise.
  723. Thus the observation model becomes</p>
  724. <div class="amsmath math notranslate nohighlight" id="equation-54905b1b-79c0-47fc-b0cc-f790a1e50550">
  725. <span class="eqno">(10)<a class="headerlink" href="#equation-54905b1b-79c0-47fc-b0cc-f790a1e50550" title="Permalink to this equation">¶</a></span>\[\begin{align}
  726. \underbrace{\begin{pmatrix} y_{1,t} \\ y_{2,t} \end{pmatrix}}_{\vy_t}
  727. &amp;=
  728. \underbrace{
  729. \begin{pmatrix}
  730. 1 &amp; 0 &amp; 0 &amp; 0 \\
  731. 0 &amp; 0 &amp; 1 &amp; 0
  732. \end{pmatrix}
  733. }_{\ldsObs}
  734. \
  735. \underbrace{\begin{pmatrix} u_t\\ \dot{u}_t \\ v_t \\ \dot{v}_t \end{pmatrix}}_{\vz_t}
  736. + \veta_t
  737. \end{align}\]</div>
  738. <p>where <span class="math notranslate nohighlight">\(\veta_t \sim \gauss(\vzero,\vR)\)</span> is the \keywordDef{observation noise}.
  739. We see that the observation matrix <span class="math notranslate nohighlight">\(\ldsObs\)</span> simply ``extracts’’ the
  740. relevant parts of the state vector.</p>
  741. <p>Suppose we sample a trajectory and corresponding set
  742. of noisy observations from this model,
  743. <span class="math notranslate nohighlight">\((\vz_{1:T}, \vy_{1:T}) \sim p(\vz,\vy|\vtheta)\)</span>.
  744. (We use diagonal observation noise,
  745. <span class="math notranslate nohighlight">\(\vR = \diag(\sigma_1^2, \sigma_2^2)\)</span>.)
  746. The results are shown below.</p>
  747. <div class="cell docutils container">
  748. <div class="cell_input docutils container">
  749. <div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="n">key</span> <span class="o">=</span> <span class="n">jax</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">PRNGKey</span><span class="p">(</span><span class="mi">314</span><span class="p">)</span>
  750. <span class="n">timesteps</span> <span class="o">=</span> <span class="mi">15</span>
  751. <span class="n">delta</span> <span class="o">=</span> <span class="mf">1.0</span>
  752. <span class="n">A</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
  753. <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">delta</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
  754. <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">delta</span><span class="p">],</span>
  755. <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
  756. <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
  757. <span class="p">])</span>
  758. <span class="n">C</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
  759. <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
  760. <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span>
  761. <span class="p">])</span>
  762. <span class="n">state_size</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">A</span><span class="o">.</span><span class="n">shape</span>
  763. <span class="n">observation_size</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">C</span><span class="o">.</span><span class="n">shape</span>
  764. <span class="n">Q</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">state_size</span><span class="p">)</span> <span class="o">*</span> <span class="mf">0.001</span>
  765. <span class="n">R</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">observation_size</span><span class="p">)</span> <span class="o">*</span> <span class="mf">1.0</span>
  766. <span class="c1"># Prior parameter distribution</span>
  767. <span class="n">mu0</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">8</span><span class="p">,</span> <span class="mi">10</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">])</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="nb">float</span><span class="p">)</span>
  768. <span class="n">Sigma0</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">state_size</span><span class="p">)</span> <span class="o">*</span> <span class="mf">1.0</span>
  769. <span class="kn">from</span> <span class="nn">jsl.lds.kalman_filter</span> <span class="kn">import</span> <span class="n">LDS</span><span class="p">,</span> <span class="n">smooth</span><span class="p">,</span> <span class="nb">filter</span>
  770. <span class="n">lds</span> <span class="o">=</span> <span class="n">LDS</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">C</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">mu0</span><span class="p">,</span> <span class="n">Sigma0</span><span class="p">)</span>
  771. <span class="nb">print</span><span class="p">(</span><span class="n">lds</span><span class="p">)</span>
  772. </pre></div>
  773. </div>
  774. </div>
  775. <div class="cell_output docutils container">
  776. <div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
  777. </pre></div>
  778. </div>
  779. <div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>LDS(A=DeviceArray([[1., 0., 1., 0.],
  780. [0., 1., 0., 1.],
  781. [0., 0., 1., 0.],
  782. [0., 0., 0., 1.]], dtype=float32), C=DeviceArray([[1, 0, 0, 0],
  783. [0, 1, 0, 0]], dtype=int32), Q=DeviceArray([[0.001, 0. , 0. , 0. ],
  784. [0. , 0.001, 0. , 0. ],
  785. [0. , 0. , 0.001, 0. ],
  786. [0. , 0. , 0. , 0.001]], dtype=float32), R=DeviceArray([[1., 0.],
  787. [0., 1.]], dtype=float32), mu=DeviceArray([ 8., 10., 1., 0.], dtype=float32), Sigma=DeviceArray([[1., 0., 0., 0.],
  788. [0., 1., 0., 0.],
  789. [0., 0., 1., 0.],
  790. [0., 0., 0., 1.]], dtype=float32), state_offset=None, obs_offset=None, nstates=4, nobs=2)
  791. </pre></div>
  792. </div>
  793. </div>
  794. </div>
  795. <div class="cell docutils container">
  796. <div class="cell_input docutils container">
  797. <div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">jsl.demos.plot_utils</span> <span class="kn">import</span> <span class="n">plot_ellipse</span>
  798. <span class="k">def</span> <span class="nf">plot_tracking_values</span><span class="p">(</span><span class="n">observed</span><span class="p">,</span> <span class="n">filtered</span><span class="p">,</span> <span class="n">cov_hist</span><span class="p">,</span> <span class="n">signal_label</span><span class="p">,</span> <span class="n">ax</span><span class="p">):</span>
  799. <span class="n">timesteps</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">observed</span><span class="o">.</span><span class="n">shape</span>
  800. <span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">observed</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">observed</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">marker</span><span class="o">=</span><span class="s2">&quot;o&quot;</span><span class="p">,</span> <span class="n">linewidth</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span>
  801. <span class="n">markerfacecolor</span><span class="o">=</span><span class="s2">&quot;none&quot;</span><span class="p">,</span> <span class="n">markeredgewidth</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">markersize</span><span class="o">=</span><span class="mi">8</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;observed&quot;</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="s2">&quot;tab:green&quot;</span><span class="p">)</span>
  802. <span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="o">*</span><span class="n">filtered</span><span class="p">[:,</span> <span class="p">:</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="n">signal_label</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="s2">&quot;tab:red&quot;</span><span class="p">,</span> <span class="n">marker</span><span class="o">=</span><span class="s2">&quot;x&quot;</span><span class="p">,</span> <span class="n">linewidth</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span>
  803. <span class="k">for</span> <span class="n">t</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">timesteps</span><span class="p">,</span> <span class="mi">1</span><span class="p">):</span>
  804. <span class="n">covn</span> <span class="o">=</span> <span class="n">cov_hist</span><span class="p">[</span><span class="n">t</span><span class="p">][:</span><span class="mi">2</span><span class="p">,</span> <span class="p">:</span><span class="mi">2</span><span class="p">]</span>
  805. <span class="n">plot_ellipse</span><span class="p">(</span><span class="n">covn</span><span class="p">,</span> <span class="n">filtered</span><span class="p">[</span><span class="n">t</span><span class="p">,</span> <span class="p">:</span><span class="mi">2</span><span class="p">],</span> <span class="n">ax</span><span class="p">,</span> <span class="n">n_std</span><span class="o">=</span><span class="mf">2.0</span><span class="p">,</span> <span class="n">plot_center</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
  806. <span class="n">ax</span><span class="o">.</span><span class="n">axis</span><span class="p">(</span><span class="s2">&quot;equal&quot;</span><span class="p">)</span>
  807. <span class="n">ax</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
  808. </pre></div>
  809. </div>
  810. </div>
  811. </div>
  812. <div class="cell docutils container">
  813. <div class="cell_input docutils container">
  814. <div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="n">z_hist</span><span class="p">,</span> <span class="n">x_hist</span> <span class="o">=</span> <span class="n">lds</span><span class="o">.</span><span class="n">sample</span><span class="p">(</span><span class="n">key</span><span class="p">,</span> <span class="n">timesteps</span><span class="p">)</span>
  815. <span class="n">fig_truth</span><span class="p">,</span> <span class="n">axs</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">()</span>
  816. <span class="n">axs</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x_hist</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">x_hist</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">],</span>
  817. <span class="n">marker</span><span class="o">=</span><span class="s2">&quot;o&quot;</span><span class="p">,</span> <span class="n">linewidth</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">markerfacecolor</span><span class="o">=</span><span class="s2">&quot;none&quot;</span><span class="p">,</span>
  818. <span class="n">markeredgewidth</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">markersize</span><span class="o">=</span><span class="mi">8</span><span class="p">,</span>
  819. <span class="n">label</span><span class="o">=</span><span class="s2">&quot;observed&quot;</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="s2">&quot;tab:green&quot;</span><span class="p">)</span>
  820. <span class="n">axs</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">z_hist</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">z_hist</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">],</span>
  821. <span class="n">linewidth</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">&quot;truth&quot;</span><span class="p">,</span>
  822. <span class="n">marker</span><span class="o">=</span><span class="s2">&quot;s&quot;</span><span class="p">,</span> <span class="n">markersize</span><span class="o">=</span><span class="mi">8</span><span class="p">)</span>
  823. <span class="n">axs</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
  824. <span class="n">axs</span><span class="o">.</span><span class="n">axis</span><span class="p">(</span><span class="s2">&quot;equal&quot;</span><span class="p">)</span>
  825. </pre></div>
  826. </div>
  827. </div>
  828. <div class="cell_output docutils container">
  829. <div class="output text_plain highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>(7.24486608505249, 23.857812213897706, 8.042076778411865, 11.636079120635987)
  830. </pre></div>
  831. </div>
  832. <img alt="../../_images/lds_7_1.png" src="../../_images/lds_7_1.png" />
  833. </div>
  834. </div>
  835. <p>The main task is to infer the hidden states given the noisy
  836. observations, i.e., <span class="math notranslate nohighlight">\(p(\vz|\vy,\vtheta)\)</span>. We discuss the topic of inference in <a class="reference internal" href="inference.html#sec-inference"><span class="std std-ref">Inferential goals</span></a>.</p>
  837. </div>
  838. </div>
  839. <script type="text/x-thebe-config">
  840. {
  841. requestKernel: true,
  842. binderOptions: {
  843. repo: "binder-examples/jupyter-stacks-datascience",
  844. ref: "master",
  845. },
  846. codeMirrorConfig: {
  847. theme: "abcdef",
  848. mode: "python"
  849. },
  850. kernelOptions: {
  851. kernelName: "python3",
  852. path: "./chapters/ssm"
  853. },
  854. predefinedOutput: true
  855. }
  856. </script>
  857. <script>kernelName = 'python3'</script>
  858. </div>
  859. <!-- Previous / next buttons -->
  860. <div class='prev-next-area'>
  861. <a class='left-prev' id="prev-link" href="hmm.html" title="previous page">
  862. <i class="fas fa-angle-left"></i>
  863. <div class="prev-next-info">
  864. <p class="prev-next-subtitle">previous</p>
  865. <p class="prev-next-title">Hidden Markov Models</p>
  866. </div>
  867. </a>
  868. <a class='right-next' id="next-link" href="nlds.html" title="next page">
  869. <div class="prev-next-info">
  870. <p class="prev-next-subtitle">next</p>
  871. <p class="prev-next-title">Nonlinear Gaussian SSMs</p>
  872. </div>
  873. <i class="fas fa-angle-right"></i>
  874. </a>
  875. </div>
  876. </div>
  877. </div>
  878. <footer class="footer">
  879. <p>
  880. By Kevin Murphy, Scott Linderman, et al.<br/>
  881. &copy; Copyright 2021.<br/>
  882. </p>
  883. </footer>
  884. </main>
  885. </div>
  886. </div>
  887. <script src="../../_static/js/index.be7d3bbb2ef33a8344ce.js"></script>
  888. </body>
  889. </html>