chap17soln


Modeling and Simulation in Python

Chapter 17

Copyright 2017 Allen Downey

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


Out[2]:
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

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 [ ]: