Microsoft Excelを用いたケモメトリクス計算(1)
−行列の扱いと重回帰−

吉村 季織, 高柳 正夫


Return

1 はじめに

実験的に得られた化学データを数学的・統計的手法により解析する「ケモメトリクス」が、近年盛んに用いられるようになってきた。これには、PCの高性能化と普及が大きく寄与していると考えられる。ケモメトリクスに関する教科書類がいくつか[1 - 4]出版されている一方で、大学等の授業ではほとんど扱われていないのが実情である。この原因として、ケモメトリクスを教えることができる教員が少ないことだけでなく、教育用に用いることのできる手軽なケモメトリクスのソフトウェアがないことが考えられる。
ケモメトリクスで最も一般的に用いられる多変量解析は、線形代数学を基本としているため、行列の計算を行わなくてはならない。基本的な行列の取り扱いは、線形代数学の教科書や授業で習得することができる。そしてケモメトリクスの教育では、計算の意味や流れの解釈に加えて、実際に計算が行えるようにする必要がある。ケモメトリクスで実用的な大きさの行列を取り扱うために、高価なMATLABなどではなく、Microsoft Excel(以後Excelと記す)を利用することができれば、教育・研究の場に広くケモメトリクスを普及させることに役立つものと思われる。
そこで本稿および続報では、ケモメトリクス計算をExcelによって行う方法について解説していく。計算の流れを明快にし、多くのExcelユーザーにとってなじみの薄い機能は使わないようにするために、Visual Basic for Applications(VBA)や計算結果のみを示す機能を使うことなく、ワークシート上で計算を展開させる。
本稿では、まずExcelを使って行列計算を行うために必要な技術、主に配列と配列数式(array formula)について解説する。続いて、ケモメトリクスの基本的な内容のうち、重回帰を用いた曲線フィッティングを行う。さらに進んだ内容に関しては、続報にて報告する。
本稿の内容は、Microsoft Windows上のExcelを想定している。そのため、Mac OS上のExcelの場合、キー操作など適宜読み替えてほしい。

2 Excelによる行列−配列−の扱い

2. 1 行列と配列

ケモメトリクスでは行列が用いられる。そこで、Excelで行列を扱う方法をまず紹介する。Excelによる行列計算法に関しては、吉村忠与志らの著書[5, 6]に他の科学技術計算とあわせて紹介されているので、参考にされたい。本稿だけでもExcelによる行列計算を習得できるようにするため、および今後のシリーズにおいて基礎となるため、以下よりExcelによる行列計算について解説する。
Excelをはじめとする表計算ソフトウェアは、2次元的に敷き詰められたセルより構成されている。セルの横方向の並びを"行"、縦方向の並びを"列"と称するように、このセルの並びそのものが行列であるといえる。Figure 1にA1からD4までの4行4列に作成した行列の例を示した。とは言っても、ただ数値を並べただけであり、特別な定義を必要とするわけではない。Excelではこのように特定の範囲のセルをまとめて扱うことができ、"配列"と呼ぶ。つまり、Excelで行列を扱うということは、配列を扱うこと(の一部)である。
Excelでの配列の表記は、左上のセルと右下のセルのアドレスを":"で結合し、Figure 1を例とすれば"A1:D4"のように表す。配列の利用の最も一般的な例の一つとして、合計を求めるSUM関数がある。Figure 1のような複数の数値の合計を求める場合、本来なら合計値を求めたいセルに"=SUM(A1,A2,A3,...,D4)"のように全てのセルを指定して記入する必要がある。しかし、ほとんどのユーザーが"=SUM(A1:D4)"のように、配列を意識せずに配列を使っているはずである。


Figure 1. A simple example of the array expressed on an Excel worksheet.

2. 2 配列数式の例

Excelでは、配列を用いて計算を行う数式や、配列を出力する数式を配列数式と呼ぶ。配列数式の最も簡単な例として、何もせずそのままの配列を返す配列数式から説明する。
Figure 1のA1からD4までの範囲内の数値を別の場所(たとえばA6からD9までの範囲)に参照する場合、A6セルに"=A1"と入力してA1セルの内容を参照した後、A6セルをコピーしてD9セルまでの4行4列に貼り付けする手順が一般的である。この場合は相対参照が有効になり、D9セルの数式は"=D4"となる。この参照を配列数式で行うには、Figure 2のような手順を用いる。まず、Figure 2(a)のように、A6からD9までのセルを選択する(例の場合、A6からD9へドラッグし選択を行った)。Figure 2(b)のように、アクティブセル(例の場合はA6セル)へ"=A1:D4"を記入する(ただし、まだ[Enter]キーは押してはならない)。この数式は、配列をそのまま記入しているだけであるが、最も簡単な配列数式の例(配列をそのまま出力するだけ)である。ここで、[Enter]キーを押す代わりに、[Ctrl]キーと[Shift]キーを押した状態で[Enter]キーを押すことで入力を確定する(この「[Ctrl]キーと[Shift]キーを押した状態で[Enter]キーを押すという入力確定法」を本稿では以降「配列入力」と呼び、「[Enter]キーのみを押す」通常の入力確定法と区別する)。そうすると、Figure 2(c)のように選択した範囲にA1:D4の配列が参照されてくる。A6からD9を移動させながら数式バーを見ると、どのセルにも同じ数式が入力されていることが分かる。これらの数式は{=A1:D4}と"{}"で囲まれており、配列入力されていることを示している。すなわち、配列入力をすることによって、1つの配列数式から返される複数の結果を配列として、ワークシート上に表示している。配列入力をした範囲は、1つの配列数式からの出力のセットであるため、その範囲の一部だけを修正・削除することはできない。数式が間違えていた場合、配列入力されている範囲のセル(どこでもよい)を選択し、数式バーをクリックするか[F2]キーを押して編集できる状態にする。編集した後、再度配列入力をしなくてはならない。


Figure 2. An example of the array formula for refering the cell range A1:D4 to A6:D9.
(a)Select the cell range A6:D9. (b)Type "=A1:D4". (c)Press [Ctrl] + [Shift] + [Enter] key.

次に配列を使った計算の例を示す。前の例で用いたA6からD9の内容を削除し、Figure 3(a)に示したようにA6セルに数値"2"を入力し、A8からD11までの4行4列の範囲を選択する。A8セルに"=A1:D4*A6"と記述し(Figure 3(a))、配列入力する。出力はFigure 3(b)に示したように、配列A1:D4の全ての要素を2倍(A6セルに記入された数値倍)した数値からなる配列になっている。これは、A8セルに"=A1*$A$6"と入力し、D11までの範囲にコピーした場合と同じ結果となる。このように、絶対・相対参照を活用した場合、配列数式を用いた場合では異なる数式でありながら、同じ計算を行うことができる。


Figure 3. Example of a calculation using the array formula. The output array (A8:D11) has elements whose values are twice those of the original array (A1:D4).
(a)Select the cell range A8:D11. Type "=A1:D4*A6" . (b)Press [Ctrl] + [Shift] + [Enter] key.

Figure 4は1行のみの配列と複数行の配列の掛け算の例である。前の例に加えて、B6からD6にも数値を入力し、A6からD6の範囲を1行4列の配列とする(Figure 4(a))。A8からD11を選択し、A8セルに"=A1:D4*A6:D6"として配列入力をする。この配列数式の出力がFigure 4(b)であり、配列A8:D11の各列に配列A6:D6の各列を掛けた数値となっている。要するに、A8セルに"=A1*A$6"と入力し、コピーしたときと同じ出力となる。この例では、2つの配列の1列目がA列となっており共通していたが、2つの配列のワークシート上での位置関係は任意であり、まったく別の場所にあっても1列目同士が計算される。また、"= A6:D6*A1:D4"のように数式内の配列の順序を逆にして計算しても、同様の結果が得られる。また、1行だけの配列との計算をしたときも同様に、行同士の計算をすることができる。


Figure 4. Product of 1 × 4 array by 4 × 4 array.
(a)Select the cell range A8:D11. Type "=A1:D4*A6:D6". (b)Press [Ctrl] + [Shift] + [Enter] key.

配列の計算の最後の例として、4行4列の配列同士の掛け算の様子をFigure 5に示す。A6からD9までの領域にFigure 5(a)のように4行4列の配列を作成し、A11からD14の範囲を選択する。A11セルに"=A1:D4*A6:D9"として配列入力する。これによって、2つの配列の行と列が同じ要素同士が掛け算されて出力される(Figure 5(b))。これを一般的な数式で表すと、A11セルに"=A1*A6"と入力してコピーした場合に相当する。


Figure 5. Calculation of product between two arrays of 4 × 4 size.
(a)Select the cell range A8:D11. Type "=A1:D4*A6:D9". (b)Press [Ctrl] + [Shift] + [Enter] key.

ここまでの計算で明らかなように、配列の計算は、行列の計算とは異なっている。また、積だけでなく、四則演算や累乗などの計算も可能である。1つ1つの例を挙げれば際限がないので、個々でいろいろと試してほしい。
配列数式の表現方法を用いると、通常の数式では多くの作業用のセルを必要とするような計算であっても、1つのセルで計算が済んでしまうことが多い。逆に、出力が非常にシンプルである一方で、数式が複雑になるため、計算の流れがつかみづらくなることもある。まずは、作業用のセルを用いて計算手順を構築してゆき、それをセルの節約などのためにまとめる段階で効果的に配列を使うと良い。前述したように、配列数式の結果は配列となるため、配列のある範囲に行や列およびセルの挿入はできない。そのような挿入を効果的に用いたい場合は、絶対・相対参照を用いた通常の数式を用いたほうが良いこともある。
最後に配列を扱いやすくするため、特定のセル範囲に名前を定義する方法を説明する。Figure 6はこれまで用いてきたA1からD4で示される配列であるが、この配列にAという名前を定義する。この範囲を選択した後、左上にある[名前ボックス]にマウスポインタを移動し、クリックしたのち"A"と入力する。こうすることで、"A"と"A1:D4"は同等となり、数式中で"A"と入力するだけで"A1:D4"と言う範囲を示すことになる。使えない文字列(C, R, A1など)があるので注意が必要である。範囲と名前がどのように定義されているかは、[Ctrl]+[F3]キーで確認できる。


Figure 6. Naming the cell range (A1:D4) "A" using name box.
(1)Select the cell range A1:D4. (2)Type "A" in [name box] and press [Enter] key.

2. 3 行列の計算

Excelには、配列を使って行列の計算を行える関数が4つ用意されている。転置行列を求める"TRANSPOSE"関数、行列の積"MMULT"関数、行列式を求める"MDETERM"関数、逆行列"MINVERSE"関数である。Figure 7にこれらの関数を用いた計算の結果を示す。これまでと同様A1からD4の配列から出発する(個々の要素の数値が異なっていることに注意。またこの範囲は、前節で"A"と定義されている)。


Figure 7. Examples of matrix calculations on an Excel worksheet.

2番目の配列(A6からD9)は、配列数式"=TRANSPOSEA(A)"を配列入力することで求めた、Aの転置行列である。この配列数式は、"=TRANSPOSE(A1:D4)"と入力しても構わないが、名前を使うことで簡素になるだけでなく、行列をより意識することができる。こうして求めた転置行列である、A6からD9までを、ATと定義した。
A11からD14は、"=MMULT(A,AT)"(行列の積なので、MMULT(AT,A)は別計算になることに注意)と配列入力して求めた2つの行列の積であり、"AAT"と定義した。注意しなくてはならない点として、MMULT関数には生成される行列の大きさに制限があるようで、行列の要素の数が5461を超える行列を生成するような計算はできない(ただしExcel 2007ならばこれ以上でも計算可能である。また、Excel 2003以前であっても、行列の積をいくつかに分けて行えばこの要素数以上でも計算可能である)。
この計算で得られた行列の行列式は、"=MDETERM(AAT)"と通常入力して求めることができる(A16セル)。ここでは行列式が1であった。行列が特異行列もしくは特異行列とみなせる場合は、行列式が0もしくは0とみなせるほど非常に小さな値になる。
今回の例では逆行列が存在することが分かったので、A18からD21(INVと定義)に"=MINVERSE(AAT)"として配列入力することで逆行列を求めた。最後に"=MMULT(AAT,INV)"にて求めた行列の積を、A23からD26に示した。対角要素以外の要素の一部に、0以外の数値が見られるが、これは逆行列を計算する際に生じたわずかな誤差によるものである。しかし、10-12程度の非常に小さな値であることを許容すれば、この行列は単位行列とみなすことができる。
ここまでの計算は、1つ1つの行列計算を順を追って行ってきた(行列式の計算を除く)。しかし、関数の出力をそのまま別の関数の入力にしてもかまわないことから、配列数式でも同様に記述の工夫を行うことによって、途中の計算結果を表示させるためのセルを節約することができる。上の例で、行列とその転置行列の積を計算する場合、一般的な式で表せば、AATとなる。これを1つの配列数式で記述すると、"=MMULT(A,TRANSPOSE(A))"となる。これによって、転置行列を展開するスペースを使うことなく、計算を行うことができる。同様に、最後の単位行列が出力されるところまで計算を行うことも可能で、式はAAT(AAT)-1であるから配列数式では
"=MMULT(MMULT(A,TRANSPOSE(A)),MINVERSE(MMULT(A,TRANSPOSE(A))))"
とすればよい。しかし、この式では"MMULT(A,TRANSPOSE(A))"という同じ計算が2回行われるため効率が悪い。"MMULT(A,TRANSPOSE(A))"を計算し、その出力した範囲(配列)をAATと名前をつけて、その後"=MMULT(AAT,MINVERSE(AAT))"とするのが効率的である。
行列の積の特殊な場合として、ベクトルの内積がある。Figure 8に示したように、2つの列ベクトルを表す範囲の名前をそれぞれVa, Vbと定義すると、Excelでこれらの積を計算するためには"=MMULT(TRANSPOSE(Va),Vb)"と配列入力すればよい(A7セル)。ベクトルの内積は2つのベクトルの対応する要素同士を掛けあわせ、さらに合計するだけであることを利用すれば、"=SUMPODUCT(Va,Vb)と入力(A8セル)するだけで、前述の式と同じ数値を得ることができる。SUMPRODUCT関数を使うことで、TRANSPOSE関数や配列入力の手間を省けるといったことだけでなく、より複雑な計算過程を含んでいるMMULT関数を用いないため、計算の効率が良いものと思われる。


Figure 8. Example of product between two vectors. "SUMPRODUCT" function is simpler expression than "MMULT" function for this calculation.

3 重回帰による曲線フィッティング

3. 1 理論

重回帰は、目的変数fを説明変数gjの線形結合で表現する方法である。ここで、j (= 1, 2,...)は説明変数の種類を表す。説明変数とその係数ajを用いて、目的変数を記述すると、

となる。これを、有限個の説明変数による線形結合の項と残差項rを用いれば、

となる。実験などによって得られた複数の目的変数を、i (= 1, 2,..., m)で区別すると、式(2)は、

と一般化できる。式(2)’を行列で表すと、

であり、fi, gij, aj, riを要素とする行列もしくはベクトルをそれぞれf, G, a. rとすれば、

と表すことができる。目的変数を説明変数のみで表現できることが理想であるから、すべてのriが0になるようなaを求めることが望まれる。n = mの場合、式(2)’で表される式の数と係数の数が等しくなるので、特定の場合を除き全てのriが0となるようなaを一意的に求めることができる。しかし、一般的な問題としてはn < mであるからすべてのriを同時に0にすることはできず、riの二乗和が最小になるようなaを最小二乗法で求めることになる。
最小二乗法では、残差項riの二乗和をEとし、このEの各係数ajによる偏微分が0になるようにする。このようにして得られる連立方程式を整理し行列で表現すると、

となる。aLSは最小二乗解として得られる係数ベクトルを明示するために、LSを下付添え字にした。この連立方程式は、式と説明変数の数が共にnとなるため、一意的に解を求めることができる。そこで、式(5)をaLSに関して解くと

を得る。
さらに、得られたaLSを用いて、

のように、計算によるfの予測値fpreを見積もることができる。

3. 2 実例

3. 2. 1 データ生成

重回帰を用いた曲線フィッティングを、Excelにて実現する方法を簡単な例を使って解説する。まず、fに当たるデータを、

を用い、x = 2980から3020まで2間隔で21点生成した(Figure 9)。B列からE列の2行目および4行目は、それぞれ次数とその係数である。この次数と係数を用いて各xにおけるfの値を、配列数式を利用して求めた。式(8)の各項の値は、[係数].x[次数]で求められる。一般的には、"=B$4*$A5^B$2"をB5セルに入力し、E25までの範囲にコピーする方法を採る。しかしここでは、慣れるという意味からも配列数式を用いることにする。配列数式で表すと、この数式は"=B4:E4*A5:A25^B2:E2"と書くことができる。B5からE25の範囲を選択し、この式を配列入力すれば、Figure 9のように各xにおける各次数の値が求まる。各xの各次数の値の合計を求めれば(たとえば、"=SUM(B5:E5)")、F列に示したfの値を得ることができる。
この計算ではB5からE25のセル領域を各次数の項を計算するために作業セルとして消費している。このような作業セルを使うことなく、配列数式を駆使して計算することも可能である。たとえば、F5セルに、"=SUM($B$4:$E$4*A5^$B$2:$E$2)"を配列入力する(出力は1つしかないが、配列入力しないとエラーになる)。これをコピーして、F6からF25までに貼り付ければよい。もう1つは、"SUMPRODUCT"関数を用いた式で、"=SUMPRODUCT($B$4:$E$4,A5^$B$2:$E$2)"をF5に通常入力する。これを、F25までコピーする。今回の例では、計算の流れをよく見えるようにするためにも、各次数の値を示す方法を採った。もし、作業用のセルを節約したいときは、配列数式を使うとよい。ただし、計算の流れがつかみづらくなってしまい、コメントなどを入れておかないと、他人が見て計算内容のわからないシートとなるばかりか、時間がたつと自分自身でもわからなくなる可能性があるので注意しなくてはならない。
Figure 9のG列のjは、平均が0で標準偏差が1000である誤差をF列の値に加えた、仮想的な偶然誤差を含む3次曲線データである。誤差の発生には、前報[7]にて紹介した"NORMINV"関数を用い、G5セルには"=F5+$G$4*NORMINV(RAND(),0,1)"という数式が入力されており、F25セルまでコピーする。発生させる誤差の標準偏差はG4セルの値で調整できるようになっている。ただし、標準偏差が0(誤差なし)のときにも対応できるよう、G4セルの値を関数に与えず(0のときはエラーとなるため)、関数の出力との積を用いた。誤差の発生に乱数を用いているので、Figure 9のG列は誤差を含んだデータの一例に過ぎない。再計算である[F9]キーを押すごとに値が変化するのを確かめてほしい。


Figure 9. Generaion of cubic curve data without and with error. *SD:standard deviation.

3. 2. 2 フィッティング

生成したデータを高次式による重回帰によってカーブフィッティングするために用いたワークシートがFigure 10である(Figure 9とは別のワークシートを用いた)。 xであるFigure 9のA5からA25の範囲をコピーし、Figure 10のB5からB25の範囲に貼り付けてある。また、fであるFigure 9のA5からA25の範囲には数式が入力されているため、この範囲をコピーした後、貼り付け先の先頭行であるFigure 10のC5セルを選択し、[Alt]→[E]→[S]→[V]→[Enter]の順にキー入力することでC5からC25の範囲に値貼り付けを行った(この手順は、Excel 97から2007に共通して用いることのできる方法であるが、ほかの方法を知っている場合は、その方法でかまわない)。Figure 10のA列のiは式(2)' のiに対応する。前節で用いたそれぞれのxの値をiで区別しxiと表す。さらにxiをまとめてベクトルxとおく(Figure 10 B列)。C列にコピーしてきたfも同様にまとめてベクトルfとし、"f"と名前定義した。
xの0次から3次まで計算して式(2)の各gijを求める必要がある。しかし、この例ではxの範囲が0を含まない上、3000付近と大きな値であるため、エクセルでの計算精度では、どの次数もこの範囲においてほとんど直線性を持つようになり、各次数間の関係が独立とみなせなくなってしまう。そのため、多重共線性の問題が現れ、フィッティングがうまく行われない。多重共線性の問題を回避するためには、各次数の項が互いに最も独立な関係を持つようにする必要がある。つまりxが0を中心とする範囲のデータとなるように変換すると良い。そこで、yi = xi - 3000と変換し、-20から20までのデータにする。この変換データがFigure 10 E列の値であり、E5からE25の範囲には、"=B5:B25-B15"もしくは"=B5:B25-3000"を配列入力した。名前ボックスを使って、E5からE25の範囲を"y"と定義した。次に、G4からJ4の範囲に次数である0から3を入力し、これを"deg"と定義した。"y"と"deg"を用いて、式(4)のGを求めた(Figure 10 G5からJ25の範囲)。この範囲に配列入力されている配列数式は"=IF(deg=0,1,y^deg)"である。IF関数の条件は、"degが0である場合、1を返す"というものである。これは、yの範囲に0が含まれるため、エラーになってしまう00の計算を回避するための処理である。条件が成り立たない場合、"y^deg"を計算する。こうして計算されたG5からJ25の範囲を"G"と定義した。
0から3次までの4項を用いているので、その係数ベクトルであるaLSは4行1列となる。これをN5からN8までの範囲に計算し、"aLS"と名前定義する。この範囲には式(6)に従った配列数式が入力されており、

"=MMULT(MINVERSE(MMULT(TRANSPOSE(G),G)),MMULT(TRANSPOSE(G),f))"      数式(a).

である。この配列数式は、式(6)をそのままExcelの数式表現に翻訳したものであり、式(6)の計算を1つの配列数式で行うことができる。しかし、前述したように計算の流れの見通しが悪くなるので、注意が必要である。
この配列数式によって、N5からN8にはそれぞれ0, -100, 0, 1が計算された。行列の積は結合法則を満たしているので、式(6)の計算のどの積から計算するかは任意である。すなわち数式(a)は式(6)を{(GTG)-1GT}fのように計算しているが、(GTG)-1(GTf)としても計算可能である。後者の式をExcelの数式で記述したのが、以下の数式(b)である。

"=MMULT(MMULT(MINVERSE(MMULT(TRANSPOSE(G),G)),TRANSPOSE(G)),f)"          数式(b).

数式(b)で計算して得られたaLSFigure 10のN11からN14)のうち、0次項の係数(N11セル)が1.14E-13となったが、これは計算の過程で発生した誤差によるものと思われ、非常に小さな値であるから0とみなすことができる。そうすれば、数式(a)と数式(b)の計算結果は同じであるといえる。


Figure 10. The least squares fitting calculation using matrix on an Excel worksheet.
Coefficient of each degree, aLS, was calculated by *1: formula (a) and *2: formula (b).

式(7)にしたがって、行列Gと最小二乗法によって求められたaLSの積を計算することで、fの予測値fpreを求めることができる。これは"=MMULT(G,aLS)"を配列入力することで計算でき、結果をFigure 10 P列に示した。
ここまでのような重回帰計算は、Excelが基本的に持っている"LINEST"関数や、"分析ツール"に含まれる"回帰分析"を用いることでも実現できる。さらに、これらの機能では、予測値fpreが実測値fにどれだけ適合しているかの目安となる重相関係数Rなどの統計データも得ることができる。重相関係数は、0から1の値をとり、1に近いほどより適合していると考えることができる。重相関係数は、式(9)に示すように、予測値の変動を実測値の変動で除し、その平方根を求めることで計算でき、以下の式で表される。

ここでffpreはそれぞれfifpreiの平均である。この式に従い、重相関係数を計算したのが、Figure 10のN17セルである。入力されている数式は、"=SQRT(DEVSQ(P5:P25)/DEVSQ(f))"である。
ffpreが同じであること、Rが1であることなどから、ほぼ完全にフィッティングされていることがわかる。このことはFigure 11ffpreのプロットがほぼ完全に重なっていることからも明らかである。これらは、3次式で生成した点を、3次式でフィッティングしていることに加え、元のデータにノイズがまったく含まれていないことによる。


Figure 11. Comparison of fittings for data with and without error. The original data without error, f (), and fitting data, fpre (+). The original data with error, j(), and fitting data, jpre (×).

しかし、実際にはノイズが含まれないデータを得ることはまずあり得ないうえ、想定した高次式で必ずフィッティングができるとも限らない。そこで、Figure 9のワークシートで誤差を含ませたデータj(G5からG25)をコピーし、Figure 10のワークシートのf(C5からC25)に値貼り付けすると、誤差を含むデータに対する3次曲線フィッティングをすることができる。この誤差込みデータとそれに対するフィッティング結果をプロットしたものが、Figure 11j)とjpre×)である。jfからずれた点ではあるものの、fに誤差を含ませたデータなので、fに沿って分散していることが分かる。そのため、フィッティング結果jprefpreと一致はしなかったが、よく似た傾向でを示している。また、Figure 10で計算される重相関係数Rは0.933...と1を下回る値となったが、比較的1に近くかなりの適合性があると思われる。
今回の例では、重回帰の手法を用いた3次曲線によるフィッティングを行った。もちろん式(1)に従い、式(6)を計算すれば、任意の次数の曲線だけでなく、任意のgjを用いて重回帰が可能である。実際に得られる実験値やスペクトルに対して、重回帰による解析を行う際には、適当なモデル・仮定に基づいて適切なgjを決めなくてはならない。その上で、実験データには誤差が含まれていること、フィッティングの式が本当に適当であるかと言ったことを考察するためにも、重相関係数などを用いて適合度合いについて熟慮する必要がある。

3. 2. 3 係数の変換

誤差の無いデータfに議論を戻し、最小二乗計算によって得られたaLSについて考察する。このaLSyの多項式の係数である。また、y = x - 3000という関係があるので、fpre

となる。これは、yで表した式であり式(8)とは異なっているが、各項を展開すると各係数を計算することができる。任意のxについてfpreを求めることが目的であるなら、式(10)を使うだけでかまわない。しかし、もしxの各次数の項の係数を求めたい場合は、この式を展開しなくてはならない。
式(10)ように簡単な式ならば展開は容易であるが、より複雑で高次の式の場合、展開するのは非常に手間がかかる。そこで、y = (x + x0)の関係があるとき、yの高次式を展開して、xの各次数の係数をExcelで求める方法を解説する。yの高次式とxの高次式の間には次のような関係が成り立たなくてはならない。

二項定理を考慮してこの関係式から、albkを用いて表すと、

となる。Figure 12は、フィッティングの例で得られたyの係数からxの係数を計算するために用いたExcelワークシートである。例では3次式であったので、次数は0から3までを用意しておく。B3セルからB6セルまでは、フィッティングによって得られた結果である。y = x - 3000であるからx0 = -3000(B9セル)である。これらの数値と式(12)を用いて、配列数式でxの各次数の係数を求める。組み合わせkClを求めるには、"COMBIN"関数を用いて、"=COMBIN(k,l)"と入力する。E3セルには以下の式が配列入力されている。


Figure 12. Coefficients conversion from y polynomial to x polynomial.

"=SUM(IF($A$3:$A$6<D3,0,$B$3:$B$6*COMBIN($A$3:$A$6,D3)*$B$9^($A$3:$A$6-D3)))"    数式(c)
ここで注意すべきことは、この数式(c)をE3セルに単独で配列入力しなくてはならないことである。E4からE6セルには、E3セルをコピーして貼り付ければよい。また、式(12)ではSの条件によりklよりも小さくなることを回避しているが、Excelの数式では"IF"関数を用いて、klよりも小さくなる場合は"0"を返すことで対応した。このようにして得られたxの各次数の係数は、-26999700000, 26999900, -9000, 1であり、式(8)と完全に一致していた。
もちろん、誤差を含んだ場合のaLSをB3からB6セルに貼り付ければ、それに応じた係数が計算される。

4 終わりに

Excelによるケモメトリクス計算の準備として、行列の扱いと重回帰計算について解説した。今回紹介した内容は、配列数式を多用している上、計算手順は計算で用いるセルを節約するために1つの式に数式を集約した。しかし、Excelは広い2次元のセル空間を持っており、その使い方は自由である。普通の数式を用い、途中計算用の作業セルを利用したワークシートの設計であってかまわない。ただし、行列の計算を関数で行うならば、どうしても配列数式を理解しなくてはならないので、この機会にマスターしてほしい。
多変量を扱う最も基本的な方法として、今回は重回帰を取り上りあげ、行列計算用の関数を用いて解いた。一般には、単回帰であれば"SLOPE"関数や"INTERCEPT"関数、重回帰であれば"LINEST"関数を用いることができる。さらに、"分析ツール"の"回帰分析"を使っても得られる。これらの方法を用いた場合と、今回紹介した方法の相違点の1つが、重回帰計算の内容をブラックボックスにするか否かである。ケモメトリクス計算の教育という観点からは、計算過程を理解すること、その過程を自身で実行してみることが重要であると考えている。そのため本稿では、計算の理論の解説とその理論に沿った実際の計算を、読者自身がExcelを用いて行えるようにした。その結果としてケモメトリクス計算の理解が深まり、ケモメトリクス計算の普及・発展につながると期待している。
続報では、重回帰の定量への応用、主成分分析(PCA)、主成分回帰(PCR)、partial least squares (PLS)回帰法、前処理で用いる平滑化、微分についてExcelの基本機能だけを用いて実行する方法を解説してゆく予定である。

参考文献

[ 1] 長谷川健, スペクトル定量分析, 講談社サイエンティフィク (2005).
[ 2] 三井利幸, ケモメトリックスの基礎と応用−分析化学と多変量解析法−, アイピーシー (2003).
[ 3] 尾崎幸洋, 宇田明史, 赤井俊雄, 化学者のための多変量解析−ケモメトリクス入門−, 講談社サイエンティフィク (2002).
[ 4] 宮下芳勝, 佐々木慎一, ケモメトリックス 化学パターン認識と多変量解析, 共立出版 (1995).
[ 5] 吉村忠与志, 厳選例題Excelで解く問題解決のための化学計算入門, 技術評論社 (2005).
[ 6] 吉村忠与志, 青山義弘, トランジスタ技術Special No.78 技術者のためのExcel活用研究, CQ出版社 (2002).
[ 7] 吉村季織, J. Comput. Chem. Jpn., 5, 29 (2006).


Return