測定値 \(\{y_i\}\) の分散を \(\sigma^2\) とすると、 最小 2 乗法で定めたパラメータ \(a\) の分散は次のように求めることができる。 まず \(a\) の表現を整理して
\begin{equation} a = \frac{N S_{xy} - S_x S_y}{N S_{xx} - S_x^2} = \sum_{i=1}^N {\left[ \left( \frac{N x_i - S_x}{N S_{xx} - S_x^2} \right) \right] y_i} \label{eq:ci1} \end{equation}
とすれば、 \(y_i\) は互いに独立なので \(a\) の分散は
\begin{equation} \var{a^2} = \sum_{i=1}^N {\left( \frac{N x_i - S_x}{N S_{xx} - S_x^2} \right)^2 \sigma^2} = \frac{\sigma^2}{(N S_{xx} - S_x^2)^2} \sum_{i=1}^N {(N x_i - S_x)^2} \label{eq:ci2} \end{equation}
ここで次の関係に注意して
\begin{equation} \sum_{i=1}^N {(N x_i - S_x)^2} = N^2 S_{xx} - 2 N S_x^2 + N S_x^2 = N (N S_{xx} - S_x^2) \label{eq:ci3} \end{equation}
勾配 \(a\) の分散として次の関係を得る (\(\var{x^2}_\mrm{s} = \overline{x^2} - \bar{x}^2\)):
\begin{equation} \var{a^2} = \frac{N}{N S_{xx} - S_x^2} \sigma^2 = \frac{1}{N \var{x^2}_\mrm{s}} \sigma^2 \label{eq:ci4} \end{equation}
同様にして
\begin{equation} \var{b^2} = \frac{S_{xx}}{N S_{xx} - S_x^2} \sigma^2 = \frac{1}{N} \left( 1 + \frac{\bar{x}^2}{\var{x^2}_\mrm{s}} \right) \sigma^2 = \frac{\sigma^2}{N} + \var{a^2} \bar{x}^2 \label{eq:ci5} \end{equation} \begin{equation} \var{ab} = -\frac{S_x}{N S_{xx} - S_x^2} \sigma^2 = - \frac{\bar{x}}{N \var{x^2}_\mrm{s}} \sigma^2 \label{eq:ci6} \end{equation}
を得ることができる。
推定パラメータの信頼区間の幅はその標準偏差に比例すると考えられ、 繰り返し実験の場合同様、測定値の標準偏差 \(\sigma\) に比例し、 データ点の数の平方根 \(N\) に反比例する。 測定値の分散 \(\sigma^2\) が分かっていないときには次式で推定することになる:
\begin{equation} \sigma^2 \approx \frac{S_\mrm{min}}{N - 2} = \frac{S_{yy} - (a S_{xy} + b S_y)}{N-2} \label{eq:ci7} \end{equation}
先に得たパラメータの表式から計算を示しましたが、 実験条件 \(x\) を
\[ x^* = x - \bar{x} \]
とおいて \(x\) の標本平均 \(\bar{x} = S_x/N\) だけシフトさせるともっと見通しが良くなります。 以下、こうしてシフトさせて得られる量には * を付けて表記します。
このようにシフトさせることで \(\bar{x^*} = S^*_x/N = 0\) となり、 得られる最小2乗法のパラメータが次のように、 測定値(ランダム変数)\(y_i\) の線形結合になっていることがはっきりします:
\[ a^* = a = \sum_{i=1}^{N} \frac{x^*_i}{S^*_{xx}} y_i, ~~~ b^* = b + a \bar{x} = \sum_{i=1}^{N} \frac{1}{N} y_i \]
\(y_i\) は独立なランダム変数で、\(\var{y^2_i} = \sigma^2\) ですから、 ここから \(a^*\) と \(b^*\) の(共)分散を評価すると次のようになります:
\[ \var{a^{*2}} = \sum_{i=1}^{N} \frac{x^{*2}_i}{S^{*2}_{xx}} \sigma^2 = \frac{S^*_{xx}}{S^{*2}_{xx}} \sigma^2 = \frac{\sigma^2}{S^*_{xx}} \] \[ \var{a^* b^*} = \sum_{i=1}^{N} \frac{x^*_i}{N S^{*}_{xx}} \sigma^2 = 0 \] \[ \var{b^{*2}} = \frac{\sigma^2}{N} \]
パラメータ \(a^*\) と \(b^*\) は無相関で独立に振舞います。
さて \(x\) をシフトさせる前の元のパラメータについての(共)分散は次のように求まり、 先の結果と一致します:
\[ \var{a^2} = \var{a^{*2}} = \frac{N}{N S_{xx} - S^2_x} \sigma^2 \] \[ \var{a b} = \var{a^* (b^* - a \bar{x})} = - \bar{x} \var{a^2} = -\frac{S_x}{N S_{xx} - S^2_x} \sigma^2 \] \[ \var{b^2} = \var{(b^* - a \bar{x})^2} = \var{b^{*2}} + \bar{x}^2 \var{a^2} = \frac{\sigma^2}{N} + \bar{x}^2 \frac{\sigma^2}{S_{xx} - N \bar{x}^2} = \frac{S_{xx}}{N S_{xx} - S^2_x} \sigma^2 \]
なお測定データの分散を評価する式 \eqref{eq:ci7} で \(N - 2\) で割っているのは、 不偏推定にするためです。 測定値が次のように揺らぐものとします:
\[ y_i = \avg{y_i} + \delta y_i \]
ここで \(\delta y_i\) は平均 0 で分散 \(\sigma^2\) のランダム変数であるとします。 また測定値の平均は、パラメーターの平均と次の関係にあると想定して最小2乗法を行っています。
\[ \avg{y_i} = \avg{a^*} x_i^* + \avg{b^*} \]
このように置くと残差 \(R_i\) は
\[ R_i = y_i - (a^* x_i^* + b^*) = \delta y_i - (\delta a^* x_i^* + \delta b^*) \]
ここで \(\delta a^*\)、\(\delta b^*\) はそれぞれ次のように表されます:
\[ \delta a^* = \sum_{i=1}^{N} \frac{x^*_i}{S^*_{xx}} \delta y_i, ~~~ \delta b^* = \sum_{i=1}^{N} \frac{1}{N} \delta y_i \]
2乗和の最小値は、次のように書くことができます;
\[ \avg{S_\mrm{min}} = \avg{\sum_{i=1}^{N} y_i R_i} = \sum_{i=1}^{N} \avg{y_i} \avg{R_i} + \sum_{i=1}^{N} \avg{\delta y_i R_i} = \sum_{i=1}^{N} \avg{\delta y_i R_i} \]
ここで測定値のゆらぎについて
\[ \avg{\delta y_i \delta y_j} = \left\{ \begin{array}{rl} & \sigma^2 & (i = j) \\ & 0 & (i \ne j) \end{array} \right. \]
という関係が成立するので
\[ \avg{\delta y_i R_i} = \avg{(\delta y_i)^2} - x_i^* \avg{\delta a^* \delta y_i} - \avg{\delta b^* \delta y_i} = \sigma^2 - \frac{x^{*2}_i}{S^*_{xx}} \sigma^2 - \frac{1}{N} \sigma^2 \]
となり、2乗和の最小値の平均値は次式で表され、
\[ \avg{S_\mrm{min}} = \sum_{i=1}^{N} \avg{\delta y_i R_i} = N \sigma^2 - \sum_{i=1}^{N} x^{*2}_i \frac{\sigma^2}{S^*_{xx}} - \sigma^2 = (N-2) \sigma^2 \]
となり、式 \eqref{eq:ci7} が不偏推定になっていることがわかります (残差2乗和の最小値は自由度 \(N - 2\) の χ2 分布に従うと考えられます)。
先にも触れたようにあくまで最ももっともらしい値にこだわる立場(最尤推定)からは、 \(N\) で割ることになりますが、 標本分散と同様、不偏分散をとるのが一般です。
先の空気の平均分子量測定の実験の例で、 計算を確認してみると次のようになります:
\[ \sigma^2 \approx \frac{S_{yy} - (a S_{xy} + b S_y)}{N-2} = \frac{1.4181 - [8.136 \x 10^{-3} \x (-75.16) + (-0.8283) \x (-2.45)]}{6-2} = 5.09 \x 10^{-5} \] \[ \sigma^2_a = \var{a^2} \approx \frac{6}{6 \x 22291 - 309.7^2} \sigma^2 = 1.59 \x 10^{-4} \x 5.09 \x 10^{-5} = 8.09 \x 10^{-9} \] \[ \sigma_a \approx 9.0 \x 10^{-5} \]
この結果が先の実験の結果を Igor で解析して得た直線の勾配 「0.0081378 ± 8.98 e-005」の ± 以下の数値 8.98 e-005(図 2-2 の係数の Textbox)になっているわけです。 市販のデータ解析ソフトで瞬時に出てくる数値の裏でどのような処理が行われているか、 お分かりいただけたでしょうか。
ここで評価しているのはあくまで分散、標準偏差であることに注意してください。 ここからパラメータの信頼区間を評価するには、 有意水準(危険率)を設定し、t 分布を用いたりといった検定の手続きを踏む必要があります。 化学では多くの場合、そういった手続きは後回しにして、 評価された標本標準偏差を信頼区間の幅であるかのように扱うことが多いです (Igor でも素知らぬ顔をして、標本標準偏差を信頼区間のように ± で表示しています)。 一般に信頼区間を \(k \sigma\) と表示し、\(k\) を包含係数と呼びます。 化学でも分野によっては(難しいことは抜きにして)、\(2 \sigma\) を信頼区間として表示することがあります。 これは有意水準をおよそ 5 % に取ったことに相当します。
測定データを \(y = ax + b\) に当てはめてパラメータ \(a\)、\(b\) を推定する時、 測定データが多い(\(N\) が大きい)ことは無論のこと、 勾配 \(a\) の評価にはできるだけ \(x\) の範囲を広くとって測定し(\(\var{x^2}_\mrm{s}\) を大きくする)、 切片 \(b\) の評価には \(x = 0\) 付近の値を取る(\(|\bar{x}|\) を小さくする)のが望ましい。
なお得られるパラメータ \(a\), \(b\) の不確かさは一般に独立でない(\(\var{ab} \ne 0\))。 特に原点から \(x\) 方向に遠く離れたデータ(\(|\bar{x}|\)が大きい)を用いてパラメータ \(b\) を推定するときには、 信頼区間の評価には注意が必要である。
最小2乗法のパラメータの不確かさに関わって、 平衡定数の温度変化から、反応にともなうエンタルピー変化、エントロピー変化を評価するのに用いられる、 ファントホッフプロットを見てみましょう。
ニトロソベンゼン類は、溶液中で単量体(M)と二量体(D)の間で平衡状態にあることが知られています。
2M ⇌ D K = [D]/[M]2
表 3-1 に示すのは、2,4,6-トリメチルニトソベンゼン(TMNB)の四塩化炭素中 CCl4 における二量化反応の平衡定数 \(K\) を、 いくつかの温度で吸光光度法を用いて測定した結果です。 平衡定数(標準状態の濃度 1 mol/L)は自然対数を取った形で示してあり、標準偏差が 0.02(平衡定数 K としては標準偏差が 2 %)、 温度の不確かさについては、ここではないものとします (取り扱い方については、後の章を見て下さい)。
| T / °C | 1000 K/T | ln K |
|---|---|---|
| 20.0 | 3.411 | 3.58 |
| 25.0 | 3.354 | 3.24 |
| 30.0 | 3.299 | 2.89 |
|
| 2,4,6-トリメチルニトロソベンゼン(TMNB)の二量化平衡 |
|
| 図 3-1a TMNB の二量化平衡定数のファントホッフプロット。 ln K の標準偏差を 0.02 として最小2乗法を行っています。 |
平衡定数の温度依存性を評価するため、 ファントホッフプロットを行い、 平衡定数の対数を熱力学温度の逆数に対してプロットした結果を、 図 3-1a に示します。 この結果を
ln K = a (1000 K/T) + b
という直線に最小2乗法で当てはめると、次のような結果を得ます;
平衡定数の対数 ln K の標準偏差が 0.02 であるのに、 それを推定するパラメーターの標準偏差が 10 倍以上もありますが、 これはパラメーターの間に強い相関があるためです。 たとえばこの結果を用いて 25 °C の平衡定数を推定したとしましょう:
ln K = 6.159 (1000 K/298.15 K) - 17.42 = 3.24
|
| 図 3-1b. 原点付近まで含めたTMNB の二量化平衡定数のファントホッフプロット。 |
表 3-1 の実験結果と同じ値になりますが、この推定の不確かさはどの程度になるでしょうか? x = 1000 K/298.15 K = 3.354 とおいて、 分散を求めてみると次のようになります:
⟨⟨(ln K)2⟩⟩
= x2 ⟨⟨a2⟩⟩
+ 2 x ⟨⟨ab⟩⟩ + ⟨⟨b2⟩⟩
= 0.7106 - 1.4214 + 0.7111 = 0.0003
大きな数同士の差し引きの結果、有効数字の最後のケタにようやく引っかかる値になります。 これでは心もとないので、もう 1 ケタ有効数字を取って計算し直すと、 0.00013、標準偏差でいうと 0.011 になります (独立な 3 つのデータの平均として分散が 1/3 になる)。
このような結果になるのは、先にも述べたようにパラメータの間の相関が強いためですが、 直観的には図 3-1b にみるように、 原点からはるか離れた地点から、切片の値を求めようとしているところにあります。 わずかな勾配のゆらぎ(σa2 = ⟨⟨a2⟩⟩)が、 大きな切片のゆらぎをもたらします。 このことは式 \eqref{eq:ci5} からも明らかで、 切片の分散は次式のように与えられます:
⟨⟨b2⟩⟩ = σ2/3 + x2 ⟨⟨a2⟩⟩ = 0.022/3 + 3.3552 × 0.06317 = 0.00013 + 0.7110 = 0.7111
つまり σb ≈ x σa で、 切片のゆらぎはほとんどが勾配のゆらぎに起因しています。 実際にこうした結果をわれわれが目にするのは、 反応のエンタルピー変化、エントロピー変化という形を通してでしょう。 反応の標準エンタルピー変化 \(\Delta_\mrm{r} H^\circ\) は、次式で与えられます:
\[ \Delta_\mrm{r} H^\circ/R = -\pdifA{\ln K}{(1/T)}{P} - \alpha T^2 = - (1000~\mrm{K}) a - \alpha T^2 \]
ここで \(R\) は気体定数、\(\alpha\) は溶媒の膨張率です(四塩化炭素では 1.2 × 10-3 K-1)。 また反応の標準エントロピー変化 \(\Delta_\mrm{r} S^\circ\) は、 \(\Delta_\mrm{r} G^\circ = -RT \ln K = \Delta_\mrm{r} H^\circ - T \Delta_\mrm{r} S^\circ\) の関係から次のようになります:
\[ \Delta_\mrm{r} S^\circ/R = b - \alpha T \]
こうした関係式を用い、 下記のような結果を得ます(± 以下は標準偏差):
反応のエンタルピー変化、エントロピー変化、それぞれに標準偏差が表示されていますが、 上に見たように、実質、 エントロピー変化の標準偏差はエンタルピー変化の標準偏差を温度(298.15 K)で割ったものになります。 ちなみにエンタルピー変化に対する溶媒の膨張に起因する項 αRT2 の寄与は 0.9 kJ mol-1 で、標準偏差よりも小さいものに止まっています(エントロピー変化についても同断です)。
実験条件が \(m\) 個の要素で与えられる場合に、パラメータ \(\vec{a}\) の共分散行列 \(\var{\vec{a} {}^t \vec{a}}\) は次式で与えられる
\begin{equation} \var{\vec{a} {}^t \vec{a}} = ({}^t X X)^{-1} {}^t X \var{\vec{Y} {}^t \vec{Y}} X ({}^t X X)^{-1} = ({}^t X X)^{-1} \sigma^2 \label{eq:ci8} \end{equation}
また実験データの分散が分かっていない時には、残差2乗和の最小値を用いて、次式で分散を推定できる:
\begin{equation} \sigma^2 \approx \frac{S_\mrm{min}}{N - m} = \frac{{}^t \vec{Y} (\vec{Y} - X \vec{a})}{N-m} \label{eq:ci9} \end{equation}
この章でパラメータの分散を評価する時に、 係数がバタバタと整理され、 結構さっぱりした形になったのには、 式 \eqref{eq:ci8} に整理される背景があったわけです。 実験条件を \(\vec{x} = (x, 1)\)、パラメータ \(\vec{a} = (a, b)\) と取れば、
\[ {}^t X X = \matx{S_{xx}}{S_x}{S_x}{N} \]
なので、共分散行列は次式のようになり、式 \eqref{eq:ci4} - \eqref{eq:ci6} が確認できます:
\[ \matx{\var{a^2}}{\var{ab}}{\var{ab}}{\var{b^2}} = [{}^t X X]^{-1} \sigma^2 = \frac{1}{N S_{xx} - S_x^2} \matx{N}{-S_x}{-S_x}{S_{xx}} \sigma^2 \]
実験値の分散の評価に用いる式 \eqref{eq:ci9} は、 以前にも示したことがありますが、 次のようにして導くことができます。 実験値がその平均値の周りに次のようにゆらいでいるとします:
\[ \vec{Y} = \avg{\vec{Y}} + \delta \vec{Y} \]
ここで \(\delta \vec{Y}\) の成分 \(Y_i\) は、互いに独立で、平均 0 で分散 \(\sigma^2\) であるとします。 するとパラメータ \(\vec{a}\) の平均について、 次の関係が成り立っていると想定していることに注意して
\[ \avg{\vec{Y}} = X \avg{\vec{a}} \]
残差は次のように表されます:
\[ \vec{R} = (\vec{Y} + \delta \vec{Y}) - X (\vec{a} + \delta \vec{a}) = \delta \vec{Y} - X \delta \vec{a} \]
したがって残差2乗和の最小値は次のように表されます
\[ S_\mrm{min} = {}^t \vec{Y} (\delta \vec{Y} - X \delta \vec{a}) \]
この平均を求めると
\[ \avg{S_\mrm{min}} = \avg{{}^t \delta \vec{Y} \delta \vec{Y}} - \avg{{}^t \delta \vec{Y} X \delta \vec{a}} \]
\(\delta y_i\) は互いに独立なので、右辺第 1 項は次のようになり
\[ \avg{{}^t \delta \vec{Y} \delta \vec{Y}} = N \sigma^2 \]
第 2 項は少し厄介ですが、 パラメーター \(\vec{a}\) の解の形をそのまま代入すると次のようになり
\[ \avg{{}^t \delta \vec{Y} X \delta \vec{a}} = = \avg{{}^t \delta \vec{Y} X ({}^t X X)^{-1} {}^t X \delta \vec{Y}} \]
これは \(\delta y_i\) に関する 2 次形式になりますが、 \(\delta y_i\) は互いに独立なので平均を取ると、 \(\avg{\delta y_i^2} = \sigma^2\) の項のみが残ります。
\[ \avg{{}^t \delta \vec{Y} X ({}^t X X)^{-1} {}^t X \delta \vec{Y}} = \mrm{Tr} [X ({}^t X X)^{-1} {}^t X] \sigma^2 \]
ここで行列 \(X ({}^t X X)^{-1} {}^t X\) の対角項の和、トレースが現れますが、 N 行 m 列と m 行 N 列の行列の積のトレースは、 順序の交換をしても変わらないので、次のように簡単になります:
\[ \mrm{Tr} [X ({}^t X X)^{-1} {}^t X] \sigma^2 = \mrm{Tr} [({}^t X X)^{-1} ({}^t X X)] \sigma^2 = m \sigma^2 \]
そして式 \eqref{eq:ci9} を得ます:
\[ \avg{S_\mrm{min}} = (N - m) \sigma^2 \]