Microsoft Word - 1

J. Comput. Chem. Jpn., Vol. 7, No. 2, pp. 71–82 (2008)

技術論文

Microsoft Excel

を用いたケモメト リクス計算

(1)

−行列の扱いと重回帰−

吉村 季織*, 高柳 正夫

国立大学法人東京農工大学農学府, 〒 183-8509 東京都府中市幸町 3-5-8 *e-mail: [email protected]

(Received: January 28, 2008; Accepted for publication: May 1, 2008; Advance publication: May 30, 2008)

近年、化学データを数学的・統計的手法により解析する「ケモメトリクス」が盛んに用いられるよ うになってきた。しかし 、日本の大学の化学教育の場ではほとんど 取り上げられていない。ケモメト リクスや数値計算の専用ソフトウェアを使うことなく、現在最も普及しているソフトウェアのひとつ である Microsoft Excel（ Excel）の基本機能を用いてケモメトリクス計算を行うことができれば 、多く の教育・研究機関で役立つものと思われる。そこで本稿では 、Excel を用いたケモ メト リクス計算シ リーズの第１弾として、Excel での行列の扱い方とその応用として多変量解析の基本である重回帰によ る曲線フィッティングを解説する。Excel で行列を扱うことは、配列を扱うことであり、配列数式に慣 れる必要がある。そこで、最も簡単な配列数式から、行列の演算、行列式、逆行列の計算を行い、基 礎的な配列の取り扱いを解説した。また、重回帰では、曲線データの生成と、生成したデータへの高 次曲線の最小二乗フィッティングを配列数式により実行した。 キーワード : Microsoft Excel, ケモメトリクス, 配列数式, 行列, 重回帰

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

行列と配列

ケモメトリクスでは行列が用いられる。そこで、Ex-celで行列を扱う方法をまず紹介する。Excel による行 列計算法に関しては、吉村忠与志らの著書 [5, 6] に他 の科学技術計算とあわせて紹介されているので、参考 にされたい。本稿だけでも Excel による行列計算を習 得できるようにするため、および今後のシリーズにお いて基礎となるため、以下より Excel による行列計算 について解説する。 Excelをはじめとする表計算ソフトウェアは、２次 元的に敷き詰められたセルより構成されている。セル の横方向の並びを “行”、縦方向の並びを “列”と称す るように、このセルの並びそのものが行列であるとい える。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] キー

(a)

(b)

(c)

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.

を 押すことで 入力を 確定する（ この「[Ctrl] キーと [Shift]キーを押した状態で [Enter] キーを押すという入 力確定法」を本稿では以降「配列入力」と呼び、「[Enter] キーのみを押す」通常の入力確定法と区別する）。そ うすると、Figure 2(c) のように選択した範囲に A1:D4 の配列が参照されてくる。A6 から D9 を移動させなが ら数式バーを見ると、どのセルにも同じ数式が入力さ れていることが分かる。これらの数式はf=A1:D4gと “fg”で囲まれており、配列入力されていることを示し ている。すなわち、配列入力をすることによって、１ つの配列数式から返される複数の結果を配列として、 ワークシート上に表示している。配列入力をした範囲 は 、1 つの配列数式からの出力のセットであるため、 その範囲の一部だけを修正・削除することはできない。 数式が間違えていた場合、配列入力されている範囲の セル（どこでもよい）を選択し 、数式バーをクリック するか [F2] キーを押して編集できる状態にする。編 集した後、再度配列入力をしなくてはならない。 次に配列を使った計算の例を示す。前の例で用いた 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 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 行 だけの配列との計算をしたときも同様に、行同士の計 算をすることができる。

(a)

(b)

Figure 3. Example of a calculation using the array for-mula. 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.

(a)

(b)

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.

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.

配列の計算の最後の例として 、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”と入力してコピーした場合 に相当する。 ここまでの計算で明らかなように、配列の計算は、 行列の計算とは異なっている。また、積だけでなく、 四則演算や累乗などの計算も可能である。１つ１つの 例を挙げれば際限がないので、個々でいろいろと試し てほしい。 配列数式の表現方法を用いると、通常の数式では多 くの作業用のセルを必要とするような計算であっても、 1つのセルで計算が済んでし まうことが多い。逆に 、 出力が非常にシンプルである一方で、数式が複雑にな るため、計算の流れがつかみづらくなることもある。 まずは、作業用のセルを用いて計算手順を構築してゆ き、それをセルの節約などのためにまとめる段階で効 果的に配列を使うと良い。前述したように、配列数式 Name box

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.

の結果は配列となるため、配列のある範囲に行や列お よびセルの挿入はできない。そのような挿入を効果的 に用いたい場合は、絶対・相対参照を用いた通常の数 式を用いたほうが良いこともある。 最後に配列を扱いやすくするため、特定のセル範囲 に名前を定義する方法を説明する。Figure 6 はこれま で用いてきた A1 から D4 で示され る配列であるが 、 この配列に A という名前を定義する。この範囲を選 択した後、左上にある [名前ボックス] にマウスポイン タを移動し 、クリックしたのち “A”と入力する。こう することで 、“A”と “A1:D4”は同等となり、数式中で “A”と入力するだけで “A1:D4”と言う範囲を示すこと になる。使えない文字列（ C, R, A1 など ）があるので 注意が必要である。範囲と名前がどのように定義され ているかは、[Ctrl]+[F3] キーで確認できる。

2.3

行列の計算

Excelには、配列を使って行列の計算を行える関数 が４つ用意されている。転置行列を求める “TRANS-POSE”関数、行列の積 “MMULT”関数、行列式を求め る “MDETERM”関数、逆行列 “MINVERSE”関数であ る。Figure 7 にこれらの関数を用いた計算の結果を示 す。これまでと同様 A1 から D4 の配列から出発する （ 個々の要素の数値が異なっていることに注意。また この範囲は、前節で “A”と定義されている）。 2番目の配列（ A6 から D9）は、配列数式 “=TRANS-POSEA(A)”を配列入力することで求めた、A の転置行 列である。この配列数式は、“=TRANSPOSE(A1:D4)” と入力しても構わないが 、名前を使うことで簡素にな るだけでなく、行列をより意識することができる。こ うして求めた転置行列である、A6 から D9 までを、AT と定義した。

Original matrix, “A”.

Transposed matrix of A, “AT”. Formula:”=TRANSPOSE(A)”

Product between A and AT, “AAT”. Formula:”=MMULT(A,AT)”

ńDetermination of AAT. Formula:”=MDETERM(AAT)”

Inverse matrix of AAT, “INV”. Formula:”=MINVERSE(AAT)”

Product between AAT and INV. Formula:”=MMULT(AAT, INV)”

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

A11から D14 は、“=MMULT(A,AT)”（ 行列の積な ので、MMULT(AT,A) は別計算になることに注意）と 配列入力して求めた 2 つの行列の積であり、“AAT”と 定義した。注意しなくてはならない点として、MMULT 関数には生成される行列の大きさに制限があるようで、 行列の要素の数が 5461 を超える行列を生成するよう な計算はできない（ただし Excel 2007 ならばこれ以上 でも計算可能である。また、Excel 2003 以前であって も、行列の積をいくつかに分けて行えばこの要素数以 上でも計算可能である）。 こ の 計 算で 得られ た 行 列 の 行 列 式は 、“=MDE-TERM(AAT)”と 通 常 入 力し て 求め るこ とが で き る （ A16 セル ）。ここでは行列式が１であった。行列が特 異行列もしくは特異行列とみなせる場合は、行列式が 0もし くは 0 とみなせるほど 非常に小さな値になる。 今回の例では逆行列が存在することが分かったので、

A18から D21（ INV と定義）に “=MINVERSE(AAT)” として配列入力することで逆行列を求めた。最後に “=MMULT(AAT,INV)”にて求めた行列の積を、A23 か ら D26 に示した。対角要素以外の要素の一部に、0 以 外の数値が見られるが 、これは逆行列を計算する際に 生じたわずかな誤差によるものである。しかし 、10-12 程度の非常に小さな値であることを許容すれば 、この 行列は単位行列とみなすことができる。 ここまでの計算は、1 つ 1 つの行列計算を順を追っ て行ってきた (行列式の計算を除く)。しかし 、関数 の出力をそのまま別の関数の入力にしてもかまわな いことから 、配列数式でも同様に記述の工夫を行う ことによって 、途中の計算結果を表示させるための セルを節約することができる。上の例で 、行列とそ の転置行列の積を計算する場合、一般的な式で表せ ば 、AA T となる。これを１つの配列数式で記述する

と、“=MMULT(A,TRANSPOSE(A))”となる。これによ って、転置行列を展開するスペースを使うことなく、 計算を行うことができる。同様に、最後の単位行列が 出力されるところまで計算を行うことも可能で、式は AA T AA T  ;1 であるから配列数式では “=MMULT(MMULT(A,TRANSPOSE(A)),MINVER SE(MMULT(A,TRANSPOSE(A))))” とすればよい。しかし 、この式では “MMULT(A,TRAN SPOSE(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 関数を用いないため、計算の効率が良 いものと思われる。

3

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

3.1

理論

重回帰は 、目的変数 f を説明変数 gjの線形結合で 表現する方法である。ここで、j (= 1, 2, )は説明変 数の種類を表す。説明変数とその係数 ajを用いて、目 的変数を記述すると、 f =a 1 g 1 +a 2 g 2 + = 1 X j=1 a j g j (1) となる。これを、有限個の説明変数による線形結合の 項と残差項 r を用いれば 、

Va

Vb

formula:”=MMULT(Va,TRANSPOSE(Vb))” and press [Ctrl] + [Shift] + [Enter] key. formula:”=SUMPRODUCT(Va,Vb)”

and press [Enter] key.

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

f = n X j=1 a j g j +r (2) となる。実験などによって得られた複数の目的変数を、 i (= 1, 2, , m)で区別すると、式 (2) は、 f i = n X j=1 a j g ij +r i (2)’ と一般化できる。式 (2) ’を行列で表すと 、 0 B B B B B B B B @ f1 f 2 . .. f

i

.. . f

m

1 C C C C C C C C A = 0 B B B B B B B B @ g11g12g1

j

g1

n

g 21 g 22 g 2

j

g 2

n

. .. ... . .. ... ... g

i

1 g

i

2 g

ij

 g

in

.. . ... ... . .. ... g

m

1g

m

2 g

mj

g

mn

1 C C C C C C C C A 0 B B B B B B B B @ a1 a 2 . .. a

j

.. . a

n

1 C C C C C C C C A + 0 B B B B B B B B @ r1 r 2 . .. r

i

.. . r

m

1 C C C C C C C C A (3) であり、fi, gij, aj, riを要素とする行列もし くはベク トルをそれぞれfGarとすれば 、 f =Ga+r (4) と表すことができる。目的変数を説明変数のみで表現 できることが理想であるから、すべての riが 0 になる ようなaを求めることが望まれる。n = m の場合、式 (2)’で表される式の数と係数の数が等しくなるので 、 特定の場合を除き全ての riが 0 となるような aを一 意的に求めることができる。しかし 、一般的な問題と しては n<mであるからすべての r iを同時に 0 にす ることはできず、riの二乗和が最小になるような aを 最小二乗法で求めることになる。 最小二乗法では 、残差項 riの二乗和を E とし 、こ の E の各係数 ajによる偏微分が 0 になるようにする。

このようにして得られる連立方程式を整理し行列で表 現すると、 G T f =G T Ga LS (5) となる。a LSは最小二乗解として得られる係数ベクト ルを明示するために、LS を下付添え字にした。この連 立方程式は、式と説明変数の数が共に n となるため、 一意的に解を求めることができる。そこで、式 (5) を a LSに関して解くと a LS =(G T G) ;1 G T f (6) を得る。 さらに、得られたa LSを用いて、 f pre =Ga LS (7) のように、計算によるf の予測値f preを見積もるこ とができる。

3.2

実例

3.2.1 データ生成 重回帰を用いた曲線フィッティングを 、Excel にて 実現する方法を簡単な例を使って解説する。まず、f に当たるデータを 、 f =;26999700000+26999900x;9000x 2 +x 3 (8) を用い、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 までに貼り付けれ ばよい。もう１つは、“SUMPRODUCT”関数を用いた 式で 、“=SUMPRODUCT($B$4:$E$4,A5ˆ$B$2:$E$2)” を F5 に通常入力する。これを、F25 までコピーする。 今回の例では、計算の流れをよく見えるようにするた めにも、各次数の値を示す方法を採った。もし 、作業 用のセルを節約したいときは、配列数式を使うとよい。 ただし 、計算の流れがつかみづらくなってしまい、コ メントなどを入れておかないと、他人が見て計算内容 のわからないシートとなるばかりか、時間がたつと自 分自身でもわからなくなる可能性があるので注意しな くてはならない。 Figure 9の G 列の'は、平均が 0 で標準偏差が 1000 である誤差を F 列の値に加えた、仮想的な偶然誤差を 含む 3 次曲線データである。誤差の発生には、前報 [7] にて紹介した “NORMINV”関数を用い、G5 セルには “=F5+$G$4*NORMINV(RAND(),0,1)”という数式が入 力されており、F25 セルまでコピーする。発生させる 誤差の標準偏差は G4 セルの値で調整できるようになっ ている。ただし 、標準偏差が 0(誤差なし) のときにも 対応できるよう、G4 セルの値を関数に与えず（ 0 のと きはエラーとなるため）、関数の出力との積を用いた。 誤差の発生に乱数を用いているので、Figure 9 の G 列 は誤差を含んだデータの一例に過ぎない。再計算であ る [F9] キーを押すごとに値が変化するのを確かめて ほしい。 3.2.2 フィッティング 生成したデータを高次式による重回帰によってカー ブフィッティングするために用いたワークシートが Fig-ure 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 9. Generaion of cubic curve data without and with error. *SD:standard deviation. 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 項を用いているので、その係数 ベクトルであるa LSは 4 行 1 列となる。これを N5 か ら N8 までの範囲に計算し 、“aLS”と名前定義する。こ の範囲には式 (6) に従った配列数式が入力されており、 “=MMULT(MINVERSE(MMULT(TRANSPOSE(G), G)),MMULT(TRANSPOSE(G),f))” 数式 (a).

Figure 10. The least squares fitting calculation using matrix on an Excel worksheet. Coefficient of each degree,a

LS, was calculated by *1: formula (a) and *2: formula (b).

である。この配列数式は、式 (6) をそのまま Excel の 数式表現に翻訳したものであり、式 (6) の計算を 1 つ の配列数式で行うことができる。しかし 、前述したよ うに計算の流れの見通しが悪くなるので、注意が必要 である。 この配列数式によって 、N5 から N8 にはそれぞ れ 0, -100, 0, 1 が 計算され た 。行列の 積は 結合法 則を 満たし ているので 、式 (6) の 計算のど の 積か ら計算するかは任意である。すなわち数式 (a) は式 (6)を n (G T G) ;1 G T o f のように計算し ているが 、 (G T G) ;1 (G T f)としても計算可能である。後者の式 を Excel の数式で記述したのが、以下の数式 (b) である。 “=MMULT(MMULT(MINVERSE(MMULT(TRANS POSE(G),G)),TRANSPOSE(G)),f)” 数式 (b). 数式 (b) で計算して得られたa LS（ Figure 10 の N11 か ら N14）のうち、0 次項の係数（ N11 セル）が 1.14E-13 となったが 、これは計算の過程で発生した誤差による ものと思われ 、非常に小さな値であるから 0 とみなす ことができる。そうすれば 、数式 (a) と数式 (b) の計 算結果は同じであるといえる。 式 (7) にしたがって、行列Gと最小二乗法によって求 められたa LSの積を計算することで、 fの予測値f pre を求めることができる。これは “=MMULT(G,aLS)”を 配列入力することで計算でき、結果を Figure 10 P 列 に示した。 ここまでのような重回帰計算は、Excel が基本的に 持っている “LINEST”関数や、“分析ツール”に含まれ る “回帰分析”を用いることでも実現できる。さらに、 これらの機能では、予測値 fpreが実測値 f にどれだけ 適合しているかの目安となる重相関係数 R などの統計 データも得ることができる。重相関係数は、0 から 1 の値をとり、1 に近いほどより適合していると考える ことができる。重相関係数は 、式 (9) に示すように 、 予測値の変動を実測値の変動で除し 、その平方根を求 めることで計算でき、以下の式で表される。 R= v u u t n X i=0 ; f prei ;f pre  2 , n X i=0 ; f i ;f  2 : (9)

-6000 -4000 -2000 0 2000 4000 6000 2980 2985 2990 2995 3000 3005 3010 3015 3020 f fpre no error x f with error ij ijpre

Figure 11. Comparison of fittings for data with and with-out error. The original data withwith-out error, f (), and fitting

data, fpre(+). The original data with error,

'( ), and fit-ting data,' pre( ). ここで f、f preはそれぞれ f i、 f prei の平均である。 こ の 式に 従い 、重 相 関 係 数を 計 算し た のが 、Fig-ure 10の N17 セルである。入力されている数式は 、 “=SQRT(DEVSQ(P5:P25)/DEVSQ(f))”である。 f と fpreが同じであること、R が 1 であることなど から、ほぼ完全にフィッティングされていることがわ かる。このことは Figure 11 で f と fpreのプロットが ほぼ完全に重なっていることからも明らかである。こ れらは、3 次式で生成した点を、3 次式でフィッティン グしていることに加え、元のデータにノイズがまった く含まれていないことによる。 しかし 、実際にはノイズが含まれないデータを得る ことはまずあり得ないうえ 、想定した高次式で必ず フィッティングができるとも限らない。そこで、Figure 9のワークシートで誤差を含ませたデータ'（ G5 から G25）をコピーし 、Figure 10 のワークシートの f（ C5 から C25）に値貼り付けすると、誤差を含むデータに 対する 3 次曲線フィッティングをすることができる。 この誤差込みデータとそれに対するフィッティング結 果をプロットしたものが 、Figure 11 の'（ ）と' pre （）である。'は f からずれた点ではあるものの、f に 誤差を含ませたデータなので 、f に沿って分散してい ることが分かる。そのため、フィッティング結果' pre も fpreと一致はしなかったが 、よく似た傾向でを示し ている。また、Figure 10 で計算される重相関係数 R は 0.933 と 1 を下回る値となったが 、比較的 1 に近く かなりの適合性があると思われる。 今回の例では、重回帰の手法を用いた 3 次曲線によ るフィッティングを行った。もちろん式 (1) に従い、式 (6)を計算すれば 、任意の次数の曲線だけでなく、任 意の gjを用いて重回帰が可能である。実際に得られ る実験値やスペクトルに対して、重回帰による解析を 行う際には、適当なモデル・仮定に基づいて適切な gj を決めなくてはならない。その上で、実験データには 誤差が含まれていること、フィッティングの式が本当 に適当であるかと言ったことを考察するためにも、重 相関係数などを用いて適合度合いについて熟慮する必 要がある。 3.2.3 係数の変換 誤差の無いデータ f に議論を戻し 、最小二乗計算に よって得られたa LSについて考察する。この a LSは y の多項式の係数である。また、y = x – 3000 という関 係があるので 、fpreは f pre = ;100y+y 3 = ;100(x;3000)+(x;3000) 3 (10) となる。これは、y で表した式であり式 (8) とは異なっ ているが 、各項を展開すると各係数を計算することが できる。任意の x について fpreを求めることが目的で あるなら、式 (10) を使うだけでかまわない。しかし 、 もし x の各次数の項の係数を求めたい場合は、この式 を展開しなくてはならない。 式 (10) ように簡単な式ならば展開は容易であるが 、 より複雑で高次の式の場合、展開するのは非常に手間 がかかる。そこで 、y = (x + x0)の関係があるとき、y の高次式を展開して、x の各次数の係数を Excel で求 める方法を解説する。y の高次式と x の高次式の間に は次のような関係が成り立たなくてはならない。 M X k =0 b k y k = M X l=0 a l x l ,ここでy=x+x 0 (11) 二項定理を考慮してこの関係式から、alを bkを用い て表すと、 a l = M X k =l b k k C l x k ;l 0 (12)

Figure 12. Coefficients conversion from y polynomial to x polynomial. となる。Figure 12 は、フィッティングの例で得られた yの係数から x の係数を計算するために用いた Excel ワークシートである。例では 3 次式であったので、次 数は 0 から 3 までを用意しておく。B3 セルから B6 セ ルまでは、フィッティングによって得られた結果であ る。y = x – 3000 であるから x0= –3000（ B9 セル ）で ある。これらの数値と式 (12) を用いて、配列数式で x の各次数の係数を求める。組み合わせkClを求めるに は 、“COMBIN”関数を用いて 、“=COMBIN(k,l)”と入 力する。E3 セルには以下の式が配列入力されている。 “=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) ではの条件により k が l よりも 小さくなることを回避しているが 、Excel の数式では “IF”関数を用いて、k が l よりも小さくなる場合は”0” を返すことで対応した。このようにして得られた x の 各次数の係数は、-26999700000, 26999900, -9000, 1 で あり、式 (8) と完全に一致していた。 もちろん 、誤差を含んだ場合のa LSを B3 から B6 セルに貼り付ければ 、それに応じた係数が計算される。

4

終わりに

Excelによるケモメトリクス計算の準備として、行 列の扱いと重回帰計算について解説した。今回紹介し た内容は、配列数式を多用している上、計算手順は計 算で用いるセルを節約するために 1 つの式に数式を集 約した。しかし 、Excel は広い 2 次元のセル空間を持っ ており、その使い方は自由である。普通の数式を用い、 途中計算用の作業セルを利用したワークシートの設計 であってかまわない。ただし 、行列の計算を関数で行 うならば 、ど うしても配列数式を理解しなくてはなら ないので、この機会にマスターしてほしい。 多変量を扱う最も基本的な方法として、今回は重回 帰を取り上りあげ 、行列計算用の関数を用いて解いた。 一般には 、単回帰であれば “SLOPE”関数や “INTER-CEPT”関数、重回帰であれば “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).

Chemometrics Calculations with Microsoft Excel (1)

–Matrix Calculations and Multi Regression–

Norio YOSHIMURA* and Masao TAKAYANAGI

Graduate School of Agriculture, Tokyo University of Agriculture and Technology 3-5-8 Saiwaicho Fuchu, Tokyo 183-8509 Japan

*e-mail: [email protected]

Although chemometrics has become widely used recently for analyzing experimental chemical data, there exist in Japan only a few instructions for the proper usage of chemometrics except for some introductory books. We have found that Microsoft Excel (Excel) is convenient, inexpensive, and popular software for chemometrics calculations especially for educational purposes. In this paper, the first of a series on chemo-metrics calculations with Excel, instructions for matrix calculations and thier applications to multivariate analysis with Excel are described in detail. Array formulas and array calculations are explained in relation to basic handlings of matrix (e.g. calculations of determination and inverse matrices). We also demonstrate applications of array formulas to calculations of polynomial curves, and to curve fittings.