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

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?

• So a,b,c are fixed and you're minimizing with respect to x? Dec 25, 2020 at 16:49
• thanks for your comment @Richard. No $x$ is given. ($\textbf{x}$ and $\textbf{y}$ are the data). And I want to find the optimal $a,b,c,d$ to minimize $\sum_{i=1}^{n}{(\text{x}_i-\text{y}_i)^2}$. It is a curve fitting problem
– ecjb
Dec 25, 2020 at 16:56
• Could you edit your question to clarify? d should be part of the argmin and the phrase "a=1 in the generated data" seems confusing. Dec 25, 2020 at 17:16
• Is there any reason why you have that weird cost function $||\mathbf{x}-\mathbf{y}||^2$ instead of least square? I mean are you sure you are not looking for $||\mathbf{y}-f(\mathbf{x})||^2$? Cause I have a hard time understanding why you want to minimize the distance between measured $\mathbf{y}$ and your features $\mathbf{x}$. Dec 25, 2020 at 19:18
• Those constraints are not basically constraints in a classical definition. One of them is simply a given parameter for $a$ which you can replace $a$ with 1 everywhere. The other one just replaces $c$ with $b$ and $d$ everywhere and then finally you just fit $\mathbf{y}$ to $f(\mathbf{x},b,d)$. Simply the scipy’s example given in their documentation works perfectly fine here. Dec 25, 2020 at 23:04

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: • Happy Holidays! I hope things are well with you. Dec 26, 2020 at 17:56
• @Richard Happy holidays and thanks for your kind message! Dec 26, 2020 at 18:24