淡路(2009)『データ同化』練習問題

淡路他(2009)「データ同化」の練習問題を解いてみました。

データ同化―観測・実験とモデルを融合するイノベーション

作者:淡路敏之,池田元美,石川洋一,蒲地政文
京都大学学術出版会

Amazon

練習問題1.1
- [i]
- [ii]
練習問題1.2
練習問題1.3
練習問題2.1

練習問題1.1

[i]

推定値 $x^{\prime}_a$ は以下の式で表され、バイアス $α_{1}η$ をもつ。

$x^{\prime}_a=α_1 x^{\prime}_1+(1-α_1)x_2=α_1 x_1+(1-α_1)x_2+α_{1}η \tag{1}$

(a)

推定値 $x^{\prime}_a$ の誤差分散は、以下のように表される。 $x_1,x_2$ の誤差相関が0の仮定を用いた。

$\begin{eqnarray} {σ_a}^{\prime2}&=&E[(x^{\prime}_a-x_t)^2] \\ &=&E[\{α_1(x_1-x_t)+(1-α_1)(x_2-x_t)+α_{1}η\}^2]\\ &=&α^2_1E[ε^2_1] + (1-α_1)^2 E[ε^2_2] +α^2_{1}η^2 +2α_1(1-α_1)E[ε_{1}ε_{2}]+2(1-α_1)α_{1}η E[ε_{2}] + α^2_{2} η E[ε_{1}] \\ &=&α^2_1E[ε^2_1] + (1-α_1)^2 E[ε^2_2] +α^2_{1}η^2 \\ &=&α^2_1 σ^2_1 + (1-α_1)^2 σ^2_2 +α^2_{1}η^2\tag{2} \end{eqnarray}$

(2)を最小にする $α_{1}$ を求める。(3)を $α_{1}$ について解くと、

$\begin{eqnarray} \frac{\partial {σ_a}^{\prime2}}{\partialα_1} = 2α_1 σ^2_1 - 2(1-α_1)σ^2_2 + 2α_1 η^2 = 0 \tag{3} \end{eqnarray}$

最適な係数 $α_1$ が求まる。

$\begin{eqnarray} α_1 = \frac{{σ_2}^{2}}{{σ_1}^{2}+{σ_2}^{2}+η^2} \tag{4} \end{eqnarray}$

(1)より $x^{\prime}_a$ のバイアスは $α_1 η$ なので、これに(4)を代入すると具体的なバイアスが求まる。

(b)
(4)を(2)に代入すると、最適推定値 $x^{\prime}_a$ の誤差分散が求まる。

$\begin{eqnarray} {σ_a}^{\prime2}　&=& \left(\frac{{σ_2}^{2}}{{σ_1}^{2}+{σ_2}^{2}+η^2}\right)^2 ({σ_1}^2 + η^2) + \left(\frac{{σ_2}^{2}}{{σ_1}^{2}+{σ_2}^{2}+η^2}\right)^2 {σ_2}^2 \\ &=& \frac{σ_2^2(σ_1^2 + η^2)\{σ_2^2 +( σ_1^2 + η^2)\}}{({σ_1}^{2}+{σ_2}^{2}+η^2)^2}\\ &=& \frac{σ_2^2(σ_1^2 + η^2)}{({σ_1}^{2}+{σ_2}^{2}+η^2)} \tag{5} \end{eqnarray}$

本文(1.11)式は

$\begin{eqnarray} {σ_a}^{2}&=&\frac{σ_1^{2} σ_2^{2}}{σ_1^2 + σ_2^2} \tag{6} \end{eqnarray}$

なので(5)と(6)の差分を計算すると、正の値になることが分かる。つまり $x_a$ はバイアスを持つ推定値を用いて作った推定値であるため、そうではない場合に比べて誤差分散が大きくなる。

[ii]

$x_1,x_2$ の誤差に相関がある（相関係数は $μ$ ）場合の最適推定値は、

(a)

$x^{μ}_a=α_1 x_1+(1-α_1)x_2 \tag{7}$

相関係数の定義より

$\begin{eqnarray} μ&=&\frac{Cov[ε_1,ε_2]}{\sqrt{V[ε_1]V[ε_2]}} \\ &=&\frac{E[ε_1,ε_2]}{\sqrt{E[ε^2_1]E[ε^2_2]}}\\ &=&\frac{E[ε_1,ε_2]}{σ_1 σ_2} \tag{8} \end{eqnarray}$

(7)の誤差分散 $(σ^μ_a)^2$ は

$\begin{eqnarray} (σ^μ_a)^2&=&E[(x^{μ}_a-x_t)^2] \\ &=&E[\{α_1(x_1-x_t)+(1-α_1)(x_2-x_t)\}^2]\\ &=&α^2_1E[ε^2_1] + (1-α_1)^2 E[ε^2_2] +2α_1(1-α_1)E[ε_{1}ε_{2}]\\ &=&α^2_1σ_1+ (1-α_1)^2 σ_2 +2α_1(1-α_1)μ σ_1 σ_2\tag{9} \end{eqnarray}$

(9)を $α_1$ で微分して0と置く。

$\begin{eqnarray} \frac{\partial (σ^μ_a)^2}{\partialα_1} = 2α_1 σ^2_1 - 2(1-α_1)σ^2_2 + 2μ σ_1 σ_2 - 4 α^2_1 μ σ_1 σ_2 = 0 \tag{10} \end{eqnarray}$

(10)を解くと

$\begin{eqnarray} α_1 = \frac{σ^2_2 - μ σ_1 σ_2}{σ^2_1 -2 μ σ_1 σ_2 + σ^2_2 } \tag{11} \end{eqnarray}$

(11)を(7)に代入すると最適推定値 $x^{μ}_a$ が得られ、

(11)を(9)に代入すると、誤差分散が得られる。

$\begin{eqnarray} (σ^μ_a)^2&=&(μ-1)^2 σ^2_1 σ^4_2 + 2(1-μ) σ^3_1 σ^3_2 + (μ-1)^2 σ^4_1 σ^2_2 \tag{12} \end{eqnarray}$

(b)
$μ=0.5$ のとき、 $(σ^μ_a)^2=13$ となり、 $μ=0.75$ のとき、 $(σ^μ_a)^2=5.25$ となりいずれも相関を考慮しない場合の誤差分散よりも値が大きくなる。

誤差に関する不偏性（バイアス0）と無相関の仮定をリラックスした場合の推定値に関する問題。一変数なので丁寧に式展開すれば特段問題ない気がする。回帰分析を勉強したときによくやったような誤差項の仮定と推定量の話とよく似ている。

練習問題1.2

(1.20)の条件付き確率は(1)のように表すことができる。

$\begin{eqnarray} p(x|x_2) p(x)= \frac{p(x_2|x)p(x)}{\int_{-\infty}^{\infty} p(x_2|x)p(x)dx} \tag{1} \end{eqnarray}$

(1.20)の分子を、2変数が正規分布に従う場合の条件付き密度関数を使って表すことができる。

分子の確率密度関数は

$\begin{eqnarray} p(x_2|x)= \frac{1}{\sqrt{2πσ^2_2}} \exp\left(-\frac{1}{2} \frac{(x_2-x)^2}{σ^2_2} \right) \tag{2} \end{eqnarray}$

$\begin{eqnarray} p(x)= \frac{1}{\sqrt{2πσ^2_1}} \exp\left(-\frac{1}{2} \frac{(x-x_1)^2}{σ^2_1} \right) \tag{3} \end{eqnarray}$

(1)の分子は(2),(3)を用いて、

$\begin{eqnarray} p(x_2|x)p(x) &=& \frac{1}{2πσ_{1}σ_{2}} \exp\left\{ -\frac{1}{2} \left(\frac{1}{σ^2_1}+\frac{1}{σ^2_2} \right) \left(x- \left( \frac{σ^2_{2} x_1 + σ^2_{1} x_2}{σ^2_2 + σ^1_2} \right)^2 \right) \right\} \exp \left\{-\frac{1}{2} \left(\frac{1}{σ^2_1}+\frac{1}{σ^2_2} \right) \frac{σ^2_{1} σ^2_{1}}{(σ^2_2 + σ^2_1)^2} (x_1 - x_2)^2 \right\} \\ &\propto& \frac{1}{2πσ_{1}σ_{2}} \exp\left\{ -\frac{1}{2} \left(\frac{1}{σ^2_1}+\frac{1}{σ^2_2} \right) \left(x- \left( \frac{σ^2_{2} x_1 + σ^2_{1} x_2}{σ^2_2 + σ^1_2} \right)^2 \right) \right\} \tag{4} \end{eqnarray}$

(4)の2つ目の指数関数は $x_1,x_2$ とパラメータのみから構成される定数であり、 $x$ の関数を考えるときは無視する。

(1),(4)より(1)の分母は

$\begin{eqnarray} \int_{-\infty}^{\infty} p(x_2|x)p(x)dx &=& \int_{-\infty}^{\infty} \frac{1}{2πσ_{1}σ_{2}} \exp\left\{ -\frac{1}{2} \left(\frac{1}{σ^2_1}+\frac{1}{σ^2_2} \right) \left(x- \left( \frac{σ^2_{2} x_1 + σ^2_{1} x_2}{σ^2_2 + σ^1_2} \right)^2 \right) \right\} dx \\ &=& \frac{1}{2πσ_{1}σ_{2}} \int_{-\infty}^{\infty} \exp\left\{ -\frac{1}{2} \left(\frac{1}{σ^2_1}+\frac{1}{σ^2_2} \right) \left(x- \left( \frac{σ^2_{2} x_1 + σ^2_{1} x_2}{σ^2_2 + σ^1_2} \right)^2 \right) \right\} dx \tag{5} \end{eqnarray}$

ガウス積分の公式(6)を用いて

$\begin{eqnarray} \int_{-\infty}^{\infty} \exp(-a x^2) = \sqrt{\frac{π}{a}} \tag{6} \end{eqnarray}$

(5)の最後の式は以下のように表せる。

$\begin{eqnarray} \frac{1}{2πσ_{1}σ_{2}} \frac{\sqrt{π}}{\sqrt{\frac{1}{2} \left(\frac{σ^2_1+σ^2_2}{σ^2_{1}σ^2_{2}}\right)}} \tag{7} \end{eqnarray}$

(4),(7)より

$\begin{eqnarray} p(x_2|x)p(x) &=& \frac{1}{\sqrt{2π\frac{σ^2_{1}σ^2_{2}}{σ^2_1+σ^2_2}}} \exp\left\{ -\frac{1}{2} \left(\frac{σ^2_1+σ^2_2}{σ^2_1 σ^2_2} \right) \left(x- \left( \frac{σ^2_{2} x_1 + σ^2_{1} x_2}{σ^2_2 + σ^1_2} \right)^2 \right) \right\} \tag{8} \end{eqnarray}$

(8)は平均 $\frac{σ^2_{2} x_1 + σ^2_{1} x_2}{σ^2_2 + σ^1_2}$ 、分散 $\frac{σ^2_{1}σ^2_{2}}{σ^2_1+σ^2_2}$ の正規分布の密度関数と同じ形になることが分かる。これで(1.20)の式を導出できる。

これは間違っているかもしれません。(4)の式変形でどうしても不要な項が出てきてしまい、等式（ $=$ ）だと成立しないため比例（ $\propto$ ）で無理矢理つじつまを合わせました。密度関数の性質が関係してそうだと考え、積分したら1になるかどうか試したのですがよく分からず、消化不良です。

以下のページを参考にしました。

inzkyk.xyz

練習問題1.3

(1.21)、(1.22)式を使って(1.23)式を導出する。

$\begin{eqnarray} R(x_a) &=& \int (x_a - \bar{x} + \bar{x} - x)^2 p(x|x_2) dx \\ &=& (x_a - \bar{x})^2 \int p(x|x_2) dx + 2(x_a - \bar{x}) \left\{ \int \bar{x} p(x|x_2) dx - \int x p(x|x_2) dx \right\} + \int ( \bar{x} - x)^2 p(x|x_2) dx \\ &=& (x_a - \bar{x})^2 + 2(x_a - \bar{x}) \left\{ \bar{x} \int p(x|x_2) dx - \bar{x} \right\} + \int ( \bar{x} - x)^2 p(x|x_2) dx\\ &=& (x_a - \bar{x})^2 + \int ( \bar{x} - x)^2 p(x|x_2) dx \tag{1} \end{eqnarray}$

最後の式展開は二項目がゼロになることを用いた。

練習問題2.1

(2.16)の導出

解析誤差共分散行列（本文(2.9)）を以下のように変形する。

$\begin{eqnarray} \mathbf{P^a} &=& \mathbf{B} - \mathbf{BH^{T}W^{T}+WRW^{T}-WHB+WHBH^{T}W^{T}} \\ &=& \mathbf{B}-\mathbf{WHB- (BH^{T} - WR - WHB^{T})W^{T}} \tag{1} \end{eqnarray}$

本文(2.12)を0とおく（ $\mathbf{P^a}$ の対角成分（解析誤差分散）の最小化問題を解く）と、

$\begin{eqnarray} \mathbf{-2BH^{T} + 2WR + 2WHBH^{T}} = 0 \tag{2} \end{eqnarray}$

(2)を(1)に代入すると(2.16)を導出できる。

$\begin{eqnarray} \mathbf{P^a} &=& \mathbf{B-WHB} \\ &=& \mathbf{(I-WH)B} \tag{3} \end{eqnarray}$

(2.15)の導出
重み行列 $\mathbf{W}$ は、本文(2.14)と付録(A.11)の逆行列補題を用い、以下のように変形する

$\begin{eqnarray} \mathbf{W} &=& \mathbf{(B^{-1}+H^{T}R^{-1}H)^{-1}H^{T}R^{-1}} \\ &=& \mathbf{BH^{T}(HBH^{T}+R)^{-1}} \tag{4} \end{eqnarray}$

(4)を(1)に代入して、付録(A.12)の逆行列補題を用いて式変形すると、

$\begin{eqnarray} \mathbf{P^a} &=& \mathbf{(I-(BH^{T}(HBH^{T}+R)^{-1})H)B} \\ &=& \mathbf{B -BH^{T}(HBH^{T}+R)^{-1}HB } \\ &=& \mathbf{(H^{T}R^{-1}H+B^{-1})^{-1}} \tag{5} \end{eqnarray}$