计算科学 - 如何使用给定的生成数据最好地编码 scipy、cvxpy 或 Convex.jl 的问题 - 吾爱随笔录

如何使用给定的生成数据最好地编码 scipy、cvxpy 或 Convex.jl 的问题

计算科学优化约束优化凸优化约束

2021-12-22 01:59:49

我有以下形式的曲线拟合问题：

y = f (x, a, b, c, d) + ε

$\textbf{y} = f(\textbf{x}, a,b,c,d) + \varepsilon$

f (x, a, b, c, d) = \frac{b}{e^{x \cdot a} + c} + d

$f(x, a,b,c,d) = \frac{b}{e^{x\cdot a}+c}+d$ 有约束

\begin{aligned} \arg min_{a, b, c, d} & \sum_{i = 1}^{n} (f (x)_{i} - y_{i})^{2} \\ s.t. & - c - \frac{b}{d} - 1 = 0 \\ a = 1 \end{aligned}

$\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}$

在哪里

$\textbf{x}, \textbf{y} \in \mathbb{R}^n$ 。 $\textbf{x}, \textbf{y}$ 是实验给出的数据
$a \in \mathbb{R}$
$b \in \mathbb{R}$
$c \in \mathbb{R}$
$d \in \mathbb{R}$ 注意：因为这个参数是由经验中的测量来测量的，因此是固定的。 $a = 1$

模拟数据和以下代码如下：

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()

我想使用优化包（例如scipy，cvxpy或Convex.jl）。

我该如何设置问题？

1个回答

和的约束，您的看起来像这样： $a = 1$ $c = -\frac{b}{d} - 1$ $f$

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

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

您基本上是在尝试解决非线性最小二乘问题：

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

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

这是一份工作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()

最终拟合如下所示：

其它你可能感兴趣的问题

上一篇SciPy 中的 nquad 集成下一篇Python中向后时间集成的问题