実際に最小２乗法を適用した例を眺めてみましょう。

空気の平均分子量を求めるために、 ２回生の実験で圧力を変化させてフラスコの重さをはかった結果を表 2-1 に示します。 この実験ではダイアフラムポンプを利用して 10 kPa 程度まで減圧しています。 実験は室温 22 °C で行い、装置の内容積は減圧栓部分などの死容積を含め 680 mL でした。 フラスコの一件の重さは 300 g 程度ありますが、 最近は 500 g まで 0.01 g までの精度ではかれるロードセルタイプの電子天秤が安価に入手できるようになり、 こうした実験が手軽にできるようになりました。

表 2-1. 圧力変化にともなうフラスコの質量の変化

P/kPa	Δw/g
101.1	0.00
80.2	-0.18
59.9	-0.35
39.6	-0.50
20.1	-0.66
8.8	-0.76

図 2-1. 空気の重さの測定。 500 mL のナス型フラスコを数 kPa まで減圧した後、 少しずつ空気を入れては圧力と重さをはかります。 ゴム管を折ってクリップで挟みこみ、圧力を止めています。 1.5 L のペットボトルを輪切りにして作った台座の上に、 フラスコは乗っています。

今日では、こうしたデータの解析は市販の科学技術用ソフト（ここでは Igor pro。excel でも余裕）を利用し、 図 2-2 のようにまずプロットして直線になることから、 理想気体の状態方程式 \(PV = nRT\) が成り立っていることを確認。 組み込みの解析機能を使って最小２乗法を適用した近似直線を得て、 空気の質量の圧力依存性を定めるまで、せいぜい数分で済む作業になりました。

図 2-2. 空気の重さの圧力依存性。

表 2-2. 圧力変化にともなうフラスコの質量の変化に関わる統計量：\(x = P / \mrm{kPa}\)、\(y = \Delta w / \mrm{g}\)

\(N\)	6	\(S_{xx}\)	22291
\(S_x\)	309.7	\(S_{xy}\)	-75.16
\(S_y\)	-2.45	\(S_{yy}\)	1.4181

\[ a = \frac{6 \x (-75.16) - 309.7 \x (-2.45)}{6 \x 22291 - 309.7^2} = 8.136 \x 10^{-3} \] \[ \frac{\rmd w}{\rmd P} = \frac{M V}{RT} = 8.136 \x 10^{-3} ~\mrm{g ~kPa^{-1}} \] \begin{eqnarray*} M &=& \frac{8.314 ~\mrm{J ~mol^{-1}~ K^{-1}} \x 295 ~\mrm{K}}{680 \x 10^{-6}~\mrm{m^3}} \x 8.136 \x 10^{-3} ~\mrm{g~ kPa^{-1}} \\ &=& 29.3 ~\mrm{g~mol^{-1}} \end{eqnarray*}

この市販のソフトの中で行われている演算を、 実際に確認してみましょう。 表 2-2 には、実験データから得られる、最小２乗法に必要な統計量をまとめました。 ここから図 2-2 の直線の勾配を、式 \eqref{eq:param5} にしたがって求めると 8.136 × 10^-3 g/kPa になり、 図 2-2 中のテキストボックスにある解析結果通りの値となります （± 以下の部分の計算については後述。 なお Igor では回帰分析でよく用いられる y = α + β x という記法にならって、y = a + b x としているようです）。 瞬時にしてこうした結果が出るのですが、 その背景でどんな計算が行われているかがお分かりいただけるでしょう。 そしてこうして求めた直線の勾配と理想気体の状態方程式から、 空気のモル質量が 29.3 g/mol と得られます（２回生の実験では、単位を付けて折り目正しく計算することを期待しています）。

最小２乗法での処理の中味を、行列の演算の形で書いたわけですが、 式の言わんとするところをイメージしやすいように、 簡単な絵にしてみましょう。 残差の式 \eqref{eq:param9a} は右図のような行列の計算になっています。

測定値 \(\vec{Y}\)、残差 \(\vec{R}\) は、\(N\) 個の要素からなるベクトルで、 それぞれに対応する \(m\) 個の実験条件は \(N\) 行 \(m\) 列の行列 \(X\) の各行 \({}^t \vec{X_i}\)（\(i = 1, 2, \ldots, N\)）で与えられ、 理論値（計算値）は、各条件と \(m\) 個の要素からなるパラメータのベクトル \(\vec{a}\) の内積の形で与えられます （統計モデルの分野などではここで ”実験条件” としている行列 \(X\) を計画行列 design matrix、 また残差の式 \eqref{eq:param9a} に相当するものを観測方程式 observation equation と呼ぶことが多いようです）。

残差２乗和は形式的に少し丁寧に書くと次のように表せます （このはなしでは 1 行 1 列の行列とスカラーは同一のものと見なします。 いろいろいい加減ですが、何分よろしく）：

\[ S = {}^t (\vec{Y} - X \vec{a}) (\vec{Y} - X \vec{a}) = \vec{Y}^2 - {}^t\vec{Y} X \vec{a} - {}^t[X \vec{a}] \vec{Y} + {}^t[X \vec{a}] X \vec{a} = \vec{Y}^2 - 2 {}^t \vec{Y} X \vec{a} + {}^t \vec{a} ~[{}^tX X] \vec{a} \]

残差２乗和の \(j\) 番目のパラメータ \(a_j\) についての偏微分を取ると、 次式のようになります（\(S\) を 1 行 1 列の行列と見なせば、 転置を取っても同じで \({}^t \vec{Y} X \vec{I_j} = {}^t [{}^t \vec{Y} X \vec{I_j}] = {}^t \vec{I_j} {}^t X \vec{Y}\) となることなどに注意）：

\[ \pdif{S}{a_j} = - 2 {}^t \vec{Y} X \vec{I_j} + {}^t \vec{I_j} ~[{}^tX X] \vec{a} + {}^t \vec{a} ~[{}^tX X] \vec{I_j} = - 2 {}^t \vec{I_j} {}^t X \vec{Y} + 2 {}^t \vec{I_j} ~[{}^tX X] \vec{a} = 0 \]

ここで \(\vec{I_j}\) は \(j\) 番目の要素だけが 1 で他は 0 のベクトルです。 \({}^t \vec{I_j}\) を左からかけるのは、行列の \(j\) 列目を取り出すことに相当します。 これをベクトルの形で整理したのが、正規方程式 \eqref{eq:param10} になります。 正規方程式もイメージしやすいように、 簡単な絵にしてみると次のようになります：

正規方程式でパラメータに掛かる \(m \x m\) の行列 \({}^tX X\) の \((k, l)\) の要素は、行列 \(X\) の \(i\)（= \(1, 2, \ldots, m\)）番目の列である、 各測定における実験条件を要素とするベクトル \(\vec{X_i}\)（= \((X_{1,i}, \ldots, X_{N,i})\)）の間の内積 \(\vec{X_k} \cdot \vec{X_l}\) になります。 先の式 \eqref{eq:param3}, \eqref{eq:param4} の正規方程式では、 実験条件を \(\vec{x} = (x, 1)\)、パラメータ \(\vec{a} = (a, b)\) と取ったことになっています。

なお残差２乗和の最小値を与える式 \eqref{eq:param12} は、 得られたパラメータ \(\vec{a}\) の式 \eqref{eq:param11} を用い

\[ {}^t \vec{a} ~[{}^tX X] \vec{a} = {}^t \vec{Y} X [({}^t X X)^{-1}] ~[{}^tX X] \vec{a} = {}^t \vec{Y} X \vec{a} \]

であることに注意すれば容易に導くことができます。

このおはなしでは扱いませんが、 パラメータの数が増えてくると、 種々の条件付きの最小２乗法を行う必要に迫られる場合が出てきます。 たとえばさまざまな組成についての実験データから、 混合液体の沸点 \(t\) について、 よく知られている純粋な成分 \(i\) の沸点 \(t_i^\bullet\) を与えるように、 沸点と組成の実験式を定めるようなケースです。 表 2-3 には３回生の実験で登場するシクロヘキサン-トルエン混合液体の 1 atm（= 101.325 kPa）における沸点のデータを示しました （化学便覧掲載のデータから一部抜粋）。

図 2-3. シクロヘキサン-トルエン系の 1 atm における沸点の組成依存性。 青い実線は純成分の沸点からモル分率についての直線関係を想定して描いたもの。

シクロヘキサン (1)-トルエン (2) 混合物の 1 atm における沸点（通常沸点）は、 シクロヘキサンのモル分率 \(x\) の増加とともに低下していきますが、 図 2-3 に見るようにモル分率に単純に比例するわけではありません。 実験値を再現するためには \(x\) の高次の項が必要で、 変動の様子からは 2 次以上、3 次ぐらいまでの項が必要なようです。

さて沸点 \(t\) をシクロヘキサンのモル分率 \(x\) の多項式で表そうとした時、 純物質の沸点のデータは信頼性が高いので、 それを再現するように、 \(t\) / °C を次のような関数 \(f(x)\) に当てはめることが考えられます：

\[ f(x) = f(1) x + f(0) (1-x) + a x (1-x) + b x^2 (1-x) = f(1) x + f(0) (1-x) + a p_1(x) + b p_2(x) \]

ここで表 2-3 から \(f(0)\) = 110.65, \(f(1)\) = 80.70 であり、また

\[ p_1(x) = x (1-x), ~~~p_2(x) = x p_1(x) \]

とおきました。 つまり混合物の沸点についてモル分率から推定される直線関係

\[ t_\mrm{lin} = x t^\bullet_1 + (1-x) t^\bullet_2 \]

からの偏倚を

\[ (t - t_\mrm{lin}) / {}^\circ \mrm{C} = a p_1(x) + b p_2(x) \]

という関数に当てはめ、 係数 \(a\)、\(b\) を実験を再現するように最小２乗法で決定しようというわけです。 混合物の沸点の直線関係からの偏倚を拡大して示すと図 2-4 のようになり、 正規方程式は式 \eqref{eq:param3}-\eqref{eq:param4} と同断です。 表 2-4 に必要な統計量をまとめ、 係数 \(a\)、\(b\) の計算も示しておきます：

図 2-4. シクロヘキサン-トルエン混合物の通常沸点の直線関係からの偏倚と実験式

表 2-4. シクロヘキサン-トルエン混合物の通常沸点の直線関係からの偏倚に対する実験式に関わる統計量。 分かりやすいように \(S_q\) を \([q]\) と表記。

\([p_1 p_1]\)	0.1634	\([y p_1]\)	-3.466
\([p_1 p_2]\)	0.0824	\([y p_2]\)	-1.648
\([p_2 p_2]\)	0.0471	\([y y]\)	75.33

\[ a = \frac{(-3.466) \x 0.0471 - (-1.648) \x 0.0824}{0.1634 \x 0.0471 - 0.0824^2} = -30.3 \] \[ b = \frac{0.1634 \x (-1.648) - 0.0824 \x (-3.466)}{0.1634 \x 0.0471 - 0.0824^2} = 18.0 \]

得られる実験式は、最終的に次のようになります （学生実験のテキスト所載の式と係数が若干異なるのは、 当てはめに用いたデータ数が異なるため）：

\begin{eqnarray*} t / {}^\circ \mrm{C} &=& 80.7 x + 110.65 (1 - x) - 30.3 x(1 - x) + 18.0 x^2 (1 - x) \\ &=& 80.7 x + 110.65 (1 - x) - 30.3 x(1 - x)(1 - 0.6 x) \end{eqnarray*}

このようにして、条件に合うように式の形を設定して最小２乗法を適用し、 シクロヘキサン-トルエン混合物の通常沸点に対する実験式を得ることができました。 ここでは一歩一歩、段階を踏んで実験式を構成する手続きを紹介しましたが、 たいていの市販の科学技術用の汎用のソフトであれば、 当てはめの関数形を与えてやれば、直ちに最小２乗法による解を与えてくれます。 もっとも計算の手法はもっと手が込んでいて、 たとえば Igor Pro であれば Curve Fitting で New Fit Function に上記の関数 \(f(x)\) を設定し、 パラメータ a、b の初期値を与え、逐次計算をする形になります（非線形最小２乗法）。

なお条件が複雑になってくると、一筋縄でいかないケースも出てきます。 そうした時には条件付きの極値問題の定石ともいえる、 ラグランジュの未定乗数法を用いるのが有効になってきます。 化学で出会う問題では、たいていの場合、そこまでする必要がないでしょうが・・・