1

Python の scipy.optimize.curve_fit 関数の数値精度に問題があります。〜15桁を希望する場合、〜8桁の精度しか得られないように思えます。次のデータ作成から作成されたデータ(この時点では人工的に作成された)があります。

ここに画像の説明を入力

ここで、項 1 ~ 10^-3、項 2 ~ 10^-6、および項 3 は ~ 10^-11 です。データでは、私はAランダムに変化します (ガウス誤差です)。次に、これをモデルに適合させようとします。

ここに画像の説明を入力

ここで、ラムダは定数であり、私は適合alphaするだけです (関数内のパラメーターです)。alphaデータ作成の項 1 と 2 もモデルに含まれているため、との間に線形関係が見られると予想されるAので、それらは完全に相殺されるはずです。

ここに画像の説明を入力

そう;

ここに画像の説明を入力

ただし、小さい場合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