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

最小二乗法と格子定数精密化
一般に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 (すなわち逆格子定数)の誤差なので、格子定数の誤差に戻すためにはさらに誤差の伝播を解いていく必要がありますが、根気さえあればできるのでここでは省略します。