分析化学の話（9） Z=XY：積の確率分布

分析化学とは直接関係ない、わき道にそれた話です。興味がなかったら、無視してください。

１．共通の疑問

1.1 共通の前提

確率変数（例：分析結果）$${Z}$$が、互いに独立な確率変数（例：測定値、データ）$${X}$$と$${Y}$$の関数になっているとします：

$$
Z = F(X, Y)　\quad \quad \quad \quad \quad \quad \quad \quad (1) 
$$

表1 各確率変数の母集団に関するパラメータ

1.2 共通の疑問

母集団に対して、誤差の伝播則[1, 2]に似た次の式が本当に成り立つでしょうか？

$$
\sigma^2= \bigg (\dfrac{\partial F}{\partial x} \bigg )_{\mu_1, \mu_2}^2 \sigma_1^2 + \bigg (\dfrac{\partial F}{\partial y} \bigg )_{\mu_1, \mu_2}^2 \sigma_2^2 \quad \quad \quad \quad  (2) 
$$

式(2)で注意したいのは、$${X}$$と$${Y}$$の確率密度関数$${f(x)}$$と$${g(y)}$$に無関係に成り立つということです。

1.3 結論

最初に結論を述べると、

【結論】
積 $${Z=XY}$$の場合、誤差（不確かさ）の伝播則がほぼ成立しそうです。

２．基礎となる式

次の式が成り立ちます[3, 4]：

$$
q(z) = \displaystyle \int_{-\infty}^\infty \int_{-\infty}^\infty\delta (z-F(x, y))f(x)g(y) \mathrm{d}x \mathrm{d}y \quad \quad \quad \quad (3)
$$

$${\delta (x)}$$：デルタ関数

３． Z=XY

3.1 具体的な問題の設定

$${F(x, y) = xy}$$だから、式(3)から$${q(z)}$$は次式で与えられます：

$$
q(z) = \displaystyle \int_{-\infty}^\infty \int_{-\infty}^\infty\delta (z-xy))f(x)g(y) \mathrm{d}x \mathrm{d}y \quad \quad \quad \quad (4)
$$

$${\delta}$$関数の性質を使います[5]。式(A5)において$${u(x)=z-xy}$$と於けば、$${u'(x)=y}$$であり、$${u(x)=0}$$の解$${\alpha}$$は$${z/y}$$だから、

$$
q(z) = \displaystyle \int_{-\infty}^\infty \dfrac{1}{|y|} f\Big( \dfrac{z}{y} \Big)g(y) \mathrm{d}y \quad \quad \quad \quad (5)
$$

$$
\sigma^2 = \displaystyle \int_{-\infty}^\infty (z-\mu)^2q(z) \mathrm{d}z \quad \quad \quad \quad (6)
$$

一方、式(2)の偏微分はというと、

$$
\bigg(\dfrac{\partial F}{\partial x} \bigg)_{\mu_1, \mu_2}= \bigg(\dfrac{\partial (xy)}{\partial x} \bigg)_{\mu_1, \mu_2}= \bigg(y\dfrac{\partial x}{\partial x} \bigg)_{\mu_1, \mu_2}= (y)_{\mu_1, \mu_2} = \mu_2
$$

$$
\bigg(\dfrac{\partial F}{\partial y} \bigg)_{\mu_1, \mu_2}= \bigg(\dfrac{\partial (xy)}{\partial y} \bigg)_{\mu_1, \mu_2}= \bigg(x\dfrac{\partial y}{\partial y} \bigg)_{\mu_1, \mu_2}= (x)_{\mu_1, \mu_2} = \mu_1
$$

だから、式(2)に該当するのは次式です：

$$
\sigma^2= \mu_2^2\sigma_1^2 + \mu_1^2\sigma_2^2 \quad \quad \quad \quad  (7) 
$$

式(6)から得られた$${\boldsymbol {\sigma}}$$と、式(7)から得られた$${\boldsymbol {\sigma}}$$とが、どの程度一致するのかが問題です。
　ちなみに、平均値は、

$$
\mu = \mu_1 \mu_2 \quad \quad \quad \quad (8)
$$

3.2 正規分布の場合

3.2.1 もとの分布

$$
f(x) = \dfrac{1}{\sqrt{2\pi} \sigma_1} \exp \bigg [ - \dfrac{(x-\mu_1)^2}{2\sigma_1^2} \bigg ]\quad \quad \quad \quad (9)
$$

$$
g(y) = \dfrac{1}{\sqrt{2\pi} \sigma_2} \exp \bigg [ - \dfrac{(y-\mu_2)^2}{2\sigma_2^2} \bigg ]\quad \quad \quad \quad (10)
$$

3.2.2 $${\boldsymbol{q(z)}}$$

式(5)から出発します。式(9)と(10)を式(5)に代入すると、

$$
q(z) = \dfrac{1}{2\pi\sigma_1\sigma_2} \displaystyle \int_{-\infty}^\infty \dfrac{1}{|y|} \exp \Big[ - \dfrac{(z/y - \mu_1)^2}{2\sigma_1^2} -\dfrac{(y-\mu_2)^2}{2\sigma_2^2} \Big]\mathrm{d}x \mathrm{d}y \quad (11)
$$

式(11)は解析的に積分を表せないので、Pythonで数値積分を行うことにしました。

3.2.3 平均値と標準偏差

$$
\mu = \displaystyle \int_{-\infty}^\infty zq(z) \mathrm{d}z
$$

$$
\sigma^2 = \displaystyle \int_{-\infty}^\infty (z-\mu)^2q(z) \mathrm{d}z
$$

もPythonで数値積分しました。

3.2.4 グラフ

i) 極端（XとYの相対標準偏差[6]が20%）な例

図1 積の分布関数q(z)（Z=XY）。標準偏差が大きい例。f(x): μ1=2, σ1=0.4、g(y): μ2=10, σ2=2。q(z)_GUM（黒）は、平均値(μ)=μ1×μ2=20で、式(7)より計算された標準偏差(σ)をもつ正規分布。ヒストグラム（正規乱数の数は40,000）はあきらかにq(z)と一致しています。Pythonで計算しました。

ii) 通常（XとYの相対標準偏差が2%）の例

図2 積の分布関数q(z)（Z=XY）。標準偏差が小さい例。f(x): μ1=2, σ1=0.04、g(y): μ2=10, σ2=0.2。q(z)_GUM（黒）は、平均値(μ)=μ1×μ2=20で、式(7)より計算された標準偏差(σ)をもつ正規分布。ヒストグラム：正規乱数の数は40,000。q(z)、q(z)_GUM、ヒストグラムが一致しています。Pythonで計算しました。

2.3.5 結果
図1と図2の結果を表にすると、

$$
\begin{array}{lcc}
\small \boldsymbol{表2} \quad 極端な例（XとYの相対標準偏差が20％）& &\\
\hline \hline
\small 方法 & \small 平均値(\mu) & \small 標準偏差(\sigma)\\
\hline \hline
\small q(z)_{\mathrm{GUM}}:式(2), (7)から計算 & \small 20.00&\small 5.66 \\
\hline \small q(z)_{\mathrm{GUM}}に式(\mathrm{A}6)の補正を加えた[7] & \small  20.00 &\small 5.71\\
\hline  \small q(z)から計算& \small 20.00&\small 5.71\\
\hline \small 乱数（ヒストグラム）& \small 20.00&\small 5.72\\
\hline
\end{array}
$$

$$
\begin{array}{lcc}
\small \boldsymbol{表3}\quad 通常の例（XとYの相対標準偏差が2％）& &\\
\hline \hline
\small 方法 & \small 平均値(\mu) & \small 標準偏差(\sigma)\\
\hline \hline
\small q(z)_{\mathrm{GUM}}:式(2), (7)から計算 & \small 20.00&\small 0.566 \\
\hline \small q(z)_{\mathrm{GUM}}に式(\mathrm{A}6)の補正を加えた[7] & \small  20.00 &\small 0.566\\
\hline  \small q(z)から計算& \small 20.00&\small 0.566\\
\hline \small 乱数（ヒストグラム）& \small 20.00&\small 0.567\\
\hline
\end{array}
$$

表2では、$${X}$$と$${Y}$$の標準偏差が大きいにもかかわらず（それぞれ相対標準偏差[6]が20%と20%）、数値の結果はほとんど同じ値を与えています。標準偏差が小さいほど（相対標準偏差が2%）$${q(z)}$$と$${q(z)_{\mathrm{GUM}}}$$とが一致してくるようです（表3）。
　上記から、誤差（不確かさ）の伝播則（式(2)、(7)）がほぼ成立していると思われます。

2.4 一様分布（矩形分布）の場合

2.4.1 もとの分布

$$
f(x) = \dfrac{1}{2a_1} \quad \quad(\mu_1 - a_1 \leq x \leq \mu_1 + a_1)\\
f(x)= 0  \quad (x<\mu_1 - a_1, \quad x> \mu_1 + a_1)
$$

$$
g(y) = \dfrac{1}{2a_2} \quad \quad(\mu_2- a_2 \leq y \leq \mu_2 + a_2)\\
g(y)= 0  \quad (x<\mu_2 - a_2, \quad y> \mu_2 + a_2)
$$

図3 f(x)とg(y)の一様分布（矩形分布）のグラフ

2.4.2 前提

$$
\mu_1 \leq \mu_2
$$

2.4.3 $${\boldsymbol{z}}$$に関するパラメータ

$$
x_1 = \mu_1 - a_1 \\
x_2 = \mu_1 + a_1 \\
y_1 = \mu_2 - a_2 \\
y_2 = \mu_2 + a_2
$$

$$
z_1 = x_1 y_1 \\
z_2 = x_2 y_1 \\
z_3 = x_1 y_2 \\
z_4 = x_2 y_2
$$

2.4.4 $${\boldsymbol{q(z)}}$$
i) $${\boldsymbol{z_2 \leq z_3}}$$の場合：

$$
z_1 < z_2 \leq z_3 < z_4
$$

$$
q(z) = 0, \quad \quad \quad (z < z_1 , \; z\geq z_4)
$$

$$
q(z) = \dfrac{1}{4a_1a_2} \displaystyle \int_{y_1}^{z/x_1} \dfrac{1}{y} \mathrm{d}y = \dfrac{1}{4a_1a_2} \Big[ \ln y \Big]_{y_1}^{z/x_1} \\ = \dfrac{1}{4a_1a_2} \Big[ \ln(z/x_1)-\ln(z_1/x_1) \Big] = \dfrac{1}{4a_1a_2}\ln \dfrac{z}{z_1}, \quad (z_1 \leq z < z_2)
$$

$$
q(z) =  \dfrac{1}{4a_1a_2} \displaystyle \int_{z/x_2}^{z/x_1}\dfrac{1}{y} \mathrm{d}y = \dfrac{1}{4a_1a_2} \ln \Big( \dfrac{x_2}{x_1} \Big), \quad (z_2 \leq z < z_3 )
$$

$$
q(z) = \dfrac{1}{4a_1a_2} \displaystyle \int_{z/x_2}^{y_2} \dfrac{1}{y} \mathrm{d}y = \dfrac{1}{4a_1a_2} \ln \dfrac{z_4}{z}, \quad (z_3 \leq z < z_4)
$$

図4 q(z)の計算の説明図(z2<z3)（f(x): μ1 = 2, a1 = 0.1、g(y): μ2 = 4, a2 = 0.5）。zが増加するにつれて、f(z/y)（灰色）が右に移動します（(1)→(2)→…→(8)→(9)）。f(z/y)がg(y)と重なる場合がq(z)が0にならない場合です（f(z-y)g(y)≠0、赤の部分）。z1 = 6.65, z2 = 7.35, z3 = 8.55, z4 = 9.45。

ii) $${\boldsymbol{z_2 > z_3}}$$の場合：

$$
z_1 < z_3 < z_2 < z_4
$$

$$
q(z) = 0, \quad \quad \quad (z < z_1 , \; z\geq z_4)
$$

$$
q(z) = \dfrac{1}{4a_1a_2} \displaystyle \int_{y_1}^{z/x_1} \dfrac{1}{y} \mathrm{d}y = \dfrac{1}{4a_1a_2} \ln \dfrac{z}{z_1}, \quad (z_1 \leq z < z_3)
$$

$$
q(z) = \dfrac{1}{4a_1a_2} \displaystyle \int_{y_1}^{y_2} \dfrac{1}{y} \mathrm{d}y =  \dfrac{1}{4a_1a_2} \ln \dfrac{y_2}{y_1}, \quad (z_3 \leq z < z_2)
$$

$$
q(z) = \dfrac{1}{4a_1a_2} \displaystyle \int_{z/x_2}^{y_2} \dfrac{1}{y} \mathrm{d}y = \dfrac{1}{4a_1a_2} \Big( \ln{y_2} - \ln{(z/x_2)} \Big) \\ =  \dfrac{1}{4a_1a_2} \Big( \ln{x_2y_2} - \ln{z} \Big) = \dfrac{1}{4a_1a_2} \ln \dfrac{z_4}{z}, \quad (z_2 \leq z < z_4)
$$

2.4.5 平均値と標準偏差

i) $${\boldsymbol{z_2 \leq z_3}}$$の場合：
平均値$${\mu}$$：

$$
\mu = \displaystyle \int_{z_1}^{z_4} zq(z) \mathrm{d}z =
\\
= \dfrac{1}{4a_1a_2} \bigg[ \dfrac{z_2^2 \ln(z_2/z_1)}{2} - \dfrac{z_2^2 - z_1^2}{4} + \dfrac{1}{2} (z_3^2 - z_2^2) \ln \Big( \dfrac{x_2}{x_1} \Big) +
\\
+\dfrac{z_4^2 - z_3^2}{4} - \dfrac{z_3^2 \ln(z_4/z_3)}{2} \bigg]
$$

標準偏差$${\sigma}$$：

$$
\sigma^2 = \displaystyle \int_{z_1}^{z_4} (z-\mu)^2q(z) \mathrm{d}z = \bar{z^2} - \mu^2
$$

$$
\bar{z^2} = \dfrac{1}{4a_1a_2} \bigg[ \dfrac{z_2^3 \ln(z_2/z_1)}{3} - \dfrac{z_2^3 - z_1^3}{9} + \dfrac{1}{3} (z_3^3 - z_2^3) \ln{\Big(\dfrac{x_2}{x_1} \Big)} -
\\
- \dfrac{z_3^3 \ln{(z_4/z_3)}}{3} + \dfrac{z_4^3 - z_3^3}{9} \bigg]
$$

ii) $${\boldsymbol{z_2 > z_3}}$$の場合：
平均値$${\mu}$$：

$$
\mu = \displaystyle \int_{z_1}^{z_4} zq(z) \mathrm{d}z =
\\
=\dfrac{1}{4a_1a_2} \bigg( \dfrac{z_3^2 \ln(z_3/z_1)}{2} - \dfrac{z_3^2 - z_1^2}{4} + \dfrac{1}{2} (z_2^2 - z_3^2) \ln { \Big( \dfrac{y_2}{y_1} \Big)} +\\
+ \dfrac{z_4^2 - z_2^2}{4} - \dfrac{z_2^2 \ln{(z_4/z_2)}}{2}
$$

標準偏差$${\sigma}$$：

$$
\sigma^2 = \displaystyle \int_{z_1}^{z_4} (z-\mu)^2q(z) \mathrm{d}z = \bar{z^2} - \mu^2
$$

$$
\bar{z^2} = \dfrac{1}{4a_1a_2} \bigg[ \dfrac{z_3^3 \ln(z_3/z_1)}{3} - \dfrac{z_3^3 - z_1^3}{9} + \dfrac{1}{3} (z_2^3 - z_3^3) \ln{\Big(\dfrac{y_2}{y_1} \Big)} -
\\
- \dfrac{z_2^3 \ln{(z_4/z_2)}}{3} + \dfrac{z_4^3 - z_2^3}{9} \bigg]
$$

2.4.6 グラフ

図5 q(z)（Z=XY）のグラフ。q(z)がヒストグラムと一致しているのがわかります。ヒストグラムは、Pythonのモジュールnumpyの一様分布に基づく乱数から作成しました（乱数の数=40,000）。Xの一様分布f(x): μ1 = 2, a1 = 0.1、Yの一様分布g(y): μ2 = 3, a2 = 0.3。z2<z3でした。 計算とヒストグラム作成はともにPythonで行いました。

2.4.7 結果
図5の結果は下の表のようでした：

$$
\begin{array}{lcc}
\small \boldsymbol{表4 }　XとYの相対標準偏差が各5％と10％& &\\
\hline \hline
\small 方法 & \small 平均値(\mu) & \small 標準偏差(\sigma)\\
\hline \hline
\small q(z)_{\mathrm{GUM}}:式(2), (7)から計算 & \small 6&\small 0.387 \\
\hline  \small q(z)から計算& \small 6.000&\small 0.387\\
\hline \small 乱数（ヒストグラム）\quad \quad & \small 5.996&\small 0.385\\
\hline
\end{array}
$$

表4と図5の結果から、式(2)または式(7)（誤差（不確かさ）の伝播則）が成立していると思われます。

文献とNote

[1] 誤差（不確かさ）の伝播則、あるいは不確かさの合成は、標準不確かさ$${u}$$を用いて次式のように表されます（$${X}$$と$${Y}$$が独立）：

$$
u^2(z)= \bigg (\dfrac{\partial F}{\partial x} \bigg )^2 u^2(x) + \bigg (\dfrac{\partial F}{\partial y} \bigg )^2 u^2(y) \quad \quad \quad \quad  (\mathrm{A}1) 
$$

標準不確かさは標準偏差と同じような意味を持っています。
[2] ISO, Guide to the Expression of Uncertainty in Measurement (1995). 略して、"GUM"。誤差（不確かさ）の伝播則に従う確率密度関数に下付きの添え字"GUM"をつけました。
[3] 緑川章一、"確率変数の和，積，商，べき乗の分布".
[4] 早川美徳、"Pythonプログラミング（ステップ７・確率密度関数とその計算）" Feb. 15, 2022.
[5] $${\delta}$$関数の性質：

$$
\delta (-x) = \delta (x) \quad \quad \quad \quad (\mathrm{A}2)
$$

$$
\displaystyle \int_{-\infty}^\infty \delta(x-a)f(x) \mathrm{d}x = f(a)\quad \quad \quad \quad (\mathrm{A}3)
$$

$$
\delta(ax) = \dfrac{1}{|a|} \delta(x)\quad \quad \quad \quad (\mathrm{A}4)
$$

$$
\delta (u(x)) = \dfrac{1}{|u'(\alpha)|} \delta(x-\alpha)　\quad \quad \quad \quad (\mathrm{A}5)
$$

ただし、$${\alpha}$$は、$${u(x)=0}$$を$${x}$$の方程式とみた解。

[6] 相対標準偏差 = 標準偏差 ÷ 平均値
[7] 収束するか個人的には不明だけど、とりあえず$${F(x,y) = xy}$$をテイラー（Talyor）展開し、2次の項までとると：

$$
F(x,y) = F(\mu_1, \mu_2) + \Big( \dfrac{\partial F}{\partial x} \Big)_{\mu_1, \mu_2}(x-\mu_1) + \Big( \dfrac{\partial F}{\partial y} \Big)_{\mu_1, \mu_2}(y-\mu_2) +
\\
+\dfrac{1}{2} \Big(\dfrac{\partial^2 F}{\partial x^2} \Big)_{\mu_1, \mu_2}(x-\mu_1)^2 +\dfrac{1}{2} \bigg[ \Big(\dfrac{\partial^2 F}{\partial x \partial y} \Big)_{\mu_1, \mu_2} + \Big(\dfrac{\partial^2 F}{\partial y \partial x} \Big)_{\mu_1, \mu_2} \bigg](x-\mu_1)(y-\mu_2)+
\\
+ \dfrac{1}{2}\Big(\dfrac{\partial^2 F}{\partial y^2} \Big)_{\mu_1, \mu_2}(y-\mu_2)^2
$$

ここで、

$$
\partial F / \partial x = y \\
\partial F / \partial y = x \\
\partial^2 F / \partial x^2 = 0 \\
\partial^2 F / \partial y^2 = 0 \\
\partial^2 F / \partial x \partial y = \partial x / \partial x = 1 \\
\partial^2 F / \partial y \partial x = \partial y / \partial y = 1
$$

を使って、

$$
F(x,y) = \mu_1 \mu_2 + \mu_2 (x-\mu_1) +
\\
+ \mu_1 (y-\mu_2) + (x-\mu_1)(y-\mu_2)
$$

$$
\mu = \displaystyle \int_{-\infty}^\infty \int_{-\infty}^\infty F(x,y)f(x)g(y) \mathrm{d}x \mathrm{d}y = \int_{-\infty}^\infty \int_{-\infty}^\infty [\mu_1 \mu_2 + \mu_2 (x-\mu_1) +
\\
+ \mu_1 (y-\mu_2) + (x-\mu_1)(y-\mu_2)] f(x)g(y) \mathrm{d}x \mathrm{d}y = \mu_1 \mu_2
$$

$$
\sigma^2 = \displaystyle \int_{-\infty}^\infty \int_{-\infty}^\infty (F(x,y) - \mu)^2 f(x)g(y) \mathrm{d}x \mathrm{d}y = \mu_2^2 \int_{-\infty}^\infty (x-\mu_1)^2 f(x)\mathrm{d}x + \\
+ \mu_1^2 \int_{-\infty}^\infty (y-\mu_2)^2 g(y)\mathrm{d}y + \int_{-\infty}^\infty (x-\mu_1)^2 f(x) \mathrm{d}x  \int_{-\infty}^\infty (y-\mu_2)^2 g(y) \mathrm{d}y \\
= \mu_2^2 \sigma_1^2 + \mu_1^2 \sigma_2^2 + \sigma_1^2 \sigma_2^2
$$

最後の式の右辺の2項は式(7)と同一。よって、補正項は、3項目：

$$
補正項 = \sigma_1^2 \sigma_2^2 \quad \quad \quad \quad \quad \quad (\mathrm{A}6)
$$

---

【免責事項】本記事は単なるメモとして書かれたもので、その正確性を必ずしも保証するものではありません。本記事によって生じたトラブル、損失、又は損害に対して一切責任を負いません。また、著者が所属する組織とは関係ありません。誤りがあればご指摘ください。クレームはご遠慮ください。