我在 python 中遇到了 scipy.optimize.curve_fit 函数的数值准确性问题。在我看来,当我想要 ~ 15 位数时,我只能获得 ~ 8 位数的准确度。我有一些从以下数据创建中获得的数据(此时是人工创建的):
其中术语 1 ~ 10^-3,术语 2 ~ 10^-6,术语 3 是 ~ 10^-11。在数据中,我A
随机变化(这是一个高斯错误)。然后我尝试将其拟合到模型中:
其中 lambda 是一个常数,我只适合alpha
(它是函数中的一个参数)。现在我希望看到 和 之间的线性关系alpha
,A
因为数据创建中的项 1 和 2 也在模型中,所以它们应该完全取消;
所以;
然而,对于小的A
(~10^-11 及以下)会发生什么,alpha
不会随 缩放A
,也就是说,随着A
变得越来越小,alpha
趋于平稳并保持不变。
作为参考,我称之为: op, pcov = scipy.optimize.curve_fit(model, xdata, ydata, p0=None, sigma=sig)
我的第一个想法是我没有使用双精度,但我很确定 python 会自动创建双精度数字。然后我认为这是文档的问题,可能会切断数字?无论如何,我可以把我的代码放在这里,但它有点复杂。有没有办法确保曲线拟合功能保存我的数字?
非常感谢你的帮助!
编辑:以下是我的代码:
# Import proper packages
import numpy as np
import numpy.random as npr
import scipy as sp
import scipy.constants as spc
import scipy.optimize as spo
from matplotlib import pyplot as plt
from numpy import ndarray as nda
from decimal import *
# Declare global variables
AU = 149597871000.0
test_lambda = 20*AU
M_Sun = (1.98855*(sp.power(10.0,30.0)))
M_Jupiter = (M_Sun/1047.3486)
test_jupiter_mass = M_Jupiter
test_sun_mass = M_Sun
rad_jup = 5.2*AU
ran = np.linspace(AU, 100*AU, num=100)
delta_a = np.power(10.0, -11.0)
chi_limit = 118.498
# Model acceleration of the spacecraft from the sun (with Yukawa term)
def model1(distance, A):
return (spc.G)*(M_Sun/(distance**2.0))*(1 +A*(np.exp(-distance/test_lambda))) + (spc.G)*(M_Jupiter*distance)/((distance**2.0 + rad_jup**2.0)**(3.0/2.0))
# Function that creates a data point for test 1
def data1(distance, dela):
return (spc.G)*(M_Sun/(distance**2.0) + (M_Jupiter*distance)/((distance**2.0 + rad_jup**2.0)**(3.0/2.0))) + dela
# Generates a list of 100 data sets varying by ~&a for test 1
def generate_data1():
data_list = []
for i in range(100):
acc_lst = []
for dist in ran:
x = data1(dist, npr.normal(0, delta_a))
acc_lst.append(x)
data_list.append(acc_lst)
return data_list
# Generates a list of standard deviations at each distance from the sun. Since &a is constant, the standard deviation of each point is constant
def generate_sig():
sig = []
for i in range(100):
sig.append(delta_a)
return sig
# Finds alpha for test 1, since we vary &a in test 1, we need to generate new data for each time we find alpha
def find_alpha1(data_list, sig):
alphas = []
for data in data_list:
op, pcov = spo.curve_fit(model1, ran, data, p0=None, sigma=sig)
alphas.append(op[0])
return alphas
# Tests the dependence of alpha on &a and plots the dependence
def test1():
global delta_a
global test_lambda
test_lambda = 20*AU
delta_a = 10.0**-20.0
alphas = []
delta_as = []
for i in range(20):
print i
data_list = generate_data1()
print np.array(data_list[0])
sig = generate_sig()
alpha = find_alpha1(data_list, sig)
delas = []
for alp in alpha:
if alp < 0:
x = 0
plt.loglog(delta_a, abs(alp), '.' 'r')
else:
x = 0
plt.loglog(delta_a, alp, '.' 'b')
delta_a *= 10
plt.xlabel('Delta A')
plt.ylabel('Alpha (at Lambda = 5 AU)')
plt.show()
def main():
test1()
if __name__ == '__main__':
main()