Calculating doubling times from data points

Question

In the code below, noisy data points with unique errors are created. From this, an exponential function is fitted to the data points, and then doubling times (10 unit windows) are calculated.

I'm uncertain how to show the unique errors in the data points in the fitted function or doubling times.

Output:

from scipy import optimize
from matplotlib import pylab as plt
import numpy as np
import pdb
from numpy import log

def exp_growth(t, x0, r):
    return x0 * ((1 + r) ** t)

def doubling_time(m, x_pts, y_pts):
    window = 10

    x1 = x_pts[m]
    y1 = y_pts[m]
    x2 = x_pts[m+window]
    y2 = y_pts[m+window]

    return (x2 - x1) * log(2) / log(y2 / y1)

# First, artificially create data points to work with
data_points = 42

# Create the x-axis
x_pts = range(0, data_points)

# Create noisy points with: y = x^2 + noise, with unique possible errors
y_pts = []
y_err = []
for i in range(data_points):
    random_scale = np.random.random()
    y_pts.append((i * i) + data_points * random_scale)
    y_err.append(random_scale * 100 + 100)

x_pts = np.array(x_pts)
y_pts = np.array(y_pts)
y_err = np.array(y_err)

# Fit to function
[x0, r], pcov  = optimize.curve_fit(exp_growth, x_pts, y_pts, p0=(0.001, 1.0))
fitted_data = exp_growth(x_pts, x0, r)

# Find doubling times
x_t2 = range(32)
t2 = []
t2_fit = []
for i in range(32):
    t2.append(doubling_time(i, x_pts, y_pts))
    t2_fit.append(doubling_time(i, x_pts, fitted_data))

# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)

ax1.plot(x_pts, y_pts, 'bo')
ax1.errorbar(x_pts, y_pts, yerr=y_err)
ax1.set_ylim([0, 2000])
ax1.set_title('Artificially created raw data points with unique errors', fontsize=8)

ax2.plot(fitted_data, 'g-')
ax2.set_ylim([0, 2000])
ax2.set_title('Fitted exponential function', fontsize=8)

ax3.plot(x_t2, t2, 'ro', label='From points')
ax3.plot(x_t2, t2_fit, 'bo', label='From fitted')
ax3.set_title('Doubling time at each point (10 unit window)', fontsize=8)
ax3.legend(fontsize='8')

plt.show()

Answer 1

Your code is doing the job, but I think there are a few ways in which it could be improved.

1. Consistent use of parameters for functions.

---

def doubling_time(m, x_pts, y_pts):
    window = 10

    x1 = x_pts[m]
    y1 = y_pts[m]
    x2 = x_pts[m+window]
    y2 = y_pts[m+window]

    return (x2 - x1) * log(2) / log(y2 / y1)

Why is window not a parameter in the function? Having it as a parameter will make it much easier to change, and more importantly, it will make future users of the code (including yourself six months/weeks/days from now) realize that in fact this function depends on a window parameter to compute its result. If you do def find_local_doubling_time(index, x_pts, y_pts, window=10):[...], then window will be a parameter with a default value of 10, so you don't have to pass it in if you don't want to.

2. I would personally find a name like index or something to be far more illustrative and informative than m. Also, the function itself could do with a slightly more informative name, such as local_doubling_time() or find_local_doubling_time().

3. The way you are generating y_pts feels very unnatural to me. Since you have x_pts already, you can just use numpy to define a y_pts without any for loops, like this:

---

x_pts = np.array(range(0, data_points))
random_scales = np.random.random(size=data_points)
y_pts = x_pts**2 + data_points*random_scales
y_err = random_scale * 100 + 100

BTW, the more common choice for the error model would be Gaussian noise instead of uniformly distributed noise. You might consider adding a comment to the code to explain your choice of the uniform model.

4. Finding the "local" or "instananeous" doubling times does not require for loops either.

---

t2 = [doubling_time(i, x_pts, y_pts) for i in x_t2]
t2_fit = [doubling_time(i, x_pts, fitted_data) for i in x_t2]

5. I didn't know about the sharex option for plt.subplots(). Very cool to learn! I like your graphs; the only improvement would be to plot the points and the curve on the same graph (i.e. combine the top panel and mid panel of the graph).

6. Scipy's optimize uses nonlinear least squares regression. In addition to comparing to the "local" results, you might also compare the NLSR results to the results of doing linear-regression on log-transformed data.

Answer 2

Especially with mathematical-like programs, you want to write docstrings with your functions and classes.

In case you don't know what a docstring is, a docstring is a short piece of documentation that describes the parameters and return value of a function. And, it can include any sort of information/summary that is needed to fully explain a function.

Here is an example:

def add(a, b):
    """
    Adds two numbers together and returns the sum
    a + b
    @param(number) -- one number to add
    @param(number) -- another number to add
    @return(number) -- the sum of both parameters

    return a + b

In the case of your code where you have a few mathematical functions, providing a docstring is greatly important because there you can write out the formula that the function is using, what each part of the formula means, and what the return of the formula does for the rest of your code.

And, in places that aren't functions but need a docstring, you can just a # comments on the side (you don't have enough of these).

---

All of your main code (aside from functions) should be put into:

if __name__ == "__main__":
    [code or main()]

See here for why.

---

Your code seems very procedural. You have a few functions here and there but other than that, everything is just one after the other.

I'm having difficulty following your code so I can't give too many specific recommendations, but you should move more code to separate functions. For example, however, you could try to make a function call create_graph that will create a graph for you in just a single call.

That way, when you put your code into if __name__ == "__main__", every single line of code that you wrote won't be in there.