如何正确猜测 LogLog 图线性回归中的初始点？答案

【问题标题】：How to correctly guess the initial points in LogLog plot linear regression?如何正确猜测 LogLog 图线性回归中的初始点？
【发布时间】：2020-10-08 23:19:48
【问题描述】：

我在下面的代码中有 5 组数据，用 5 个不同颜色的错误条表示（我没有显示大写字母）。 errorbar plot 在两个轴上都以对数刻度显示。使用curvefit，我试图找到通过这些误差线的最佳线性回归。但是，我定义的幂律方程似乎不容易找到 5 条线的最佳拟合斜率。我的期望是所有 5 条彩色线都应该是直的，带有负斜率。我很难弄清楚应该在曲线拟合过程中指定哪个起点p0。即使有了我最初难以猜测的值，我仍然没有得到所有的直线，其中一些与我的观点相差太远。这里有什么问题？

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

x_mean = [2.81838293e+20, 5.62341325e+20, 1.12201845e+21, 2.23872114e+21, 4.46683592e+21, 8.91250938e+21, 1.77827941e+22]

mean_1 = [52., 21.33333333, 4., 1., 0., 0.,  0.]
mean_2 = [57., 16.66666667, 5.66666667, 2.33333333, 0.66666667, 0., 0.33333333]
mean_3 = [67.33333333, 20., 8.66666667, 3., 0.66666667, 1., 0.33333333]
mean_4 = [79.66666667, 25., 8.33333333, 3., 1., 0., 0.]
mean_5 = [54.66666667, 16.66666667, 8.33333333, 2., 2., 1., 0.]

error_1 = [4.163332, 2.66666667, 1.15470054, 0.57735027, 0., 0., 0.]
error_2 = [4.35889894, 2.3570226, 1.37436854, 0.8819171, 0.47140452, 0., 0.33333333]
error_3 = [4.7375568, 2.5819889, 1.69967317, 1., 0.47140452, 0.57735027, 0.33333333]
error_4 = [5.15320828, 2.88675135, 1.66666667, 1., 0.57735027, 0., 0.]
error_5 = [4.26874949, 2.3570226, 1.66666667, 0.81649658, 0.81649658, 0.57735027, 0.]


newX = np.logspace(20, 22.3)
def myExpFunc(x, a, b):
    return a*np.power(x, b)

popt_1, pcov_1 = curve_fit(myExpFunc, x_mean, mean_1, sigma=error_1, absolute_sigma=True, p0=(4e31,-1.5))
popt_2, pcov_2 = curve_fit(myExpFunc, x_mean, mean_2, sigma=error_2, absolute_sigma=True, p0=(4e31,-1.5))
popt_3, pcov_3 = curve_fit(myExpFunc, x_mean, mean_3, sigma=error_3, absolute_sigma=True, p0=(4e31,-1.5))
popt_4, pcov_4 = curve_fit(myExpFunc, x_mean, mean_4, sigma=error_4, absolute_sigma=True, p0=(4e31,-1.5))
popt_5, pcov_5 = curve_fit(myExpFunc, x_mean, mean_5, sigma=error_5, absolute_sigma=True, p0=(4e31,-1.5))


fig, ax1 = plt.subplots(figsize=(3,5))
ax1.errorbar(x_mean, mean_1, yerr=error_1, ecolor = 'magenta', fmt= 'mo', ms=0, elinewidth = 1, capsize = 0, capthick=0)
ax1.errorbar(x_mean, mean_2, yerr=error_2, ecolor = 'red', fmt= 'ro', ms=0, elinewidth = 1, capsize = 0, capthick=0)
ax1.errorbar(x_mean, mean_3, yerr=error_3, ecolor = 'orange', fmt= 'yo', ms=0, elinewidth = 1, capsize = 0, capthick=0)
ax1.errorbar(x_mean, mean_4, yerr=error_4, ecolor = 'green', fmt= 'go', ms=0, elinewidth = 1, capsize = 0, capthick=0)
ax1.errorbar(x_mean, mean_5, yerr=error_5, ecolor = 'blue', fmt= 'bo', ms=0, elinewidth = 1, capsize = 0, capthick=0)

ax1.plot(newX, myExpFunc(newX, *popt_1), 'm-', label='{:.2f} \u00B1 {:.2f}'.format(popt_1[1], pcov_1[1,1]**0.5))
ax1.plot(newX, myExpFunc(newX, *popt_2), 'r-', label='{:.2f} \u00B1 {:.2f}'.format(popt_2[1], pcov_2[1,1]**0.5))
ax1.plot(newX, myExpFunc(newX, *popt_3), 'y-', label='{:.2f} \u00B1 {:.2f}'.format(popt_3[1], pcov_3[1,1]**0.5))
ax1.plot(newX, myExpFunc(newX, *popt_4), 'g-', label='{:.2f} \u00B1 {:.2f}'.format(popt_4[1], pcov_4[1,1]**0.5))
ax1.plot(newX, myExpFunc(newX, *popt_5), 'b-', label='{:.2f} \u00B1 {:.2f}'.format(popt_5[1], pcov_5[1,1]**0.5))
ax1.legend(handlelength=0, loc='upper right', ncol=1, fontsize=10)

ax1.set_xlim([2e20, 3e22])
ax1.set_ylim([2e-1, 1e2])
ax1.set_xscale("log")
ax1.set_yscale("log")
plt.show()

【问题讨论】：

首先：尝试扩展您的数据。将x 除以1e21，您基本上不再需要初始猜测了。真正的a,b 稍后通过简单的错误传播得到。第二：如果处理错误，有一些为零，函数实际上可以通过输入a = 0 实现完美匹配不是一个好主意。如果这是一次测量，我很确定测量存在相对误差和绝对误差，因此零信号当然具有零相对误差，但总误差肯定大于零。那会是什么？
最后提示：制作mean, popt 等字典。像这样，您可以循环它并避免在拟合和绘图中重复代码。
关于您的第一个建议，看来我仍然需要指定 p0，但我想我不太理解您对错误栏的讨论。哪里有错误栏的 NaN，确实我没有从实验中获得的数据，但我认为你是说由于缺乏约束，我应该将错误栏的大小增加到 (-inf, inf)。正确的？你介意分享一个你正在应用你所说的话的代码吗？
以下是应用您对约束的建议之前的症状： py:734: RuntimeWarning: 除以零在 true_divide transform = 1.0 / sigma /home/username/anaconda3/lib/python3.7/site 中遇到-packages/scipy/optimize/minpack.py:808: OptimizeWarning: 无法估计参数的协方差 category=OptimizeWarning)
您好，在我在这里复制的内容中，如果不考虑任何错误，一切正常。如果你输入一个零错误的零值是不是你的意思是“没有数据”？这有很大的不同。由于误差作为权重，零误差为无限权重；这没有意义。没有错误的点将是一个约束。那么，是否应该将零值作为“无可用数据”从列表中删除？

标签： python-3.x optimization scipy curve-fitting

【解决方案1】：

你的 X 数字太大了。也许您可以尝试记录双方的日志并进行拟合？如：

log Y = log(a) + b*log(X)

那时你甚至不需要curve_fit，它是一个标准的线性回归。

编辑

请查看我的粗略且未经过很好检查的实现（注意：我只有 Python 2，因此请调整以适应）：

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

x_mean = [2.81838293e+20, 5.62341325e+20, 1.12201845e+21, 2.23872114e+21, 4.46683592e+21, 8.91250938e+21, 1.77827941e+22]

mean_1 = [52., 21.33333333, 4., 1., 0., 0.,  0.]
mean_2 = [57., 16.66666667, 5.66666667, 2.33333333, 0.66666667, 0., 0.33333333]
mean_3 = [67.33333333, 20., 8.66666667, 3., 0.66666667, 1., 0.33333333]
mean_4 = [79.66666667, 25., 8.33333333, 3., 1., 0., 0.]
mean_5 = [54.66666667, 16.66666667, 8.33333333, 2., 2., 1., 0.]

error_1 = [4.163332, 2.66666667, 1.15470054, 0.57735027, 0., 0., 0.]
error_2 = [4.35889894, 2.3570226, 1.37436854, 0.8819171, 0.47140452, 0., 0.33333333]
error_3 = [4.7375568, 2.5819889, 1.69967317, 1., 0.47140452, 0.57735027, 0.33333333]
error_4 = [5.15320828, 2.88675135, 1.66666667, 1., 0.57735027, 0., 0.]
error_5 = [4.26874949, 2.3570226, 1.66666667, 0.81649658, 0.81649658, 0.57735027, 0.]


def powerlaw(x, amp, index):
    return amp * (x**index)
    
# define our (line) fitting function
def fitfunc(p, x):
    return p[0] + p[1] * x   

def errfunc(p, x, y, err):
    out = (y - fitfunc(p, x)) / err
    out[~np.isfinite(out)] = 0.0
    return out

pinit = [1.0, -1.0]
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)

for indx in range(1, 6):

    mean = eval('mean_%d'%indx)
    error = eval('error_%d'%indx)
    logx = np.log10(x_mean)
    logy = np.log10(mean)
    logy[~np.isfinite(logy)] = 0.0
    logyerr = np.array(error) / np.array(mean)
    logyerr[~np.isfinite(logyerr)] = 0.0

    out = optimize.leastsq(errfunc, pinit, args=(logx, logy, logyerr), full_output=1)

    pfinal = out[0]
    covar = out[1]
    
    index = pfinal[1]
    amp = 10.0**pfinal[0]

    indexErr = np.sqrt(covar[0][0] ) 
    ampErr = np.sqrt(covar[1][1] ) * amp

    ##########
    # Plotting data
    ##########

    ax1.plot(x_mean, powerlaw(x_mean, amp, index), label=u'{:.2f} \u00B1 {:.2f}'.format(pfinal[1], covar[1,1]**0.5))     # Fit
    ax1.errorbar(x_mean, mean, yerr=error, fmt='k.', label='__no_legend__')  # Data
    ax1.set_title('Best Fit Power Law', fontsize=18, fontweight='bold')
    ax1.set_xlabel('X', fontsize=14, fontweight='bold')
    ax1.set_ylabel('Y', fontsize=14, fontweight='bold')
    ax1.grid()

    ax2.loglog(x_mean, powerlaw(x_mean, amp, index), label=u'{:.2f} \u00B1 {:.2f}'.format(pfinal[1], covar[1,1]**0.5))
    ax2.errorbar(x_mean, mean, yerr=error, fmt='k.', label='__no_legend__')  # Data
    ax2.set_xlabel('X (log scale)', fontsize=14, fontweight='bold')
    ax2.set_ylabel('Y (log scale)', fontsize=14, fontweight='bold')
    ax2.grid(b=True, which='major', linestyle='--', color='darkgrey')
    ax2.grid(b=True, which='minor', linestyle=':', color='grey')


ax1.legend()
ax2.legend()
plt.show()

图片：

【讨论】：

谢谢，考虑到约束的重要性，即使这样做似乎也没有考虑数据的存在或不存在。您介意实施您的建议以了解我的意思吗？
我不确定我是否理解：我们在谈论什么样的数据存在/不存在？有哪些限制？
如何在python3中修改这两行？ mean = eval('mean_%d'%indx) error = eval('error_%d'%indx) 我不确定你的意思！
这是一种快速获取mean_1、mean_2等...以及error_1、error_2等值的方法。你不需要我修改这些行，这是基本的 Python，有很多方法可以做到（也许制作一个 2D numpy 数组来存储错误和方法？）。我的方法是一种快速而肮脏的解决方案，使用非正统的工具（eval）。而且我不在乎，它可以工作，而且可以轻松修改。
据我所知，“eval”也适用于 Python 3...