python的几种非线性回归

前一阵子有人和我吐槽过matlab内置的几个线性数据拟合工具满足不了需求，今天正好看到了一个关于使用scipy进行非线性回归的工具使用方法，写下来备忘。

Theory

Given pairs:{(ti,yi),i=1,...,n}\{(t_i,y_i), i=1,...,n\}{(ti,yi),i=1,...,n} estimate parameters x defining a nonlinear function ϕ(t;x)\phi(t;x)ϕ(t;x) estimate parameters xxx defining a nonlinear function phi(t;x)phi(t;x)phi(t;x), assuming the models:
yi=ϕ(ti;x)+ϵiy_i = \phi(t_i;x) + \epsilon_i yi=ϕ(ti;x)+ϵi
where ϵi\epsilon_iϵi is the measurement errors. Then use the least-square estimation to search xxx as the solution of the optimization problem:
12∑i=1n(ϕ(ti;x)−yi)2→min(x)\frac{1}{2}\sum_{i = 1}^{n}(\phi(t_i;x)-y_i)^2 \rightarrow min(x) 21i=1∑n(ϕ(ti;x)−yi)2→min(x)
Then, we introduce a sublinear function ρ(z)\rho(z)ρ(z) and formulate a new least-square like optimization problem:
12∑i=1nρ((ϕ(ti;x)−yi)2)→min(x)\frac{1}{2}\sum_{i = 1}^{n}\rho((\phi(t_i;x)-y_i)^2) \rightarrow min(x) 21i=1∑nρ((ϕ(ti;x)−yi)2)→min(x)

Standard least-square method (linear) ρ(z)=z\rho(z) = zρ(z)=z.
Huber loss:ρ(z)=z(z<=1)z(z>1)\rho(z) = z (z<=1) \sqrt{z} (z>1)ρ(z)=z(z<=1)z(z>1).
Smooth approximation to absolute value loss:ρ(z)=2(1+z−1)\rho(z) = 2(\sqrt{1+z}-1)ρ(z)=2(1+z−1).
Cauchy loss: ρ(z)=ln(1+z)\rho(z) = ln(1+z)ρ(z)=ln(1+z)
Loss by arctan:ρ(z)=arctan(z)\rho(z) = arctan(z)ρ(z)=arctan(z)

Code

首先生成一些数据：
f(t;A,σ,ω)=Ae−σtsin⁡(ωt)f(t; A, \sigma, \omega) = A e^{-\sigma t} \sin (\omega t)f(t;A,σ,ω)=Ae−σtsin(ωt)

def generate_data(t, A, sigma, omega, noise=0, n_outliers=0, random_state=0):y = A * np.exp(-sigma * t) * np.sin(omega * t)rnd = np.random.RandomState(random_state)error = noise * rnd.randn(t.size)outliers = rnd.randint(0, t.size, n_outliers)error[outliers] *= 35return y + error
########model parameters
A = 2
sigma = 0.1
omega = 0.1 * 2 * np.pi
x_true = np.array([A, sigma, omega])
noise = 0.1
t_min = 0
t_max = 30
t_train = np.linspace(t_min, t_max, 30)
y_train = generate_data(t_train, A, sigma, omega, noise=noise, n_outliers=4)

定义残差

def fun(x, t, y):return x[0] * np.exp(-x[1] * t) * np.sin(x[2] * t) - y

开始进行回归：

x0 = np.ones(3)
#Run standard least squares:
from scipy.optimize import least_squares
res_lsq = least_squares(fun, x0, args=(t_train, y_train))
#Run robust least squares with loss='soft_l1', set f_scale to 0.1 which means that inlier residuals are approximately lower than 0.1.
res_robust = least_squares(fun, x0, loss='soft_l1', f_scale=0.1, args=(t_train, y_train))
#post
t_test = np.linspace(t_min, t_max, 300)
y_test = generate_data(t_test, A, sigma, omega)
y_lsq = generate_data(t_test, *res_lsq.x)
y_robust = generate_data(t_test, *res_robust.x)plt.plot(t_train, y_train, 'o', label='data')
plt.plot(t_test, y_test, label='true')
plt.plot(t_test, y_lsq, label='lsq')
plt.plot(t_test, y_robust, label='robust lsq')
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend();

可以看到非线性回归的效果要好于标准最小二乘法