Scilab を学ぶ

作成日 : 2022-05-14
最終更新日 :

Scilab とは

Scilab は「サイラボ」と読む。Scilab はオープンソースの数値計算システムである。 商用の数式処理システム MATLAB との類似機能をもつ。制御システムの設計に使われることも多い。 私は、主に数値解析の観点から調べていこうと思う。

Scilab ヘルプ

2023-07-26 現在最新のバージョンは、Scilab Version 2023.1.0 このバージョンのヘルプは、Scilabヘルプ であり、 大分類は次のようになっている。カッコ内は私の注釈だ。バージョン 6.1.1 との違いは少ないが、6.1.1 にあったARnoldi PACKage(大規模行列固有値問題を解くパッケージ、ARPACK)が本バージョンでは省かれている。

最初にインストールしたのは Scilab Version 6.1.1 である。6.1.1 のヘルプは、 Scilabヘルプ (help.scilab.org) にある。幸いなことに、かなりの部分で日本語訳がある。このヘルプの大分類は次のようになっている。 カッコ内は私の注釈だ。

全部訳すつもりはなかったが、意地になって訳してしまった。Atom をアトムとしただけなど、 訳にもなっていないものもある。この大分類の下に、下位分類が何段階もある。

行列の操作は上記の Linear Algebra だけではなく、 Elementary Functions や CACSD にもあることに注意が必要だ。

初等関数 - 初等行列

初等関数 - 初等行列

フランク行列

Scilabヘルプ >> Elementary Functions >> elementarymatrices に、 testmatrix という関数がある。これは、「ヒルバート,フランクのような特殊行列を作成」 するものだ。[y] = testmatrix(name, n) という形式で、 name には 'magi', 'frk', 'hilb' しか選べない。'magi' は魔方陣を作る。 'frk' は Scilab のページによればフランク行列(Franck matrix)を返すのだが、 フランク行列は聞いたことがない。試してみよう。 

--> testmatrix('frk',5)
ans  =

  5.   4.   3.   2.   1.
  4.   4.   3.   2.   1.
  0.   3.   3.   2.   1.
  0.   0.   2.   2.   1.
  0.   0.   0.   1.   1.

フランク行列とは、どうやら、次の式で表される`n` 次正方行列のことをいうらしい。

`A= (a_(ij)) = {(n+1-j, i le j),(n-j,i = j - 1),(0, otherwise):}`
私の調べた限り、この行列は、Franck 行列ではなく、 Frank 行列という名前のはずだ（スペリングに注意）。これ以上は各自で調べてほしい。 なお、「ヒルバート」とはヒルベルト行列のことだ。

固有値と固有ベクトル

固有値・固有ベクトル・特異値分解

`n` 次の正方行列 `A` の固有値(eigenvalue) `lambda` と固有ベクトル(eigenvector) `x` は次の等式

`Ax = lambdax`
を満たす `lambda` と `x` として定義される。この固有値と固有ベクトルを求める関数として、 Scilab では、eigs と spec が定義されている。どちらも固有値と固有ベクトルを求めることができるが、 ヘルプの説明を見ていると、 eigs の説明は「行列の最大固有値と固有ベクトルを計算」であり、spec は「行列とペンシルの固有値」である。 両者はどう違うのか、どちらを使えばいいのかがわからない。少し調べてみた。

まず、eigs について調べる。一般には `n` 次正方行列の固有値は全部で重複をこめて `n` 個ある。このうちの絶対値が最大となる固有値を eigs の最大固有値といっている。 実際には特に指定のない場合は絶対値の大きな順から最大計 6 個の固有値と固有ベクトルを求める。

一方で、spec について、ペンシルというのが聞き慣れないことばである。ペンシルとは行列束 (matrix pencil)のことである。伊理正夫「線形代数Ⅱ」p.279 から引用すると次のとおりである。

`K= CC` の上の同じ型 (`m times n` 型) の二つの行列 `A, B` に対して， 行列の族 `{muA - lambdaB | lambda, mu in CC; lamda = mu = 0` ではない `}`のことを `A` と `B` から作られる行列束（matrix pencil; pencil of matrices）と呼ぶ （後略）

なお、Scilab では `A - sB` の型の行列をペンシルと呼んでいる。

さて、今考えている行列 `A` の固有値と固有ベクトルを求めるにはどちらを使えばいいのだろうか。この答は、 https://help.scilab.org/docs/6.1.1/en_US/eigs.html にある次のコメントでわかる。これを抄訳する。

この eigs のページには参照「spec を見よ」とあるが、eigs と spec との違いがわからないかもしれない。 私の知る限り eigs は ARnoldi パッケージモジュール (ARPACK) を使っている。eigs は疎行列の扱いに優れるが、 spec は密行列を対象にしている。ARPACK は疎行列を扱う工夫がこらされていて、 行列の特性に応じたアルゴリズムをうまく選んでいる。 コメント求む。

核と像

kernel では核（カーネル・ヌル空間）が求められる。 なっとくする演習・行列ベクトルの演習問題 7.7[1](i) から引用する。

`A=((3,2,8),(-1,1,-1),(0,1,-1))` で定義される線形写像 `f_A:RR^3 -> RR^3` を考える。 `"Ker" f_A` と `"Im" f_A` の 1 組の基底を求めよ。 

--> A=[
3 2 8
-1 1 -1 
0 1 1
];

--> W=kernel(A);

--> W
W  = 

-0.8164966
-0.4082483
 0.4082483

こんな値が出てしまったがいいのだろうか。書籍の答は次のように出ている。 p.265 には紙面の都合で行ベクトルで表すと出ている。

`lt (-2, 1, 1)gt`

よく見ると、書籍の答のベクトルを長さが 1 になるよう正規化したものが Scilab で得られたベクトルである。 ではこんどは像を求めてみよう。そう思ったが、像を求める関数が見当たらない。 私が調べた結果では、rowcomp という関数が使える。[W,rk]=rowcomp(A) は、 X の行ベクトルの最初の dim 行に像を表すベクトルを格納するという。

この Scilab では像 (image) ということばは出てこない。 代わりに 範囲 (range) ということばを使っているらしい。 ともかく、求めてみよう。 同書の解答では、`lt (3, -1, 1), (2, 1, 1) gt` なのだが、はたしてどうか。 数値計算の結果だからということだが、わからない。少し計算させてみた。 

--> A=[
3 2 8
-1 1 -1 
0 1 1
];

--> [X,dim]=range(A,1);
X  = 

-0.9850394   0.1161092  -0.1273424
-0.0335558  -0.854035   -0.5191322
-0.1690309  -0.5070926   0.8451543
dim  = 

2.

本当に、`((-0.9850394), (0.1161092), (-0.1273424))` と `((-0.0335558),( -0.854035),( -0.5191322))` が像なのだろうか。少し計算すると、 誤差の範囲内で、 `((-0.9850394), (0.1161092), (-0.1273424)) = -0.2434516 ((3),(-1),(0)) -0.1273424((2),(1),(1))` および `((-0.0335558),( -0.854035),(-0.5191322)) = 0.3349028((3),(-1),(0)) -0.5191322((2),(1),(1))` が成り立つから、きっと像なのだろう。それにしても、まったくどのようなアルゴリズムを使っているかわからないから困る。 少なくとも宿題の解答には使えないだろう。

なお、`A` が正方行列ならば、そしてたまたまここでは `A` が正方行列なのだが、range が使える。 [X,dim]=range(A,k) は、正方行列 `A^k` の像を、 X の行ベクトルの最初の dim 行に像を表すベクトルを格納するという。結果は同じだ。

数値微分

数値微分に関しては、 Differential Equations, Integration の numderivative を使う。diff も使える。 なお、記号微分は多項式に関しては使える。Polynomials にある derivat — 有理行列の微分を参照してほしい。 しかし、一般の超越関数が入る記号微分は見当たらない。

統計

大分類「統計」の下の小分類「Cumulated Distribution Functions（累積分布関数）」に、 cdft （ステューデントの `t` 分布）がある。これには、次のようなインターフェースがある。

[P,Q]=cdft("PQ",T,Df)
[T]=cdft("T",Df,P,Q)
[Df]=cdft("Df",P,Q,T)

自由度 `Df` の `t` 分布の確率密度関数 `t(x; Df)` とすると、P (または Q), T, Df の３変数は次の式を満たす：

`Q = int_-oo^T t(x; Df) dx , P = 1 - Q`

したがって、P (または Q), T, Df のうちの 2 つの値がわかれば、残りの 1 つの値がわかるということを意味する。 実際に、t 分布表の左上にある、`n = 1、alpha = 0.100` に相当する t = 6.314 を例にとろう。

まず、第１の式は、`t` 分布の値 `t` と自由度 `Df` がわかれば、下側確率 P と 上側確率 Q が求められることを意味する。 Scilab の入出力は次のようになる。

   --> cdft("PQ",6.314,1)
    ans  =
   
      0.9500019
   
   --> [P,Q]=cdft("PQ",6.314,1)
    P  = 
   
      0.9500019
    Q  = 
   
      0.0499981

P = 0.90 ではなく、Q = 0.10 でもないのは、私が見た表が、`int_t^oo f(x;df)dx = alpha//2` となる `t` を求める表だからで、 表の値も Scilab の値も正しい。

次に第２の式は 上側確率 P 、下側確率 Q と自由度 `Df` がわかれば、`t` 分布の値 `t` がわかることを意味する。 Scilab の入出力は次のようになる。

    --> cdft("T",1,0.95,0.05)
        ans  =
        
          6.3137515
    --> cdft("T",1,0.95)

      cdft: 入力引数の数に誤りがあります: 4 個の引数を指定してください.
          
    --> cdft("T",1,0.95,0.15)
          
      cdft: P + Q ≠ 1

この場合、引数は 4 つ必要だ。どうせ Q = 1 - P だから第4 変数を省略してもいいだろう、 と考えて引数を 3 つにすると叱られる。また `P + Q != 1` となる P と Q を入れても叱られる。

最後に第３の式は 上側確率 P 、下側確率 Q と`t` 分布の値 `t` がわかれば、自由度 `Df` の値が求められることを意味する。 Scilab の入出力は次のようになる。

    --> cdft("Df",0.95,0.05,6.314)
        ans  =
        
          0.9999768

これで統計表の値を検算できるほか、自由度が非整数の場合の t 値から上側確率や下側確率を計算できるので、 Welch の t 検定ができるようになる。

リンク集

まりんきょ学問所 ＞ コンピュータの部屋 ＞ Scilab を学ぶ

MARUYAMA Satosi