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

Question

在下面的代码中，我有 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()

Answer 1

你的 X 数字太大了。也许你可以尝试取两边的对数并拟合？例如：

\n

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

\n

此时您甚至不需要 curve_fit，它是标准线性回归。

\n

编辑

\n

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

\n

import numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.optimize as optimize\n\nx_mean = [2.81838293e+20, 5.62341325e+20, 1.12201845e+21, 2.23872114e+21, 4.46683592e+21, 8.91250938e+21, 1.77827941e+22]\n\nmean_1 = [52., 21.33333333, 4., 1., 0., 0.,  0.]\nmean_2 = [57., 16.66666667, 5.66666667, 2.33333333, 0.66666667, 0., 0.33333333]\nmean_3 = [67.33333333, 20., 8.66666667, 3., 0.66666667, 1., 0.33333333]\nmean_4 = [79.66666667, 25., 8.33333333, 3., 1., 0., 0.]\nmean_5 = [54.66666667, 16.66666667, 8.33333333, 2., 2., 1., 0.]\n\nerror_1 = [4.163332, 2.66666667, 1.15470054, 0.57735027, 0., 0., 0.]\nerror_2 = [4.35889894, 2.3570226, 1.37436854, 0.8819171, 0.47140452, 0., 0.33333333]\nerror_3 = [4.7375568, 2.5819889, 1.69967317, 1., 0.47140452, 0.57735027, 0.33333333]\nerror_4 = [5.15320828, 2.88675135, 1.66666667, 1., 0.57735027, 0., 0.]\nerror_5 = [4.26874949, 2.3570226, 1.66666667, 0.81649658, 0.81649658, 0.57735027, 0.]\n\n\ndef powerlaw(x, amp, index):\n    return amp * (x**index)\n    \n# define our (line) fitting function\ndef fitfunc(p, x):\n    return p[0] + p[1] * x   \n\ndef errfunc(p, x, y, err):\n    out = (y - fitfunc(p, x)) / err\n    out[~np.isfinite(out)] = 0.0\n    return out\n\npinit = [1.0, -1.0]\nfig = plt.figure()\nax1 = fig.add_subplot(2, 1, 1)\nax2 = fig.add_subplot(2, 1, 2)\n\nfor indx in range(1, 6):\n\n    mean = eval(\'mean_%d\'%indx)\n    error = eval(\'error_%d\'%indx)\n    logx = np.log10(x_mean)\n    logy = np.log10(mean)\n    logy[~np.isfinite(logy)] = 0.0\n    logyerr = np.array(error) / np.array(mean)\n    logyerr[~np.isfinite(logyerr)] = 0.0\n\n    out = optimize.leastsq(errfunc, pinit, args=(logx, logy, logyerr), full_output=1)\n\n    pfinal = out[0]\n    covar = out[1]\n    \n    index = pfinal[1]\n    amp = 10.0**pfinal[0]\n\n    indexErr = np.sqrt(covar[0][0] ) \n    ampErr = np.sqrt(covar[1][1] ) * amp\n\n    ##########\n    # Plotting data\n    ##########\n\n    ax1.plot(x_mean, powerlaw(x_mean, amp, index), label=u\'{:.2f} \\u00B1 {:.2f}\'.format(pfinal[1], covar[1,1]**0.5))     # Fit\n    ax1.errorbar(x_mean, mean, yerr=error, fmt=\'k.\', label=\'__no_legend__\')  # Data\n    ax1.set_title(\'Best Fit Power Law\', fontsize=18, fontweight=\'bold\')\n    ax1.set_xlabel(\'X\', fontsize=14, fontweight=\'bold\')\n    ax1.set_ylabel(\'Y\', fontsize=14, fontweight=\'bold\')\n    ax1.grid()\n\n    ax2.loglog(x_mean, powerlaw(x_mean, amp, index), label=u\'{:.2f} \\u00B1 {:.2f}\'.format(pfinal[1], covar[1,1]**0.5))\n    ax2.errorbar(x_mean, mean, yerr=error, fmt=\'k.\', label=\'__no_legend__\')  # Data\n    ax2.set_xlabel(\'X (log scale)\', fontsize=14, fontweight=\'bold\')\n    ax2.set_ylabel(\'Y (log scale)\', fontsize=14, fontweight=\'bold\')\n    ax2.grid(b=True, which=\'major\', linestyle=\'--\', color=\'darkgrey\')\n    ax2.grid(b=True, which=\'minor\', linestyle=\':\', color=\'grey\')\n\n\nax1.legend()\nax2.legend()\nplt.show()\n

\n

图片：

\n

\n