python - Python 中 scipy.optimize.curve_fit 的数值精度

Question

我对 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(模型, xdata, ydata, p0=None, sigma=sig)

我的第一个想法是我没有使用 double ，但我很确定Python会自动创建 double 字。然后我想可能是文档的问题导致了数字被切断？不管怎样，我可以把我的代码放在这里，但有点复杂。有没有办法确保曲线拟合功能保存我的数字？

非常感谢您的帮助!

编辑:以下是我的代码:

# 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()

Answer 1

我相信这与此处使用的最小化算法以及可获得的最大精度有关。

我记得几年前在数字食谱中读到过它，我看看是否可以为您找到引用。

编辑:

链接到数字食谱here - 跳至第 394 页，然后阅读该章。注意404页第三段:

"Indulge us a final reminder that tol should generally be no smaller than the square root of your machine’s floating-point precision."

和mathematica提到如果你想要准确性，那么你需要采用不同的方法，并且他们实际上不会使用 LMA 除非问题被认为是平方和问题。

鉴于您只是进行一维拟合，尝试实现他们在该章中提到的拟合算法之一可能是一个很好的练习。

您实际上想要实现什么目标？据我了解，您本质上是在尝试计算出添加到曲线中的随机噪声量。但这并不是你真正在做的事情 - 除非我理解错误......

编辑2:

因此，在了解如何生成数据后，您发现数据和您所应用的模型存在问题。

你基本上是在拟合这个的两侧:

您实际上是在尝试将高斯的高度与随机数相匹配。您没有将高斯拟合到这些数字的频率。

看看你的代码，从你所说的来看，这不是你的最终目标，你只是想习惯优化方法？

如果你随机调整距太阳的距离，然后拟合数据，看看是否可以最小化找到生成数据集的距离，这会更有意义吗？