radu
/
SSMBook
mirror da https://github.com/probml/ssm-book.git


			
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856
							
<!DOCTYPE html>

<html>
  <head>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0" />
    <title>&lt;no title&gt; &#8212; State Space Models: A Modern Approach</title>
    
  <link href="../_static/css/theme.css" rel="stylesheet">
  <link href="../_static/css/index.ff1ffe594081f20da1ef19478df9384b.css" rel="stylesheet">

    
  <link rel="stylesheet"
    href="../_static/vendor/fontawesome/5.13.0/css/all.min.css">
  <link rel="preload" as="font" type="font/woff2" crossorigin
    href="../_static/vendor/fontawesome/5.13.0/webfonts/fa-solid-900.woff2">
  <link rel="preload" as="font" type="font/woff2" crossorigin
    href="../_static/vendor/fontawesome/5.13.0/webfonts/fa-brands-400.woff2">

    
    <link rel="stylesheet" type="text/css" href="../_static/pygments.css" />
    <link rel="stylesheet" type="text/css" href="../_static/sphinx-book-theme.css?digest=c3fdc42140077d1ad13ad2f1588a4309" />
    <link rel="stylesheet" type="text/css" href="../_static/togglebutton.css" />
    <link rel="stylesheet" type="text/css" href="../_static/copybutton.css" />
    <link rel="stylesheet" type="text/css" href="../_static/mystnb.css" />
    <link rel="stylesheet" type="text/css" href="../_static/sphinx-thebe.css" />
    <link rel="stylesheet" type="text/css" href="../_static/panels-main.c949a650a448cc0ae9fd3441c0e17fb0.css" />
    <link rel="stylesheet" type="text/css" href="../_static/panels-variables.06eb56fa6e07937060861dad626602ad.css" />
    
  <link rel="preload" as="script" href="../_static/js/index.be7d3bbb2ef33a8344ce.js">

    <script data-url_root="../" id="documentation_options" src="../_static/documentation_options.js"></script>
    <script src="../_static/jquery.js"></script>
    <script src="../_static/underscore.js"></script>
    <script src="../_static/doctools.js"></script>
    <script src="../_static/clipboard.min.js"></script>
    <script src="../_static/copybutton.js"></script>
    <script>let toggleHintShow = 'Click to show';</script>
    <script>let toggleHintHide = 'Click to hide';</script>
    <script>let toggleOpenOnPrint = 'true';</script>
    <script src="../_static/togglebutton.js"></script>
    <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>
    <script src="../_static/sphinx-book-theme.d59cb220de22ca1c485ebbdc042f0030.js"></script>
    <script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"
const thebe_selector = ".thebe,.cell"
const thebe_selector_input = "pre"
const thebe_selector_output = ".output, .cell_output"
</script>
    <script async="async" src="../_static/sphinx-thebe.js"></script>
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
    <meta name="viewport" content="width=device-width, initial-scale=1" />
    <meta name="docsearch:language" content="None">
    

    <!-- Google Analytics -->
    
  </head>
  <body data-spy="scroll" data-target="#bd-toc-nav" data-offset="80">
    
    <div class="container-fluid" id="banner"></div>

    
    <div class="container-xl">
      <div class="row">
          
<div class="col-12 col-md-3 bd-sidebar site-navigation show" id="site-navigation">
    
        <div class="navbar-brand-box">
    <a class="navbar-brand text-wrap" href="../index.html">
      
      
      <h1 class="site-logo" id="site-title">State Space Models: A Modern Approach</h1>
      
    </a>
</div><form class="bd-search d-flex align-items-center" action="../search.html" method="get">
  <i class="icon fas fa-search"></i>
  <input type="search" class="form-control" name="q" id="search-input" placeholder="Search this book..." aria-label="Search this book..." autocomplete="off" >
</form><nav class="bd-links" id="bd-docs-nav" aria-label="Main">
    <div class="bd-toc-item active">
        <ul class="nav bd-sidenav">
 <li class="toctree-l1">
  <a class="reference internal" href="../root.html">
   State Space Models: A Modern Approach
  </a>
 </li>
</ul>
<ul class="nav bd-sidenav">
 <li class="toctree-l1 has-children">
  <a class="reference internal" href="ssm/ssm_index.html">
   State Space Models
  </a>
  <input class="toctree-checkbox" id="toctree-checkbox-1" name="toctree-checkbox-1" type="checkbox"/>
  <label for="toctree-checkbox-1">
   <i class="fas fa-chevron-down">
   </i>
  </label>
  <ul>
   <li class="toctree-l2">
    <a class="reference internal" href="ssm/ssm_intro.html">
     What are State Space Models?
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="ssm/hmm.html">
     Hidden Markov Models
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="ssm/lds.html">
     Linear Gaussian SSMs
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="ssm/nlds.html">
     Nonlinear Gaussian SSMs
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="ssm/inference.html">
     States estimation (inference)
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="ssm/learning.html">
     Parameter estimation (learning)
    </a>
   </li>
  </ul>
 </li>
 <li class="toctree-l1 has-children">
  <a class="reference internal" href="hmm/hmm_index.html">
   Hidden Markov Models
  </a>
  <input class="toctree-checkbox" id="toctree-checkbox-2" name="toctree-checkbox-2" type="checkbox"/>
  <label for="toctree-checkbox-2">
   <i class="fas fa-chevron-down">
   </i>
  </label>
  <ul>
   <li class="toctree-l2">
    <a class="reference internal" href="hmm/hmm_filter.html">
     HMM filtering (forwards algorithm)
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="hmm/hmm_smoother.html">
     HMM smoothing (forwards-backwards algorithm)
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="hmm/hmm_viterbi.html">
     Viterbi algorithm
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="hmm/hmm_parallel.html">
     Parallel HMM  smoothing
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="hmm/hmm_sampling.html">
     Forwards-filtering backwards-sampling algorithm
    </a>
   </li>
  </ul>
 </li>
 <li class="toctree-l1 has-children">
  <a class="reference internal" href="lgssm/lgssm_index.html">
   Linear-Gaussian SSMs
  </a>
  <input class="toctree-checkbox" id="toctree-checkbox-3" name="toctree-checkbox-3" type="checkbox"/>
  <label for="toctree-checkbox-3">
   <i class="fas fa-chevron-down">
   </i>
  </label>
  <ul>
   <li class="toctree-l2">
    <a class="reference internal" href="lgssm/kalman_filter.html">
     Kalman filtering
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="lgssm/kalman_smoother.html">
     Kalman (RTS) smoother
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="lgssm/kalman_parallel.html">
     Parallel Kalman Smoother
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="lgssm/kalman_sampling.html">
     Forwards-filtering backwards sampling
    </a>
   </li>
  </ul>
 </li>
 <li class="toctree-l1 has-children">
  <a class="reference internal" href="extended/extended_index.html">
   Extended (linearized) methods
  </a>
  <input class="toctree-checkbox" id="toctree-checkbox-4" name="toctree-checkbox-4" type="checkbox"/>
  <label for="toctree-checkbox-4">
   <i class="fas fa-chevron-down">
   </i>
  </label>
  <ul>
   <li class="toctree-l2">
    <a class="reference internal" href="extended/extended_filter.html">
     Extended Kalman filtering
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="extended/extended_smoother.html">
     Extended Kalman smoother
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="extended/extended_parallel.html">
     Parallel extended Kalman smoothing
    </a>
   </li>
  </ul>
 </li>
 <li class="toctree-l1 has-children">
  <a class="reference internal" href="unscented/unscented_index.html">
   Unscented methods
  </a>
  <input class="toctree-checkbox" id="toctree-checkbox-5" name="toctree-checkbox-5" type="checkbox"/>
  <label for="toctree-checkbox-5">
   <i class="fas fa-chevron-down">
   </i>
  </label>
  <ul>
   <li class="toctree-l2">
    <a class="reference internal" href="unscented/unscented_filter.html">
     Unscented filtering
    </a>
   </li>
   <li class="toctree-l2">
    <a class="reference internal" href="unscented/unscented_smoother.html">
     Unscented smoothing
    </a>
   </li>
  </ul>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="quadrature/quadrature_index.html">
   Quadrature and cubature methods
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="postlin/postlin_index.html">
   Posterior linearization
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="adf/adf_index.html">
   Assumed Density Filtering
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="vi/vi_index.html">
   Variational inference
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="pf/pf_index.html">
   Particle filtering
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="smc/smc_index.html">
   Sequential Monte Carlo
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="learning/learning_index.html">
   Offline parameter estimation (learning)
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="tracking/tracking_index.html">
   Multi-target tracking
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="ensemble/ensemble_index.html">
   Data assimilation using Ensemble Kalman filter
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="bnp/bnp_index.html">
   Bayesian non-parametric SSMs
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="changepoint/changepoint_index.html">
   Changepoint detection
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="timeseries/timeseries_index.html">
   Timeseries forecasting
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="gp/gp_index.html">
   Markovian Gaussian processes
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="ode/ode_index.html">
   Differential equations and SSMs
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="control/control_index.html">
   Optimal control
  </a>
 </li>
 <li class="toctree-l1">
  <a class="reference internal" href="../bib.html">
   Bibliography
  </a>
 </li>
</ul>

    </div>
</nav> <!-- To handle the deprecated key -->

<div class="navbar_extra_footer">
  Powered by <a href="https://jupyterbook.org">Jupyter Book</a>
</div>

</div>


<main class="col py-md-3 pl-md-4 bd-content overflow-auto" role="main">
    
    <div class="topbar container-xl fixed-top">
    <div class="topbar-contents row">
        <div class="col-12 col-md-3 bd-topbar-whitespace site-navigation show"></div>
        <div class="col pl-md-4 topbar-main">
            
            <button id="navbar-toggler" class="navbar-toggler ml-0" type="button" data-toggle="collapse"
                data-toggle="tooltip" data-placement="bottom" data-target=".site-navigation" aria-controls="navbar-menu"
                aria-expanded="true" aria-label="Toggle navigation" aria-controls="site-navigation"
                title="Toggle navigation" data-toggle="tooltip" data-placement="left">
                <i class="fas fa-bars"></i>
                <i class="fas fa-arrow-left"></i>
                <i class="fas fa-arrow-up"></i>
            </button>
            
            
<div class="dropdown-buttons-trigger">
    <button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn" aria-label="Download this page"><i
            class="fas fa-download"></i></button>

    <div class="dropdown-buttons">
        <!-- ipynb file if we had a myst markdown file -->
        
        <!-- Download raw file -->
        <a class="dropdown-buttons" href="../_sources/chapters/scratchpad.ipynb"><button type="button"
                class="btn btn-secondary topbarbtn" title="Download source file" data-toggle="tooltip"
                data-placement="left">.ipynb</button></a>
        <!-- Download PDF via print -->
        <button type="button" id="download-print" class="btn btn-secondary topbarbtn" title="Print to PDF"
                onclick="printPdf(this)" data-toggle="tooltip" data-placement="left">.pdf</button>
    </div>
</div>

            <!-- Source interaction buttons -->

<div class="dropdown-buttons-trigger">
    <button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn"
        aria-label="Connect with source repository"><i class="fab fa-github"></i></button>
    <div class="dropdown-buttons sourcebuttons">
        <a class="repository-button"
            href="https://github.com/probml/ssm-book"><button type="button" class="btn btn-secondary topbarbtn"
                data-toggle="tooltip" data-placement="left" title="Source repository"><i
                    class="fab fa-github"></i>repository</button></a>
        <a class="issues-button"
            href="https://github.com/probml/ssm-book/issues/new?title=Issue%20on%20page%20%2Fchapters/scratchpad.html&body=Your%20issue%20content%20here."><button
                type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip" data-placement="left"
                title="Open an issue"><i class="fas fa-lightbulb"></i>open issue</button></a>
        
    </div>
</div>

            <!-- Full screen (wrap in <a> to have style consistency -->

<a class="full-screen-button"><button type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip"
        data-placement="bottom" onclick="toggleFullScreen()" aria-label="Fullscreen mode"
        title="Fullscreen mode"><i
            class="fas fa-expand"></i></button></a>

            <!-- Launch buttons -->

<div class="dropdown-buttons-trigger">
    <button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn"
        aria-label="Launch interactive content"><i class="fas fa-rocket"></i></button>
    <div class="dropdown-buttons">
        
        <a class="binder-button" href="https://mybinder.org/v2/gh/probml/ssm-book/main?urlpath=tree/chapters/scratchpad.ipynb"><button type="button"
                class="btn btn-secondary topbarbtn" title="Launch Binder" data-toggle="tooltip"
                data-placement="left"><img class="binder-button-logo"
                    src="../_static/images/logo_binder.svg"
                    alt="Interact on binder">Binder</button></a>
        
        
        <a class="colab-button" href="https://colab.research.google.com/github/probml/ssm-book/blob/main/chapters/scratchpad.ipynb"><button type="button" class="btn btn-secondary topbarbtn"
                title="Launch Colab" data-toggle="tooltip" data-placement="left"><img class="colab-button-logo"
                    src="../_static/images/logo_colab.png"
                    alt="Interact on Colab">Colab</button></a>
        
        
    </div>
</div>

        </div>

        <!-- Table of contents -->
        <div class="d-none d-md-block col-md-2 bd-toc show noprint">
            
            <div class="tocsection onthispage pt-5 pb-3">
                <i class="fas fa-list"></i> Contents
            </div>
            <nav id="bd-toc-nav" aria-label="Page">
                <ul class="simple visible nav section-nav flex-column">
</ul>

            </nav>
        </div>
    </div>
</div>
    <div id="main-content" class="row">
        <div class="col-12 col-md-9 pl-md-3 pr-md-0">
            <!-- Table of contents that is only displayed when printing the page -->
            <div id="jb-print-docs-body" class="onlyprint">
                <h1><no title></h1>
                <!-- Table of contents -->
                <div id="print-main-content">
                    <div id="jb-print-toc">
                        
                        <div>
                            <h2> Contents </h2>
                        </div>
                        <nav aria-label="Page">
                            <ul class="simple visible nav section-nav flex-column">
</ul>

                        </nav>
                    </div>
                </div>
            </div>
            
              <div>
                
  <div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1">### Import standard libraries</span>

<span class="kn">import</span> <span class="nn">abc</span>
<span class="kn">from</span> <span class="nn">dataclasses</span> <span class="kn">import</span> <span class="n">dataclass</span>
<span class="kn">import</span> <span class="nn">functools</span>
<span class="kn">from</span> <span class="nn">functools</span> <span class="kn">import</span> <span class="n">partial</span>
<span class="kn">import</span> <span class="nn">itertools</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<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>

<span class="kn">import</span> <span class="nn">jax</span>
<span class="kn">import</span> <span class="nn">jax.numpy</span> <span class="k">as</span> <span class="nn">jnp</span>
<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>
<span class="c1">#from jax.scipy.special import logit</span>
<span class="c1">#from jax.nn import softmax</span>
<span class="kn">import</span> <span class="nn">jax.random</span> <span class="k">as</span> <span class="nn">jr</span>


<span class="kn">import</span> <span class="nn">distrax</span>
<span class="kn">import</span> <span class="nn">optax</span>

<span class="kn">import</span> <span class="nn">jsl</span>
<span class="kn">import</span> <span class="nn">ssm_jax</span>
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">inspect</span>
<span class="kn">import</span> <span class="nn">inspect</span> <span class="k">as</span> <span class="nn">py_inspect</span>
<span class="kn">import</span> <span class="nn">rich</span>
<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>
<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>

<span class="k">def</span> <span class="nf">print_source</span><span class="p">(</span><span class="n">fname</span><span class="p">):</span>
    <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>
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># meta-data does not work yet in VScode</span>
<span class="c1"># https://github.com/microsoft/vscode-jupyter/issues/1121</span>

<span class="p">{</span>
    <span class="s2">&quot;tags&quot;</span><span class="p">:</span> <span class="p">[</span>
        <span class="s2">&quot;hide-cell&quot;</span>
    <span class="p">]</span>
<span class="p">}</span>


<span class="c1">### Install necessary libraries</span>

<span class="k">try</span><span class="p">:</span>
    <span class="kn">import</span> <span class="nn">jax</span>
<span class="k">except</span><span class="p">:</span>
    <span class="c1"># For cuda version, see https://github.com/google/jax#installation</span>
    <span class="o">%</span><span class="k">pip</span> install --upgrade &quot;jax[cpu]&quot; 
    <span class="kn">import</span> <span class="nn">jax</span>

<span class="k">try</span><span class="p">:</span>
    <span class="kn">import</span> <span class="nn">distrax</span>
<span class="k">except</span><span class="p">:</span>
    <span class="o">%</span><span class="k">pip</span> install --upgrade  distrax
    <span class="kn">import</span> <span class="nn">distrax</span>

<span class="k">try</span><span class="p">:</span>
    <span class="kn">import</span> <span class="nn">jsl</span>
<span class="k">except</span><span class="p">:</span>
    <span class="o">%</span><span class="k">pip</span> install git+https://github.com/probml/jsl
    <span class="kn">import</span> <span class="nn">jsl</span>

<span class="k">try</span><span class="p">:</span>
    <span class="kn">import</span> <span class="nn">rich</span>
<span class="k">except</span><span class="p">:</span>
    <span class="o">%</span><span class="k">pip</span> install rich
    <span class="kn">import</span> <span class="nn">rich</span>
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="p">{</span>
    <span class="s2">&quot;tags&quot;</span><span class="p">:</span> <span class="p">[</span>
        <span class="s2">&quot;hide-cell&quot;</span>
    <span class="p">]</span>
<span class="p">}</span>


<span class="c1">### Import standard libraries</span>

<span class="kn">import</span> <span class="nn">abc</span>
<span class="kn">from</span> <span class="nn">dataclasses</span> <span class="kn">import</span> <span class="n">dataclass</span>
<span class="kn">import</span> <span class="nn">functools</span>
<span class="kn">import</span> <span class="nn">itertools</span>

<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>

<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>


<span class="kn">import</span> <span class="nn">jax</span>
<span class="kn">import</span> <span class="nn">jax.numpy</span> <span class="k">as</span> <span class="nn">jnp</span>
<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>
<span class="kn">from</span> <span class="nn">jax.scipy.special</span> <span class="kn">import</span> <span class="n">logit</span>
<span class="kn">from</span> <span class="nn">jax.nn</span> <span class="kn">import</span> <span class="n">softmax</span>
<span class="kn">from</span> <span class="nn">functools</span> <span class="kn">import</span> <span class="n">partial</span>
<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>

<span class="kn">import</span> <span class="nn">inspect</span>
<span class="kn">import</span> <span class="nn">inspect</span> <span class="k">as</span> <span class="nn">py_inspect</span>
<span class="kn">import</span> <span class="nn">rich</span>
<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>
<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>

<span class="k">def</span> <span class="nf">print_source</span><span class="p">(</span><span class="n">fname</span><span class="p">):</span>
    <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>
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">ssm_jax</span>
<span class="kn">from</span> <span class="nn">ssm_jax.hmm.models</span> <span class="kn">import</span> <span class="n">GaussianHMM</span>

<span class="n">print_source</span><span class="p">(</span><span class="n">GaussianHMM</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_html"><pre style="white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace">class GaussianHMM<span style="font-weight: bold">(</span>BaseHMM<span style="font-weight: bold">)</span>:
    def __init__<span style="font-weight: bold">(</span>self,
                 initial_probabilities,
                 transition_matrix,
                 emission_means,
                 emission_covariance_matrices<span style="font-weight: bold">)</span>:
        <span style="color: #008000; text-decoration-color: #008000">""</span>"_summary_

        Args:
            initial_probabilities <span style="font-weight: bold">(</span>_type_<span style="font-weight: bold">)</span>: _description_
            transition_matrix <span style="font-weight: bold">(</span>_type_<span style="font-weight: bold">)</span>: _description_
            emission_means <span style="font-weight: bold">(</span>_type_<span style="font-weight: bold">)</span>: _description_
            emission_covariance_matrices <span style="font-weight: bold">(</span>_type_<span style="font-weight: bold">)</span>: _description_
        <span style="color: #008000; text-decoration-color: #008000">""</span>"
        super<span style="font-weight: bold">()</span>.__init__<span style="font-weight: bold">(</span>initial_probabilities,
                         transition_matrix<span style="font-weight: bold">)</span>

        self._emission_distribution = tfd.MultivariateNormalFullCovariance<span style="font-weight: bold">(</span>
            emission_means, emission_covariance_matrices<span style="font-weight: bold">)</span>

    @classmethod
    def random_initialization<span style="font-weight: bold">(</span>cls, key, num_states, emission_dim<span style="font-weight: bold">)</span>:
        key1, key2, key3 = jr.split<span style="font-weight: bold">(</span>key, <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">3</span><span style="font-weight: bold">)</span>
        initial_probs = jr.dirichlet<span style="font-weight: bold">(</span>key1, jnp.ones<span style="font-weight: bold">(</span>num_states<span style="font-weight: bold">))</span>
        transition_matrix = jr.dirichlet<span style="font-weight: bold">(</span>key2, jnp.ones<span style="font-weight: bold">(</span>num_states<span style="font-weight: bold">)</span>, <span style="font-weight: bold">(</span>num_states,<span style="font-weight: bold">))</span>
        emission_means = jr.normal<span style="font-weight: bold">(</span>key3, <span style="font-weight: bold">(</span>num_states, emission_dim<span style="font-weight: bold">))</span>
        emission_covs = jnp.tile<span style="font-weight: bold">(</span>jnp.eye<span style="font-weight: bold">(</span>emission_dim<span style="font-weight: bold">)</span>, <span style="font-weight: bold">(</span>num_states, <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">1</span>, <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">1</span><span style="font-weight: bold">))</span>
        return cls<span style="font-weight: bold">(</span>initial_probs, transition_matrix, emission_means, emission_covs<span style="font-weight: bold">)</span>

    # Properties to get various parameters of the model
    @property
    def emission_distribution<span style="font-weight: bold">(</span>self<span style="font-weight: bold">)</span>:
        return self._emission_distribution

    @property
    def emission_means<span style="font-weight: bold">(</span>self<span style="font-weight: bold">)</span>:
        return self.emission_distribution.mean<span style="font-weight: bold">()</span>

    @property
    def emission_covariance_matrices<span style="font-weight: bold">(</span>self<span style="font-weight: bold">)</span>:
        return self.emission_distribution.covariance<span style="font-weight: bold">()</span>

    @property
    def unconstrained_params<span style="font-weight: bold">(</span>self<span style="font-weight: bold">)</span>:
        <span style="color: #008000; text-decoration-color: #008000">""</span>"Helper property to get a PyTree of unconstrained parameters.
        <span style="color: #008000; text-decoration-color: #008000">""</span>"
        return tfb.SoftmaxCentered<span style="font-weight: bold">()</span>.inverse<span style="font-weight: bold">(</span>self.initial_probabilities<span style="font-weight: bold">)</span>, \
               tfb.SoftmaxCentered<span style="font-weight: bold">()</span>.inverse<span style="font-weight: bold">(</span>self.transition_matrix<span style="font-weight: bold">)</span>, \
               self.emission_means, \
               PSDToRealBijector.forward<span style="font-weight: bold">(</span>self.emission_covariance_matrices<span style="font-weight: bold">)</span>

    @classmethod
    def from_unconstrained_params<span style="font-weight: bold">(</span>cls, unconstrained_params, hypers<span style="font-weight: bold">)</span>:
        initial_probabilities = tfb.SoftmaxCentered<span style="font-weight: bold">()</span>.forward<span style="font-weight: bold">(</span>unconstrained_params<span style="font-weight: bold">[</span><span style="color: #000080; text-decoration-color: #000080; font-weight: bold">0</span><span style="font-weight: bold">])</span>
        transition_matrix = tfb.SoftmaxCentered<span style="font-weight: bold">()</span>.forward<span style="font-weight: bold">(</span>unconstrained_params<span style="font-weight: bold">[</span><span style="color: #000080; text-decoration-color: #000080; font-weight: bold">1</span><span style="font-weight: bold">])</span>
        emission_means = unconstrained_params<span style="font-weight: bold">[</span><span style="color: #000080; text-decoration-color: #000080; font-weight: bold">2</span><span style="font-weight: bold">]</span>
        emission_covs = PSDToRealBijector.inverse<span style="font-weight: bold">(</span>unconstrained_params<span style="font-weight: bold">[</span><span style="color: #000080; text-decoration-color: #000080; font-weight: bold">3</span><span style="font-weight: bold">])</span>
        return cls<span style="font-weight: bold">(</span>initial_probabilities, transition_matrix, emission_means, emission_covs, 
*hypers<span style="font-weight: bold">)</span>

</pre>
</div></div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Set dimensions</span>
<span class="n">num_states</span> <span class="o">=</span> <span class="mi">5</span>
<span class="n">emission_dim</span> <span class="o">=</span> <span class="mi">2</span>

<span class="c1"># Specify parameters of the HMM</span>
<span class="n">initial_probs</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="n">num_states</span><span class="p">)</span> <span class="o">/</span> <span class="n">num_states</span>
<span class="n">transition_matrix</span> <span class="o">=</span> <span class="mf">0.95</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">num_states</span><span class="p">)</span> <span class="o">+</span> <span class="mf">0.05</span> <span class="o">*</span> <span class="n">jnp</span><span class="o">.</span><span class="n">roll</span><span class="p">(</span><span class="n">jnp</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">num_states</span><span class="p">),</span> <span class="mi">1</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="n">emission_means</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">column_stack</span><span class="p">([</span>
    <span class="n">jnp</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">jnp</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span> <span class="o">*</span> <span class="n">jnp</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">num_states</span><span class="o">+</span><span class="mi">1</span><span class="p">))[:</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span>
    <span class="n">jnp</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">jnp</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span> <span class="o">*</span> <span class="n">jnp</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">num_states</span><span class="o">+</span><span class="mi">1</span><span class="p">))[:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
<span class="p">])</span>
<span class="n">emission_covs</span> <span class="o">=</span> <span class="n">jnp</span><span class="o">.</span><span class="n">tile</span><span class="p">(</span><span class="mf">0.1</span><span class="o">**</span><span class="mi">2</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">emission_dim</span><span class="p">),</span> <span class="p">(</span><span class="n">num_states</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">))</span>

<span class="n">hmm</span> <span class="o">=</span> <span class="n">GaussianHMM</span><span class="p">(</span><span class="n">initial_probs</span><span class="p">,</span>
                       <span class="n">transition_matrix</span><span class="p">,</span>
                       <span class="n">emission_means</span><span class="p">,</span> 
                       <span class="n">emission_covs</span><span class="p">)</span>

<span class="n">print_source</span><span class="p">(</span><span class="n">hmm</span><span class="o">.</span><span class="n">sample</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<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.)
</pre></div>
</div>
<div class="output text_html"><pre style="white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace">    def sample<span style="font-weight: bold">(</span>self, key, num_timesteps<span style="font-weight: bold">)</span>:
        <span style="color: #008000; text-decoration-color: #008000">""</span>"Sample a sequence of latent states and emissions.

        Args:
            key <span style="font-weight: bold">(</span>_type_<span style="font-weight: bold">)</span>: _description_
            num_timesteps <span style="font-weight: bold">(</span>_type_<span style="font-weight: bold">)</span>: _description_
        <span style="color: #008000; text-decoration-color: #008000">""</span>"
        def _step<span style="font-weight: bold">(</span>state, key<span style="font-weight: bold">)</span>:
            key1, key2 = jr.split<span style="font-weight: bold">(</span>key, <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">2</span><span style="font-weight: bold">)</span>
            emission = self.emission_distribution.sample<span style="font-weight: bold">(</span><span style="color: #808000; text-decoration-color: #808000">seed</span>=<span style="color: #800080; text-decoration-color: #800080">key1</span><span style="font-weight: bold">)</span>
            next_state = self.transition_distribution.sample<span style="font-weight: bold">(</span><span style="color: #808000; text-decoration-color: #808000">seed</span>=<span style="color: #800080; text-decoration-color: #800080">key2</span><span style="font-weight: bold">)</span>
            return next_state, <span style="font-weight: bold">(</span>state, emission<span style="font-weight: bold">)</span>

        # Sample the initial state
        key1, key = jr.split<span style="font-weight: bold">(</span>key, <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">2</span><span style="font-weight: bold">)</span>
        initial_state = self.initial_distribution.sample<span style="font-weight: bold">(</span><span style="color: #808000; text-decoration-color: #808000">seed</span>=<span style="color: #800080; text-decoration-color: #800080">key1</span><span style="font-weight: bold">)</span>

        # Sample the remaining emissions and states
        keys = jr.split<span style="font-weight: bold">(</span>key, num_timesteps<span style="font-weight: bold">)</span>
        _, <span style="font-weight: bold">(</span>states, emissions<span style="font-weight: bold">)</span> = lax.scan<span style="font-weight: bold">(</span>_step, initial_state, keys<span style="font-weight: bold">)</span>
        return states, emissions

</pre>
</div></div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">distrax</span>
<span class="kn">from</span> <span class="nn">distrax</span> <span class="kn">import</span> <span class="n">HMM</span>

<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
    <span class="p">[</span><span class="mf">0.95</span><span class="p">,</span> <span class="mf">0.05</span><span class="p">],</span>
    <span class="p">[</span><span class="mf">0.10</span><span class="p">,</span> <span class="mf">0.90</span><span class="p">]</span>
<span class="p">])</span>

<span class="c1"># observation matrix</span>
<span class="n">B</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
    <span class="p">[</span><span class="mi">1</span><span class="o">/</span><span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">6</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">6</span><span class="p">],</span> <span class="c1"># fair die</span>
    <span class="p">[</span><span class="mi">1</span><span class="o">/</span><span class="mi">10</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">10</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">10</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">10</span><span class="p">,</span> <span class="mi">1</span><span class="o">/</span><span class="mi">10</span><span class="p">,</span> <span class="mi">5</span><span class="o">/</span><span class="mi">10</span><span class="p">]</span> <span class="c1"># loaded die</span>
<span class="p">])</span>

<span class="n">pi</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">])</span>

<span class="p">(</span><span class="n">nstates</span><span class="p">,</span> <span class="n">nobs</span><span class="p">)</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">shape</span><span class="p">(</span><span class="n">B</span><span class="p">)</span>

<span class="n">hmm</span> <span class="o">=</span> <span class="n">HMM</span><span class="p">(</span><span class="n">trans_dist</span><span class="o">=</span><span class="n">distrax</span><span class="o">.</span><span class="n">Categorical</span><span class="p">(</span><span class="n">probs</span><span class="o">=</span><span class="n">A</span><span class="p">),</span>
            <span class="n">init_dist</span><span class="o">=</span><span class="n">distrax</span><span class="o">.</span><span class="n">Categorical</span><span class="p">(</span><span class="n">probs</span><span class="o">=</span><span class="n">pi</span><span class="p">),</span>
            <span class="n">obs_dist</span><span class="o">=</span><span class="n">distrax</span><span class="o">.</span><span class="n">Categorical</span><span class="p">(</span><span class="n">probs</span><span class="o">=</span><span class="n">B</span><span class="p">))</span>

<span class="nb">print</span><span class="p">(</span><span class="n">hmm</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>&lt;distrax._src.utils.hmm.HMM object at 0x7fde82c856d0&gt;
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="n">print_source</span><span class="p">(</span><span class="n">hmm</span><span class="o">.</span><span class="n">sample</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_html"><pre style="white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace">  def sample<span style="font-weight: bold">(</span>self,
             *,
             seed: chex.PRNGKey,
             seq_len: chex.Array<span style="font-weight: bold">)</span> -&gt; Tuple:
    <span style="color: #008000; text-decoration-color: #008000">""</span>"Sample from this HMM.

    Samples an observation of given length according to this
    Hidden Markov Model and gives the sequence of the hidden states
    as well as the observation.

    Args:
      seed: Random key of shape <span style="font-weight: bold">(</span><span style="color: #000080; text-decoration-color: #000080; font-weight: bold">2</span>,<span style="font-weight: bold">)</span> and dtype uint32.
      seq_len: The length of the observation sequence.

    Returns:
      Tuple of hidden state sequence, and observation sequence.
    <span style="color: #008000; text-decoration-color: #008000">""</span>"
    rng_key, rng_init = jax.random.split<span style="font-weight: bold">(</span>seed<span style="font-weight: bold">)</span>
    initial_state = self._init_dist.sample<span style="font-weight: bold">(</span><span style="color: #808000; text-decoration-color: #808000">seed</span>=<span style="color: #800080; text-decoration-color: #800080">rng_init</span><span style="font-weight: bold">)</span>

    def draw_state<span style="font-weight: bold">(</span>prev_state, key<span style="font-weight: bold">)</span>:
      state = self._trans_dist.sample<span style="font-weight: bold">(</span><span style="color: #808000; text-decoration-color: #808000">seed</span>=<span style="color: #800080; text-decoration-color: #800080">key</span><span style="font-weight: bold">)</span>
      return state, state

    rng_state, rng_obs = jax.random.split<span style="font-weight: bold">(</span>rng_key<span style="font-weight: bold">)</span>
    keys = jax.random.split<span style="font-weight: bold">(</span>rng_state, seq_len - <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">1</span><span style="font-weight: bold">)</span>
    _, states = jax.lax.scan<span style="font-weight: bold">(</span>draw_state, initial_state, keys<span style="font-weight: bold">)</span>
    states = jnp.append<span style="font-weight: bold">(</span>initial_state, states<span style="font-weight: bold">)</span>

    def draw_obs<span style="font-weight: bold">(</span>state, key<span style="font-weight: bold">)</span>:
      return self._obs_dist.sample<span style="font-weight: bold">(</span><span style="color: #808000; text-decoration-color: #808000">seed</span>=<span style="color: #800080; text-decoration-color: #800080">key</span><span style="font-weight: bold">)</span>

    keys = jax.random.split<span style="font-weight: bold">(</span>rng_obs, seq_len<span style="font-weight: bold">)</span>
    obs_seq = jax.vmap<span style="font-weight: bold">(</span>draw_obs, <span style="color: #808000; text-decoration-color: #808000">in_axes</span>=<span style="font-weight: bold">(</span><span style="color: #000080; text-decoration-color: #000080; font-weight: bold">0</span>, <span style="color: #000080; text-decoration-color: #000080; font-weight: bold">0</span><span style="font-weight: bold">))(</span>states, keys<span style="font-weight: bold">)</span>

    return states, obs_seq

</pre>
</div></div>
</div>

    <script type="text/x-thebe-config">
    {
        requestKernel: true,
        binderOptions: {
            repo: "binder-examples/jupyter-stacks-datascience",
            ref: "master",
        },
        codeMirrorConfig: {
            theme: "abcdef",
            mode: "python"
        },
        kernelOptions: {
            kernelName: "python3",
            path: "./chapters"
        },
        predefinedOutput: true
    }
    </script>
    <script>kernelName = 'python3'</script>

              </div>
              
            
                <!-- Previous / next buttons -->
<div class='prev-next-area'>
</div>
            
        </div>
    </div>
    <footer class="footer">
  <p>
    
      By Kevin Murphy, Scott Linderman, et al.<br/>
    
        &copy; Copyright 2021.<br/>
  </p>
</footer>
</main>


      </div>
    </div>
  
  <script src="../_static/js/index.be7d3bbb2ef33a8344ce.js"></script>

  </body>
</html>