誤差伝搬の法則をシミュレーションで理解しよう

はじめに

測定などを行う際、その測定結果には多かれ少なかれ誤差が含まれます。実験などをしていると、そのような誤差が含まれた測定結果を複数用いて計算した結果には、どの程度の誤差になるのかを評価したい場合があります。その時に用いられるのが誤差伝搬の法則です。（イメージが湧きにくいと思うので例を挙げると、力学の実験で力$F$を直接測れない場合、質量$m$と加速度$a$を測定して、運動方程式$ma=F$より、$F$を間接的に求めることがあります。しかし、測定には誤差（$\sigma_\ast$）が含まれているので、$m±\sigma_m$と$a±\sigma_a$となります。もちろん、$F$にもこの誤差が伝搬し$F±\sigma_F$となり、この$\sigma_F$を求めるのに誤差伝搬の法則が使われます。）
また、一般に精密な計測器には、精度が仕様として標準偏差を用いて記載されています。（このような計測器を用いた測定の場合、同じものを複数回測ると同じ結果が得られることはまれで、多少のばらつきが生じます。このばらつき度合いは標準偏差で表すことができ、計測器の精度となります。）

本記事では、誤差伝搬の法則をシミュレーションを通してイメージを掴んでもらい、数式でさらに理解を深めることを目標としています。（ちなみに著者は、徹夜で書いた実験レポートの誤差伝搬の計算をミスして、再レポートになった苦い思い出があります。誤差の取り扱いには気をつけましょう。）

誤差伝搬の詳細な導出についても記事を書いているのでよろしければ

もご覧いただけると幸いである。

誤差伝搬の法則

測定量を$x_i$、その標準偏差（精度）を$\sigma_{x_i}$とすると、$n$個の測定量を用いて計算される$y$は、

y = f(x_1, x_2, \cdots, x_n)

と表すことができます。このとき、$y$の誤差$\sigma_y$は誤差伝搬の法則より、

\begin{align}
\sigma_y &= \sqrt{\left(\frac{\partial f}{\partial x_1} \right)^2\sigma_{x_1}^2 + \left(\frac{\partial f}{\partial x_2} \right)^2\sigma_{x_2}^2 + \cdots + \left(\frac{\partial f}{\partial x_n} \right)^2\sigma_{x_n}^2}\\
&= \sqrt{\sum\limits_{i=1}^n \left(\frac{\partial f}{\partial x_i} \right)^2\sigma_{x_i}^2}
\end{align}

となります。注意点として、この法則は測定量$x_1, x_2, \cdots x_n$が互いに独立なときにのみ成り立ちます。本記事では扱いませんが、独立でない場合は共分散などを考慮する必要があります。

測定量が独立でない場合の誤差伝搬の法則

測定量を$x_i$、その標準偏差（精度）を$\sigma_{x_i}$とすると、$n$個の測定量を用いて計算される$y$は、

y = f(x_1, x_2, \cdots, x_n)

と表すことができます。各測定量が独立でない時、全ての測定量間の相関を$x_i, x_j$の共分散$\sigma_{x_i,x_j}$を用いて$y$の誤差$\sigma_y$を表すと、

\sigma_y = \sqrt{\sum\limits_{i=1}^n \left(\frac{\partial f}{\partial x_i} \right)^2\sigma_{x_i}^2 + \sum\limits_{i=1}^n \sum\limits_{j=1(j\neq i)}^n \frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\sigma_{x_i,x_j}}

となります。特に、全ての測定量が独立の場合は、共分散$\sigma_{x_i,x_j}$が0になり、先ほどの独立した誤差伝搬の法則と一致します。

実装コード

Google Colabで作成したコードは、こちらにあります。

各種インポート

各種インポート

import numpy as np
import matplotlib.pyplot as plt

実行時の各種バージョン

Name	Version
numpy	1.21.6
matplotlib	3.2.2

シミュレーション

誤差伝搬の法則は複数の測定量で成り立つが、簡単のために2つの測定量を仮定してシミュレーションを行います。

シミュレーションの流れです。

1. 信号1（$s_1$）、信号2（$s_2$）にそれぞれ正規分布に従う乱数を生成させます。（本記事では、$s_1$を$\mu=100, \sigma=10$、$s_2$を$\mu=200, \sigma=10$としました。）
2. それらの$s_1$, $s_2$を用いて任意の演算した信号（$s$）を求めます。
3. 1~2の処理を十分な回数行います。（本記事では100万回行いました。）
4. 十分な回数行ったら、$s$の平均と不偏標準偏差を求めます。
5. 4の結果と誤差伝搬の法則によって得られる誤差を比較します。

このシミュレーションにより、本来$s_1=100$、$s_2=200$と測定したいところ測定精度により$s_1=100±5$、$s_2=200±10$となる場合に、$s_1$, $s_2$を使って演算した$s$がどの程度のばらつくかを確認することができます。十分な回数行うので、本当に誤差伝搬の法則に従っているか見ていきましょう。下図右側がシミュレーションのイメージ図です。

信号1+信号2

以下手順です。
まず、独立した乱数により生成された100万個の$s_1$（$\mu_{s_1}=100, \sigma_{s_1}=5$の正規分布$N_1(100, 25)$）、$s_2$（$\mu_{s_2}=200, \sigma_{s_2}=10$の正規分布$N_2(200, 100)$）を使って、

s=s_1+s_2

を計算します。その分布$s$の平均（$mean\_s$）、不偏標準偏差（$sigma\_s$）を求めます。
次に、誤差伝搬の法則から導かれる平均（厳密には誤差伝搬の法則からではなく普通に理論的に導かれる）、誤差を求めます。なお、区別するために先頭に$law$をつけています。

law\_mean\_s = \mu_{s_1} + \mu_{s_2}

\begin{align}
law\_sigma\_s &= \sqrt{\left(\frac{\partial s}{\partial s_1} \right)^2\sigma_{s_1}^2 + \left(\frac{\partial s}{\partial s_2} \right)^2\sigma_{s_2}^2}\\
&= \sqrt{\sigma_{s_1}^2 + \sigma_{s_2}^2}
\end{align}

最後に、シミュレーションの結果（$mean\_s$, $sigma\_s$）と誤差伝搬の法則の結果（$law\_mean\_s$, $law\_sigma\_s$）を比較します。

信号1+信号2のコード

信号1+信号2

# s=s1+s2の場合
np.random.seed(2022)
n = 1000000

# 信号1(s1)
mu_s1 = 100
sigma_s1 = 5
s1 = np.random.normal(mu_s1, sigma_s1, n)

# 信号2(s2)
mu_s2 = 200
sigma_s2 = 10
s2 = np.random.normal(mu_s2, sigma_s2, n)

# 信号1+信号2
s = s1 + s2

# シミュレーションの結果の平均値と不偏標準偏差
mean_s = np.mean(s)
sigma_s = np.std(s, ddof=1)

# 誤差伝搬の法則に従う平均値と誤差
law_mean_s = mu_s1 + mu_s2
law_sigma_s = np.sqrt(sigma_s1**2 + sigma_s2**2)

print(f'シミュレーションの平均値: {mean_s:.3f}, 誤差伝搬の法則で求めた平均値: {law_mean_s:.3f}')
print(f'シミュレーションの不偏標準偏差: {sigma_s:.3f}, 誤差伝搬の法則で求めた誤差: {law_sigma_s:.3f}')

# グラフプロット
color_list = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(2, 2, figsize=(14, 14))
axs[0, 0].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[0, 0].axvline(mu_s1 - sigma_s1, c=color_list[0], ls='--', label='mu_s1±sigma_s1')
axs[0, 0].axvline(mu_s1 + sigma_s1, c=color_list[0], ls='--')
axs[0, 0].set_xlabel('value')
axs[0, 0].set_ylabel('counts')
axs[0, 0].legend()

axs[0, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[0, 1].axvline(mu_s2 - sigma_s2, c=color_list[1], ls='--', label='mu_s2±sigma_s2')
axs[0, 1].axvline(mu_s2 + sigma_s2, c=color_list[1], ls='--')
axs[0, 1].set_xlabel('value')
axs[0, 1].set_ylabel('counts')
axs[0, 1].legend()

axs[1, 0].hist(s, bins=20, alpha=0.6, color=color_list[2], label='s=s1+s2')
axs[1, 0].axvline(mean_s - sigma_s, c='r', ls='--', alpha=0.5, label='Simulation mean_s±sigma_s')
axs[1, 0].axvline(mean_s + sigma_s, c='r', ls='--', alpha=0.5)
axs[1, 0].axvline(law_mean_s - law_sigma_s, c='k', ls='--', alpha=0.5, label='Law law_mean_s±law_sigma_s')
axs[1, 0].axvline(law_mean_s + law_sigma_s, c='k', ls='--', alpha=0.5)
axs[1, 0].set_xlabel('value')
axs[1, 0].set_ylabel('counts')
axs[1, 0].legend()

axs[1, 1].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[1, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[1, 1].hist(s,  bins=20, alpha=0.6, color=color_list[2], label='s=s1+s2')
axs[1, 1].set_xlabel('value')
axs[1, 1].set_ylabel('counts')
axs[1, 1].legend(loc='upper center')

plt.show()

シミュレーションの結果

シミュレーションの平均値: 300.007, 誤差伝搬の法則で求めた平均値: 300.000
シミュレーションの不偏標準偏差: 11.171, 誤差伝搬の法則で求めた誤差: 11.180

シミュレーションの結果

上段の図は、それぞれの信号をヒストグラムで表したもので、正規分布に従っています。
下段左図の緑の分布が$s=s_1+s_2$のシミュレーション結果で、赤線がその結果から求めた平均±不偏標準偏差で、黒線が誤差伝搬の法則から導かれる平均±誤差になります。このようにほとんど一致していることが図から読み取れると思います。
また、下段右側はそれぞれのヒストグラムを並べて表示したものです。$s=s_1+s_2$の場合は、緑の分布が一番広がって見えることが図から読み取れると思います。

信号1-信号2

信号1-信号2もシミュレーションしてみます。基本的に設定は信号1+信号2と同じです。
信号($s$)を

s=s_1-s_2

とします。誤差伝搬の法則より、

law\_mean\_s = \mu_{s_1} - \mu_{s_2}

\begin{align}
law\_sigma\_s &= \sqrt{\left(\frac{\partial s}{\partial s_1} \right)^2\sigma_{s_1}^2 + \left(\frac{\partial s}{\partial s_2} \right)^2\sigma_{s_2}^2}\\
&= \sqrt{\sigma_{s_1}^2 + \sigma_{s_2}^2}
\end{align}

となります。誤差伝搬の法則と従っているかシミュレーションで確かめます。

信号1-信号2のコード

信号1-信号2

# s=s1-s2の場合
np.random.seed(2022)
n = 1000000

# 信号1(s1)
mu_s1 = 100
sigma_s1 = 5
s1 = np.random.normal(mu_s1, sigma_s1, n)

# 信号2(s2)
mu_s2 = 200
sigma_s2 = 10
s2 = np.random.normal(mu_s2, sigma_s2, n)

# 信号1-信号2
s = s1 - s2

# シミュレーションの結果の平均値と不偏標準偏差
mean_s = np.mean(s)
sigma_s = np.std(s, ddof=1)

# 誤差伝搬の法則に従う平均値と誤差
law_mean_s = mu_s1 - mu_s2
law_sigma_s = np.sqrt(sigma_s1**2 + sigma_s2**2)

print(f'シミュレーションの平均値: {mean_s:.3f}, 誤差伝搬の法則で求めた平均値: {law_mean_s:.3f}')
print(f'シミュレーションの不偏標準偏差: {sigma_s:.3f}, 誤差伝搬の法則で求めた誤差: {law_sigma_s:.3f}')

# グラフプロット
color_list = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(2, 2, figsize=(14, 14))
axs[0, 0].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[0, 0].axvline(mu_s1 - sigma_s1, c=color_list[0], ls='--', label='mu_s1±sigma_s1')
axs[0, 0].axvline(mu_s1 + sigma_s1, c=color_list[0], ls='--')
axs[0, 0].set_xlabel('value')
axs[0, 0].set_ylabel('counts')
axs[0, 0].legend()

axs[0, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[0, 1].axvline(mu_s2 - sigma_s2, c=color_list[1], ls='--', label='mu_s2±sigma_s2')
axs[0, 1].axvline(mu_s2 + sigma_s2, c=color_list[1], ls='--')
axs[0, 1].set_xlabel('value')
axs[0, 1].set_ylabel('counts')
axs[0, 1].legend()

axs[1, 0].hist(s, bins=20, alpha=0.6, color=color_list[2], label='s=s1-s2')
axs[1, 0].axvline(mean_s - sigma_s, c='r', ls='--', alpha=0.5, label='Simulation mean_s±sigma_s')
axs[1, 0].axvline(mean_s + sigma_s, c='r', ls='--', alpha=0.5)
axs[1, 0].axvline(law_mean_s - law_sigma_s, c='k', ls='--', alpha=0.5, label='Law law_mean_s±law_sigma_s')
axs[1, 0].axvline(law_mean_s + law_sigma_s, c='k', ls='--', alpha=0.5)
axs[1, 0].set_xlabel('value')
axs[1, 0].set_ylabel('counts')
axs[1, 0].legend()

axs[1, 1].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[1, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[1, 1].hist(s,  bins=20, alpha=0.6, color=color_list[2], label='s=s1-s2')
axs[1, 1].set_xlabel('value')
axs[1, 1].set_ylabel('counts')
axs[1, 1].legend(loc='upper center')

plt.show()

シミュレーションの結果

シミュレーションの平均値: -100.006, 誤差伝搬の法則で求めた平均値: -100.000
シミュレーションの不偏標準偏差: 11.177, 誤差伝搬の法則で求めた誤差: 11.180

シミュレーションの結果

左下の図でほとんど赤線と黒線が被っていますね。誤差伝搬の法則とシミュレーションの結果が良く一致しています。

信号1*信号2

信号1*信号2もシミュレーションしてみます。基本的に設定は信号1+信号2と同じです。
信号($s$)を

s=s_1\cdot s_2

とします。誤差伝搬の法則より、

law\_mean\_s = \mu_{s_1} \cdot \mu_{s_2}

\begin{align}
law\_sigma\_s &= \sqrt{\left(\frac{\partial s}{\partial s_1} \right)^2\sigma_{s_1}^2 + \left(\frac{\partial s}{\partial s_2} \right)^2\sigma_{s_2}^2}\\
&= \left|\mu_{s_1} \cdot \mu_{s_2}\right|\sqrt{\left(\frac{\sigma_{s_1}}{\mu_{s_1}}\right)^2 + \left(\frac{\sigma_{s_2}}{\mu_{s_2}}\right)^2}
\end{align}

となります。誤差伝搬の法則と従っているかシミュレーションで確かめます。

信号1*信号2のコード

信号1*信号2

# s=s1*s2の場合
np.random.seed(2022)
n = 1000000

# 信号1(s1)
mu_s1 = 100
sigma_s1 = 5
s1 = np.random.normal(mu_s1, sigma_s1, n)

# 信号2(s2)
mu_s2 = 200
sigma_s2 = 10
s2 = np.random.normal(mu_s2, sigma_s2, n)

# 信号1*信号2
s = s1 * s2

# シミュレーションの結果の平均値と不偏標準偏差
mean_s = np.mean(s)
sigma_s = np.std(s, ddof=1)

# 誤差伝搬の法則に従う平均値と誤差
law_mean_s = mu_s1 * mu_s2
law_sigma_s = abs(mu_s1 * mu_s2) * np.sqrt((sigma_s1/mu_s1)**2 + (sigma_s2/mu_s2)**2)
print(f'シミュレーションの平均値: {mean_s:.3f}, 誤差伝搬の法則で求めた平均値: {law_mean_s:.3f}')
print(f'シミュレーションの不偏標準偏差: {sigma_s:.3f}, 誤差伝搬の法則で求めた誤差: {law_sigma_s:.3f}')

# グラフプロット
color_list = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(2, 2, figsize=(14, 14))
axs[0, 0].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[0, 0].axvline(mu_s1 - sigma_s1, c=color_list[0], ls='--', label='mu_s1±sigma_s1')
axs[0, 0].axvline(mu_s1 + sigma_s1, c=color_list[0], ls='--')
axs[0, 0].set_xlabel('value')
axs[0, 0].set_ylabel('counts')
axs[0, 0].legend()

axs[0, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[0, 1].axvline(mu_s2 - sigma_s2, c=color_list[1], ls='--', label='mu_s2±sigma_s2')
axs[0, 1].axvline(mu_s2 + sigma_s2, c=color_list[1], ls='--')
axs[0, 1].set_xlabel('value')
axs[0, 1].set_ylabel('counts')
axs[0, 1].legend()

axs[1, 0].hist(s, bins=20, alpha=0.6, color=color_list[2], label='s=s1*s2')
axs[1, 0].axvline(mean_s - sigma_s, c='r', ls='--', alpha=0.5, label='Simulation mean_s±sigma_s')
axs[1, 0].axvline(mean_s + sigma_s, c='r', ls='--', alpha=0.5)
axs[1, 0].axvline(law_mean_s - law_sigma_s, c='k', ls='--', alpha=0.5, label='Law law_mean_s±law_sigma_s')
axs[1, 0].axvline(law_mean_s + law_sigma_s, c='k', ls='--', alpha=0.5)
axs[1, 0].set_xlabel('value')
axs[1, 0].set_ylabel('counts')
axs[1, 0].legend()

axs[1, 1].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[1, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[1, 1].hist(s,  bins=20, alpha=0.6, color=color_list[2], label='s=s1*s2')
axs[1, 1].set_xscale('log')
axs[1, 1].set_xlabel('value')
axs[1, 1].set_ylabel('counts')
axs[1, 1].legend(loc='upper center')

plt.show()

シミュレーションの結果

シミュレーションの平均値: 20000.683, 誤差伝搬の法則で求めた平均値: 20000.000
シミュレーションの不偏標準偏差: 1413.869, 誤差伝搬の法則で求めた誤差: 1414.214

シミュレーションの結果

やはり、左下の図でほとんど赤線と黒線が被っています。誤差伝搬の法則万能ですね。
ところで、左下の図をよく見ると、右に裾が長い非対称なグラフになっているのがわかると思います。この分布は何に従っているのでしょうか。補足で少し触れるのでよろしければご覧ください。

信号1/信号2

信号1/信号2もシミュレーションしてみます。基本的に設定は信号1+信号2と同じです。
信号($s$)を

s=\frac{s_1}{s_2}

とします。誤差伝搬の法則より、

law\_mean\_s = \frac{\mu_{s_1}}{\mu_{s_2}}

\begin{align}
law\_sigma\_s &= \sqrt{\left(\frac{\partial s}{\partial s_1} \right)^2\sigma_{s_1}^2 + \left(\frac{\partial s}{\partial s_2} \right)^2\sigma_{s_2}^2}\\
&= \left|\frac{\mu_{s_1}}{\mu_{s_2}}\right|\sqrt{\left(\frac{\sigma_{s_1}}{\mu_{s_1}}\right)^2 + \left(\frac{\sigma_{s_2}}{\mu_{s_2}}\right)^2}
\end{align}

となります。誤差伝搬の法則と従っているかシミュレーションで確かめます。

信号1/信号2のコード

信号1/信号2

# s=s1/s2の場合
np.random.seed(2022)
n = 1000000

# 信号1(s1)
mu_s1 = 100
sigma_s1 = 5
s1 = np.random.normal(mu_s1, sigma_s1, n)

# 信号2(s2)
mu_s2 = 200
sigma_s2 = 10
s2 = np.random.normal(mu_s2, sigma_s2, n)

# 信号1/信号2
s = s1 / s2

# シミュレーションの結果の平均値と不偏標準偏差
mean_s = np.mean(s)
sigma_s = np.std(s, ddof=1)

# 誤差伝搬の法則に従う平均値と誤差
law_mean_s = mu_s1 / mu_s2
law_sigma_s = (mu_s1 / mu_s2) * np.sqrt((sigma_s1/mu_s1)**2 + (sigma_s2/mu_s2)**2)
print(f'シミュレーションの平均値: {mean_s:.3f}, 誤差伝搬の法則で求めた平均値: {law_mean_s:.3f}')
print(f'シミュレーションの不偏標準偏差: {sigma_s:.3f}, 誤差伝搬の法則で求めた誤差: {law_sigma_s:.3f}')

# グラフプロット
color_list = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(2, 2, figsize=(14, 14))
axs[0, 0].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[0, 0].axvline(mu_s1 - sigma_s1, c=color_list[0], ls='--', label='mu_s1±sigma_s1')
axs[0, 0].axvline(mu_s1 + sigma_s1, c=color_list[0], ls='--')
axs[0, 0].set_xlabel('value')
axs[0, 0].set_ylabel('counts')
axs[0, 0].legend()

axs[0, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[0, 1].axvline(mu_s2 - sigma_s2, c=color_list[1], ls='--', label='mu_s2±sigma_s2')
axs[0, 1].axvline(mu_s2 + sigma_s2, c=color_list[1], ls='--')
axs[0, 1].set_xlabel('value')
axs[0, 1].set_ylabel('counts')
axs[0, 1].legend()

axs[1, 0].hist(s, bins=20, alpha=0.6, color=color_list[2], label='s=s1/s2')
axs[1, 0].axvline(mean_s - sigma_s, c='r', ls='--', alpha=0.5, label='Simulation mean_s±sigma_s')
axs[1, 0].axvline(mean_s + sigma_s, c='r', ls='--', alpha=0.5)
axs[1, 0].axvline(law_mean_s - law_sigma_s, c='k', ls='--', alpha=0.5, label='Law law_mean_s±law_sigma_s')
axs[1, 0].axvline(law_mean_s + law_sigma_s, c='k', ls='--', alpha=0.5)
axs[1, 0].set_xlabel('value')
axs[1, 0].set_ylabel('counts')
axs[1, 0].legend()

axs[1, 1].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[1, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[1, 1].hist(s,  bins=20, alpha=0.6, color=color_list[2], label='s=s1/s2')
axs[1, 1].set_xscale('log')
axs[1, 1].set_xlabel('value')
axs[1, 1].set_ylabel('counts')
axs[1, 1].legend(loc='upper center')

plt.show()

シミュレーションの結果

シミュレーションの平均値: 0.501, 誤差伝搬の法則で求めた平均値: 0.500
シミュレーションの不偏標準偏差: 0.036, 誤差伝搬の法則で求めた誤差: 0.035

シミュレーションの結果

やはり、左下の図でほとんど赤線と黒線が被っています。
ところで、左下の図のように、右に裾が長い非対称なグラフになっています。補足で少し触れるのでよろしければご覧ください。

信号1*ln(信号2)

最後に、もう少し複雑な信号1*ln(信号2)でシミュレーションしてみます。基本的に設定は信号1+信号2と同じです。
信号($s$)を

s=s_1\ln (s_2)

とします。誤差伝搬の法則より、

law\_mean\_s = \mu_{s_1}\ln (\mu_{s_2})

\begin{align}
law\_sigma\_s &= \sqrt{\left(\frac{\partial s}{\partial s_1} \right)^2\sigma_{s_1}^2 + \left(\frac{\partial s}{\partial s_2} \right)^2\sigma_{s_2}^2}\\
&= \sqrt{\left(\ln (\mu_{s_2})\right)^2\sigma_{s_1}^2+\left(\frac{\mu_{s_1}}{\mu_{s_2}}\right)^2\sigma_{s_2}^2}
\end{align}

となります。誤差伝搬の法則と従っているかシミュレーションで確かめます。

信号1*ln(信号2)のコード

信号1*ln(信号2)

# s=s1*ln(s2)の場合
np.random.seed(2022)
n = 1000000

# 信号1(s1)
mu_s1 = 100
sigma_s1 = 5
s1 = np.random.normal(mu_s1, sigma_s1, n)

# 信号2(s2)
mu_s2 = 200
sigma_s2 = 10
s2 = np.random.normal(mu_s2, sigma_s2, n)

# 信号1*ln(信号2)
s = s1 * np.log(s2)

# シミュレーションの結果の平均値と不偏標準偏差
mean_s = np.mean(s)
sigma_s = np.std(s, ddof=1)

# 誤差伝搬の法則に従う平均値と誤差
law_mean_s = mu_s1 * np.log(mu_s2)
law_sigma_s = np.sqrt((np.log(mu_s2)*sigma_s1)**2 + (mu_s1/mu_s2 * sigma_s2)**2)
print(f'シミュレーションの平均値: {mean_s:.3f}, 誤差伝搬の法則で求めた平均値: {law_mean_s:.3f}')
print(f'シミュレーションの不偏標準偏差: {sigma_s:.3f}, 誤差伝搬の法則で求めた誤差: {law_sigma_s:.3f}')

# グラフプロット
color_list = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(2, 2, figsize=(14, 14))
axs[0, 0].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[0, 0].axvline(mu_s1 - sigma_s1, c=color_list[0], ls='--', label='mu_s1±sigma_s1')
axs[0, 0].axvline(mu_s1 + sigma_s1, c=color_list[0], ls='--')
axs[0, 0].set_xlabel('value')
axs[0, 0].set_ylabel('counts')
axs[0, 0].legend()

axs[0, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[0, 1].axvline(mu_s2 - sigma_s2, c=color_list[1], ls='--', label='mu_s2±sigma_s2')
axs[0, 1].axvline(mu_s2 + sigma_s2, c=color_list[1], ls='--')
axs[0, 1].set_xlabel('value')
axs[0, 1].set_ylabel('counts')
axs[0, 1].legend()

axs[1, 0].hist(s, bins=20, alpha=0.6, color=color_list[2], label='s=s1*ln(s2)')
axs[1, 0].axvline(mean_s - sigma_s, c='r', ls='--', label='Simulation mean_s±sigma_s')
axs[1, 0].axvline(mean_s + sigma_s, c='r', ls='--')
axs[1, 0].axvline(law_mean_s - law_sigma_s, c='k', ls='--', label='Law law_mean_s±law_sigma_s')
axs[1, 0].axvline(law_mean_s + law_sigma_s, c='k', ls='--')
axs[1, 0].set_xlabel('value')
axs[1, 0].set_ylabel('counts')
axs[1, 0].legend()

axs[1, 1].hist(s1, bins=20, alpha=0.6, color=color_list[0], label='s1')
axs[1, 1].hist(s2, bins=20, alpha=0.6, color=color_list[1], label='s2')
axs[1, 1].hist(s,  bins=20, alpha=0.6, color=color_list[2], label='s=s1*ln(s2)')
axs[1, 1].set_xlabel('value')
axs[1, 1].set_ylabel('counts')
axs[1, 1].legend(loc='upper center')

plt.show()

シミュレーションの結果

シミュレーションの平均値: 529.712, 誤差伝搬の法則で求めた平均値: 529.832
シミュレーションの不偏標準偏差: 26.940, 誤差伝搬の法則で求めた誤差: 26.959

シミュレーションの結果

やはり、左下の図でほとんど赤線と黒線が被っています。少し複雑な演算ですが、誤差伝搬の法則に従っているのがわかります。

補足

独立した正規分布に従った信号2つの和、差、積、商がどんな分布になるか気になったので調べました。

s=s1+s2, s=s1-s2の分布について

確率変数$X_1, X_2$が、それぞれ$N(\mu_1, \sigma_1^2), N(\mu_2, \sigma_2^2)$に独立に従うとき、定数$a_1,~a_2 \in \mathbb{R}$に対し、

a_1X_1+a_2X_2 \sim N(a_1\mu_1+a_2\mu_2,a_1^2\sigma_1^2+a_2^2\sigma_2^2)

が成り立ちます。このことを正規分布の再生性と呼ばれています。確かに、$s=s_1+s_2,~s=s_1-s_2$のシミュレーションの結果も、正規分布に従っているのているのがわかると思います。

正規分布の再生性の導出などは、

が参考になります。

s=s1*s2の分布について

結論を先に述べると、正規分布ではなさそうです。分布は

のP17の$\rho=0$（2つの分布が独立）に従っていると思われます。式が複雑で検証できていないので、間違えていたらご指摘いただけると幸いです。

ちなみに、

を参考にすると、標準正規分布$N(0,1)$に従う独立した$x,~y$の乱数の積$Z=xy$は、第2種変形ベッセル関数に従い

p_{Z}(z)={\frac {K_{0}(|z|)}{\pi }}

となるそうです。

s=s1/s2の分布について

結論を先に述べると、正規分布ではなさそうです。分布は

に従っていると思われます。式が複雑で検証できていないので、間違えていたらご指摘いただけると幸いです。

ちなみに、

を参考にすると、標準正規分布$N(0,1)$に従う独立した$x,~y$の乱数の尚$Z=x/y$は、コーシー分布に従い

p_{Z}(z)={\frac {1/\pi }{1+z^{2}}}

となるそうです。

まとめ

正規分布に従う誤差を入れてシミュレーションで実装した場合に、誤差伝搬の法則に従うか調べました。その結果、複雑な式でも誤差伝搬の法則に従うことがわかりました。
また、補足では誤差伝搬の結果の分布について触れましたが、積や商の場合は複雑な分布になり奥が深いことがわかりました。ただ、このような演算した結果が複雑な分布の場合でも、誤差伝搬の法則によりしっかりと誤差を見積もることができるのが素晴らしいですね。

最後に、著者もよく使うのですが、式だけ見てもイメージできない場合は、一旦立ち止まりシミュレーションしてみると理解が捗るかも知れません。本記事を通して、誤差伝搬の法則を式とシミュレーションの両面で理解してもらえたら幸いです。

参考資料

• 誤差伝搬の法則について

• 正規分布の再生性について

• 独立した2つの正規分布X*Yの分布

• 独立した2つの正規分布X/Yの分布