fit_negative_exponential.py

from scipy.optimize import curve_fit
from scipy import stats
import matplotlib.pyplot as plt


def fit_exponential_neg(x, a, b, c):
    return a * np.exp(-b * x) + c

X = np.array(rpkm_log['mean'])
Y = np.array(rpkm_log['qv2'])


ci = 0.99
# Convert to percentile point of the normal distribution.
# See: https://en.wikipedia.org/wiki/Standard_score
pp = (1. + ci) / 2.

nstd = stats.norm.ppf(pp)

# Find best fit.
parameters, covariance_matrix = curve_fit(fit_exponential_neg, X, Y)
# Standard deviation errors on the parameters.
perr = np.sqrt(np.diag(pcov))
# Add nstd standard deviations to parameters to obtain the upper confidence
# interval.
popt_up = parameters + (nstd * perr)
popt_dw = parameters - (nstd * perr)

fig, axis = plt.subplots(2, sharex=True)
# Plot data and best fit curve.
axis[0].scatter(X, Y)
x = np.linspace(0, 6.5, 100)
axis[0].plot(x, fit_exponential_neg(x, *parameters), c='g', lw=2.)
axis[0].plot(x, fit_exponential_neg(x, *popt_up), c='r', lw=2.)
axis[0].plot(x, fit_exponential_neg(x, *popt_dw), c='r', lw=2.)
axis[0].text(12, 0.5, '{}% confidence interval'.format(ci * 100.))
axis[0].set_title("fit")

# plot residuals
residuals = Y - fit_exponential_neg(X, *parameters)
fres = sum(residuals ** 2)

axis[1].scatter(X, residuals)
axis[1].set_title("residuals")
plt.show()