標準誤差に誤差伝搬式を適用する

標準誤差に誤差伝搬の式を適用

誤差伝搬式と呼ばれるものがあります。下記に誤差伝搬の式の説明を記載します。
尚、COV(a,b)は共分散なのでa,bの誤差が独立ならばゼロになります。

y(a,b)とするとyの誤差σ_{y}は,aの誤差σ_{a}及びbの誤差σ_{b}を使用して下記となる.

σ_{y}^2 = (\frac{∂y}{∂a})^2 σ_{a}^2 +  (\frac{∂y}{∂b})^2 σ_{b}^2+ 2COV(a,b)\frac{∂y}{∂a}\frac{∂y}{∂b}

今回は標準誤差に誤差伝搬の式を適用できるか考えてみました。

##　標準誤差とは？
そもそも標準誤差とは何か、統計学の教科書的な説明を下記に書きます。

ある標本データから、母分布のパラメータを推定するとします。
(例として、母分布が正規分布であれば、パラメータは平均と標準偏差です。)
ここでさらに、パラメータ推定時に導出される標準誤差を使って、信頼区間というものも考えます。

ちなみに標準誤差の導出方法に関して詳細は割愛しますが、
最尤法でパラメータ推定した場合は対数尤度関数を色々といじくると(微分して、期待値をとって～等)
標準誤差を導出することができます。

95%信頼区間(=推定したパラメータ+/-2×標準誤差)を考えると
上記の95%信頼区間の推定を、母分布からサンプリングされた新しい標本データで100回繰り返します。

100回試行したら,95回は95%信頼区間の中に母分布のパラメータ真値を含み、5回は含まない。
標準誤差とは上記のような意味合いで使用します。
ざっくりいうと、推定したパラメータがどれくらい信頼できるかの指標のようなものです。

似た言葉として標準偏差σがありますが、これは実際にデータが取り合える分布のことをさします。

標準誤差と標準偏差の違い

標準誤差はデータがどれくらい信頼できるかの指標で、標準偏差はデータのとりえる範囲を表します。
ここで最初の表題に戻ります。
誤差伝搬の式は標準偏差(というか誤差)に対してはよく使っていましたが、今回疑問に思ったことは、
「標準誤差にも適用できるのだろうか？」ということです。
そのまま適用しちゃって問題ないとは思っていたのですが、念には念を入れて数値シミュレーションしてみました。

数値シミュレーションしてみた

流れとしては、
①下記に従うデータを発生させて最尤法によりパラメータa,b,cを推定します。
　ここでεは正規分布型のランダムノイズになります。

y = a+bx^2+cx^{-2} + ε

②b,cが推定できたら,d=(c/b)^(1/4)というパラメータを計算します。
③dの標準誤差をb,cの標準誤差から誤差伝搬の式で計算します。
　
この①②③を複数回繰り返し、②で計算したdのバラつきがdの標準誤差と考えます。
これが③の計算値と合っていれば、標準誤差にも誤差伝搬の式は適用できます。
③の計算も複数実施すので平均値を採用しました。

下記に実際のpythonコードと出力結果です
おおよそ②のバラつきと③の計算結果が合っているかと思います。

main.py

import numpy as np
from numpy.random import *
from scipy.optimize import curve_fit
import statistics
import math


def create_toydata(x_list):
    """
    x_listを引数として,誤差付きのy_listを作成する
    :param x_list: xのリスト
    :return y_list: xの関数であるyのリスト
    :return true_d: 真のd
    """

    # パラメータ真値 数字は適当
    a = 2e2
    b = 2e-5
    c = 2e7

    y_list = [None] * len(x_list)

    for i, x in enumerate(x_list):
        # データに正規分布の誤差を加える。2という数字は適当に決めた
        noize = 2 * randn()
        y_list[i] = func(x, a, b, c) + noize

    # c,bの関数であるdを計算する
    true_d = (c/b)**(1/4)

    return y_list, true_d


def func(x, a, b, c):
    """fittingする関数を定義する。関数自体は適当に決めた"""
    return a + b * x**2 + c * x**(-2)


def cal_se(C, B, C_se, B_se, COV):
    """
    相関がある場合の誤差伝搬則から(C/B)^(1/4)の誤差を求める
    :param C: パラメータC
    :param B: パラメータB
    :param C_se: Cの標準誤差
    :param B_se: Bの標準誤差
    :param COV: C,Bの共分散
    :return: (C/B)^(1/4)の標準誤差
    """
    del_c = (1/4) * C**(-3/4) * B**(-1/4)
    del_b = (-1/4) * C**(1/4) * B**(-5/4)

    return math.sqrt((del_c* C_se)**2 + (del_b * B_se)**2 + 2 * COV * del_c * del_b)


if __name__ == '__main__':

    # 試行回数
    n = 1000

    a_list = [None] * len(range(n))
    b_list = [None] * len(range(n))
    c_list = [None] * len(range(n))
    d_list = [None] * len(range(n))
    d_se_list = [None] * len(range(n))

    # cnt:95%信頼区間内に真値が入った数
    cnt = 0

    # n回下記試行を繰り返す
    for i in range(n):
        # toydataから、最尤法で「y= a + x^2 + c^-2」の係数a,b,cを求める
        # b,cの関数であるdを計算する。
        # b,cの標準誤差からdの標準誤差を計算する。

        # 数字は適当
        x_st = 500
        x_end = 2500
        x_num = 200

        x_list = np.linspace(x_st, x_end, x_num).tolist()
        y_list, true_d = create_toydata(x_list=x_list)

        popt, pcov = curve_fit(func, x_list, y_list)

        # # dを計算する
        d = (popt[2]/popt[1])**(1/4)
        d_list[i] = d

        # dの標準誤差を,b,cの標準誤差から計算する
        d_se = cal_se(popt[2], popt[1], math.sqrt(pcov[2][2]), math.sqrt(pcov[1][1]), pcov[1][2])
        d_se_list[i] = d_se

        # 95%信頼区間を求める
        d_low, d_up = d - 2 * d_se, d + 2 * d_se

        # 95%信頼区間内に真値が入った回数を調べる
        if d_low < true_d and true_d < d_up:
            cnt += 1


    # n回試行したdの標準偏差を求める(これがdの標準誤差とみなす)
    true_d_se = statistics.pstdev(d_list)
    print("dの標準誤差:", true_d_se)

    # 試行回数毎にdの標準誤差を計算するので、最終的な結果は平均とする
    cal_d_se = statistics.mean(d_se_list)
    print("計算したdの標準誤差の平均:", cal_d_se)

    # 95%信頼区間内に標準誤差が入る割合を調べる.95%信頼区間の定義から0.95になるはず
    print("95%信頼区間内に標準誤差が入った割合:", cnt/n)

下記が実行結果です。1行目が②のバラつきで2行目が③の計算結果です。

dの標準誤差: 2.2655785060979126
計算したdの標準誤差の平均: 2.2443534560700673
95%信頼区間内に標準誤差が入った割合: 0.947