異なる精度を持った測定値 \(\{y_i\}\) を取り扱う場合、それぞれの測定値の分散を \(\sigma_i^2\) とすると、 \(1/\sigma_i^2\) の重みを付けた残差2乗和 \(S^*\) を最小にすることを考えればよい:
\begin{equation} S^* = \sum_{i=1}^N {\frac{[y_i - (ax_i + b)]^2}{\sigma_i^2}} \label{eq:weight1} \end{equation}
化学で出会う典型的な例としては、有効数字3ケタのデータの線形の式へのあてはめ、 相対誤差が一定と見なせるデータのあてはめの問題がある。 こうした場合、その対数を取ったものへの線形の式のあてはめは、 先の分散一定とした取り扱いで十分であることに注意する:
\begin{equation} \ln [y (1 \pm \delta y)] = \ln y \pm \delta y \label{eq:weight2} \end{equation}
したがって例えばサーミスターの抵抗 \(R\) を各温度で有効数字 3 ケタ程度で測定し、 そこからサーミスターの抵抗の温度依存性を与える関係式 \(R_0 \exp(B/T)\) の \(B\) パラメータを決める場合には、 \(\ln (R/\Omega)\) と \(1/T\) について、分散を一定と仮定して前節の扱いに従えばよい。
凝固点降下による分子量の決定の実験に関わって、 サーミスターの抵抗の温度依存性を調べた学生さんたちの結果から紹介しましょう。
|
図 4-1a. ln (R / kΩ) を 1000 K/ T に対してプロットしたもの。 青線はデータの重みを一定に取って得た線形最小2乗法で得た直線。 テキストボックスの結果は、 得られた線形最小2乗法のパラメーター。 |
|
| 図 4-1b. R / kΩ を 1000 K/ T に対してプロットしたもの。 青線はデータの重みを一定に取って、\(f(T) = R_{25} \exp{1000 B (1/T - 1/298.15 \mrm{~K})}\) という関数にデータの重みを一定に取った非線形最小2乗法で当てはめて得た曲線。 テキストボックスの結果は、 得られた非線形最小2乗法のパラメーター。 |
得られたサーミスターの電気抵抗 \(R\) について、 熱力学温度の逆数に対して(Igor Pro のデフォルトの設定のまま、) 図 4-1a では電気抵抗 \(R\) の対数に対して線形最小2乗法を適用し、 図 4-1b では電気抵抗 \(R\) そのものに対して関数
\[ f(T) = R_{25} \exp{\left[ 1000 B (1/T - 1/298.15 \mrm{~K}) \right ]} \]
に非線形最小2乗法を適用しています。 両者から得られる B 係数は異なっていて、 次のようになっています:
カタログ値は 3250 K なので、非線形最小2乗法の結果はいささか偏倚が大きすぎます。 これは非線形最小2乗法でデータの重みを誤って取っているためです。 特にここで扱っているサーミスターの抵抗の実験データは、2 から 20 kΩ ぐらいまでおよそ 10 倍違い、 単純に処理すると、最大 100 倍くらいの重みの違いが出てきます。
こうした事態に対処できるように、 たいていの科学技術用のデータ解析ソフトには、重みの異なるデータへの対応措置が備えられています。 Igor Pro であれば、「Curve Fitting」の中の「Data Options」の項に「Weighting」 として、wave を指定することができるようになっています(デフォルトでは _none_)。 たとえばこの場合は相対誤差が一定だと見なせるので、 データの重みとして抵抗値に比例する wave を作成して、標準偏差の形で「Weighting」として指定すれば、 非線形最小2乗法でも、対数を取って重みを一定とした図 4-1a と同じ B = 3229.2 K という結果を得ます。
同様のことは反応速度の解析の際、 グッゲンハイム法などを適用する場合に起きてきます。 反応速度の場合には、サーミスターの電気抵抗とはちがって、 相対精度が一定ではなく、絶対精度を一定と見なせるケースが多く、 こことは話が逆になりますが注意が必要です。 たいていの場合、反応のせいぜい 3 半減期程度までを扱い、 変化量が半減する程度なのであまり気にする必要はありません。 けれどもノイズレベルが高い場合には大きな問題になり(下手をすると負の数の対数を取ることになる!)、 非線形の最小2乗法に従うのが薦められます。
なお多次元のデータを扱う一般的な場合には、重み因子の行列 \(W\) を考えると(\(W\)は対角要素が \(1/\sigma_i^2\) の対角行列)、 残差2乗和は \(S = {}^t (Y - X \vec{a}) W (Y - X \vec{a})\) で表され、正規方程式は次の形で表示できる:
\begin{equation} {}^t X W X \vec{a} = {}^t X W \vec{Y} \label{eq:weight3} \end{equation}
したがってパラメータ \(\vec{a}\) は、形式的には次式で与えられる:
\begin{equation} \vec{a} = ({}^t X W X)^{-1} {}^t X W \vec{Y} \label{eq:weight4} \end{equation}
実験条件 \(x\) には不確かさがないとしてきたが、 実際には実験条件にも測定値 \(y\) と同様に不確かさが見込まれる場合も多い。 こうした場合の関係式のパラメータ推定について考えてみよう。
実験条件 \(x\) の不確かさ \(\epsilon_x\) が実験条件と独立であれば、 測定値の分散を調整すれば、前章までの結果はほとんどそのまま成立する (\(\epsilon_x\) も平均 0 の正規分布に従うとする)。 つまり最初に想定する実験条件 \(x\) と測定結果 \(y\) のモデルが次式のように表されるので (測定結果の不確かさを \(\epsilon_y\) とする)
\begin{equation} y = a (x + \epsilon_x) + b + \epsilon_y = a x + b + (a \epsilon_x + \epsilon_y) \label{eq:cond1} \end{equation}
測定結果の分散 \(\sigma^2\) として、次式を取ればよい:
\begin{equation} \sigma^2 = \var{\epsilon^2} = \var{(a \epsilon_x + \epsilon_y)^2} = a^2 \var{\epsilon_x^2} + \var{\epsilon_y^2} \label{eq:cond2} \end{equation}
実験条件の側に不確かさがあっても、 それが独立に振舞うなら、パラメーターの値には影響を及ぼさないというのは、 線形最小2乗法の大事なポイントです。
ただしこの「不確かさが独立に振舞う」というのがポイントで、 化学でよく問題になるのは、多項式への当てはめを行う場合です。 たとえば次のような2次式に当てはめる時
\[ y = a_0 + a_1 x + a_2 x^2 + \epsilon_y \]
\(x\) に不確かさが存在すると
\[ y = a_0 + a_1 (x + \epsilon_x) + a_2 (x + \epsilon_x)^2 + \epsilon_y \approx a_0 + a_1 x + a_2 x^2 + (a_1 + 2 a_2 x) \epsilon_x + \epsilon_y \]
となって、実験データの不確かさ(データの重み)が実験条件に依存するようになります。
\[ \sigma^2 \approx (a_1 + 2 a_2 x)^2 \var{\epsilon_x^2} + \var{\epsilon_y^2} \]
多項式への当てはめには慎重でないといけないということがしばしばいわれますが、 その根拠の一つはこうしたところにあります。