1

我在 python 中遇到了 scipy.optimize.curve_fit 函数的数值准确性问题。在我看来,当我想要 ~ 15 位数时,我只能获得 ~ 8 位数的准确度。我有一些从以下数据创建中获得的数据(此时是人工创建的):

在此处输入图像描述

其中术语 1 ~ 10^-3,术语 2 ~ 10^-6,术语 3 是 ~ 10^-11。在数据中,我A随机变化(这是一个高斯错误)。然后我尝试将其拟合到模型中:

在此处输入图像描述

其中 lambda 是一个常数,我只适合alpha(它是函数中的一个参数)。现在我希望看到 和 之间的线性关系alphaA因为数据创建中的项 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()
4

1 回答 1

2

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

我记得几年前在数字食谱中读过它,我会看看我是否可以为你挖掘参考。

编辑:

此处链接到数字食谱- 跳至第 394 页,然后阅读该章节。请注意第 404 页的第三段:

“请给我们一个最后的提醒,tol通常不应小于您机器浮点精度的平方根。”

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

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

你到底想达到什么目的?据我了解,您实际上是在尝试计算添加到曲线中的随机噪声量。但这并不是你真正在做的事情——除非我理解错了……

编辑2:

因此,在阅读了您如何生成数据之后,您正在应用的数据和模型存在问题。

你基本上适合这个的两个方面:

在此处输入图像描述

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

查看您的代码,根据您所说的判断,这不是您的最终目标,您只是想习惯优化方法吗?

如果您随机调整与太阳的距离,然后拟合数据并查看是否可以最小化以找到生成数据集的距离,那会更有意义?

于 2013-07-08T21:23:24.547 回答