Числовая точность с scipy.optimize.curve_fit в Python

Question

У меня возникают проблемы с численной точностью функции scipy.optimize.curve_fit в python. Мне кажется, что я могу получить ~ 8 цифр точности, когда я желаю ~ 15 цифр. У меня есть некоторые данные (в этой точке искусственно созданной) изготовлены из следующих данных: создание

, где термин 1 ~ 10^-3, срок 2 ~ 10^-6, а срок 3 ~ 10^-11. В данных я меняю A случайно (это ошибка Гаусса). Затем я стараюсь, чтобы соответствовать этим к модели:

где лямбда является постоянной, и я только подхожу alpha (он является параметром в функции). Теперь я ожидаю увидеть линейную связь между alpha и A, потому что термины 1 и 2 при создании данных также находятся в модели, поэтому они должны отлично отменить;

Итак,

Однако то, что происходит, для малых A (~ 10^-11 и ниже), alpha не масштабируется с A, то есть, как и A становится все меньше и меньше, alpha уровней, и остается постоянным.

Для справки, я вызываю следующее: оп, pcov = scipy.optimize.curve_fit (модель, XData, ydata, p0 = None, сигма = сиг)

Моя первая мысль была, что я не использовал двойная точность, но я уверен, что python автоматически создает числа в двойной точности. Тогда я подумал, что это проблема с документацией, возможно, которая отсекает цифры? В любом случае, я мог бы поставить свой код здесь, но это довольно сложно. Есть ли способ гарантировать, что функция подгонки кривой сохранит мои цифры?

Большое вам спасибо за помощь!

EDIT: Ниже мой код:

# 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:

«Побалуйте нам NAL напоминание о том, что фи tol обычно должно быть не меньше, чем квадратный корень из плавающей точкой точности фл вашей машины."

И mathematica оговориться, что если вы хотите точности, то вам нужно перейти на другой метод, и что они не Infact использовать LMA, если проблема не будет признана как сумма проблемы квадратов.

Учитывая, что вы просто выполняете одномерную подгонку, это может быть хорошим упражнением, чтобы попытаться реализовать один из подходящих алгоритмов, упомянутых в этой главе.

Что вы на самом деле пытаетесь достичь? я понимаю об этом, вы, по сути, пытаетесь определить количество случайных шумов, которые вы добавили к t он криво. Но это не совсем то, что вы делаете, - если я не понял неправильно ...

Edit2:

Итак, после прочтения, как вы генерировать данные, есть проблема с данными и моделью вас» повторно применяя.

Вы по существу подгонки две стороны этого:

Вы по существу пытаетесь соответствовать высоты гауссовских к случайным числам. Вы не подходите гауссову к частоте этих чисел.

Посмотрите на свой код и, судя по тому, что вы сказали, это не конечная цель, и вы просто хотите привыкнуть к методу оптимизации?

Было бы разумнее, если бы вы случайно отрегулировали расстояние от солнца, а затем подходите к данным и видите, можете ли вы свести к минимуму, чтобы найти расстояние, которое сгенерировало набор данных?