曲线拟合 - 函数失败

Question

曲线是：

import numpy as np
import scipy.stats as sp
from scipy.optimize import curve_fit
from lmfit import minimize, Parameters, Parameter, report_fit#
import xlwings as xw
import os
import pandas as pd

我尝试了运行来自 scipy 的简单曲线拟合：这个returns

Out[156]: 
(array([ 1.,  1.]), array([[ inf,  inf],
        [ inf,  inf]]))

好的，所以我想我需要开始限制我的参数，这是求解器所做的。我用 lmfit 来做这个：

params = Parameters()

params.add

我试图通过更改起始参数将 python 推向正确的方向：

或使用曲线拟合：

Answer 1

我认为部分问题在于您只有 5 个观测值，其中 2 个与 x 的值相同，并且该模型不能完美地代表您的数据。我还建议尝试将模型日志与数据日志相匹配。而且，如果您希望 n2 为 ~10，则应将其用作起始值。

任意将 min 和 max 应用于计算模型，尤其是值非常接近您的数据范围时，意义不大。这样做会阻止 fit 方法探索更改参数值对模型函数的影响。

通过对您的代码进行一些修改以省略第一个数据点（这似乎与模型脱节），我发现：

import numpy as np
import scipy.stats as sp
from lmfit import Parameters, minimize, fit_report

import matplotlib.pyplot as plt

x_data = np.array((1e-04,9e-01,9.5835e-01,9.8e-01,9.9e-01,9.9e-01))
y_data = np.array((250e3,1e6,2.5e6,5e6,7.5e6,10e6))

x_data = x_data[1:]
y_data = y_data[1:]

def func(params, x, data, m1=250e3,m2=10e6):
    n1 = params['n1'].value
    n2 = params['n2'].value

    model = n2 + n1*(sp.norm.ppf(x))
    return np.log(data) - model

params = Parameters()

params.add('n1', value= 1, min=0.01, max=20)
params.add('n2', value= 10, min=0, max=20)

result = minimize(func, params, args=(x_data[1:-1], y_data[1:-1]))

print(fit_report(result))
n1 = result.params['n1'].value
n2 = result.params['n2'].value 
ym = np.exp(n2 + n1*(sp.norm.ppf(x_data)))

plt.plot(x_data, y_data, 'o')
plt.plot(x_data, ym, '-')
plt.legend(['data', 'fit'])
plt.show()

给出

的报告

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 15
    # data points      = 3
    # variables        = 2
    chi-square         = 0.006
    reduced chi-square = 0.006
    Akaike info crit   = -14.438
    Bayesian info crit = -16.241
[[Variables]]
    n1:   1.85709072 +/- 0.190473 (10.26%) (init= 1)
    n2:   11.5455736 +/- 0.390805 (3.38%) (init= 10)
[[Correlations]] (unreported correlations are <  0.100)
    C(n1, n2)                    = -0.993

和一个情节

Answer 2

当我们提供一个好的起点时，曲线拟合运行顺利。我们可以得到一个

通过 sp.norm.ppf(x_data) 和 np.log(y_data)
或先拟合自由（未裁剪）模型

或者，如果您希望计算机在没有 "help"

的情况下找到解决方案

使用像 basin hopping 这样的随机算法（不幸的是它不是 100% 自主的，我不得不从默认值增加步长）
暴力破解（这需要用户提供搜索网格）

所有四种方法产生相同的结果，并且比 excel 结果更好。

import numpy as np
import scipy.stats as sp
from scipy.optimize import curve_fit, basinhopping, brute

def func1(x, n1, n2):
    return np.clip(np.exp(n2 + n1*(sp.norm.ppf(x))),25e4,10e6)

def func_free(x, n1, n2):
    return np.exp(n2 + n1*(sp.norm.ppf(x)))

def sqerr(n12, x, y):
    return ((func1(x, *n12) - y)**2).sum()

x_data = np.array((1e-04,9e-01,9.5835e-01,9.8e-01,9.9e-01,9.9e-01))
y_data = np.array((250e3,1e6,2.5e6,5e6,7.5e6,10e6))

# get a good starting point

# either by linear regression
lin = sp.linregress(sp.norm.ppf(x_data), np.log(y_data))

# or by using the free (non-clipped) version of the formula
(n1f, n2f), infof = curve_fit(func_free, x_data, y_data, (1, 1))

# use those on the original problem 
(n1, n2), info = curve_fit(func1, x_data, y_data, (lin.slope, lin.intercept))
(n12, n22), info2 = curve_fit(func1, x_data, y_data, (n1f, n2f))

# OR

# use basin hopping
hop = basinhopping(sqerr, (1, 1), minimizer_kwargs=dict(args=(x_data, y_data)), stepsize=10)

# OR 

# brute force it
brt = brute(sqerr, ((-100, 100), (-100, 100)), (x_data, y_data), 201, full_output=True)

# all four solutions are essentially the same:
assert np.allclose((n1, n2), (n12, n22))
assert np.allclose((n1, n2), hop.x)
assert np.allclose((n1, n2), brt[0])

# we are actually a bit better than excel
n1excel, n2excel = 1.7925, 11.6771

print('solution', n1, n2)
print('error', ((func1(x_data, n1, n2) - y_data)**2).sum())
print('excel', ((func1(x_data, n1excel, n2excel) - y_data)**2).sum())

输出：

solution 2.08286042997 11.1397332743
error 3.12796761241e+12
excel 5.80088578059e+12

备注：一个简单的优化——为了简单起见我省略了它，因为无论如何事情已经足够快了——本来可以将 sp.norm.ppf 从模型函数中拉出来。这是可能的，因为它不依赖于拟合参数。因此，当我们的任何求解器调用该函数时，它总是首先执行完全相同的计算 - sp.norm.ppf(x_data) - 所以我们不妨预先计算它。

这个观察也是我们在线性回归中使用sp.norm.ppf(x_data)的原因。

Answer 3

这是使用 scipy.optimize.curve_fit 执行回归的简单方法：

import matplotlib.pyplot as plt
import scipy.optimize as opt
import scipy.stats as stats
import numpy as np

% matplotlib inline


# Objective
def model(x, n1, n2):
    return np.exp(n2 + n1*(stats.norm.ppf(x)))

# Data
x_samp = np.array((1e-04, 9e-01, 9.5835e-01, 9.8e-01, 9.9e-01, 9.9e-01))
y_samp = np.array((250e3, 1e6, 2.5e6, 5e6 ,7.5e6, 10e6))

x_lin = np.linspace(min(x_samp), max(x_samp), 50)  # for fitting

# Regression
p0 = [5, 5]                                        # guessed params
w, cov = opt.curve_fit(model, x_samp, y_samp, p0=p0)     
print("Estimated Parameters", w)  
y_fit = model(x_lin, *w)

# Visualization
plt.plot(x_samp, y_samp, "ko", label="Data")
plt.plot(x_lin, y_fit, "k--", label="Fit")
plt.title("Curve Fitting")
plt.legend(loc="upper left")

输出

Estimated Parameters [  2.08285048  11.13975585]

详情

您的数据是按原样绘制的，没有进行转换。如果模型和 初始参数 、p0 合理地适合您的数据，则最容易执行此类回归。当这些项目被提供给 scipy.optimize.curve_fit 时，返回权重或优化的 估计参数 、w 的元组，连同 covariance matrix。我们可以从矩阵的对角线计算一个标准差：

p_stdev = np.sqrt(np.diag(cov)) 
print("Std. Dev. of Params:", p_stdev)
# Std. Dev. of Params: [ 0.42281111  0.95945127]

我们通过再次向模型提供这些估计参数并绘制拟合线来直观地评估拟合优度。

曲线拟合 - 函数失败

Curve fit - function failing

python

excel

curve-fitting