［Octave］固有値、固有ベクトルを求める

eig関数を使う。固有値、固有ベクトルを求めたい行列を引数に与えると、固有ベクトルと固有値が一緒に返ってくる。

以下に求めた例を示す。見てのとおり、対称行列でなくても求めることができる。

>> mxaa = [3, -2; 1, 0;]
mxaa =
   3  -2
   1   0
>> [mxxx, l] = eig(mxaa)
mxxx =
   0.8944   0.7071
   0.4472   0.7071
l =
Diagonal Matrix
   2   0
   0   1

1つ目の戻り値が固有ベクトル、2つ目が固有値。例えば1つ目の固有値と固有ベクトルを取り出すには、以下のようにする。

>> l(1, 1)
ans = 2
>> mxxx(:, 1)
ans =
   0.8944
   0.4472

2つ目の固有値と固有ベクトル。

>> l(2, 2)
ans = 1
>> mxxx(:, 2)
ans =
   0.7071
   0.7071

以下、2つの行列の固有値と固有ベクトルを求めた例。

>> mxaa = [6, -3, 5; -1, 4, -5; -3, 3, -4]
mxaa =
   6  -3   5
  -1   4  -5
  -3   3  -4
>> [mxxx, l] = eig(mxaa)
mxxx =
   7.0711e-01  -3.0151e-01  -3.4940e-15
   7.0711e-01  -9.0453e-01   8.5749e-01
   2.7081e-15  -3.0151e-01   5.1450e-01
l =
Diagonal Matrix
   3   0   0
   0   2   0
   0   0   1
>> mxaa = [3, -1, 1; 0, 2, 1; 0, 0, 3]
mxaa =
   3  -1   1
   0   2   1
   0   0   3
>> [mxxx, l] = eig(mxaa)
mxxx =
   1.0000   0.7071        0
        0   0.7071   0.7071
        0        0   0.7071
l =
Diagonal Matrix
   3   0   0
   0   2   0
   0   0   3

［Octave］警告メッセージ「warning: Invalid UTF-8 byte sequences have been replaced.」

ソースコードが書かれたファイルを適切な文字コードで保存をしていない。Windows版Octaveは初期状態ではUTF-8を扱う設定になっているため、ソースファイルをUTF-8で保存し直す。以下、実験。

以下の2行をテキストエディタに貼り付けてシフトJISで保存をする。

% 日本語
fprintf("ABCあいうXYZ\n")

GUI版で実行してみる。

>> nihongo
warning: Invalid UTF-8 byte sequences have been replaced.
warning: Invalid UTF-8 byte sequences have been replaced.
ABC■■■XYZ

エラーメッセージが表示される。ソースファイルをUTF-8で保存し直して同じように実行する

>> nihongo
ABCああいいううXYZ

警告メッセージは表示されなくなったが、表示がおかしい。現在の最新版のOctaveはマルチバイト文字の処理がうまく行えず、日本語をうまく表示することはできない。

なお、テキストファイルはシフトJISのままで、

1. メニュー「編集」→［設定」
2. タブ「エディタ」の下の方にある「読込，保存に使用するテキストのエンコード」が「UTF-8」になっているはず。これを「SHIFT-JIS」にして、「OK」ボタンをクリック。

と設定し、同じように動作させると、

>> nihongo
ABCああいいううXYZ¥n>>

警告メッセージは表示されなくなったが、依然表示がおかしくマルチバイト文字の処理がうまく行えていないことがわかる。これらはCLI版でも同様である。

[Octave］t分布におけるパーセント点を求める

tinv関数を使う。以下は自由度5のt分布において、上側5%点（=2.015048…）、両側5%点（=2.570582…）を求めた例。

>> pkg load statistics
>> tinv(1 - 0.05, 5)
ans = 2.0150
>> tinv(1 - 0.05 / 2, 5)
ans = 2.5706

tinv関数は下側パーセント点の値しか得られないため、上側パーセント点の値を得るには1からその確率の値を引くことで得られる。両側パーセント点は下側パーセント点の2倍のため、それを2倍している。

［Octave］t分布における分布関数の値を求める

tcdf関数を使う。以下は自由度5のt分布において、確率変数Xが2における上側（=0.05096974…）、両側（=0.1019395…）の分布関数の値を求めた例。

>> pkg load statistics
>> 1 - tcdf(2, 5)
ans = 0.050970
>> 2 * (1 - tcdf(2, 5))
ans = 0.1019

tcdf関数は下側の値しか得られないため、上側の値を得るには1からその値を引くことで得られる。両側は下側の2倍のため、それを2倍している。

［Octave］代数方程式の根を求める

roots関数を使う。以下は2x3+3x2+8x-5=0の3つの根（0.5,-1+2i,-1-2i）を求めた例。虚部も求めることができる。引数には、次数の高いほうから係数を行列でまとめて与える（この場合は[2 3 8 -5]とする）。

>> roots([2 3 8 -5])
ans =
  -1.0000 + 2.0000i
  -1.0000 - 2.0000i
   0.5000 +      0i

戻り値は複素数型。複素数型から実部だけを取り出すにはreal関数を、虚部だけを取り出すにはimag関数を使う。

>> ri = roots([2 3 8 -5]);
>> iscomplex(ri)
ans = 1
>> real(ri)
ans =
  -1.0000
  -1.0000
   0.5000
>> imag(ri)
ans =
   2
  -2
   0

［Octave］固有値と固有ベクトルを求める

eig関数を使う。固有値と固有ベクトルが行列でそれぞれ戻り値になるため、戻り値は2つ指定する。

>> mxaa = [3, 1, 1; 1, 2, 0; 1, 0, 2]
mxaa =
   3   1   1
   1   2   0
   1   0   2
>> [mxv, mxl] = eig(mxaa)
mxv =
  -5.7735e-01  -1.7174e-16   8.1650e-01
   5.7735e-01  -7.0711e-01   4.0825e-01
   5.7735e-01   7.0711e-01   4.0825e-01
mxl =
Diagonal Matrix
   1   0   0
   0   2   0
   0   0   4

この例では固有値と固有ベクトルは3つあり、例えば3番目の固有値と固有ベクトルを取り出すには、以下のようにする。

>> mxl(3, 3)
ans = 4.0000
>> mxv(:, 3)
ans =
   0.8165
   0.4082
   0.4082

［Octave］t分布における確率密度関数の値を求める

tpdf関数を使う。以下は確率変数が-5～5の範囲において、自由度1と5のt分布における確率密度関数の値を求めてグラフにした例。

>> pkg load statistics
>> xx = -5:0.1:5;
>> hold off
>> plot(xx, tpdf(xx, 5), 'b')
>> hold on
>> plot(xx, tpdf(xx, 1), 'k')
>> text(1, 0.3, "n=5")
>> text(0.8, 0.1, "n=1")

［Octave］QR分解を行う

行列AをQR分解で行列Qと行列Rに分解することを考える。

A = Q R

行列QはQ^TQ=I（Iは単位行列）となる直交行列。行列Rは上三角行列（対角成分はすべて0以外で、対角成分より下の要素はすべて0）。

qr関数を使うと、それぞれQR分解した際の行列Qと行列Rが得られる。

>> a = [1, 2, 3; 4, 5, 6; 9, 9, 8]
a =
   1   2   3
   4   5   6
   9   9   8
>> [q, r] = qr(a)
q =
  -0.10102  -0.71841  -0.68825
  -0.40406  -0.60254   0.68825
  -0.90914   0.34762  -0.22942
r =
   -9.89949  -10.40457  -10.00051
    0.00000   -1.32095   -2.98951
    0.00000    0.00000    0.22942

Q^TQ = I となるか試す。

>> q' * q
ans =
   1.0000e+00  -9.5469e-17   3.7850e-17
  -9.5469e-17   1.0000e+00   7.4583e-17
   3.7850e-17   7.4583e-17   1.0000e+00

［Octave］シリアル日付値について

Octave（MATLAB）では、日付型変数の扱い方の1つとして、シリアル日付値というものがある。これは、個々の日を西暦元年1月0日からの日数として表したものである。datenum関数に年、月、日を与えると、このシリアル日付値を返す。逆にdatevec関数にシリアル日付値を与えると、それが示す日付をベクトルで返す。

代表的なシリアル日付値は以下のとおり。

西暦元年1月1日　→　1
西暦元年1月2日　→　2
西暦2年1月1日　→　367
西暦2019年1月1日　→　737426

以下は、Octaveでの実行例。

>> datenum(0, 1, 1)
ans =  1
>> datenum(0, 1, 2)
ans =  2
>> datenum(1, 1, 1)
ans =  367
>> datenum(2019, 1, 1)
ans =  737426
>> datevec(1)
ans =
   0   1   1   0   0   0
>> datevec(2)
ans =
   0   1   2   0   0   0
>> datevec(367)
ans =
   1   1   1   0   0   0
>> datevec(737426)
ans =
   2019      1      1      0      0      0

［Octave］標準正規分布に従う確率変数（乱数）を得る

randn関数を使う。引数には、発生させる乱数の数を与える。

> randn(1)
ans =  2.9256
> randn(3)
ans =
   0.14849  -1.65450  -0.63153
   1.00336   0.56755   0.23638
   0.13845  -0.40872   1.38998

引数が1つの場合は、乱数からなるその次数の正方行列を返す。引数を2つ与えると、任意の大きさの行列を返す。

> randn(3, 2)
ans =
   1.229068  -0.953506
   0.579371  -0.849324
  -0.082636   0.257047