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