Surrogating a function with a machine learning estimator
========================================================

System dynamics generally represents the relationships between model
elements as either analytical equations, or lookup tables. However, in
some situations we may be presented with relationships that are not well
estimated by equations, but involve more than a single input leading to
a single output. When confrontied with this situation, other paradigms

.. image:: ../../../source/models/Manufacturing_Defects/Defects.png
   :width: 400 px

.. code:: ipython3

    %pylab inline
    import pysd
    import numpy as np
    import pandas as pd


.. parsed-literal::

    Populating the interactive namespace from numpy and matplotlib


.. code:: ipython3

    model = pysd.read_vensim('../../models/Manufacturing_Defects/Defects.mdl')

.. code:: ipython3

    data = pd.read_csv('../../data/Defects_Synthetic/Manufacturing_Defects_Synthetic_Data.csv')
    data.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>Unnamed: 0</th>
          <th>Workday</th>
          <th>Time per Task</th>
          <th>Defect Rate</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0</td>
          <td>0.303114</td>
          <td>0.060810</td>
          <td>0.023022</td>
        </tr>
        <tr>
          <th>1</th>
          <td>1</td>
          <td>0.263133</td>
          <td>0.052325</td>
          <td>0.023017</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2</td>
          <td>0.230397</td>
          <td>0.065387</td>
          <td>0.015868</td>
        </tr>
        <tr>
          <th>3</th>
          <td>3</td>
          <td>0.265632</td>
          <td>0.044866</td>
          <td>0.032806</td>
        </tr>
        <tr>
          <th>4</th>
          <td>4</td>
          <td>0.298651</td>
          <td>0.038648</td>
          <td>0.035234</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    plt.scatter(data['Workday'], data['Time per Task'], c=data['Defect Rate'], linewidth=0, alpha=.6)
    plt.ylabel('Time per Task')
    plt.xlabel('Length of Workday')
    plt.xlim(0.15, .45)
    plt.ylim(.01, .09)
    plt.box('off')
    plt.colorbar()
    plt.title('Defect Rate Measurements')
    plt.figtext(.88, .5, 'Defect Rate', rotation=90, verticalalignment='center');



.. image:: Surrogating_with_regression_files/Surrogating_with_regression_5_0.png


.. code:: ipython3

    from sklearn.svm import SVR
    
    Factors = data[['Workday','Time per Task']].values
    Outcome = data['Defect Rate'].values
    regression = SVR()
    regression.fit(Factors, Outcome)




.. parsed-literal::

    SVR()



.. code:: ipython3

    def new_defect_function():
        """ Replaces the original defects equation with a regression model"""
        workday = model.components.length_of_workday()
        
        time_per_task = model.components.time_allocated_per_unit()
        return regression.predict([[workday, time_per_task]])[0]
    
    model.components.defect_rate = new_defect_function

.. code:: ipython3

    model.components.defect_rate()




.. parsed-literal::

    0.0566754576475818



.. code:: ipython3

    model.run().plot();



.. image:: Surrogating_with_regression_files/Surrogating_with_regression_9_0.png