Demo of a KDE plot beside timeseries set
========================================

.. code:: ipython3

    %pylab inline
    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>207.734860</td>
        </tr>
        <tr>
          <th>1</th>
          <td>7109.125870</td>
        </tr>
        <tr>
          <th>2</th>
          <td>4768.041210</td>
        </tr>
        <tr>
          <th>3</th>
          <td>2191.302602</td>
        </tr>
        <tr>
          <th>4</th>
          <td>16930.863016</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>207.734860</td>
          <td>7109.125870</td>
          <td>4768.041210</td>
          <td>2191.302602</td>
          <td>16930.863016</td>
          <td>14676.223891</td>
          <td>2058.091604</td>
          <td>9229.459073</td>
          <td>1215.333008</td>
          <td>22711.067942</td>
          <td>...</td>
          <td>5438.319103</td>
          <td>13029.384661</td>
          <td>28268.791836</td>
          <td>7543.867757</td>
          <td>12566.422839</td>
          <td>1565.757571</td>
          <td>1029.656855</td>
          <td>2802.757320</td>
          <td>18335.083353</td>
          <td>3311.616813</td>
        </tr>
        <tr>
          <th>2</th>
          <td>413.392371</td>
          <td>14147.160482</td>
          <td>9488.402008</td>
          <td>4360.692179</td>
          <td>33692.417401</td>
          <td>29205.685543</td>
          <td>4095.602292</td>
          <td>18366.623556</td>
          <td>2418.512686</td>
          <td>45195.025205</td>
          <td>...</td>
          <td>10822.255014</td>
          <td>25928.475475</td>
          <td>56254.895754</td>
          <td>15012.296837</td>
          <td>25007.181449</td>
          <td>3115.857566</td>
          <td>2049.017141</td>
          <td>5577.487067</td>
          <td>36486.815873</td>
          <td>6590.117458</td>
        </tr>
        <tr>
          <th>3</th>
          <td>616.993307</td>
          <td>21114.814748</td>
          <td>14161.559198</td>
          <td>6508.387860</td>
          <td>50286.356243</td>
          <td>43589.852579</td>
          <td>6112.737873</td>
          <td>27412.416393</td>
          <td>3609.660567</td>
          <td>67454.142895</td>
          <td>...</td>
          <td>16152.351567</td>
          <td>38698.575381</td>
          <td>83961.138633</td>
          <td>22406.041626</td>
          <td>37323.532473</td>
          <td>4650.456561</td>
          <td>3058.183825</td>
          <td>8324.469516</td>
          <td>54457.031068</td>
          <td>9835.833097</td>
        </tr>
        <tr>
          <th>4</th>
          <td>818.558233</td>
          <td>28012.792471</td>
          <td>18787.984817</td>
          <td>8634.606583</td>
          <td>66714.355696</td>
          <td>57830.177944</td>
          <td>8109.702098</td>
          <td>36367.751302</td>
          <td>4788.896970</td>
          <td>89490.669408</td>
          <td>...</td>
          <td>21429.147154</td>
          <td>51340.974288</td>
          <td>111390.319083</td>
          <td>29725.848967</td>
          <td>49516.719987</td>
          <td>6169.709566</td>
          <td>4057.258842</td>
          <td>11043.982140</td>
          <td>72247.544111</td>
          <td>13049.091579</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 100 columns</p>
    </div>



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

.. code:: ipython3

    # left side should have all traces plotted
    plt.subplot2grid((1,4), loc=(0,0), colspan=3)
    [plt.plot(result.index, result[i], 'b', alpha=.02) for i in result.columns]
    plt.ylim(0, max(result.iloc[-1]))
    
    # right side has gaussian KDE on last timestamp
    plt.subplot2grid((1,4), loc=(0,3))
    seaborn.kdeplot(y=result.iloc[-1])
    plt.ylim(0, max(result.iloc[-1]));
    plt.yticks([])
    plt.xticks([])
    
    plt.suptitle('Emissions scenarios under uncertainty', fontsize=16);



.. image:: marginal_density_plot_files/marginal_density_plot_9_0.png