python - 使用最小二乘法(Levenberg-Marquardt 算法)将双曲线拟合到线性数据的程序未按预期工作
问题描述
我有一个一维数组数据,我试图使用三个参数将其建模为双曲线。我正在尝试使用 scipy.optimize 库中的 leastsq 函数来实现 Levenberg Marquardt 算法。但是,我的程序卡在一个数字被零除的迭代中,我不明白为什么。
一些背景知识:一维数组数据基本上是不同盒子大小的空隙值。我已经从一些声音文件中生成了漏洞数据,可以在此处找到其上下文。
在算法中,最小二乘函数采用三个输入:
(a)三个参数的初始猜测
(b)最小二乘问题的 x 坐标 - 在我的问题中,这基本上是一个从 1 到 100 的整数的一维数组
(c)最小二乘问题的 y 坐标 - 这是存储空隙值的一维数组。因此,腔隙值是 x 的函数,其中 x 在 1 到 100 之间变化。
双曲线是使用三个参数 a、b 和 c 建模的
代码给出以下错误:
“溢出错误:无法将浮点无穷大转换为整数”
编码:
#import
from scipy import *
from scipy.optimize import leastsq
import matplotlib.pylab as plt
import numpy as np
import codecs, json
from math import *
# Define your function to calculate the residuals.
#The fitting function holds your parameter values.
def residuals(p, y, x):
err = y-pval(x,p)
return err
def pval(x, p):
z = x
for i in range(100):
print(x)
print(x[i]**p[1])
z[i] = p[0]/(x[i]**p[1])+p[2]
return z
#read in your data
obj_text = codecs.open('textfiles\CC1.json', 'r', encoding='utf-8').read()
b_new = json.loads(obj_text)
data = np.array(b_new)
x = np.arange(1,101)
y = data[1:101]
#guess at initial parameters
A1_0=1.0
A2_0=1.0
A3_0=0.5
#leastsq package calls the Levenberg-Marquardt algorithm
pname = (['A1','A2','A3'])
p0 = array([A1_0 , A2_0, A3_0])
plsq = leastsq(residuals, p0, args=(y, x), maxfev=2000)
# Now, plot your data
plt.plot(x,y,'xo',x,pval(x,plsq[0]),'x')
title('Least-squares fit to data')
xlabel('x')
ylabel('y')
legend(['Data', 'Fit'],loc=4)
# Your best-fit paramters are kept within plsq[0].
print(plsq[0])
根据错误,x 的值在迭代中的某个点变为 0,并且第一个参数 a 最终被零除,从而产生错误。
为了排除故障,我在执行代码时打印了值 x[i]^b 和数组 x,您可以在此处查看这些值。我看到数组 x 正在被修改,这不应该发生。x 应该保持从 1 到 100 的一维自然数数组,并且不会在迭代中被修改。我无法确定修改数组 x 的代码到底在哪里。
我希望数组 x 保持不变,并且打印参数 a、b 和 c 的最后三个值的代码。
编辑:我对我的代码进行了一些更改,之后它成功运行。以下是任何人感兴趣的编辑:
没有将 z 定义为 z = x,而只是将其定义为 z = np.arange(1,101)。结果是数组 x 不再改变,这是预期的。
将数组 x 和 y 的数据类型更改为使用浮点数
x = np.array(x, dtype=np.float64)
- 我又一次陷入了绘制数据的代码中。我得到了错误“'标题'未定义。xlabel,ylabel的类似错误。所以我只是删除了这些行并坚持使用
plt.plot(x,y,'red',x,pval(x,plsq[0]),'blue')
plt.show()
解决方案
不是对您的问题的直接回答,但由于您使用的是指数 ( **
),因此我强烈建议您Decimal
预先将所有数字转换为,以避免在大值上的浮点运算中固有的精度损失。
例如:
import decimal
decimal.getcontext().prec = 100
A1_0=Decimal("1.0")
A2_0=Decimal("1.0")
A3_0=Decimal("0.5")
x = [Decimal(f) for f in x]
y = [Decimal(f) for f in y]
也许你的零会“变成”一个接近零的小值......
推荐阅读
- javascript - React 422 中的 Axios.post,错误:(数据验证失败。)
- python - 熊猫根据“邻居”删除行
- java - 如何仅从我感兴趣的 OpenCV 中提取特征并为这些特征创建 .so?
- npm - 使用 np 包发布到 npm 或 yarn
- javascript - Karma 不运行 Jasmine 实现回调
- php - 正则表达式查找所有“?” 没有“ ”
- android - 如何在 RecyclerView 中设计播放列表项,如 youtube kids 应用程序
- javascript - 从输入文本中删除的只是数字的值与我的简单数组不匹配
- python - 如何根据另一个小部件的位置放置小部件?
- javascript - 将所有 JS 代码和变量视为执行上下文的属性是否正确?