Finding x for a bi-exponential curve fit?

Question

I am using python to batch process some data and plot it. I can fit it quite well using scipy.curve_fit, a bi-exponential function and some sensible initial guesses. Here is a code snippet:

def biexpfunc(x, a, b, c, d, e):
    y_new = []
    for i in range(len(x)):
        y = (a * np.exp(b*x[i])) + (c * np.exp(d*x[i])) + e
        y_new.append(y)
    return y_new

x = np.linspace(0, 160, 100)
y = biexpfunc(x, 50, -0.2, 50, -0.1, 10)
jitter_y = y + 0.5 *np.random.rand(len(y)) - 0.1
plt.scatter(x, jitter_y)


sigma = np.ones(len(x))
sigma[[0, -1]] = 0.01
popt, pcov = curve_fit(biexpfunc, x, jitter_y, p0 = (50, -0.2, 50, -0.1, 10), 
sigma = sigma)
x_fit = np.linspace(0, x[-1])
y_fit = biexpfunc(x_fit, *popt)
plt.plot(x_fit, y_fit, 'r--')

plt.show()

I know how to interpolate this to find y for a given value of x (by putting it back into the function), but how can I find x for a given value of y? I feel like there must be a sensible method that doesn't require re-arrangement and defining a new function (partially because maths is not my strong suit and I don't know how to!). If the curve fits the data well is there a way to simply read off a value? Any assistance would be greatly appreciated!

Answer 1

Turns out, your question has nothing to do with curve fitting but is actually about root finding. Scipy.optimize has a whole arsenal of functions for this task. Choosing and configuring the right one is sometimes difficult. I might not be the best guide here, but since no-one else stepped up...
Root finding tries to determine x-values for which f(x) is zero. To find an x0 where f(x0) is a certain y0-value, we simply transform the function into g(x) = f(x)-y0. Since your function is monotonous, not more than one root is to be expected for a given y-value. We also know the x-interval in which to search, so bisect seems to be a reasonable strategy:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit, bisect

def biexpfunc(x, a, b, c, d, e):
    return (a * np.exp(b*x)) + (c * np.exp(d*x)) + e

np.random.seed(123)
x = np.linspace(0, 160, 100)
y = biexpfunc(x, 50, -0.2, 50, -0.1, 10)
jitter_y = y + 0.5 *np.random.rand(len(y)) - 0.1

fig, ax = plt.subplots(figsize=(10, 8))
ax.scatter(x, jitter_y, marker="x", color="blue", label="raw data")

#your curve fit routine
sigma = np.ones(len(x))
sigma[[0, -1]] = 0.01
popt, pcov = curve_fit(biexpfunc, x, jitter_y, p0 = (50, -0.2, 50, -0.1, 10), sigma = sigma)
x_fit = np.linspace(x.min(), x.max(), 100)
y_fit = biexpfunc(x_fit, *popt)
ax.plot(x_fit, y_fit, 'r--', label="fit")

#y-value for which we want to determine the x-value(s)
y_test=55
test_popt = popt.copy()
test_popt[-1] -= y_test

#here, the bisect method tries to establish the x for which f(x)=0
x_test=bisect(biexpfunc, x.min(), x.max(), args=tuple(test_popt))
#we calculate the deviation from the expected y-value
tol_test, = np.abs(y_test - biexpfunc(np.asarray([x_test]), *popt))

#and mark the determined point in the graph
ax.axhline(y_test, ls="--", color="grey")
ax.axvline(x_test, ls="--", color="grey")
ax.plot(x_test, y_test, c="tab:orange", marker="o", markersize=15, alpha=0.5)
ax.annotate(f"X: {x_test:.2f}, Y: {y_test:.2f}\ntol: {tol_test:.4f}", 
            xy=(x_test, y_test), xytext=(50, 50), textcoords="offset points", 
            arrowprops=dict(facecolor="tab:orange", shrink=0.05),)

ax.legend(title="root finding: bisect")
plt.show()

Sample output:

Another way to determine roots for more complex functions is, surprise, root. The script is mainly identical, only the root-routine is slightly different, for instance, we can choose the root-finding method:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit, root

def biexpfunc(x, a, b, c, d, e):
    return (a * np.exp(b*x)) + (c * np.exp(d*x)) + e

np.random.seed(123)
x = np.linspace(0, 160, 100)
y = biexpfunc(x, 50, -0.2, 50, -0.1, 10)
jitter_y = y + 0.5 *np.random.rand(len(y)) - 0.1

fig, ax = plt.subplots(figsize=(10, 8))
ax.scatter(x, jitter_y, marker="x", color="blue", label="raw data")

#your curve fit routine
sigma = np.ones(len(x))
sigma[[0, -1]] = 0.01
popt, pcov = curve_fit(biexpfunc, x, jitter_y, p0 = (50, -0.2, 50, -0.1, 10), sigma = sigma)
x_fit = np.linspace(x.min(), x.max(), 100)
y_fit = biexpfunc(x_fit, *popt)
ax.plot(x_fit, y_fit, 'r--', label="fit")

#y-value for which we want to determine the x-value(s)
y_test=55
test_popt = popt.copy()
test_popt[-1] -= y_test

#calculate corresponding x-value with root finding
r=root(biexpfunc, x.mean(), args=tuple(test_popt), method="lm")
x_test, = r.x 
tol_test, = np.abs(y_test - biexpfunc(r.x, *popt))

#mark point in graph
ax.axhline(y_test, ls="--", color="grey")
ax.axvline(x_test, ls="--", color="grey")
ax.plot(x_test, y_test, c="tab:orange", marker="o", markersize=15, alpha=0.5)
ax.annotate(f"X: {x_test:.2f}, Y: {y_test:.2f}\ntol: {tol_test:.4f}", 
            xy=(x_test, y_test), xytext=(50, 50), textcoords="offset points", 
            arrowprops=dict(facecolor="tab:orange", shrink=0.05))

ax.legend(title="root finding: lm")
plt.show()

Sample output:

The graphs look in this case identical. This is not necessarily for every function so; just like for curve fitting, the correct approach can dramatically improve the outcome.