More than 5 years have passed since last update.
はじめに
実験結果の処理などで計算をしていると，誤差を考慮しなければならないことが多々あります。
誤差を含む値を使用して計算した関数値の誤差は
$$e=\sqrt{\sum_i^ne_i^2\left(\frac{\partial f}{\partial x_i}\right)^2}$$
となります(ここでこの式の導出をすることはしません)。
ご覧の通り，これを計算するのは面倒です。
割り算でも偏微分が面倒なのに，指数関数とか対数関数が出てきたらとんでもないことになります。
そこで，「退屈なことをPythonにやらせる」ことにしました。
実行環境
- Windows 10
- Python 3.8.3
f-stringを利用しているので3.6以降で動作します。
ライブラリ
一応numpyを使っていますが，実は必要ないです。
numpy.sqrtしか使っていませんが，これはx**(1/2)で代替可能です。
なぜわざわざimportしたかと言えば，科学計算でよく出てくる指数関数・対数関数を使えるようにするためです。
copyは標準で入っていると思います。
ソースコード
import numpy
from copy import copy
from math import inf
def calc_eps():
a = 4 / 3
b = a - 1
c = b + b + b
return 1 - c
EPS = calc_eps()
def eval_func(func, beta):
if callable(func):
return func(beta)
for i, b in enumerate(beta):
func = func.replace(f'${{{i+1}}}', str(b))
return eval(func)
def diff(func, beta, index):
last_d = 0.0
d = inf
err = inf
last_e = inf
eps = 1
step = 1.1 # Step for eps
v = eval_func(func, beta)
while abs(last_d - d) > EPS:
if eps <= EPS:
break
beta2 = copy(beta)
beta2[index] += eps
tmp = eval_func(func, beta2) - v
if tmp == 0.0:
break
err = abs((last_d - d) / (step - 1)) # Estimate error
if err > last_e:
break
last_e = err
last_d = d
d = tmp / eps
eps /= step
return last_d
def parse_data(data, err=None):
beta = [x[0] for x in data] if err is None else data
err = [x[-1] for x in data] if err is None else err
return beta, err
def calc_error(func, beta, err):
E = 0.0
for i, e in enumerate(err):
E += (diff(func, beta, i)**2) * (e**2)
return numpy.sqrt(E)
def calc_with_error(func, data, err=None):
beta, err = parse_data(data, err)
return eval_func(func, beta), calc_error(func, beta, err)
def get_digit(value):
return int(numpy.log10(value) // 1)
def separate_exp(value):
d = get_digit(value)
m = value / (10**d)
return m, d
def roundat(val, digit):
d = get_digit(val)
sgn = val / numpy.absolute(val)
coeff = 10 ** (d - digit)
val /= coeff
val = int(val + 0.5 * sgn) * coeff
return val
def format_value(val, err, show_sign=False):
if val < 0:
val = -val
sgn = '-'
else:
sgn = '' if not show_sign else '+'
s_err, d_err = separate_exp(err)
s_err = int(round(s_err))
if s_err >= 10:
s_err //= 10
d_err += 1
d_val = get_digit(val)
d = d_val - d_err
s_val = roundat(val, d)
m = f'{s_val:e}'.split('e')[0] + '0' * d
width = 1 if d <= 0 else 2 + d
m = m[:width]
err = str(s_err)
if not d_val >= d_err:
d_val = d_err
return f'{sgn}{m}({err})E{d_val}'
def main():
f = '${1} * ${2} / ${3}'
data = [(25.4, 0.5), (1.56, 0.01), (13.5, 0.2)]
value, error = calc_with_error(f, data)
print(f'Value: {value}')
print(f'Error: {error}')
print(format_value(value, error))
if __name__ == '__main__':
main()
解説
ダラダラと長いので畳んでおきます。
実行例
上記のコードで示した25.4(5)×1.56(1)/13.5(2)では
Value: 2.9351111111111114
Error: 0.07471981498926973
2.94(7)E0
と表示されます。
検証
実際に手計算で頑張ると，
$$e=\sqrt{\left(\frac{1.56}{13.5}\right)^2\times0.5^2+\left(\frac{25.4}{13.5}\right)^2\times0.01^2+\left(-\frac{25.4\times1.56}{13.5^2}\right)^2\times0.2^2}=0.0747\cdots$$
となり，Pythonでも正確に計算できていることが確認できました。
最後に
ゴリゴリと力業でなんとかしている部分が多くなってしまったので，よりよい書き方等がありましたら教えていただければ幸いです。
Register as a new user and use Qiita more conveniently
- You get articles that match your needs
- You can efficiently read back useful information
- You can use dark theme