Modeling and Simulation in Python

Chapter 17

License: Creative Commons Attribution 4.0 International



In [1]:

    
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

Data

We have data from Pacini and Bergman (1986), "MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test", Computer Methods and Programs in Biomedicine, 23: 113-122..



In [2]:

    
data = pd.read_csv('data/glucose_insulin.csv', index_col='time')

Here's what the glucose time series looks like.



In [3]:

    
plot(data.glucose, 'bo', label='glucose')
decorate(xlabel='Time (min)',
         ylabel='Concentration (mg/dL)')

And the insulin time series.



In [4]:

    
plot(data.insulin, 'go', label='insulin')
decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')

For the book, I put them in a single figure, using subplot



In [5]:

    
subplot(2, 1, 1)
plot(data.glucose, 'bo', label='glucose')
decorate(ylabel='Concentration (mg/dL)')

subplot(2, 1, 2)
plot(data.insulin, 'go', label='insulin')
decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')

savefig('figs/chap17-fig01.pdf')









    



Saving figure to file figs/chap17-fig01.pdf

Interpolation

We have measurements of insulin concentration at discrete points in time, but we need to estimate it at intervening points. We'll use interpolate, which takes a Series and returns a function:

The return value from interpolate is a function.



In [6]:

    
I = interpolate(data.insulin)









    Out[6]:





<function modsim.modsim.interpolate.<locals>.wrapper(x)>

We can use the result, I, to estimate the insulin level at any point in time.



In [7]:

    
I(7)









    Out[7]:





68.0

I can also take an array of time and return an array of estimates:



In [8]:

    
t_0 = get_first_label(data)
t_end = get_last_label(data)
ts = linrange(t_0, t_end, endpoint=True)
I(ts)
type(ts)









    Out[8]:





numpy.ndarray

Here's what the interpolated values look like.



In [9]:

    
plot(data.insulin, 'go', label='insulin data')
plot(ts, I(ts), color='green', label='interpolated')

decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')

savefig('figs/chap17-fig02.pdf')









    



Saving figure to file figs/chap17-fig02.pdf

Exercise: Read the documentation of scipy.interpolate.interp1d. Pass a keyword argument to interpolate to specify one of the other kinds of interpolation, and run the code again to see what it looks like.



In [10]:

    
# Solution

I = interpolate(data.insulin, kind='cubic')
plot(data.insulin, 'go', label='insulin data')
plot(ts, I(ts), color='green', label='interpolated')

decorate(xlabel='Time (min)',
         ylabel='Concentration ($\mu$U/mL)')

Exercise: Interpolate the glucose data and generate a plot, similar to the previous one, that shows the data points and the interpolated curve evaluated at the time values in ts.



In [11]:

    
# Solution

G = interpolate(data.glucose)
plot(data.glucose, 'bo', label='gluciose data')
plot(ts, G(ts), color='blue', label='interpolated')

decorate(xlabel='Time (min)',
         ylabel='Concentration (mg/dL)')

Under the hood



In [12]:

    
source_code(interpolate)









    



def interpolate(series, **options):
    """Creates an interpolation function.

    series: Series object
    options: any legal options to scipy.interpolate.interp1d

    returns: function that maps from the index of the series to values
    """
    if has_nan(series.index):
        msg = """The Series you passed to interpolate contains
                 NaN values in the index, which would result in
                 undefined behavior.  So I'm putting a stop to that."""
        raise ValueError(msg)

    if not is_strictly_increasing(series.index):
        msg = """The Series you passed to interpolate has an index
                 that is not strictly increasing, which would result in
                 undefined behavior.  So I'm putting a stop to that."""
        raise ValueError(msg)

    # make the interpolate function extrapolate past the ends of
    # the range, unless `options` already specifies a value for `fill_value`
    underride(options, fill_value='extrapolate')

    # call interp1d, which returns a new function object
    x = magnitudes(series.index)
    y = magnitudes(series.values)
    interp_func = interp1d(x, y, **options)
    units = get_units(series.values[0])

    def wrapper(x):
        return interp_func(magnitudes(x)) * units
    return wrapper



In [ ]:

	glucose	insulin
time
0	92	11
2	350	26
4	287	130
6	251	85
8	240	51
10	216	49
12	211	45
14	205	41
16	196	35
19	192	30
22	172	30
27	163	27
32	142	30
42	124	22
52	105	15
62	92	15
72	84	11
82	77	10
92	82	8
102	81	11
122	82	7
142	82	8
162	85	8
182	90	7