【问题标题】:Power law data fitting is not correct幂律数据拟合不正确
【发布时间】:2020-05-07 16:51:27
【问题描述】:

我在尝试使用对数 y 轴拟合以下数据时做错了什么。代码和生成的图表如下。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

fig, ax = plt.subplots()
x = np.array([88.08064516, 264.24193548, 440.40322581, 616.56451613, 792.72580645, 968.88709677, 1145.0483871, 1321.20967742, 1497.37096774, 1673.53225806, 1849.69354839, 2025.85483871, 2202.01612903, 2378.17741935, 2554.33870968, 2730.5, 2906.66129032, 3082.82258065, 3258.98387097, 3435.14516129, 3611.30645161, 3787.46774194, 3963.62903226, 4139.79032258, 4315.9516129, 4492.11290323, 4668.27419355, 4844.43548387, 5020.59677419, 5196.75806452, 5372.91935484, 5549.08064516])
y = np.array([210737, 2175, 514, 158, 90, 46, 27, 22, 10, 11, 3, 7, 3, 2, 0, 1, 1, 1, 0, 0, 1, 0, 0,0, 0, 0, 1, 0, 0, 0, 0,1])
popt, pcov = curve_fit(func, x, y)
ax.plot(x, func(x, *popt), 'g--')
ax.plot(x,  y, 'ro', label='data')

ax.set_yscale('log')  # I need to have the y-axis logarithmic
plt.show()

【问题讨论】:

  • curve_fit 使用什么函数?它是来自一个包还是你自己写的?
  • 对不起,它来自 from scipy.optimize import curve_fit
  • 所以我猜问题是执行拟合后的popt 值。你能打印出来吗?
  • 这很有帮助,因为curve_fit 以某种方式喜欢 在标准化范围内进行操作。但是我已经删除了我的评论,因为现在我看到您的数据看起来像是 在对数图中 的负参数的指数!!!!如果您尝试适合自己的不同功能,您会有更好的运气
  • 例如,从ln(y) = a exp(-(b x)) + c 我们有(即Emacs 的calc 有)y = exp(c + a exp(-(b x))),你可以尝试拟合这个双指数,不是吗?

标签: python matplotlib curve-fitting power-law


【解决方案1】:

所有代码实际上都在工作。请参见下图,我运行了您的代码,但在绘图之前没有将缩放比例放在 y 轴上。你可以看到这条线实际上是合适的,除了第一个点它是一个异常值。解决这个问题的方法是在拟合函数之前先缩放 y 值。试试看,如果您需要更多帮助,请告诉我们。

【讨论】:

  • 是的,第一个点肯定是异常值,但为了使 y 轴可读,我想要一个对数刻度
  • @JamesPhillips 你误解了我想说的话,我只是表明问题中的代码产生了一个“体面”的拟合,因为 y 轴是对数缩放的,所以看起来更糟。我没有建议一个新的解决方案,所以很明显我知道它不适合数据。另外,看看作者如何同意第一点是异常值。当然,它确实符合对数缩放后创建的曲线,但很明显,在这张图片中,一个点看起来不像其他点,请参见异常值的定义(远离主体或系统或与主体或系统分离的人或事物)。
  • 如果您愿意,可以将其删除。我认为这确实回答了这个问题。问题是“我做错了什么?”参考标题“幂律数据拟合不正确”。我的回答是幂律数据拟合是正确的,我指出作者可能做错的是应用缩放,这使得幂律缩放看起来不正确。因此,我回答了这个问题。我的意思是作者从其他人那里得到了答案,所以一定要删除它。
【解决方案2】:

如果没有初步猜测,您的数据很难拟合。因此,在将 curve_fit 称为 p0 时添加提供猜测(如文档中所述):

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

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

x = np.array([88.08064516, 264.24193548, 440.40322581, 616.56451613, 792.72580645, 968.88709677, 1145.0483871, 1321.20967742, 1497.37096774, 1673.53225806, 1849.69354839, 2025.85483871, 2202.01612903, 2378.17741935, 2554.33870968, 2730.5, 2906.66129032,3082.82258065, 3258.98387097, 3435.14516129, 3611.30645161, 3787.46774194, 3963.62903226, 4139.79032258, 4315.9516129, 4492.11290323, 4668.27419355, 4844.43548387, 5020.59677419, 5196.75806452, 5372.91935484, 5549.08064516])
y = np.array([210737, 2175, 514, 158, 90, 46, 27, 22, 10, 11, 3, 7, 3, 2, 0, 1, 1, 1, 0, 0, 1, 0, 0,0, 0, 0, 1, 0, 0, 0, 0,1])

p0 = [20000,0.003,1]
popt, pcov = curve_fit(func, x, y, p0=p0)

fig, ax = plt.subplots()
ax.plot(x, func(x, *popt), 'g--', label = 'fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
ax.plot(x,  y, 'ro', label='data')

输出:

请注意,这有助于您解决问题。您仍然必须实现对数轴。对于日志。适合我建议修剪数据,例如,

y_pruned = np.where(y<1, 1, y)
popt, pcov = curve_fit(func, x, np.log(y_pruned), p0=p0)
ax.plot(x, func(x, *popt), 'g--', label = 'fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
ax.plot(x,  np.log(y_pruned), 'ro', label='data')

这会产生:

【讨论】:

  • 完美,我试过了,输出是想要的。谢谢!
  • @AlexandraBradan 但请注意,原始数据可能是双指数。初始衰减非常快(您可以通过查看图中的 1/e 并查看 x 和 y 是否匹配来快速检查这一点)。
  • @JamesPhillips ac 可以通过视觉检查从 y 值中提取。 b我通过查看图表的第一个点和第二个点之间的急剧衰减已经猜到了。
  • @JamesPhillips 谢谢,顺便看了你出色的在线曲线拟合工具,继续努力!
  • @AlexandraBradan 不要过度拟合这个问题,如果数据来自一个没有两个不同衰减率的过程,那么选择单指数。
【解决方案3】:

我可以用数据做的最好的是通过对数缩放两组数据值,然后尝试使用该功能。我包括以下代码和图表。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

fig, ax = plt.subplots()
x = np.array([88.08064516, 264.24193548, 440.40322581, 616.56451613, 792.72580645, 968.88709677, 1145.0483871, 1321.20967742, 1497.37096774, 1673.53225806, 1849.69354839, 2025.85483871, 2202.01612903, 2378.17741935, 2554.33870968, 2730.5, 2906.66129032, 3082.82258065, 3258.98387097, 3435.14516129, 3611.30645161, 3787.46774194, 3963.62903226, 4139.79032258, 4315.9516129, 4492.11290323, 4668.27419355, 4844.43548387, 5020.59677419, 5196.75806452, 5372.91935484, 5549.08064516])
y = np.array([210737, 2175, 514, 158, 90, 46, 27, 22, 10, 11, 3, 7, 3, 2, 0, 1, 1, 1, 0, 0, 1, 0, 0,0, 0, 0, 1, 0, 0, 0, 0,1])

x = np.log(x)
y = np.log(y + 1) # Need to add something to make log work

popt, pcov = curve_fit(func, x, y)
ax.plot(x, func(x, *popt), 'g--')
ax.plot(x,  y, 'ro', label='data')

plt.show()

【讨论】:

  • 谢谢你的解决方案! span>
猜你喜欢
  • 2021-06-04
  • 1970-01-01
  • 1970-01-01
  • 2016-06-24
  • 2019-06-14
  • 1970-01-01
  • 2017-04-07
  • 2020-05-15
  • 1970-01-01
相关资源
最近更新 更多