%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
#!/usr/bin/env python
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
# Model to fit the data with
def hill(x, baseline, amplitude, tconstant, hillcoef):
"""
Return the hill equation curve for the x values.
"""
return baseline+amplitude*(x**hillcoef)/(x**hillcoef+tconstant**hillcoef)
# Bootstrap resampling for fit function
class DataFit:
def __init__(self, name, function):
self.name = name
data = np.loadtxt(self.name, delimiter=',')
self.xdata = data[:, 0]
self.ydata = data[:, 1]
self.xfit = np.arange(self.xdata[0], self.xdata.max()+0.1, 0.1)
# set the fit function/parameters
self.set_fitfunc(function)
self.compute_yfit()
def set_fitfunc(self, namefunc):
import inspect
self.fitfunc = namefunc
f = inspect.getargspec(self.fitfunc)
n_params = len(f.args)-1 # minus the x argument
self.p0 = [1]*n_params
self.p0_args = f.args[1:]
def compute_yfit(self):
from scipy.optimize import curve_fit
self.popt, self.pcov = curve_fit(self.fitfunc, self.xdata,
self.ydata, p0=self.p0)
params = []
params.extend(self.popt)
params.insert(0, self.xfit)
yf = apply(self.fitfunc, params)
self.yfit = yf
def bootstrap_fitfunc(self, n, pos=1, alpha=0.05):
"""
Bootstrap n times the fit of the data and store the n fits in a matrix,
as well as the parameters of interest (position = pos) we want the
1-alpha confidence interval
"""
from scipy.optimize import curve_fit
# calculate residual from fit
params = list(self.popt)
params.insert(0, self.xdata)
ydatafit = apply(self.fitfunc, params)
r = self.ydata - ydatafit
# create matrix to store result of bootstrap
self.matrix_bootstrap = np.zeros((self.xfit.size, n))
boot_param = np.zeros(n)
for i in xrange(n):
new_popt = curve_fit(self.fitfunc, self.xdata,
ydatafit+r[sp.random.random_integers(0, len(r)-1,
len(r))], p0=self.popt)[0]
p_newfit = list(new_popt)
p_newfit.insert(0, self.xfit)
self.matrix_bootstrap[:, i] = apply(self.fitfunc, p_newfit)
boot_param[i] = new_popt[pos]
# rank rows of matrix and extract confidence band limits
ranked_matrix = self.matrix_bootstrap.copy()
ranked_matrix[:].sort()
self.ci_fitup = ranked_matrix[:, round(n*(1-alpha/2))-1]
self.ci_fitlow = ranked_matrix[:, round(n*(alpha/2))-1]
# rank parameter array and extract confidence interval
self.bootstrap_param = boot_param
boot_param.sort()
self.ci_param = list([boot_param[round(n*(alpha/2)-1)],
boot_param[n*(1-alpha/2)-1]])
# alpha to be fixed (% confidence interval) and number of resampling
alpha = 0.05
n = 10000
# Data file (2 columns, time and fluorescence value)
adult_hk2y_0glu = "all-data-adult/a-HK2-0glu/all_hk2yfp-0glu-adult.txt"
!cat "all-data-adult/a-HK2-0glu/all_hk2yfp-0glu-adult.txt"
0.000000000000000000e+00,1.000000000000000000e+02 5.000000000000000000e+00,9.487799673991361260e+01 1.000000000000000000e+01,8.529382442343876392e+01 1.500000000000000000e+01,5.301675032696944356e+01 2.000000000000000000e+01,5.474612750955677143e+01 2.500000000000000000e+01,5.494646319740675011e+01 3.000000000000000000e+01,5.581134093699490251e+01 3.500000000000000000e+01,5.301347056102324018e+01 4.000000000000000000e+01,5.248478635942989712e+01 0.000000000000000000e+00,1.000000000000000142e+02 5.000000000000000000e+00,9.766301445964390382e+01 1.000000000000000000e+01,8.335914624421469910e+01 1.500000000000000000e+01,6.994231683234843899e+01 2.000000000000000000e+01,6.632247437197489148e+01 2.500000000000000000e+01,6.785414306916727867e+01 3.000000000000000000e+01,6.452127332368395685e+01 3.500000000000000000e+01,6.702088088186738446e+01 4.000000000000000000e+01,6.327554485869725909e+01 0.000000000000000000e+00,1.000000000000000000e+02 5.000000000000000000e+00,9.381905735028574611e+01 1.000000000000000000e+01,7.672965305938846825e+01 1.500000000000000000e+01,6.353102124855913502e+01 2.000000000000000000e+01,6.463960801025592673e+01 2.500000000000000000e+01,5.752326109231942297e+01 3.000000000000000000e+01,5.584883754110204990e+01 3.500000000000000000e+01,5.420924022965913736e+01 4.000000000000000000e+01,5.092442135611415210e+01 0.000000000000000000e+00,1.000000000000000000e+02 5.000000000000000000e+00,1.000315123434486111e+02 1.000000000000000000e+01,8.759491156544339674e+01 1.500000000000000000e+01,7.022145854521663466e+01 2.000000000000000000e+01,6.690111141256915062e+01 2.500000000000000000e+01,6.691069117897242791e+01 3.000000000000000000e+01,5.895301768982740498e+01 3.500000000000000000e+01,6.006103221194981501e+01 4.000000000000000000e+01,5.839247381770547918e+01 0.000000000000000000e+00,1.000000000000000000e+02 5.000000000000000000e+00,8.853891736627421949e+01 1.000000000000000000e+01,7.096656069056011518e+01 1.500000000000000000e+01,4.899466124998252070e+01 2.000000000000000000e+01,4.912350342529322234e+01 2.500000000000000000e+01,4.626165214214815080e+01 3.000000000000000000e+01,4.587950925468438612e+01 3.500000000000000000e+01,4.432741178055862008e+01 4.000000000000000000e+01,4.092405335163461899e+01 0.000000000000000000e+00,1.000000000000000000e+02 1.000000000000000000e+00,1.066917762666061407e+02 5.000000000000000000e+00,9.397218797878801411e+01 1.000000000000000000e+01,7.455478915224961156e+01 1.500000000000000000e+01,6.440424014716889189e+01 2.000000000000000000e+01,5.395052979538768056e+01 2.500000000000000000e+01,5.406268855406460005e+01 3.500000000000000000e+01,6.284732061368104894e+01 4.500000000000000000e+01,5.602127623926315891e+01 0.000000000000000000e+00,1.000000000000000000e+02 1.000000000000000000e+00,1.026075315477147711e+02 5.000000000000000000e+00,7.827116956837613770e+01 1.000000000000000000e+01,7.854658848740508859e+01 1.500000000000000000e+01,5.604802668054081494e+01 2.000000000000000000e+01,6.276481603318313773e+01 2.500000000000000000e+01,5.640087882146369225e+01 3.500000000000000000e+01,5.537822214396339859e+01 4.500000000000000000e+01,6.802845539964606303e+01
# Fit of the data
d2y = DataFit(adult_hk2y_0glu, hill)
d2y.p0 = [100., -40., 10., 2.]
d2y.compute_yfit()
-c:6: RuntimeWarning: invalid value encountered in double_scalars
# Bootstrap of the residual of the fit
d2y.bootstrap_fitfunc(n, pos=2, alpha=alpha)
# tweak the plots
def format_axe(ax):
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.tick_params(axis='both', direction='out', width=3, length=7,
labelsize=16, pad=8)
ax.spines['left'].set_linewidth(3)
ax.spines['bottom'].set_linewidth(3)
fig = plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(221) # Data and original fit
ax2 = fig.add_subplot(222) # Individual bootstrapped fits (10000)
ax3 = fig.add_subplot(223) # Histogram of the time constants and 95% confidence intervals
ax4 = fig.add_subplot(224) # Original fit + gray shaded area for 95% confidence intervals of the fit
format_axe(ax1)
format_axe(ax2)
format_axe(ax3)
format_axe(ax4)
ax1.set_xlim(-1, 46)
ax1.set_ylim(20, 120)
ax2.set_xlim(-1, 46)
ax2.set_ylim(20, 120)
ax3.set_xlim(6, 14)
ax4.set_xlim(-1, 46)
ax4.set_ylim(20, 120)
ax1.set_xlabel('Time (min)', size=20)
ax2.set_xlabel('Time (min)', size=20)
ax3.set_xlabel('t$_{1/2}$ (min)', size=20)
ax4.set_xlabel('Time (min)', size=20)
ax1.set_ylabel('fluo (a.u.)', size=20)
ax2.set_ylabel('fluo (a.u.)', size=20)
ax3.set_ylabel('#', size=20)
ax4.set_ylabel('fluo (a.u.)', size=20)
ax1.plot(d2y.xdata, d2y.ydata, 'o', markersize=11, markerfacecolor='w',
markeredgecolor='k', markeredgewidth=2)
ax1.plot(d2y.xfit, d2y.yfit, '-r', lw=3, alpha=1)
ax2.plot(d2y.xfit, d2y.matrix_bootstrap, '-', color='0.5', lw=1)
ax2.plot(d2y.xfit, d2y.yfit, '-y', lw=3)
ax3.hist(d2y.bootstrap_param, bins=100, facecolor='0.5', edgecolor='w')
ax3.axvline(x=d2y.popt[2], linewidth=3, color='r')
ax3.axvline(d2y.ci_param[0], ymin=0, ymax=0.3, linewidth=2, color='r')
ax3.axvline(d2y.ci_param[1], ymin=0, ymax=0.3, linewidth=2, color='r')
ax4.plot(d2y.xdata, d2y.ydata, 'o', markersize=11, markerfacecolor='w',
markeredgecolor='k', markeredgewidth=2)
ax4.plot(d2y.xfit, d2y.yfit, '-r', lw=3, alpha=1)
ax4.fill_between(d2y.xfit, d2y.ci_fitlow, d2y.ci_fitup, color='k',
alpha=0.5, edgecolor='none')
ax4.plot(d2y.xfit, d2y.ci_fitup, '-r', lw=2, alpha=0.5)
ax4.plot(d2y.xfit, d2y.ci_fitlow, '-r', lw=2, alpha=0.5)
plt.tight_layout()