最小二乗法と格子定数精密化

最小二乗法と格子定数精密化
一般にn個の観測値の組

$(X^1_1,X^2_1,\cdots,X^m_1,Z_1),(X^1_2,X^2_2,\cdots,X^m_2,Z_2),\cdots,(X^1_n,X^2_n,\cdots,X^m_n,Z_n)$

が与えられたとき、パラメータについて線形的な観測方程式にフィッティングさせるためには、二乗残差を最小にすればいいということは、よく知られています。(上肩添え字は累乗ではなくパラメータの個数です。念のため。)

簡単に解き方を示すと、

$a = \left(\begin{array}{ccc} a^1 \\ a^2 \\ \vdots \\ a^m \end{array}\right) , X= \left( \begin{array}{cccc} X^1_1 & X^2_1 & \ldots & X^m_1 \\ X^1_2 & X^2_2 & \ldots & X^m_2 \\ \vdots & \vdots & \ddots & \vdots \\ X^1_n & X^2_n & \ldots & X^m_n \\ \end{array}\right) , Y=\left(\begin{array}{c}Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{array} \right)$
$W= \left( \begin{array}{cccc} W_1 & 0 & \ldots & 0 \\ 0 & W_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & W_n \\ \end{array}\right)$

としたとき (ただし $W$ は重みを対角成分に持つ対角行列) 、

$\Delta^2=(\tilde{a} \tilde{X} - \tilde{Y})W(X a -Y)$

を最小にすればいいから
$\delta \Delta^2 = \delta \tilde{a} \tilde{X} W (X a - Y) + (\tilde{a} \tilde{X} - \tilde{Y}) W X \delta a = 2\delta\tilde{a}(\tilde{X} W X a - \tilde{X} W Y) = 0$
よって
$a = (\tilde{X} W X)^{-1} \tilde{X} W Y$
となります。
格子定数の計算に際しては、結晶系によって観測方程式は変わりますが、たとえば三斜晶系のときは
$a= \left(\begin{array}{c}a^1 \\ a^2 \\ \vdots \\ a^m \end{array} \right) = \left(\begin{array}{c}a^*a^* \\ b^*b^* \\ c^*c^* \\ b^*c^*\cos\alpha^* \\ c^*a^*\cos\beta^* \\ a^*b^*\cos\gamma^* \end{array} \right)$

のようになります。ここで、 $a^*, b^*$ などは逆格子定数を示しています。

重みは角度に依存していると考えられます。 $\theta$ の測定にのみ誤差が生じていると考えると、 $\theta$ の変化に対する $1/d^2$ の応答は、ブラッグ条件から
$G = 1/d^2 = 4\sin^2\theta/\lambda$
としたとき
$\partial{G}/\partial{\theta} = 4\sin 2\theta/\lambda$
となります。つまり、　 $\theta$ 　が　 $\partial{\theta}$ 　変化すると　 $1/d^2$ 　は　 $4\sin 2\theta/\lambda$ 　だけ変化するということです。ということは　 $1/d^2$ 　に対する重み行列対角成分は誤差の2乗の逆数のである　 $1/\sin^{2}2\theta$ 　が適当であると考えられます。
$W= \left( \begin{array}{cccc} 1/\sin^{2}2\theta_1 & 0 & \ldots & 0 \\ 0 & 1/\sin^{2}2\theta_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & 1/\sin^{2}2\theta_n \\ \end{array}\right)$
ただし $1/\sin^{2}2\theta$ は各 $1/d^2$ の分散の逆数の絶対値ではなくその比率をあらわしています。
にもかかわらず上述の $a = (\tilde{X} W X)^{-1} \tilde{X} W Y$ という関係式で最適値を求めることができます。なぜなら、重み行列 $W$ が2回作用することによってキャンセルアウトするからです。

これで一応パラメータ $a$ は求まりました。

次は誤差についてのお話です。誤差の伝播式を考えると、もし正しい $W$ を使っていれば、行列 $(\tilde{X} W X)^{-1}$ の対角成分が各パラメータの分散になります。しかしながら、上で述べたように $W$ は割合しか決めていないのでその絶対値が分かりません。そこで、少し目先を変えて、分散の定義を考えてみます。すなわち
$\displaystyle N - P = \sum^{N}_{i=0}\cfrac{\delta_i^2}{s_i}$
(ただし $N$ : データ個数, $P$ : パラメータ数, $\delta$ : データの偏差, $s$ : データの分散)
という関係です。この式を用いて、求まったパラメータから以下の量を計算すれば各パラメータの分散の絶対値がわかるということになります。
$s_i = \cfrac{\displaystyle\sum^{N}_{i=0}\delta_i^2w_i}{w_i (N-P)}$
この分散の平方根が誤差になります。これらはあくまで $a$ (すなわち逆格子定数)の誤差なので、格子定数の誤差に戻すためにはさらに誤差の伝播を解いていく必要がありますが、根気さえあればできるのでここでは省略します。