Octave の詳しい説明は省きますが、この数値計算プログラムを行列計算電卓として使う方法を説明します。線型代数についての基本的な知識があるものとして説明しますので説明不足の点があるかも知れません。Octave を使うとややこしい計算をコンピュータに任せることができるので線型代数の定理の意味を直観的に把握することができます。
Octave を起動するにはコンソールから octave と入力します。終了するときは、quit または exit です。
変数に横ベクトルを代入するには次のようにします。
octave:> a = [1, 2, 3, 4] a = 1 2 3 4
a = [1, 2, 3, 4]のように行末に ; をつけないときは代入の結果が表示されますが、次のように行末に ; を付けると結果の表示をさせないようにできます。
octave:> a = [1, 2, 3, 4]; octave:>
縦ベクトルの代入は次のようにします。
octave:> a = [1; 2; 3] a = 1 2 3
変数の後に ' をつけると転置行列(行列の要素が複素数のときは共役転置行列)になります。したがって、a' は横ベクトルになります。
octave:> b = a' b = 1 2 3
ベクトルの加算は + でできます。
octave:> a = [1; 2; 3]; octave:> b = [4; 5; 6]; octave:> c = a + b c = 5 7 9
ベクトルとスカラーの積は * でできます。
octave:> 2 * a ans = 2 4 6
結果を代入する変数を指定していないときは ans という変数に代入されます。変数の内容は変数名をタイプするだけで見ることができます。
octave:> ans ans = 2 4 6
縦ベクトルと横ベクトルの積は順序によって内積になったり、行列になったりします。
octave:> a = [1, 2, 3]; octave:> b = [4; 5; 6]; octave:> a * b ans = 32 octave:> b * a ans = 4 8 12 5 10 15 6 12 18
縦ベクトルどうしの内積は次のようにして計算します。
octave:> a = [1; 2; 3]; octave:> b = [4; 5; 6]; octave:> a' * b ans = 32
横ベクトルどうしの内積は次のようになります。
octave:> a = [1, 2, 3]; octave:> b = [3, 4, 5]; octave:> a * b' ans = 26
行列の代入は次のようになります
octave:> A = [1, 2, 3; 4, 5, 6; 7, 9, 8] A = 1 2 3 4 5 6 7 9 8
行列の要素はインデックスで指定できます。i 行 j 列の要素は A(i, j) です。: 記号は全てのインデックスを表すので、A(i, :) は i 行のベクトル、A(:, j) は j 行のベクトルになります。
行列の転置は ' または .' で行います。実数の行列の場合両方共同じ意味ですが複素数の行列では意味が違いますのでマニュアルで確認してください。
octave:> A' ans = 1 4 7 2 5 9 3 6 8
逆行列は inv() 関数で求めます。
octave:> inv(A) ans = -1.55556 1.22222 -0.33333 1.11111 -1.44444 0.66667 0.11111 0.55556 -0.33333
列ベクトルと行列の積は * で計算できます。
octave:> x = [1; 2; 3]; octave:> y = A * x y = 14 32 49
連立方程式の解は A \ y で計算できます。
octave:> A \ y ans = 1.0000 2.0000 3.0000
行列の演算子にはこの外にも要素毎の掛算をする .* や要素毎の割り算をする ./ など色々ありますが詳しくはマニュアルを見てください。
固有ベクトルと固有値は次のように計算できます。
octave:> [V, D] = eig( A )
V =
0.23489 -0.72518 0.63305
0.53276 0.67762 0.27572
0.81302 -0.12224 -0.72334
D =
15.91994 0.00000 0.00000
0.00000 -0.36313 0.00000
0.00000 0.00000 -1.55681
きちんと計算できているか確認してみましょう
octave:> V * D * inv( V ) ans = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 9.0000 8.0000
行列のランクは rank() 関数で計算できます。
octave:> rank( A ) ans = 3
行列式は det() 関数で計算できます。
octave:> det( A ) ans = 9.0000
行列の対角線成分は diag() で取りだせます。
octave:> diag( A ) ans = 1 5 8
要素が全て 0 の m 行 n 列行列は zeros(m, n) で作れます。
octave:> zeros(2, 3) ans = 0 0 0 0 0 0
要素が全て 1 の場合は次のようにします。
octave:> ones(2, 3) ans = 1 1 1 1 1 1
n行n列の単位行列は eye( n ) で作ります。
octave:> eye( 3 ) ans = 1 0 0 0 1 0 0 0 1
m行n列の乱数行列は rand(m, n) で作れます。
octave:> rand(2, 3) ans = 0.275375 0.820867 0.189354 0.472215 0.920324 0.027263
行列の結合は次のようにできます。
octave:> A = [1, 2; 3, 4]; octave:> B = [5, 6; 7, 8]; octave:> C = [A, B] C = 1 2 5 6 3 4 7 8 octave:> D = [A; B] D = 1 2 3 4 5 6 7 8
行列の各要素はインデックスを使って指定することができます。: は全てのインデックスを表します。従って A(1,:) は第１行の行ベクトルになります。
octave:> A(1,1) ans = 1 octave:> A(1,:) ans = 1 2 octave:> A(:,1) ans = 1 3
行列計算の応用として多変量解析の主成分分析を計算してみましょう。次の表は、「Excel で学ぶ多変量解析入門」(管民郎著、オーム社)に載っていた例で生徒の三教科のテストの成績です。
| 数学 | 国語 | 英語 | |
| A | 90 | 65 | 55 |
| B | 46 | 50 | 50 |
| C | 63 | 63 | 58 |
| D | 50 | 40 | 48 |
| E | 55 | 32 | 40 |
先ず生のデータの行列を作成します。
octave:> A = [90, 65, 55; > 46, 50, 50; > 63, 63, 58; > 50, 40, 48; > 55, 32, 40] A = 90 65 55 46 50 50 63 63 58 50 40 48 55 32 40
次にこのデータの学課の成績の間の相関係数を計算します。
octave:> C = corrcoef( A ) C = 1.00000 0.67781 0.50006 0.67781 1.00000 0.95435 0.50006 0.95435 1.00000
相関係数行列の固有値と固有ベクトルを計算します
octave:> [V, D] = eig( C ) V = -0.19761 0.84166 0.50256 0.75072 -0.19974 0.62970 -0.63038 -0.50171 0.59238 D = 0.02021 0.00000 0.00000 0.00000 0.54106 0.00000 0.00000 0.00000 2.43873
行列 D をみると固有値が最も大きいのは 2.43873 でその固有ベクトルは行列 V の第3列の縦ベクトルで全て正であることが分かります。したがって、この縦ベクトルの要素が生徒の総合学力を計算するための各教科の得点の係数になります。そこでこれを生データの行列に掛けてやって評価点を計算します。まず、行列 V の第3列を取りだします。
octave:> v1 = V(:, 3) v1 = 0.50256 0.62970 0.59238
また生データの各得点からその平均値を引いて標準偏差で割った基準値の行列 Ascore を次のように作ります。下の手順のうち std() は各列の標準偏差を計算する関数。kron() は行列のクロネッカー積で行ベクトルをくり返す行列を作成するのに利用します。center() 関数は各要素の価から列の平均値を引いた行列をつくる関数。./ 演算子は行列の対応する要素どうしの割り算です。
octave:> sd = std( A ) sd = 17.5129 14.3003 6.9426 octave:> Asd = kron(sd, ones(5,1)) Asd = 17.5129 14.3003 6.9426 17.5129 14.3003 6.9426 17.5129 14.3003 6.9426 17.5129 14.3003 6.9426 17.5129 14.3003 6.9426 octave:> Ascore = center( A ) ./ Asd Ascore = 1.66735 1.04893 0.69138 -0.84509 0.00000 -0.02881 0.12562 0.90907 1.12349 -0.61669 -0.69928 -0.31688 -0.33119 -1.25871 -1.46919
この Ascore と v1 を掛けると生徒の第１主成分得点が計算されます。第１主成分の関係式の係数は全て正なので、第１主成分の得点は総合学力を表していると考えられます。
octave:> Ascore * v1 ans = 1.90801 -0.44177 1.30111 -0.93798 -1.82937
この結果から、総合学力は A が最も高いことが分かります。また二番目に固有値の大きい固有ベクトル v2 は次のようになります。
octave:> v2 = V(:, 2) v2 = 0.84166 -0.19974 -0.50171
この固有ベクトルで示される主成分分析の関係式の係数は数学が正で、国語と英語は負です。従ってこの係数で計算される主成分得点は正ならその生徒が理系であることを示し、負なら文系であることを示していると考えられます。各生徒の第２主成分得点は次のようになります。
octave:> Ascore * v2 ans = 0.84695 -0.69683 -0.63952 -0.22038 0.70978
これを見ると A と E の生徒は理系で、B、C、D の生徒は文系であることが分かります。
データマイニングで有名になった多変量解析ですが、ほとんどの手法が相関係数行列の固有値の計算を用いています。Octave を使うとこのように簡単に相関係数行列の固有値と固有ベクトルを手に入れることができます。