How to best code a problem with scipy, cvxpy or Convex.jl with given generated data

Question

I have a curve fitting problem of the form: $$ \textbf{y} = f(\textbf{x}, a,b,c,d) + \varepsilon $$ $$ f(x, a,b,c,d) = \frac{b}{e^{x\cdot a}+c}+d $$ with the constraint \begin{equation} \begin{aligned} \arg\min_{a,b,c,d} \quad & \sum_{i=1}^{n}{(f(\text{x})_i-\text{y}_i)^2}\\ \textrm{s.t.} \quad & -c-\frac{b}{d}-1 = 0\\ \quad & a = 1 \\ \end{aligned} \end{equation}

where

• $\textbf{x}, \textbf{y} \in \mathbb{R}^n$. $\textbf{x}, \textbf{y}$ are the data given by an experiment
• $a \in \mathbb{R}$
• $b \in \mathbb{R}$
• $c \in \mathbb{R}$
• $d \in \mathbb{R}$ Note: $a = 1$ because this parameter is being measured by a measurement in an experience and is, therefore fixed.

The simulated data and the following code are as follows:

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
n = 20
np.random.seed(1)
def sigm01(x,a,b,c,d):
  return b/(np.exp(x*a)+c)+d

start = -10
end = -start
x = np.arange(start,end,0.01)
a = 1
b = -2
c = 1
d = 1
y = sigm01(x,a,b,c,d) + np.random.randn(len(x))/10
#plt.plot(x,y)
#plt.show()

I would like to use an optimization package (such as scipy, cvxpy or Convex.jl).

How could I setup the problem?

Answer 1

With your constraints of $a = 1$ and $c = -\frac{b}{d} - 1$, your $f$ looks like this:

$$\hat{\mathbf{y}} = f(\mathbf{x},b,d) = \frac{bd}{d\exp{(\mathbf{x})}-b - d} + d$$

You're trying to solve a nonlinear least square problem basically:

$$b,d = \arg \min \sum_{i=1}^{N} (f(x_{i})-y_{i})^{2}$$

It's a job for scipy.optimize.curve_fit:

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

def func(x,b,d):
    return ((b*d)/(d*np.exp(x)-b-d))+d

xdata = np.linspace(-10, 10, 10000)
y = func(xdata, -2, 1)
np.random.seed(1729)
y_noise = 0.2 * np.random.normal(size=xdata.size)
ydata = y + y_noise

popt, pcov = curve_fit(func, xdata, ydata, bounds=([-3, 0], [-1, 2.]))

plt.plot(xdata, ydata, 'b-', label='data')
plt.plot(xdata, func(xdata, *popt), 'r-',
         label='fit: b=%5.3f, d=%5.3f' % tuple(popt))
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('curve_fit.png')
plt.show()

The final fit looks like this: