Demo of a KDE plot beside timeseries set (Interactive)
======================================================

.. code:: ipython3

    %pylab inline
    #%matplotlib inline
    #import matplotlib.pyplot as plt
    import mpld3
    from mpld3 import plugins
    #mpld3.enable_notebook()
    #%matplotlib notebook
    #
    import pysd
    import numpy as np
    import pandas as pd
    import seaborn


.. parsed-literal::

    Populating the interactive namespace from numpy and matplotlib


Load the model using PySD
~~~~~~~~~~~~~~~~~~~~~~~~~

The model is a basic, 1-stock carbon bathtub model

.. code:: ipython3

    model = pysd.read_vensim('../../models/Climate/Atmospheric_Bathtub.mdl')
    print(model.doc)


.. parsed-literal::

                       Real Name                    Py Name Subscripts  Units  \
    0                  Emissions                  emissions       None   None   
    1  Excess Atmospheric Carbon  excess_atmospheric_carbon       None   None   
    2                 FINAL TIME                 final_time       None  Month   
    3               INITIAL TIME               initial_time       None  Month   
    4            Natural Removal            natural_removal       None   None   
    5           Removal Constant           removal_constant       None   None   
    6                    SAVEPER                    saveper       None  Month   
    7                  TIME STEP                  time_step       None  Month   
    8                       Time                       time       None   None   
    
           Limits       Type Subtype                                     Comment  
    0  (nan, nan)   Constant  Normal                                        None  
    1  (nan, nan)   Stateful   Integ                                        None  
    2  (nan, nan)   Constant  Normal          The final time for the simulation.  
    3  (nan, nan)   Constant  Normal        The initial time for the simulation.  
    4  (nan, nan)  Auxiliary  Normal                                        None  
    5  (nan, nan)   Constant  Normal                                        None  
    6  (0.0, nan)  Auxiliary  Normal  The frequency with which output is stored.  
    7  (0.0, nan)   Constant  Normal           The time step for the simulation.  
    8  (nan, nan)       None    None                  Current time of the model.  


Generate a set of parameters to use as input
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Here, drawing 100 constant values for the ``Emissions`` parameter from
an exponential distribution

.. code:: ipython3

    n_runs = 100
    runs = pd.DataFrame({'Emissions': np.random.exponential(scale=10000, size=n_runs)})
    runs.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Emissions</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>13139.288132</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10414.796999</td>
        </tr>
        <tr>
          <th>2</th>
          <td>8505.605002</td>
        </tr>
        <tr>
          <th>3</th>
          <td>55251.743520</td>
        </tr>
        <tr>
          <th>4</th>
          <td>3650.497881</td>
        </tr>
      </tbody>
    </table>
    </div>



Run the model with the various parameters
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    result = runs.apply(lambda p: model.run(params=dict(p))['Excess Atmospheric Carbon'],
                        axis=1).T
    result.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>0</th>
          <th>1</th>
          <th>2</th>
          <th>3</th>
          <th>4</th>
          <th>5</th>
          <th>6</th>
          <th>7</th>
          <th>8</th>
          <th>9</th>
          <th>...</th>
          <th>90</th>
          <th>91</th>
          <th>92</th>
          <th>93</th>
          <th>94</th>
          <th>95</th>
          <th>96</th>
          <th>97</th>
          <th>98</th>
          <th>99</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>...</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
        </tr>
        <tr>
          <th>1</th>
          <td>13139.288132</td>
          <td>10414.796999</td>
          <td>8505.605002</td>
          <td>55251.743520</td>
          <td>3650.497881</td>
          <td>7490.206883</td>
          <td>5908.429778</td>
          <td>4526.126740</td>
          <td>3969.536964</td>
          <td>42777.645014</td>
          <td>...</td>
          <td>9654.733184</td>
          <td>16864.050543</td>
          <td>15725.059944</td>
          <td>7426.330819</td>
          <td>10563.284279</td>
          <td>5586.831254</td>
          <td>24433.431913</td>
          <td>13290.730858</td>
          <td>6245.476841</td>
          <td>11057.440643</td>
        </tr>
        <tr>
          <th>2</th>
          <td>26147.183383</td>
          <td>20725.446028</td>
          <td>16926.153954</td>
          <td>109950.969604</td>
          <td>7264.490783</td>
          <td>14905.511697</td>
          <td>11757.775259</td>
          <td>9006.992213</td>
          <td>7899.378559</td>
          <td>85127.513578</td>
          <td>...</td>
          <td>19212.919037</td>
          <td>33559.460581</td>
          <td>31292.869289</td>
          <td>14778.398329</td>
          <td>21020.935715</td>
          <td>11117.794196</td>
          <td>48622.529506</td>
          <td>26448.554408</td>
          <td>12428.498913</td>
          <td>22004.306879</td>
        </tr>
        <tr>
          <th>3</th>
          <td>39024.999681</td>
          <td>30932.988567</td>
          <td>25262.497417</td>
          <td>164103.203428</td>
          <td>10842.343756</td>
          <td>22246.663463</td>
          <td>17548.627285</td>
          <td>13443.049031</td>
          <td>11789.921738</td>
          <td>127053.883456</td>
          <td>...</td>
          <td>28675.523031</td>
          <td>50087.916518</td>
          <td>46705.000540</td>
          <td>22056.945165</td>
          <td>31374.010637</td>
          <td>16593.447508</td>
          <td>72569.736124</td>
          <td>39474.799721</td>
          <td>18549.690765</td>
          <td>32841.704452</td>
        </tr>
        <tr>
          <th>4</th>
          <td>51774.037816</td>
          <td>41038.455680</td>
          <td>33515.477445</td>
          <td>217713.914913</td>
          <td>14384.418199</td>
          <td>29514.403711</td>
          <td>23281.570790</td>
          <td>17834.745281</td>
          <td>15641.559485</td>
          <td>168560.989636</td>
          <td>...</td>
          <td>38043.500985</td>
          <td>66451.087897</td>
          <td>61963.010478</td>
          <td>29262.706532</td>
          <td>41623.554810</td>
          <td>22014.344287</td>
          <td>96277.470675</td>
          <td>52370.782582</td>
          <td>24609.670698</td>
          <td>43570.728050</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 100 columns</p>
    </div>



Draw a static plot showing the results, and a marginal density plot
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This would be for making graphics for a publication, where you don’t
want an interactive view, but you want fine control over what the figure
looks like.

This is relatively simple, because we rely on the plotting library
``seaborn`` to generate the KDE plot, instead of working out the
densities ourselves.

.. code:: ipython3

    import matplotlib.pylab as plt
    # define when to show the density
    density_time = 85
    
    # left side: plot all traces, slightly transparent
    plt.subplot2grid((1,4), loc=(0,0), colspan=3)
    [plt.plot(result.index, result[i], 'b', alpha=.02) for i in result.columns]
    ymax = result.max().max()
    plt.ylim(0, ymax)
    
    # left side: add marker of density location
    plt.vlines(density_time, 0, ymax, 'k')
    plt.text(density_time, ymax, 'Density Calculation', ha='right', va='top', rotation=90)
    
    # right side: gaussian KDE on selected timestamp
    plt.subplot2grid((1,4), loc=(0,3))
    seaborn.kdeplot(y=result.loc[density_time])
    plt.ylim(0, ymax)
    plt.yticks([])
    plt.xticks([])
    plt.xlabel('Density')
    
    plt.suptitle('Emissions scenarios under uncertainty', fontsize=16);



.. image:: marginal_density_plot_interactive_files/marginal_density_plot_interactive_9_0.png


Interactive plot using python backend
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The following would be for lightweight exploration

.. code:: ipython3

    import matplotlib as mpl
    from ipywidgets import interact, IntSlider
    sim_time = 200
    slider_time = IntSlider(description = 'Select Time for plotting Density',
                            min=0, max=result.index[-1], value=1)

.. code:: ipython3

    @interact(density_time=slider_time)
    def update(density_time): 
        ax1 = plt.subplot2grid((1,4), loc=(0,0), colspan=3)
        [ax1.plot(result.index, result[i], 'b', alpha=.02) for i in result.columns]
        ymax = result.max().max()
        ax1.set_ylim(0, ymax)
    
        # left side: add marker of density location
        ax1.vlines(density_time, 0, ymax, 'k')
        ax1.text(density_time, ymax, 'Density Calculation', ha='right', va='top', rotation=90)
    
        # right side: gaussian KDE on selected timestamp
        ax2 = plt.subplot2grid((1,4), loc=(0,3))
        seaborn.kdeplot(y=result.loc[density_time], ax=ax2)
        ax2.set_ylim(0, ymax)
        ax2.set_yticks([])
        ax2.set_xticks([])
        ax2.set_xlabel('Density')
    
        plt.suptitle('Emissions scenarios under uncertainty', fontsize=16);
        



.. parsed-literal::

    interactive(children=(IntSlider(value=1, description='Select Time for plotting Density'), Output()), _dom_clas…


Interactive figure with javascript background
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This script would prepare interactive graphics to share on a webpage,
independent of the python backend.

.. code:: ipython3

    fig, ax = plt.subplots(3, 3, figsize=(6, 6))
    fig.subplots_adjust(hspace=0.1, wspace=0.1)
    ax = ax[::-1]
    
    X = np.random.normal(size=(3, 100))
    for i in range(3):
        for j in range(3):
            ax[i, j].xaxis.set_major_formatter(plt.NullFormatter())
            ax[i, j].yaxis.set_major_formatter(plt.NullFormatter())
            points = ax[i, j].scatter(X[j], X[i])
            
    plugins.connect(fig, plugins.LinkedBrush(points))



.. image:: marginal_density_plot_interactive_files/marginal_density_plot_interactive_14_0.png