誤差伝搬

このページではGitHubにあげたcalculate_error_propagation.ipynbの中身をマークダウン形式で書き下します。

情報

以下ではipynbでのoutputをこの中に表示します。

ヒント

このようにした理由としてはGitHubのviwerの共同が安定していないためです。 これまではみれていたソースが急に見れなくなったりしているようです。

代数計算で誤差伝搬を計算する。

BasicAlgebraicCalculationClassを継承した自作代数計算クラスErrorPropagationClass使って、誤差伝搬の代数計算を行ったのちに誤差を見積もってみます。

また描画するためのライブラリを使って計算結果の可視化を行ってみます。

BasicAlgebraicCalculationClassについては別の場所で説明したので省略します。

まずは必要なライブラリをimportしてきます。

from fumi_python_libraries import ErrorPropagationClass
import matplotlib.pyplot as plt
import numpy
import sympy
使い方

まず誤差伝搬をしたい関数と変数を定義します。

この記法はBasicAlgebraicCalculationClassと全く一緒です。

err_pro = ErrorPropagationClass()
var = err_pro.create_variables( ["\\alpha","\\beta"] )
fun = sympy.sqrt( var[0] / var[1] )
err_pro.set_left_hand_side( 'f(\\alpha, \\beta)' )
err_pro.set_right_hand_side( fun )
err_pro.set_variable_symbols()
err_pro.print_formula()
情報
$$ f(\alpha, \beta) = \sqrt{\frac{\alpha}{\beta}} $$

set_variable_symbols()は引数で変数のリストを渡しますが、何も与えない場合は直近でcreate_variables()で宣言した変数が登録されます。

partial_derivative( bool )で登録した代数を変数として、それぞれの誤差伝搬を計算します。

引数は結果を出力するかを選択できます。

loop処理の中でこのような計算をする場合、出力したくないので選べるようになっています。

err_pro.partial_derivative( True )
情報
$$ {\delta}f=\sqrt{\sum_{i=0}^{n} \left( \frac{{\partial} f}{{\partial} x_i} {\delta}x_i \right)^2} $$ $$ x_0 = \alpha $$ $$ x_1 = \beta $$ $$ \frac{\partial f(\alpha, \beta) }{ \partial \alpha } = \frac{\sqrt{\frac{\alpha}{\beta}}}{2 \alpha} $$ $$ \frac{\partial f(\alpha, \beta) }{ \partial \beta } = - \frac{\sqrt{\frac{\alpha}{\beta}}}{2 \beta} $$

今回の場合$\alpha$ , $\beta$ が測定量です。

例えば以下のような精度で測定できたとしましょう。

$$ \alpha = 17120 \pm \sqrt{17120} $$ $$ \beta = 18571 \pm \sqrt{18571} $$

測定量の設定と誤差の設定はそれぞれset_value_series( value list )とset_error_series( error list )で渡すことができます。

そして、calculate_erorr( bool, round digit )で誤差伝搬を計算し、丸めて出力してくれます。

val = [17120., 18571.]
err = [numpy.sqrt(x) for x in val]
err_pro.set_value_series( val )
err_pro.set_error_series( err )
result = err_pro.calculate_erorr( True, 6 )
情報
$$ f(\alpha, \beta) \pm \delta f(\alpha, \beta) = 0.960139 \pm 0.005086 $$
応用例
Woods–Saxon Potentialとその微分を計算し、誤差つきで描画してみる。

少し高級な例を計算してみます。

まずWoods–Saxon Potentialを定義してみましょう。

err_pro = ErrorPropagationClass()
var = err_pro.create_variables( ["r","R_i","a_i"] )
fun = 1 / ( 1 + sympy.exp( ( var[0] - var[1] ) / var[2] ) )
err_pro.set_left_hand_side( 'f(r, R_i, a_i)' )
err_pro.set_right_hand_side( fun )
err_pro.set_variable_symbols()
err_pro.print_formula()
情報
$$ f(r, R_i, a_i) = \frac{1}{e^{\frac{- R_{i} + r}{a_{i}}} + 1} $$

標準的な原子核の教科書では分母の順番は逆になっている気がしますが、これはsympyが勝手に並べ替えています。

並べ直さない方法があると助かるのですが…

そして誤差伝搬を計算します。

ここでは$r_i$ , $R_i$, $a_i$ 全ての要素の誤差伝搬を考えます。

結果は以下のとおりです。

dum = err_pro.partial_derivative( True )
情報
$$ {\delta}f=\sqrt{\sum_{i=0}^{n} \left( \frac{{\partial} f}{{\partial} x_i} {\delta}x_i \right)^2} $$ $$ x_0 = r $$ $$ x_1 = R_i $$ $$ x_2 = a_i $$ $$ \frac{\partial f(r, R_i, a_i) }{ \partial r } = - \frac{e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} $$ $$ \frac{\partial f(r, R_i, a_i) }{ \partial R_i } = \frac{e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} $$ $$ \frac{\partial f(r, R_i, a_i) }{ \partial a_i } = \frac{\left(- R_{i} + r\right) e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i}^{2} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} $$

さて次にWoods Saxonの微分を計算してみます。

dif_err_pro = err_pro.copy_class()
func = dif_err_pro.differentiate( 'r' , True )
dif_err_pro.print_formula()
情報
$$ \frac{d}{dr}f(r, R_i, a_i) = - \frac{e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} $$

その時の誤差伝搬は以下のとおりです。

dum = dif_err_pro.partial_derivative( True )
情報
$$ {\delta}f=\sqrt{\sum_{i=0}^{n} \left( \frac{{\partial} f}{{\partial} x_i} {\delta}x_i \right)^2} $$ $$ x_0 = r $$ $$ x_1 = R_i $$ $$ x_2 = a_i $$ $$ \frac{\partial \frac{d}{dr}f(r, R_i, a_i) }{ \partial r } = - \frac{e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i}^{2} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} + \frac{2 e^{\frac{2 \left(- R_{i} + r\right)}{a_{i}}}}{a_{i}^{2} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{3}} $$ $$ \frac{\partial \frac{d}{dr}f(r, R_i, a_i) }{ \partial R_i } = \frac{e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i}^{2} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} - \frac{2 e^{\frac{2 \left(- R_{i} + r\right)}{a_{i}}}}{a_{i}^{2} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{3}} $$ $$ \frac{\partial \frac{d}{dr}f(r, R_i, a_i) }{ \partial a_i } = \frac{e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i}^{2} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} + \frac{\left(- R_{i} + r\right) e^{\frac{- R_{i} + r}{a_{i}}}}{a_{i}^{3} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{2}} - \frac{2 \left(- R_{i} + r\right) e^{\frac{2 \left(- R_{i} + r\right)}{a_{i}}}}{a_{i}^{3} \left(e^{\frac{- R_{i} + r}{a_{i}}} + 1\right)^{3}} $$
描画してみる
def get_function_value_lists(num=50):
    ws_r0 = 1.15 # fm
    ws_M = 208.
    ws_R = ws_r0 * ws_M**(1./3.)
    ws_a = 0.54
    
    fac_calc = ErrorPropagationClass()
    var = fac_calc.create_variables( ["r"] )
    fun = 1 / ( 1+ sympy.exp( ( var[0] - ws_R ) / ws_a ) )
    fac_calc.set_formula( 'f', fun )
    factor = fac_calc.definite_integral('r', 0, 30, True)
    
    ws_x  = []
    ws_y  = []
    ws_ye = []
    dws_y  = []
    dws_ye = []

    err_pro = ErrorPropagationClass()
    var = err_pro.create_variables( ["r","R_i","a_i"] )
    fun = 1 / ( 1 + sympy.exp( ( var[0] - var[1] ) / var[2] ) )
    dum = err_pro.set_left_hand_side( 'f(r, R_i, a_i)' )
    err_pro.set_right_hand_side( fun )
    err_pro.set_variable_symbols()
    err_pro.partial_derivative()
    
    dif_err_pro = err_pro.copy_class()
    dif_err_pro.differentiate('r', True)

    dif_err_pro.partial_derivative()
    
    for i in range(num):
        r_val = 15. * i / num
        ws_x.append( r_val )
        
        val = [r_val, ws_R, ws_a ]
        err = [0, 0.05*ws_M**(1./3.), 0.05]

        err_pro.set_value_series( val )
        err_pro.set_error_series( err )
        result = err_pro.calculate_erorr()
        
        ws_ye.append( result[0] / float( factor ) )
        ws_y.append( result[1] / float( factor ) )
        
        dif_err_pro.set_value_series( val )
        dif_err_pro.set_error_series( err )
        dif_result = dif_err_pro.calculate_erorr()
        dws_ye.append( dif_result[0] / float( factor ) )
        dws_y.append( dif_result[1] / float( factor ) )
        
        
    x = numpy.array(ws_x)
    y1 = numpy.array(ws_y)
    y1e = numpy.array(ws_ye)
    y2 = numpy.array(dws_y)
    y2e = numpy.array(dws_ye)

    return x, y1, y1e, y2, y2e
x, y1, y1e, y2, y2e = get_function_value_lists(100)

plt.figure(figsize=(10, 6))

plt.fill_between( x, y1 - y1e, y1 + y1e, color='red', edgecolor='gray', alpha=0.1 )
plt.fill_between( x, y2 - y2e, y2 + y2e, color='blue', edgecolor='gray', alpha=0.1 )
plt.plot(x, y1, color='red', label='Woods–Saxon')
plt.plot(x, y2, color='blue', label='Differentiation of Woods–Saxon')

plt.xlabel('Radius [fm]')
plt.ylabel('Value')
plt.title('Plot with Error Hatch')
plt.legend()
plt.show()
情報

fig1 fig1