第1章 何ができるか:MATLAB

専修大学人間科学部 情報処理心理学実習 で使用している教科書 『MATLAB プログラミング入門』 の補足説明を行います.

注意:数式が適切に表示されるまで,少し時間がかかるかもしれません.

1.1 ここから始めよう

初回の授業ではこのセクションの内容について学習します. MATLAB の起動,授業用フォルダの作成,簡単なプログラムの実行, を行います.以下のページを利用します.

1.2 数の計算

1.2.1 式の通りに書けばよい:+, -, *, /

四則演算の式の書き方は, みなさんが使い慣れているであろう Excel と同じです. テキストの例題 ex1201 を行ってください. 短いので,コマンドウィンドウでキーインすればよいでしょう. 以下のように,演算記号の前後にスペースを空けてもいいです. この方が,式が読みやすいかもしれません.


2.0 + 3.0 * (5.0 - 16.0 / 5.0)

テキストの旧版では,掛け算と割り算にはそれぞれ .* と ./ が用いられていました. 通常の演算記号の前にドット(.)がついています. これは,行列(matrix)やベクトル(vector)において, 要素ごとの演算を行うための演算記号です. テキストの 1.3.3 節,および,2.3.3 節に説明があります. 数の掛け算と割り算に .* および ./ を用いてもよいのですが, ドットがついていない演算記号の方がなじみがあるでしょうから, ここでは * および / を使っておきましょう.

1.2.2 1行に収まらないときは:...

テキストに書いてあるとおり, MATLAB ではひとつの実行文を1行に書くことを原則としています. 改行すると,MATLAB はその1行に書いてあることを解釈して, 処理を実行します

実行文が長くて2行以上に分けて書きたいときには. 行の終わりにピリオド(ドット)を3つ打って改行します. こうすることで,MATLAB はまだ文が続くということを理解し, その行の処理をすぐに実行しないで文の終わりを待ってくれます.

例題 ex1202 を行ってください. 短いので,コマンドウィンドウでキーインすればよいでしょう. 行末の連続ピリオドの前に半角のスペースを入れた方が, 式が見やすくなります. 改行のための連続ピリオドの直前が整数だと, 最初のピリオドが小数点と解釈されてエラーになってしまいます. 改行の連続ピリオドの前には半角スペースを入れる習慣にしておきましょう

テキスト(改訂版第1刷)では掛け算と割り算の演算記号にドットがついていますが, これは不要です.通常の演算記号を使ってください. テキストの旧版の記述ではドットつきの演算記号でよかったのですが, 改訂版で変更すべき演算記号が,いくつかそのままになっているようです. 今後の増刷で修正されるかもしれません.


12345.6789 + 9876.54321 * (12345.6789 - 98765.4321) / ...
12345.6789)

1.2.3 結果を表示しない

ひとつの文を複数行にわたって書くのとは逆に, 複数の文を1行に書くこともできます.

同じ行で文を区切るのは「,」または「;」です. 「;」をつけると値をコマンドウィンドウに表示しません. (p.10)

たくさんの結果が得られる文(たとえば, とても多くの乱数を発生させる命令)では, 「;」を付け忘れるとコマンドウィンドウが出力であふれてしまいます. 注意しましょう. 特に忘れやすいのは行末です. 文末には「;」をつける,と覚えておいてもよいかもしれません. 結果を表示させたいときだけ「;」を外します.

どのプログラミング言語でも, プログラムの中にコメントを書いておくための記号があります. MATLAB ではそれが「%」です.

「%」で始まる行はコメントとして扱われ,計算には関係ありません. (p.11)

コメントは何らかのメモを書いておくものです. たとえば, 複雑ですぐには意味がわからないコードに対して, 何を計算するためのコードなのかをメモしておくことができます. 入門の段階では,プログラムは短くて簡単でしょうから, コメントを書いておく必要はあまりないかもしれません. 第4章で関数を書くときに使うぐらいでしょう.

MATLAB は,% 記号があると, そこから行末までをコメントと解釈して無視します. 例題 ex1203 では独立した1行にコメントが書かれていますが, 実行したい計算が書いてあるすぐ右側にコメントを書くこともできます. たとえば,コマンドウィンドウで次のようにキーインしてみてください.


3 + 5 % コメント

% からのコメントは無視され,3 + 5 という計算が実行されます.

1.2.4 べき乗:^

べき乗とは,同じ数を複数回掛け算することです. たとえば,「2のべき乗」とは, 2 2 , 2 3 , 2 4 , のことをを意味します.MATLAB での書き方では, 2^2, 2^3, 2^4, ... のことです. 数学ではよく使う言葉ですが, 中学と高校の数学ではこの言葉を習っていないかもしれません. もちろん,言葉を習っていなくても, 中学数学で扱う計算です.

例題 ex1204b に使われている realmax および realmin はそれぞれ, 次のセクション 1.2.5 で出てくる prod 関数と同様に, MATLAB の関数のひとつです. しかし,数学で円周率を π という記号で表すのと同様に, realmax および realmin は定数を表す記号だと考えても差し支えありません. テキストには,

realmax と realmin は使っているコンピュータで MATLAB が扱える数の最大と最小を与える定数です.(p.11)

と書かれています. なお,MATLAB での円周率を表す記号は pi です. コマンドウィンドウで pi とキーインしてみてください. 円周率が表示されます.

例題 ex1204b の出力では, 指数表記という数の表記法が用いられています. Excel でも使われている表記法ですが,知らない人もいるかもしれません. テキストにも書かれているように, 1.7977e+308 は 1.7977 × 10 308 を意味します. これは 1.7977 の小数点を右側に308ずらした数ですね. 2.2251e-308 は 2.2251 × 10 -308 を意味します. マイナス乗の意味を知らない(忘れた)かもしれませんが, 難しくはありません. これは 2.2251 の小数点を左側に308ずらした数です.

1.2.5 階乗:prod

みなさんはおそらく Excel で関数を使った経験があるでしょう. MATLAB の関数も Excel の関数と同じ役割を果たします. すなわち,入力を与えると,決められた処理を実行し, 処理の結果を返します

Excel で階乗を求めるには,関数 FACT を使います. セルに =FACT(10) と入力すれば, 10! が計算され,結果が表示されます.

MATLAB には,階乗を求めるための専用の関数はありません. しかし,prod 関数を用いれば階乗を計算することができます. prod 関数の "prod" は product(積)を意味しています. この関数は,数値の集まりを与えられたとき, それら数値をすべて掛け算した結果を返します. コマンドウィンドウで,


prod([1 3 5 10])

とキーインしてみてください. これら4つの数の積である 150 が返されます.

prod 関数については,1.4.3 節(p.21)でもう一度取り上げます.

Excel での関数と同様に,MATLAB での関数への入力は, 関数名の右側にある ( ) の中に書きます. 上の例では, 1 3 5 10 という ベクトル(vector)が入力です. これは数を並べたものとして理解しておけば十分です. MATLAB でベクトルを作るには, 並べたい要素を半角のスペースで区切って並べ, 全体を [ ] で囲みます. スペースの代わりに,コンマ(,)で数を区切って, 1 3 5 10 と書くこともできます.

例題 ex1205 での prod 関数への入力である 1:10 は, 「1から10までの数を並べたベクトルを作れ」 という命令です.コマンドウィンドウで


1:10

とキーインして,確認してください.

ベクトルについては第2章で学習します.

1.2.6 複素数を表すには:i, j

数学で2次方程式の解の公式を習ったとき, 2乗するとマイナス1になる数をアルファベット小文字の i で表しました.すなわち,

i = - 1

です.テキストに書いてあるように,この数を虚数単位と呼びます.

MATLAB では,ij が虚数単位を表す記号です. コマンドウィンドウで,i^2 および j^2 とキーインしてみてください. いずれも -1 という値が返されます.

ab を実数として, a + b i と表される数を複素数と呼びます. b i は虚数単位 ib 倍です. 最近の高校数学のカリキュラムでは複素数を学習しています.

例題 ex1206a および ex1206b をやってみてください. コマンドウィンドウでキーインすればよいでしょう. テキスト(改訂版第1刷)では掛け算の演算記号にドットがついていますが, これは不要です.例題 ex1206b には誤植があります. ans = 5.0000 - 2.0000i ですね.

コマンドウィンドウで i あるいは j とキーインすると, 0.0000 + 1.0000i という結果が返されます.

MATLAB には2つの実数 a, b から複素数 a + b i を構成する関数 complex が用意されています(複素数のことを,英語で complex number と書きます). コマンドウィンドウで complex(5, -2) とキーインしてみてください. 例題 ex1206a および ex1206b と同じ結果が得られます.

1.2.7 ふつうの式の通りに書けばよい

複素数の四則演算規則は以下のようになります. 乗法および除法の計算では i 2 = - 1 に注意してください. 除法の計算では,最終形を α + β i という形にするために,分子と分母に a - b i を掛け算しています.

加法
a + b i + c + d i = a + c + b + d i
減法
a + b i - c + d i = a - c + b - d i
乗法
a + b i c + d i = a c + a d i + b c i + b d i 2 = a c - b d + a d + b c i
除法
c + d i a + b i = c + d i a - b i a + b i a - b i = a c - b c i + a d i - b d i 2 a 2 - b 2 i 2 = a c + b d a 2 + b 2 + a d - b c a 2 + b 2 i

例題 ex1207 をやってみてください. コマンドウィンドウでキーインすればよいでしょう. テキスト(改訂版第1刷)では掛け算および割り算の演算記号にドットがついていますが, これは不要です.

1.2.8 複素共役を作るには:'

複素数 a + b i に対して, a - b i を「共役(きょうやく)複素数」あるいは「複素共役」と呼びます. a - b i の共役複素数は a + b i です.

複素数 z = a - b i を共役複素数に変換するには, z にアポストロフィ(')をつけてやります. これは,+ や * と同じく,特定の演算を行う記号です. 複素数 z に対して,共役複素数に変換するという演算を実行します. 特定の演算を行う記号を, かたい言葉で 演算子(operator)と呼びます.

例題 ex1208 をやってください. コマンドウィンドウでキーインすればよいでしょう. 例題での変換とは逆に,


Y = 10 - 6*i;
X = Y'

もやってみましょう.

共役複素数を求める演算子「'」は, 本来は行列の転置を行うものです. これについては 2.3.4 節に解説があります.

1.3 変数とベクトルを使う

1.3.1 変数の名前

プログラミングにおいて, 変数は扱いたいものをしまっておく箱にたとえられます. たとえば,例題 ex1207 での X = 10 + 6*i という式は, 10 + 6 i という複素数を作り, それを X という名前の箱にしまっていると考えられます.

MATLAB での変数名には,アルファベットの大文字, アルファベットの小文字,数字, 下線( _ )を用いることができます. 変数名はアルファベットで始まる必要があります (大文字でも小文字でもよい). 大文字と小文字は区別されます. A と a は別の変数です. テキスト p.14 に書かれているように, MATLAB で用意されている関数名や変数名は小文字で始まりますので, ユーザーが自分で作成する変数名は大文字で始めるとよいでしょう

作成した変数はワークスペースに保存されます. コマンドウィンドウの右側にある 「ワークスペースブラウザ」に, ワークスペースに保存されている変数の一覧が示されています. 下図は,例題 ex1302a を実行した直後のワークスペースブラウザです. この例題で用いられた5つの変数(A, B, C, X1, X2)が示されています.

workspace browser

コマンドウィンドウで who あるいは whos とキーインすると, ワークスペースにある変数のリストが示されます. 試してみてください.

ここから先の例題では,同じ名前の変数が繰り返し使用されます. 前の例題で作成した変数がワークスペースに残っていると, 混乱のもとになりかねません. たとえば,変数を作成し忘れたまま, 前の例題で作成した同じ名前の変数を使ってしまうかもしれません. そこで,新しい例題に取り組むときには, ワークスペースを空にしておくとよいでしょう. ワークスペースを空にするには, コマンドウィンドウで clear とキーインします.

1.3.2 2次方程式を解く

例題 ex1302a で解いているのは,二次方程式 x 2 - 3 x + 4 = 0 です.公式を使って得られる2つの解を, それぞれ変数 X1 および X2 に代入しています. この代入を行っているプログラムの2行目と3行目の最後には, セミコロン(;)がつけられていませんので, これら変数に代入されている値(解)が表示されます.

プログラムの1行目では, A, B, C の値を小数点以下第1位の 0 まで入力していますが, これは単に


A = 1; B = -3; C = 4;

としてかまいません.

プログラムの2行目と3行目での,


(B^2 - 4*A*C)^0.5

という部分は,平方根を求める sqrt 関数を用いて,


sqrt(B^2 - 4*A*C)

と書くこともできます. ある数の 0.5 乗とは,その数の正の平方根のことです.

例題 ex1302a で,変数 X1 および X2 に入れられた値が, 2次方程式 x 2 - 3 x + 4 = 0 の解であることを確認しておきましょう. X1 および X2 を,この2次方程式の左辺に代入して, その値が 0 になればいいですね. コマンドウィンドウで,


A * X1^2 + B * X1 + C

および


A * X2^2 + B * X2 + C

とキーインしてください.いずれも 0 が返されます.

例題 ex1302b で使われている roots 関数は方程式の解を求める関数です. 方程式の解のことを根(こん)と呼びます.英語で root(s) です. この関数への入力はベクトルです. n + 1 個の数字を並べたベクトルを与えると, n 次方程式の解が返されます.例題 ex1302b では, 入力が 1.0 -3.0 4.0 というベクトルであり,roots 関数は2次方程式 x 2 - 3 x + 4 = 0 の解を返します.

1.3.3 たくさんの2次方程式を解く

テキスト p.16 の 1.3 式

x 2 - 3 x + 4 = 0 3.5 x 2 + 5.7 x + 8 = 0 5 x 2 - 5 x + 8 = 0

は連立方程式ではありません.3つの2次方程式であり, それぞれ解を持ちます.

このセクションで説明していることのポイントは, 2次方程式の解き方ではなく, 同じ計算を繰り返すのならベクトル(あるいは行列)を効果的に使うとよい ということです. ベクトルを使えば, 3つの2次方程式をひとつずつ解くのではなく, まとめて解くことができます.

例題 ex1303 は 1.3.2 節の例題 ex1302a とほとんど同じです.ただし, 変数 A, B, C それそれに格納されているのはスカラー(ひとつの数)ではなく, 3つの数をひとまとめにしたベクトル(1.4 式)であることに注意してください.

この例題での掛け算,割り算,べき乗はそれぞれ, 演算記号の前にドット(.)がついています. 1.2.1 節で簡単に述べたように, これらは行列やベクトルにおいて要素ごとの演算を行うための演算記号です. プログラムの2行目と3行目において, 4 .* A および 2 .* A という演算はドットを書かなくてもいいですが, その他はすべてドットが必要です. 3つのベクトル A, B, C(ex1303 の1行目)をキーインした後,

がどのような結果になるか,キーインして試してみるとよいでしょう. B^2 と A * C はエラーが返ってきます.

ドットつきの掛け算,割り算,べき乗がどのような計算を実行するのか, テキスト p.17 の (1.5) 式から (1.19) 式をよく検討してください.

1.3.4 もっとたくさんの2次方程式を解く

「もっとたくさんの2次方程式」という点は特に重要ではありません. このセクションで説明したいのは,等間隔に数が並んだベクトルの作り方です. これについては 2.1.2 節および 2.2.6 節であらためて説明されます.

テキスト (1.21) 式,および,(1.22) 式で説明されているように, a:d:b という式は 「a から始まって,d きざみで,b までの数列(ベクトル)を作る」 という意味になります. たとえば,0:0.1:1 という式は,例題 1304b に示されているように, 1 0.1 0.2 1.0 というベクトルを作ります. 数列の第 k 項は, a + k - 1 d と表されます.

例題 ex1304b のベクトル B = 0:0.11:1 で示されているように, 式 a:d:b での b は, a から始めて d を順に加えていってできる数になっていなくてもかまいません. この場合,終点である d を超えない数までが作られます. たとえば,例題 ex1304b のベクトル B = 0:0.11:1 では, 数列の最終項は 0.99(第10項)となり, 次の 1.10 は終点1を超えてしまうので数列に含まれません. テキストに書かれているように(p.19), 式 a:d:b で作られるベクトルの要素の数は, a + k - 1 d b を満たす自然数 k となります.

1.4 組み込み関数を使う

1.4.1 級数を計算する:for

プログラミング言語には,一般に, 同じ処理を繰り返すための仕組みが用意されています. 「for ループ」と呼ばれるひとまとまりの記述は, 繰り返し処理のための最も基本的な仕組みです. あらかじめ決められた回数の繰り返し処理を実行するために使われます

for ループの構文は,例題 ex1401a に示されているように,


for index = values
   statements
end

というものです.for index = values と end の間に書かれた statements を, 決められた回数だけ実行します.statements はいくつあってもかまいませんが, 例題 ex1401a では S = S + I.^2 だけです. for index = values,および,end の行の最後には,セミコロンはつけません. セミコロンがついていなくても,これらの行に関しては何も表示されません.

例題 ex1401a のプログラムでの, S = S + I.^2 という式(2乗のドットは省略可)は, 数学での等式ではないので注意してください. 等号(=)は,右辺の計算を実行して, その結果を左辺に代入するという意味です. MATLAB に限らず,多くのプログラミング言語で, 等号はこうした代入を意味します.

上記 for ループの構文での index は繰り返しを制御するための変数です. 例題 ex1401a では I です. 数字を横方向に並べたベクトルが values として与えられると, まず最初の要素が index にセットされます. 例題 ex1401a では,values として与えられたのは 1 2 100 というベクトルですから,I の最初の値は 1 となります. index の値は,statements が繰り返されるたびに, ベクトルの次の値に更新されます. index の値がベクトルの最終要素の値(例題 ex1401a では 100)となり, その値のもとで statements が実行されると, for ループは終了します.

for ループの繰り返し回数は values として与えられるベクトルの要素の数です. I = 1:100 とすれば100回,I = 2:2:100 とすれば50回になります. ここで 2:2:100 は 2 4 6 100 というベクトルです.

values としては,数字を横方向に並べたベクトルだけでなく, 数字を縦横に並べた配列(行列)を与えることもできます. これについては第3章で学習します(例題 3102b).

例題 ex1401a のプログラムでは,S の値が次々と更新され, 最終的に 1 から 100 までの2乗和が求められます. そのプロセスを追ってみましょう. 1行目に S = 0 と書かれていますから,変数 S には最初, 0 という値が入っています. 一般に,変数には最初に何らかの値(初期値)が入れられます. これを変数の 初期化(initialization)と呼びます. 2行目から for ループに入ります. I の値は最初 1 にセットされます. 3行目の指示は「右辺の S + I.^2 という計算を実行して, その結果を左辺の S に代入せよ」です. S + I.^2 = 0 + 1^2 ですから, S(左辺)には 1 2 という値(= 1)が代入されます. 4行目に end があるので,1回目のループはこれで終了し, 2回目のループに入ります. I の値は更新されて 2 になります. S + I.^2 = 1^2 + 2^2 ですから, S(左辺)には 1 2 + 2 2 という値(= 5)が代入されます. このようにして, S の値は 1 2 + 2 2 + 3 2 + と更新されていきます. 最後のループの直前では,I = 99, S = 1 2 + 2 2 + + 99 2 となっています. 最後のループにおいて,I の値は更新されて100となります. S + I.^2 = (1^2 + 2^2 + ... + 99^2) + 100^2 ですので, S(左辺)には 1 2 + 2 2 + + 100 2 が代入されます.これで for ループの繰り返し処理は終わり, ループを抜けて次に進みます.

プログラムの最終行である5行目には,S だけが書かれています. セミコロンが付けられていませんので, ここで変数 S の値が画面に出力されます. この値はもちろん,1 から 100 までの自然数の2乗和です.

for ループの記述は複数行にわたるので, コマンドウィンドウで直接にキーインするのではなく, エディターを開いてプログラムを書きます. 例題 ex1401b のようにプログラムを1行に書くことはおすすめしません. 読みにくいだけです. エディターを開くには,MATLAB の [ホーム] タブにある, [新規スクリプト] をマウスで左クリックします.

テキストの例題 1401a を見ると, for と end ではさまれた部分(3行目)は字下げされているのがわかります. このような字下げ(インデント)を適切に行うと, プログラムが読みやすくなります. エディターで for I = 1:100(プログラムの2行目)と書くと, 次の行から自動的に字下げがなされます(下図).

begining of a for-loop

繰り返し行う処理の記述を終え,end と書くと, 下図のように for ループの記述がわかりやすく表示されます. 繰り返し処理の記述が長くなっても,その部分は字下げされますから, for ループの記述がどこから始まり, どこで終わっているのかがすぐにわかります.

end of a for-loop

プログラムを書き終えたら,ex1401a.m というファイル名で保存し, コマンドウィンドウで ex1401a とキーインしてプログラムを実行してください.

1.4.2 組み込み関数 sum を使う

for ループを使った繰り返し処理の実行は,MATLAB に限らず, プログラミングの基本技術です. しかしMATLAB では, for ループを使うことが最善の方法ではない場合がよくあります. ベクトルあるいは行列でデータを用意し, それに MATLAB の組み込み関数を適用した方が, よいプログラムとなります.たとえば, 前節(1.4.1)での例題 ex1401a および ex1401b は, 例題 ex1402a あるいは ex1402b のように書き換えられます. 一般に, for ループの使用を避けて MATLAB の組み込み関数を使う方が, 実行速度が速くなります.これについては第3章で学習します.

例題 ex1402a および ex1402b で使用している sum 関数は, 数の並びであるベクトルが入力として与えられると, それらの数の総和を返します.Excel の SUM 関数と同じですね. MATLAB の sum 関数は,ベクトルだけでなく, 数字を縦横に並べた配列(行列)を入力としてとることもできます(第3章参照).

例題 ex1402a および ex1402b では,ベクトル 1 2 100 の要素それぞれを2乗してベクトル 1 2 2 2 100 2 を作り,sum 関数ですべての要素の和を求めています. ベクトルの要素ごとの2乗を計算するので, .^2 のドット(.)は省略できません.

1.4.3 乗積を計算する:prod

prod 関数は 1.2.5 節(テキスト p.12)ですでに学習しています. この関数は,数値の集まり(ベクトルです)を与えられたとき, それら数値をすべて掛け算した結果を返します. コマンドウィンドウで prod([1 3 5 10]) とキーインしてみてください. これら4つの数の積である 150 が返されます

テキスト p.21 の (1.24) 式に出てくる記号 は掛け算を意味します.たとえば,(1.24) 式では, i = 1 100 1 + 1 i = 1 + 1 1 1 + 1 2 1 + 1 100 となります. 記号 i を 1 から 100 まで変化させてできる,すべての項の積を計算します. 足し算を意味するシグマ記号 は知っていますね.これの掛け算バージョンが記号 です.これはギリシャ文字 π (パイ)の大文字です.

この計算を for ループを利用して行うプログラムが, 例題 ex1403a です.すべてを1行に書いていますが, 1.4.1 節の解説でも述べたように,これはおすすめしません. 改行してプログラムを書いてください.

1.4.2 節で述べられているように, 一般に MATLAB では, for ループを回避してプログラムを書いた方がよいプログラムとなります. ベクトルあるいは行列でデータを用意し, それに MATLAB の組み込み関数を適用します. この方法で計算を実行するプログラムが例題 ex1403b です.

ones 関数はこの例題で初めて出てきた関数ですね. この関数は,指定された行数と列数の, すべての要素が 1 である行列(matrix)を生成します(2.2.2 節参照). 行列は数を縦および横の両方向に並べたものです. 並べるものは数とは限らないのですが,当面はこの理解で問題ありません. 横方向の要素の並びを「行」, 縦方向の要素の並びを「列」と呼びます. ones 関数での行数と列数の指定は,ones(行数, 列数) とします. コマンドウィンドウで


ones(2,3)

とキーインしてみてください. すべての要素が 1 の,2行3列の行列が表示されます.

例題 ex1403b では ones(1, 100) としていますから, 1行100列の行列が生成されます. これはすなわち,1 を横方向に100個並べたベクトルです. 要素を横方向に並べたベクトルを 行ベクトル(row vector)と呼びます.

例題 ex1403b の解説をしておきます. 最初に,


X = 1:100;

として,1 から 100 までの数を横方向に並べたベクトルを作ります. これを変数 X に代入します. 式の末尾にセミコロン(;)がついていますから, この式の結果は表示されません. 1 を100個表示されても邪魔なだけなので,その方がよいでしょう.

次の式,


X = 1 + ones(1,100) ./ X;

は,少し複雑ですね.ones 関数を使って作ったベクトルと, ベクトルを代入した変数 X を使った計算を実行して, その結果を変数 X に代入しなおしています.

ones(1, 100) ./ X という部分の意味を考えます. ones(1, 100) は,数字の 1 を横方向に100個並べたベクトルを作るのでした. 変数 X には,1 から 100 までの数を横方向に並べたベクトルが代入されています. ドットのついている割り算(./)は, ベクトルでの対応する要素ごとに割り算を実行するのでした(1.3.3 節を参照). つまり,ones(1, 100) ./ X という式は,


[1 1 1 ... 1] ./ [1 2 3 ... 100]

という演算を実行することになるので,その結果は

1 1 1 2 1 3 1 100

というベクトルになります.このベクトルの各要素に 1 を加える演算


1 + ones(1,100) ./ X;

により,ベクトル

1 + 1 1 1 + 1 2 1 + 1 3 1 + 1 100

が作られ,これが X に代入されます.最後に,


P = prod(X)

としてベクトル X に prod 関数を適用すれば,(1.24) 式の結果が得られます. この結果は変数 P に代入されています. 式の最後にセミコロンがついていませんので, P に代入された計算結果が表示されます.

1.5 グラフでみる

1.5.1 複利を計算する

元金100万円を年利1%で銀行に預けると, 1年後には元金が1.01倍されて101万円になります. 2年後には101万円の1.01倍で102万100円になります. 1年たつたびに,そのときの金額が1.01倍されるわけです. 一般に,元金を a 円, 利率を r (年利1%なら, r = 0.01 ), 預ける年数を n とすると, n 年後の元利合計は (1.25) 式の

f a r n = a ( 1 + r ) n

となります.左辺の

f a r n

は,元利合計が arn の関数であるということ,すなわち, 元金(a),利率(r), 預入期間(n)の3つによって決まるということを表しています.

この 1.5 節では,元利合計の変化をグラフで表します. 最初の金額(100万円)から始まって,1年後, 2年後,最後は10年後までの元利合計を示します. グラフを描く準備として, 1.5.1 節では1年ごとの元利合計を並べたベクトルを作ります. そのプログラムが例題 ex1501 に示されています.

このプログラムの2行目で, N乗の演算がドットつき(.^)であることに注意してください. これは行列やベクトルにおいて要素ごとの演算を行うための演算記号でしたね(1.3.3 節を参照). Nは数値ではなく,ベクトル 0 1 10 です.N乗される対象は (1 + R) ですから,1.01 という数値です. つまり,この (1 + R).^N という式は,


  (1.01).^[0 1 2 3 4 5 6 7 8 9 10]

と書いていることになります. これがどのような計算を実行することになるのかは, (1.12) 式(テキスト p.17)で示されています. 計算を実行すると,

1.01 0 1.01 1 1.01 2 1.01 10

というベクトルが得られます.このベクトルの最初の要素は 1.01 0 = 1 であることに注意してください( 1.01 0 = 0 ではありません).

N がベクトルのとき,数値のN乗の演算でドットは省略できません. 下図に示すように,ドットなしの 2^N はエラーとなります.

power

例題 ex1501 を実行すると,元金と, 1年おきに10年後までの元利合計が表示されます. テキストで説明されているように,指数形式が用いられてることに注意してください. これは Excel でも使われている表示形式で, とても大きな数,および,とても小さな数(ゼロに近い数)を表すのに便利です. 元金は 1.0e+006 * 1.0000 と表示されています. この数は,Excel では 1.E+06 と表示されます.通常の数式で書けば, 1.0000 × 10 6 となります.これは,1.0000 の小数点を右側に6桁ずらした,あるいは, 1 のあとに 0 を6個つけた数なので,元金は100万円ですね.

元金および元利合計の, 10(円)の位と1(円)の位は表示されていないということに注意してください. たとえば,元金の表示での,1.0000 の最後の 0 は,金額の100の位です. 3年後の元利合計は 1.0303 × 10 6 と表示されており,これは103万300円と読むことができます. しかし,本当の元利合計は103万301円です.

テキストにも書いてあるように,元利合計を1円の桁まで表示するには, 指数形式の表示での小数点以下を,もっと下の桁まで表示する必要があります. コマンドウィンドウでの数値の表示は,デフォルトで short という形式に設定されています. この形式では,小数点以下の表示は4桁です. これを long に変更することで,小数点以下15桁まで表示されるようになります. コマンドウィンドウで,format long とキーインします. 数値形式を short に戻すには,コマンドウィンドウで format short とキーインします. あるいは,単に format とキーインすれば, デフォルトの設定にリセットされます.

数値形式は, [設定] ダイアログボックスという GUI を使って変更することもできます. このダイアログボックスを表示するには, コマンドウィンドウで preferences と入力して [Enter] キーを押します.あるいは,下図に示すように, [ホーム] タブの [環境] セクションで, [設定] をマウスで左クリックします.

[ホーム]タブ[環境]セクション[設定]

[設定] ダイアログボックスの左側の枠(「ペイン(pane)」と呼ばれます)から, [コマンドウィンドウ] をマウスで左クリックします. 右側の [MATLAB コマンドウィンドウの設定] の上部を見ると, 数値形式が short になっていることがわかります.

[設定]ダイアログボックス[コマンドウィンドウ]

数値形式を long に変更し,[OK] ボタンを押してダイアログボックスを閉じます.

format long

[設定] ダイアログボックスでの設定変更は,MATLAB を終了させても有効なままです. この例題が終わったら,数値形式をもとの short に戻しておきましょう.

1.5.2 グラフに表示する:plot

例題 ex1501 で求めた, 時間 N (N = 0, 1, ..., 10)に対する元利合計 F をグラフに表すには, plot 関数を用います. この関数は 1.1.5 節(p.7)で一度使っています. 数値の並びであるベクトル X および Y について,


plot(X, Y)

とすると,X の各要素に対する Y の値がプロットされ,線で結ばれます.

グラフのタイトルおよび軸ラベルには日本語を使うこともできます. 例題 ex1502b でのタイトルおよび軸ラベルの設定に日本語を使ってみましょう. たとえば,

グラフタイトルおよび軸ラベルを日本語で書く

のように設定すると,次のグラフが作成されます.

タイトルおよび軸ラベルが日本語で表示されたグラフ

1.6 1次方程式を解く

1.6.1 逆行列で解く

セクション 1.6 のタイトルは「1次方程式を解く」となっていますが, 正確には(1次の)連立方程式を解きます. 行列とベクトルを用いると,連立方程式を,ひとつの1次方程式 a x = b のようにとらえて解くことができます. ここで,ab は具体的な整数であり, a は 0 ではないとします. x は未知数で, この方程式を満たす x を見つけることが, 方程式を解くということを意味します. たとえば, 2 x = 6 を解くと, x = 3 となります.

文系の人は行列とベクトルの演算になじみがないと思います. MATLAB は行列とベクトルの演算を簡単に実行できることが特徴なので, これについてある程度は理解しておくことが必要になります. 以下では,テキストでの例題 ex1601 を理解するために必要な, 最小限の説明をします.

ベクトルは数を並べたもので,これまでも使ってきました. ただし,これまでは数を横に並べた行ベクトルでしたが, 例題 ex1601 では数値を縦に並べたベクトルを使います. 数を縦に並べたベクトルを 列ベクトル(column vector)と呼びます. 行ベクトルと列ベクトルはきちんと使い分ける必要があるのですが, ここではあまり深入りしないでおきます. このセクション(1.6.1)の最後に, 行列とベクトルの掛け算について説明します. そこでは,行ベクトルと列ベクトルを使い分ける理由が少しわかると思います.

1.4.3 節ですでに説明したように, 行列は,数を縦および横の両方向に並べたものです. 並べるものは数とは限らないのですが,当面はこの理解で問題ありません. MATLAB で行列

A = 1 2 3 4 5 6 7 8 9

を作るには,以下のコードのように, 行ベクトルを作る途中で改行を指示します. 改行の指示はセミコロン(;)です.


A = [1 2 3; 4 5 6; 7 8 9]

行列において,横方向の要素の並びを「行」, 縦方向の要素の並びを「列」と呼びます. 上の行列は3行3列の行列です. 1行目は [1 2 3],2行目は [4 5 6],3行目は [7 8 9] です.

1次の連立方程式(テキスト 1.26 式)

x 1 + 2 x 2 + 3 x 3 = - 1 2 x 1 + 3 x 2 + x 3 = 0 3 x 1 + x 2 + 2 x 3 = 1

を,行列とベクトルを使って表すと,

1 2 3 2 3 1 3 1 2 x 1 x 2 x 3 = - 1 0 1

となります.この表現のルールは後(セクション1.6.1 の説明の最後)で説明します. まずは,これを1次方程式のようにとらえて解いてみます.

テキスト (1.28) 式の文字

A = 1 2 3 2 3 1 3 1 2
b = - 1 0 1
x = x 1 x 2 x 3

を用いると,行列とベクトルを使った上の連立方程式は,

A x = b

という形をしていることがわかります(テキスト 1.27 式). 一般に,数学のテキストでは, 行列はアルファベットの大文字で, ベクトルは小文字(多くの数学書では小文字の太文字)で表します. この式は,見慣れた1次方程式 a x = b と同じ形ですね.実際,行列とベクトルで表された連立方程式は, 1次方程式のように解くことができます.

通常の1次方程式 a x = b を解く方法のひとつは, 両辺に a の逆数を掛け算することです (ab は整数で, a 0 であるとします).すなわち,

1 a × a x = 1 a × b

より,方程式の解

x = b a

が得られます.両辺を a で割ってもいいのですが, 行列とベクトルの「1次方程式」 A x = b では,実は両辺を A で割るということができないのです. しかし,両辺に A の「逆数」を掛け算することはできます.

行列 A は数ではありませんから, その「逆数」もやはり数ではありません. この,行列での逆数を逆行列(inverse matrix)と呼び, 記号 A - 1 で表します.読み方は,「A インバース」です.

整数 a にその逆数 1 a を掛けると 1 になります. 行列 A に逆行列 A - 1 を掛けると,行列の「1」になります. この,行列の 1 を,単位行列(unit matrix)と呼び, 記号 I あるいは E で表します.

A A - 1 = A - 1 A = I

となります.行列 A が3行3列のとき, その逆行列 A - 1 との掛け算によって得られる単位行列 I は,

I = 1 0 0 0 1 0 0 0 1

という行列です.

整数 ab については交換法則が成り立ちます(実数でも同じ).すなわち, a b = b a です.しかし, 行列 AB については,一般に交換法則は成り立ちません.すなわち, A B B A です.しかし. 行列 A とその逆行列 A - 1 については,掛け算の順序に関係なく, 結果は単位行列となります.

行列とベクトルで表現した連立方程式(見かけは1次方程式)を解く準備ができました. 方程式

A x = b

の両辺に,左側から A の逆行列 A - 1 を掛け算してやると,

A - 1 A x = A - 1 b

となります.これは,通常の1次方程式 a x = b では,

1 a × a x = 1 a × b

という操作を行ったことに相当します. 通常の1次方程式では,逆数を右側から掛けて

a x × 1 a = b × 1 a

としてもいいのですが, 行列の方程式ではこれができません(説明は省略します). そこでわざわざ,逆行列を「左から掛ける」と述べています.

計算を続けます.等式,

A - 1 A x = A - 1 b

の左辺 A - 1 A x では, A - 1 A から計算することができます.これは単位行列 I となるのでした.単位行列は数字の 1 に相当する行列ですから, I x = x です.よって,われわれが解いている行列の方程式は,

x = A - 1 b

となります.これで,文字の式としては最終的な結果が得られました. あとは,この等式の右辺を具体的に計算してやれば, もとの連立方程式の解が得られます.

ここからは MATLAB の助けを借りましょう. 逆行列 A - 1 を計算するには,inv 関数を使います. 以下のコードをキーインしてみてください. 行列 A を作って, その逆行列を計算するコードです.


A = [1 2 3; 2 3 1; 3 1 2]
inv(A)

下図のように,行列 A の逆行列が得られます.

A inverse

得られた逆行列をベクトル b に左から掛けてやれば, もとの連立方程式の解が得られます. ベクトル b は3つの数を縦に並べたベクトルです. 数を横に並べた行ベクトルをこれまでのように作っておいて, 最後にアポストロフィ(')をつけると, 数を縦に並べた列ベクトルを作ることができます. 次のコードをキーインして実行してください. 1行目の最後にあるアポストロフィを忘れないでください.


b = [-1 0 1]'
x = inv(A) * b

すると,次のように,連立方程式の解が得られます.

answer

方程式の解は小数を使って表されていますね. MATLAB では分数の表示はされません. この解を分数で表現すると次のようになります.

x = 2 3 - 1 3 - 1 3

つまり,連立方程式

x 1 + 2 x 2 + 3 x 3 = - 1 2 x 1 + 3 x 2 + x 3 = 0 3 x 1 + x 2 + 2 x 3 = 1

の解は,

x 1 x 2 x 3 = 2 3 - 1 3 - 1 3

となります.得られた解をもとの連立方程式に代入して, 解が正しく得られていることを確認してください. たとえば,最初の等式の左辺は,

x 1 + 2 x 2 + 3 x 3 = 2 3 + 2 × ( - 1 3 ) + 3 × ( - 1 3 ) = - 1

となって,等式の右辺と一致します.

ここまでの説明で,例題 ex1601 のプログラムは理解できるでしょう. 行列だけでなく,ベクトルも大文字で表していることに注意してください.

このプログラムの3行目は,得られた解が正しいかどうかの確認です. 得られた解

X = x 1 x 2 x 3 = 2 3 - 1 3 - 1 3

が正しいならば, A X - B の値は「0」になるはずです.ただし,これは数の 0 ではなく, 0 を縦に3つ並べたベクトルです.すなわち,

A X - B = 0 0 0

となることを確認します. 連立方程式は3つの等式から構成されていますから, それぞれの等式において,左辺から右辺を引くと 0 になることを確認するのです.

プログラムを実行してみると, テキスト p.25 に示されたような結果が出力されます. 得られた解が正しいかどうか確認した結果 (誤差 E)が,期待と違うように見えます. 3つの 0 が並んでいません. 下図は MATLAB R2017a でこのプログラムを実行したときの出力です. テキストの出力と数値が異なりますが, 3つの 0 が並んでいないということは同じです.

output of ex1601

誤差 E の出力では, 指数形式が用いられていることに注意してください.上の出力図での


1.0e-015 *

という表記は,その下に並んだ3つの数(0, 0.0555, -0.1110)に, 10のマイナス15乗(分数で表すと,分子が 1,分母が10の15乗) を掛け算するということを意味しています. すなわち,誤差 E は,

E = 0 × 10 - 15 0.0555 × 10 - 15 - 0.1110 × 10 - 15

と出力されているのです.最初の数はもちろん 0 ですね. 2番目と3番目の数は 0 ではありませんが, 非常に小さな数値です. 実際のところ,これは 0 なのです. MATLAB に限らず,コンピュータ・プログラムを使った計算では, 本来はぴったり 0 のはずの数値が, このように非常に小さい(0 に近い)数値として出力されることがあります.

このセクション( 1.6.1)の最後に, 1次の連立方程式(テキスト 1.26 式)

x 1 + 2 x 2 + 3 x 3 = - 1 2 x 1 + 3 x 2 + x 3 = 0 3 x 1 + x 2 + 2 x 3 = 1

を,行列とベクトルを用いて表現すると,

1 2 3 2 3 1 3 1 2 x 1 x 2 x 3 = - 1 0 1

となる理由を説明します. 右辺のベクトルが元の連立方程式の右辺と対応していることはすぐにわかります. あとは, 左辺での行列とベクトルの掛け算がどのように行われるかを説明すればいいですね.

実数を要素として持つ,任意の3行3列の行列を,

a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33

と表すことにします.2桁の添え字の左側は行番号, 右側は列番号を表しています.たとえば, a 12 は行列の第1行第2列にある要素です.

一般に,3行3列の行列を左側から掛け算することのできるベクトルは, 要素数が3つの列ベクトルだけです. 行列と同じようにサイズを表現すれば,3行1列のベクトルです. 実数を要素として持つ,任意の3行1列の列ベクトルを,

x y z

と表すことにします.

3行3列の行列と,3行1列の列ベクトルの掛け算は, 次のように定義されます.

a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 x y z = a 11 x + a 12 y + a 13 z a 21 x + a 22 y + a 23 z a 31 x + a 32 y + a 33 z

演算の結果として得られる列ベクトルの, 最初の要素(第1行第1列の要素)に注目してください. この要素は,行列の第1行と, ベクトルの第1列の要素の積和になっています. すなわち, a 11 x a 12 y a 13 z の積(掛け算)を作り, 次にそれらすべての和(足し算)をとっています. もちろん,掛け算の結果が負の数になれば, 最終的に得られる式では引き算が現れることになります。

演算の結果として得られる列ベクトルの, 次の要素(第2行第1列の要素)に注目してください. この要素は,行列の第2行と, ベクトルの第1列の要素の積和になっています.

演算の結果として得られる列ベクトルの, 最後の要素(第3行第1列の要素)に注目してください. この要素は,行列の第3行と, ベクトルの第1列の要素の積和になっています.

一般に,nn 列の行列と, 列ベクトルの積(掛け算)は,

と定義されます.

行列とベクトルの積についての以上の説明で, 1次の連立方程式(テキスト 1.26 式)

x 1 + 2 x 2 + 3 x 3 = - 1 2 x 1 + 3 x 2 + x 3 = 0 3 x 1 + x 2 + 2 x 3 = 1

を,行列とベクトルを用いて表現すると,

1 2 3 2 3 1 3 1 2 x 1 x 2 x 3 = - 1 0 1

となることがわかりましたね,左辺の行列とベクトルの積を実行すれば, 元の連立方程式の左辺が現れます. この行列の要素はすべて正の数ですが, もとの方程式の中に引き算があれば, 行列の要素には負の数が現れることになります.

テキストの (1.30) 式では,行列と行列の積の定義が書かれています. これについての説明はテキスト第2章(2.3.2 節)で行うことにします. 行列と行列の積は,その特別な場合として, ここまでに説明してきた行列とベクトルの積を含みます. そのことだけ述べておきます. 3行3列の行列と,3行1列のベクトルの掛け算の定義は, (1.30) 式の記号を使って次のように書くことができます.

a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 b 11 b 21 b 31 = i = 1 3 a 1 i b i 1 i = 1 3 a 2 i b i 1 i = 1 3 a 3 i b i 1

シグマ記号()を使った式の意味は知っていると思います. 上の式では,記号 i を 1 から 3 まで変化させてできる, すべての項の和を計算します.たとえば,

i = 1 3 a 1 i b i 1 = a 11 b 11 + a 12 b 21 + a 13 b 31

です.

テキストの (1.31) 式および (1.32) 式で説明されているドット(.)つきの演算は, 1.3.3 節で説明したので理解できるでしょう. テキスト第2章(2.3.3 節)でも説明しています.

1.6.2 演算子で解く

前のセクション(1.6.1)で行列とベクトルを用いた等式の変形を説明したとき, 行列で両辺を割ることはできないと述べました. そのため,整数での逆数に相当する逆行列を両辺に(左から)掛け算しました.

ところが,MATLAB では行列での割り算を行う式を書くことができてしまいます. もちろん,数学で行列の割り算というものは定義されていないのですから, MATLAB で何か特別な「割り算」が可能なわけではありません. 割り算の式を書くと,逆行列を掛け算するのとは異なった数学的な方法で, 等式の解を求めてくれます.

MATLAB での行列の「割り算」を行う演算子は, 行列の左除算演算子(matrix left division operator)と呼ばれる \ と, 行列の右除算演算子(matrix right division operator)と呼ばれる / です. 日本語環境では,左除算演算子はバックスラッシュでなく, 「円マーク」であることに注意してください.

前のセクションでの計算では,方程式 A x = b を解くために,両辺において, 行列 A の逆行列を左から掛けました. コードでは,


X = inv(A) * B

としました(例題 ex1601).

これで得られたのと同じ結果が,例題 ex1602 のプログラムで得られます. すなわち,上のコードの代わりに,


X = A \ B

と書きます.ここで,A \ という表記が, A の逆行列を左から掛けるという意味だと考えるとわかりやすいでしょう.

少し注意しておきたいのは,X = inv(A) * B という式と, X = A \ B という式では,MATLAB で実行される計算が異なるということです. 後者の式では,A の逆行列を計算することなく,方程式の解 X を求めます. 計算方法の詳細についての説明は省略します(線形代数の知識があるなら, MATLAB の inv 関数のリファレンスページを参照してください). MATLAB では,逆行列を使わない, 除算演算子(\)を用いた式を用いることが推奨されています.

行列の右除算演算子(/)については, / A という表記が, A の逆行列を右から掛けるという意味だと考えるとわかりやすいでしょう.

1.6.3 危ない方程式を解くと

逆行列を求めることができるのは,行数と列数が同じ行列だけです. 行数と列数が同じ行列のことを 正方行列(square matrix)と呼びます.

テキストで注意されているように, 正方行列がいつでも逆行列を持つとは限りません. 逆行列を持たない正方行列があります.たとえば,

1 2 3 4 5 6 7 8 9

という正方行列は逆行列を持ちません. (なぜ逆行列を持たないことがわかるのかについての説明は省略します.) 逆行列を持つ行列を正則行列(non-singular matrix あるいは regular matrix), 逆行列を持たない行列を特異行列(singular matrix) と呼びます.

行列 A が特異行列である(逆行列を持たない)とき, 行列とベクトルで表現した連立方程式 A x = b は,解をひとつも持たないか,あるいは無数の解をもちます. たとえば,連立方程式

1 2 3 4 5 6 7 8 9 x 1 x 2 x 3 = - 1 0 1

の解は無数にあり, t を任意の実数として,

x 1 x 2 x 3 = t + 5 3 - 2 t - 4 3 t

と表すことができます.

行列で表現された連立方程式 A x = b において,A が特異行列であるなら, A の逆行列を両辺に左から掛け算するという解法は実行できないことになります.

このような場合に,例題 ex1601 のように逆行列を使って連立方程式を解こうとすると, MATLAB は警告のメッセージを返します. たとえば,

x 1 + 2 x 2 = 3 2 x 1 + 4 x 2 = 6

という方程式(解は無数にあります)を MATLAB R2017a で解こうとすると, 次のような警告メッセージが出されて, 解は求められませんでした.

warning about singular matrix

行列が特異であるときに,警告メッセージとともに, 方程式の解が出力されることもあります. しかし,この解は信頼できません. たとえば,連立方程式

1 2 3 4 5 6 7 8 9 x 1 x 2 x 3 = - 1 0 1

を,逆行列を使って解こうとすると, 次のような結果が返されました(MATLAB R2017a での出力). 出力された解は誤りです.

warning about singular matrix

この問題を,セクション 1.6.2 で学習した除算演算子を用いて解くと, 警告は出ましたが,無数にある正しい解のひとつが得られました.

warning about singular matrix

出力された解が正しいか誤っているかという結果よりも, 計算に問題があると警告してくれていることが重要です. このような警告が出される場合には,何らかの対処をすることになります.

行列が特異行列でなくても, 特異行列とほとんど同じ行列を用いた「危ない方程式」では, 得られる解はやはり信頼できないものとなります. 例題 ex1603 では,

A = 1 2 3 4 5 6 7 8 8.99999999999999

という行列を使っています.上で例示した特異行列と比べると, 第3行3列の要素がほんの少し異なっている(9 のかわりに 8.99999999999999 となっている)だけですね.

こうした,特異行列に非常に近い行列を用いた方程式を解く場合にも, MATLAB は警告メッセージを返してくれます.例題 ex1603 をやってみてください.

例題 ex1603 のプログラム3行目で用いている det 関数は,与えられた行列の 行列式(determinant)を計算する関数です. 行列式は行列から計算される数値で, 特異行列では 0 となります. 行列式が 0 ならば,その行列は逆行列を持たないということですね. 行列式についての詳しい説明は省略します.

例題 ex1603 を実行すると,出力される警告の最後に, RCOND という数値があります. これについての簡単な説明はテキストにあります. ただし,(1.33) 式と (1.34) 式には誤りがあります(改訂版第1刷).

(1.33) 式として提示されている

A × A - 1

は,行列の 条件数(condition number)と呼ばれます. 条件数は 1 以上の実数値をとります. 連立方程式の表現に使われる行列の条件数が大きいほど, 「危ない方程式」となります. 例題 ex1603 で返されている RCOND の値は, 行列 A の条件数ではなく, 条件数の逆数です.したがって,(1.33) 式は, 正しくは

1 A × A - 1

となります.条件数の逆数ですから, この値が 0 に近いほど「危ない方程式」となります.最大値は 1 です. 例題 ex1603 では RCOND = 7.709882 × 10 - 17 ですから,非常に小さい(0 に近い)値です.

条件数の定義に現れる A は,行列の 1ノルム(1-norm)と呼ばれるものです. その定義は 1.34 式ですが,誤植があります. シグマ記号の下側が j = 1 となっていますが,正しくは i = 1 です.文系の人にとっては,(1.34) 式は恐ろしい形かもしれませんが, 意味は簡単です. 行列を構成する要素それぞれの絶対値をとって, 縦方向に足し算します.その最大値が1ノルムです. たとえば,

A = 1 2 3 4 5 6 7 8 9

という行列の場合,第1列の合計は12,第2列は15,第3列は18ですから, A = 18 となります.

1.6.4 解がない方程式を解く:最小2乗法

このセクションでは,変数 xy の間に (1.35) 式の関係 y = a x + b が成り立っていると仮定し, データからこの直線を推定します. 統計学で 回帰分析(regression analysis) と呼ばれるデータ分析の手法です. データから推定される,2変数の関係を表す直線のことを, 回帰直線(regression line)と呼びます. 例題 ex1604a で扱うデータに対して回帰分析を実行し, データと回帰直線を図示したのがテキスト p.30 の図です. このセクションの最初に,回帰分析について簡単に説明しておきます. すでに回帰分析についての知識がある人は, 例題 ex1604a あるいは ex1604b の説明に進んでください.

データは,n 個の測定対象から, 変数 xy それぞれについて得ます. 第1の測定対象からは測定値の組 x 1 y 1 が,第2の測定対象からは測定値の組 x 2 y 2 が得られます. データの例として, n 人の学生から, 学習時間(x) とテスト成績(y)のデータが得られたという状況が考えられます.

得られた n 組の測定値は, 平面上に n 個の点としてプロットすることができます. データを表す点をプロットした図を 散布図(scattergram)と呼びます. テキスト p.30 の図は散布図に回帰直線を書き込んだものです.

データを,n 個の点としてではなく, ベクトルとしてとらえることもできます, 変数 xy について, n 個の測定対象からデータが得られているとすると, データは,

X = x 1 x 2 x n
Y = y 1 y 2 y n

と表現することができます. データを n 個の点の集合としてとらえることも, ベクトルとしてとらえることも, いずれも重要な見方です. MATLAB ではデータをベクトルとしてとらえることの方が多いと思いますが, 柔軟に解釈を切り替えてください.

データを散布図で表し.2変数の関係を表す直線を推定することを考えます. 測定対象がひとつだけのときには, 2変数間の関係を表す直線を推定することはできません. 平面においてデータ x 1 y 1 を通る直線は無数に考えられ,ひとつに決めることができないからです.

測定対象が2つの場合には,直線は簡単に決まります. 平面上で x 1 y 1 x 2 y 2 の2点を通る直線は1本しかないからです. 直線の式を具体的に求めるには,連立方程式

y 1 = a x 1 + b y 2 = a x 2 + b

を解きます.この連立方程式で, 未知数は直線の傾き a と切片 b であることに注意してください.

測定対象が3つ以上になると,テキスト p.28 でも解説されているように, n 組の測定値すべてにあてはまる直線は, 一般には存在しません. 2変数 xy の散布図において3点以上がプロットされると, これらすべてを通る直線を引くことはできないからです (すべての点がたまたま1本の直線上にある場合は例外です). そこで,データを表す点の全体にわたって, 直線と点との距離がなるべく小さくなるようにします.

ここで,このような直線が引けたと想定しましょう. まだ引くことのできていない直線を引けたと想定するのは奇妙に思えますが, これで数学的な議論をうまく進めることができます.

変数 x についての測定値

X = x 1 x 2 x n

に対応する直線上の値を, 推定値であることを意味するハット(^)を使って,

Y ^ = y ^ 1 y ^ 2 y ^ n

と表すことにします.すなわち,

y ^ 1 y ^ 2 y ^ n = a x 1 x 2 x n + b b b = a x 1 + b a x 2 + b a x n + b

という関係です. ベクトルの足し算やベクトルの定数倍については説明していませんが, 上の式で演算規則は理解できると思います.ここで, y ^ i は直線上の値なので,下図に示すように実際のデータ y i とは異なります.直線上の値と実際のデータとの差を e i と表すことにします.

regression analysis

データに最もよくあてはまる直線は, データ全体にわたって誤差 e i がなるべく小さくなるような直線であると考えられます. 誤差には正の値も負の値もあるので, すべて2乗して合計し,この合計の値が最小になるように, 直線の傾き a と切片 b を決めます. この方法を最小2乗法(least squares method) と呼びます.誤差の2乗和を最小にするという意味です. 数式で書くと,

i = 1 n e i 2 = i = 1 n ( y ^ i - y i ) 2 = i = 1 n { ( a x i + b ) - y i } 2

を最小にするよう, 直線の傾き a と切片 b を決めます. この式はテキストの (1.38) 式です. テキストの (1.38) 式には 1/2 という数字がついていますが, これは後で行う計算の都合です.

(1.38) 式を最小にする傾き a と切片 b を求めるには, 偏微分という数学的手法を用います. これはそれほど難しくありませんが,ここでは説明を省略します. 結果として,(1.39) および (1.40) 式から構成される連立方程式を解くと, ab を求めることができます.

以上で回帰分析の簡単な説明を終え, 続いて例題 ex1604a について説明します.この例題はスキップして, 例題 ex1406b に進んでもかまいません.

例題 ex1604a のプログラム1行目で乱数発生のための設定をしていますが, この設定の書き方は現在の MATLAB では推奨されていません. データ X および Y はテキストでの値と異なったものになりますが, プログラムの1行目を


rng(1234);

としてください.rng 関数は乱数発生器に (たね,seed)を設定します. この関数の名前は random number generation から来ています. 種には負でない整数を用います.この例題では 1234 です. MATLAB は,こうした種をもとに, ある決められたアルゴリズム(乱数発生器)を用いて乱数を生成します. 種が同じであれば,同じ乱数が生成されます. したがって,この例題の実行結果はいつも同じになります.

MATLAB は,起動されるたびに, 乱数発生器を同じ状態にリセットします. 種も常に同じ値です.このため, MATLAB では起動後に生成される乱数はいつも同じになります. それでは困る場合には,乱数を使う前に,


rng('shuffle');

というコマンドを実行します.

プログラムの2行目では,回帰分析に用いるデータを人工的に作っています. 変数 X に代入されているのは,

X = 1 2 3 4 5

という列ベクトルです. 数を横に並べた行ベクトルを作っておいて, 最後にアポストロフィ(')をつけると, 数を縦に並べた列ベクトルを作ることができるのでした.

変数 Y は,変数 X の値に乱数を加えて作っています. テキストで解説されているように(p.29), randn 関数は, 平均 0,分散 1 の正規分布に従う乱数(正規乱数)を発生させます. 変数 Y の値を変数 X と同じにすると, 5組のデータは (1, 1), (2, 2), (3, 3), (4, 4), (5,5) となりますから,これらすべての点は直線 y = x の上にあります.これで回帰分析を行ってもつまらないので, 乱数を加えて直線 y = x の上下に点を散らしています. 作られたデータをプロットしてみましょう. プログラムの2行目までを書いて実行した後(もっと先まで, あるいはプログラム全体を実行した後でもいいです),


plot(X, Y, 'x')

というコマンドをキーインしてください. すると,以下のようなグラフが現れます.

ex1604a Data Plot

回帰分析では, 誤差の2乗和を最小にするという意味で, これら5つの点にもっともよくあてはまる直線の式を求めます.

ここでデータのプロットに用いた plot 関数は,1.5.2 節で学習しました. 数値の並びであるベクトル X および Y について,


plot(X, Y)

とすると,X の各要素に対する Y の値がプロットされ, 直線で結ばれるのでした.

上図および例題 ex1604a では,データ X および Y は プロットだけを行い, プロットされた点を線で結ばないようにしています. plot 関数への入力で,X と Y に続いて,'x' という入力があります. これは点をプロットするマーカーをクロス印にする指定です. 詳しくは第5章で学習します.

プログラムの3行目では,(1.41) 式で表されている, 2行2列の行列 Xm を構成しています. テキストでの印刷にミスがあります(改訂版第1刷). sum(X.^2) と sum(X) の間に半角のスペースを入れてください. これまでは数値を直接に入力して行列を構成していましたが, ここでは sum(X.^2) などの式が用いられていますね. もちろん,これらの式はいずれも数値を返しますから, 結果としては数値を並べた行列

X m = i = 1 5 x i 2 i = 1 5 x i i = 1 5 x i i = 1 5 1 = 55 15 15 5

が得られます.

プログラムの4行目では,(1.41) 式で表されている, 2行1列の列ベクトル Ym を構成しています. この列ベクトルは,行ベクトルにアポストロフィをつける方法ではなく, 2行1列の行列として構成されています. 要素は Xm と同様に式で表現されています, 第1行の要素である sum(X .* Y) は,X と Y の要素ごとに掛け算を行った後, その結果をすべて合計することを意味しています. 結果として得られるベクトルは,

Y m = i = 1 5 x i y i i = 1 5 y i = 65.8243 16.8823

となります(MATLAB R2017a での実行結果).

例題 ex1604a のプログラム中で作られているベクトルや行列がよくわからなければ, その内容を表示させてください.プログラム中のセミコロン(;)を外すか, プログラムの実行後にベクトルあるいは行列を表す変数名をキーインすればいいですね. プログラムの2行目から4行目で作られている Y, Xm, Ym の内容は, 以下のようになります(MATLAB R2017a での実行結果).

Y, Xm and Ym in ex1604a

プログラムの5行目では, (1.39) 式と (1.40) 式から構成される連立方程式を解いて, 直線の傾き a と切片 b を求めています. この連立方程式を行列とベクトルで表現すると,

X m = i = 1 5 x i 2 i = 1 5 x i i = 1 5 x i i = 1 5 1 = 55 15 15 5
Y m = i = 1 5 x i y i i = 1 5 y i = 65.8243 16.8823

として,

X m a b = Y m

となります.この方程式を,1.6.2 節で学習した除算演算子で解いています. コード中の AB は単一の変数名です. A と B の掛け算という意味ではないので注意してください.

プログラムの6行目以降では,データと回帰直線の図示を行っています. 6行目では.X の値

X = 1 2 3 4 5

に対する回帰直線上の値を求めて, それらを変数 Y2 に代入しています. コード中の XX は単一の変数名で,

X X = 1 1 2 1 3 1 4 1 5 1

という行列です. X と X の掛け算という意味ではないので注意してください. ここで用いている ones 関数は,1.4.3 節で学習しました. 行数と列数を指定して(ここでは5行1列), すべての要素が 1 のベクトルあるいは行列を作ります. 行列 XX は5行1列の列ベクトルを2つ並べて構成しています. 列ベクトルを横に並べるか,行ベクトルを縦に並べると, 行列を構成することができます.

直線の傾き a と切片 b はすでに求められていますから,

Y 2 = 1 1 2 1 3 1 4 1 5 1 a b = 1 a + b 2 a + b 3 a + b 4 a + b 5 a + b

という計算によって,Y2 が具体的に求められます. 正方行列でない行列とベクトルとの掛け算については説明していませんが, ここでの計算規則はわかると思います. Y2 を求めるこの計算は (1.37) 式の左辺

x 1 1 x 2 1 x 3 1 x 4 1 x 5 1 a b

です.変数 X の値に対する回帰直線上の値を得ます. これはもちろん,実際のデータ Y とは一致しません.

変数 X と Y2 の値の組をプロットして線で結べば, 回帰直線を描くことができます.

データ (X, Y) と回帰直線 (X, Y2) を同一の平面に図示するには, plot 関数の入力にこれらの組を順に並べます. データは直線で結ばずに,クロス印を使ってプロットだけを行います. これを行っているのがプログラムの7行目


plot(X, Y, 'x', X, Y2)

です.plot 関数の使い方は,第5章で詳しく学びます.

プログラム8行目の title 関数は, 現在表示されているグラフにタイトルをつける関数です. 9行目の xlabel 関数および ylabel 関数はそれぞれ, 現在の座標軸の x 軸および y 軸にラベルをつける関数です.

例題 ex1604a の解説は以上です.

データから回帰直線(の傾きと切片)を求めるには, (1.39) および (1.40) 式から構成される連立方程式を解くことが正攻法です. 例題 ex1604a はまさにこの方法で回帰直線を求めました.

しかし MATLAB では,データがあたかもすべて直線上にあるかのような方程式を書き, 1.6.2 節で学習した除算演算子を用いて方程式を解くと, 回帰直線を得ることができます. この方法で回帰分析を実行しているのが例題 ex1604b です.

データ x i y i を表す n 個の点が, 散布図においてすべて直線上にあるとします. 行列とベクトルを用いて表すと.

y 1 y 2 y n = x 1 1 x 2 1 x n 1 a b = a x 1 + b a x 2 + b a x n + b

となります. もちろん,現実にはこのような関係は成立していません. ところが,行列の除算演算子を用いてこの方程式を解くと, 回帰直線の傾きと切片が得られるのです.

例題 ex1604b のプログラム2行目までは,例題 ex1604a と同じです. データ

X = x 1 x 2 x 3 x 4 x 5 = 1 2 3 4 5

に乱数を加えてデータ

Y = y 1 y 2 y 3 y 4 y 5

を作っています. 変数 XX には行列が代入されており,

X X = x 1 1 x 2 1 x 3 1 x 4 1 x 5 1 = 1 1 2 1 3 1 4 1 5 1

です.プログラムの3行目で.方程式

y 1 y 2 y n = x 1 1 x 2 1 x n 1 a b

を,除算演算子を用いて解いています. 本来は成立していない方程式ですから強引なようですが, これで回帰直線の傾きと切片を求めることができます.

プログラムの4行目以降は,例題 ex1604a でも行った, データと回帰直線の図示です.

1.7 グラフィカルユーザインタフェイス

このセクションで書くプログラム(例題ex1701)は, これまでに書いたプログラムに比べれば比較的長いです. こういった長いプログラムは, テキストに書かれたプログラムを全部入力してから実行するのではなく, 部分的に動作をテストすることを繰り返して完成させます. 長いプログラムを入力してから実行すると, エラーが出たときに原因を突き止めて対処することが難しくなります.

このセクションでは,すべての処理をひとつのプログラムに書くのではなく, 処理の一部を切り離して別のプログラムとする方法も経験します, 切り離された別プログラムは, サブルーチン(subroutine), あるいは,関数(function)と呼ばれます. 関数については第4章で学習します.

1.7.1 ボタンを付ける:uicontrol

例題 ex1701 のプログラムを書く前に, 自作の関数 ex1701sub を完成させます. この関数がきちんと動作することを確認してから, 例題 ex1701 のプログラムに取り掛かってください. 自作の関数は,MATLAB に初めから組み込まれている関数と同じように, 他のプログラムから呼び出して利用することができます.

関数を作成するには,プログラムを書く時と同様に, MATLAB のエディターを使用します. これまではエディターで「スクリプト」を書いてきました. ここで書くのは「関数」です. 関数もスクリプトの一形式だと考えることもできますが, MATLAB では区別されています. MATLAB の [ホーム] タブの [新規作成]から, [関数] をマウスで左クリックします.

Create Function (menu)

すると,関数を定義するためのフォーマットでエディターが起動します. これまでスクリプトを書くためにエディターを起動したときは, 何も書かれていない空白のフォーマットでした.

Create Function (editor)

関数の定義は function というキーワードで始まります. テキストの関数定義 ex1701sub を見てください. 1行目は


function ex1701sub(action)

となっています.ここで,ex1701sub が関数名, action は関数への入力引数(input argument)です. 「引数」は「ひきすう」と読みます. 入力引数とは関数への入力のことです. 例題 ex1701 では,この関数に具体的に入力されるのは, 'helix' か,'sinc' か,'exit' の,いずれかの文字列です. 具体的な入力に複数の可能性があるので, さまざまな値をとる変数のように, 入力引数を action と表記しています. 文字列(strings)とは文字の並びです. MATLAB では文字列を引用符で囲みます. テキスト第2章の 2.5 節で文字列の扱いについて学びます.

テキスト(改訂版)の第1刷では関数名に誤植がありました. 関数名が sub10801 となっているテキストを持っている人は, ex1701sub に訂正してください. 入力引数に対して, 関数の出力を出力引数(output argument)と呼びます. この関数 ex1701sub では出力引数は扱いません.

関数を書くためのフォーマットで MATLAB のエディターを起動すると, 1行目は


function [ output_args ] = Untitled( input_args )

となっています. テキストの関数定義 ex1701sub の1行目と同じになるように修正してください. キーワード function は残し, その後の [ output_args ] = は消去します. 次の Untitled は関数名 ex1701sub に,input_args は action に書き換えます.

エディターであらかじめ書かれている


%UNTITLED この関数の概要をここに記述
%   詳細説明をここに記述

は注釈です.ここに何を書くのかについては第4章で学習します. 今はこのままにしておいてかまいません.あるいは, テキストの関数定義 ex1701sub の2行目のように書き換えてもいいです.

定義される関数 ex1701sub は, 関数を呼び出して使用するときに与えられる入力引数によって, 異なった動作をします. 入力引数が 'helix' という文字列である場合には, テキスト p.34 の図(グラフ)を描きます. 入力引数が 'sinc' という文字列である場合には, テキスト p.35 の図を描きます. 入力引数がその他の文字列(例題 ex1701 では, 'exit')である場合には, グラフ用のウィンドウ(Figure ウィンドウ と呼ばれます)が閉じられます.

このような条件分岐を可能にしているのは, 第3章 3.5 節で学習する, switch 文と呼ばれる次の構文です.


switch switch_expression
   case case_expression
      statements
   case case_expression
      statements
    ...
   otherwise
      statements
end

switch 文は「式 switch_expression の値に注目して, 特定の値 case_expression に応じた実行文 statements を実行せよ」 という意味だと考えてください. 1行目の switch_expression は,数値または文字列を与える式です. 関数 ex1701sub では,ここに入力引数 action が使われています. この action には 'helix' か,'sinc' か,'exit' の,いずれかの文字列が入るのでした. 関数の実行時にどの文字列が action として指定されるかによって, 2行目以降で指定される異なった動作をします.

テキストでは switch 文の1行目が switch(action) となっていますが, action に括弧をつけずに


switch action

としてください.switch の後に半角のスペースを入れます. この方が switch 文の構文とあっています.

2行目以降は switch_expression の値に応じた場合分けです. 関数 ex1701sub での case 'helix' は, 式 switch_expression の値が文字列 'helix' である場合を意味します. 同様に,case 'sinc' は, 式 switch_expression の値が文字列 'sinc' である場合を意味します. それぞれの場合に行う動作が, 上述の構文での statements の位置に記述されます.

式 switch_expression の値が文字列 'helix' である場合の実行文は,


t = 0:pi/50:10*pi;
plot3(sin(t), cos(t), t)

です.最初の文の意味はわかりますね,ベクトル t を作っています. ベクトル t の最初の要素は 0 です.そこから π/50 きざみで, 10 * π まで要素を並べます. 次の文では plot3 関数が使われています. これは,2次元のグラフの描画に用いた plot 関数の, 3次元バージョンです.ベクトル t の各要素について, sin(t) と cos(t) の値を計算すると,3次元の座標 sin 0 cos 0 0 sin π 50 cos π 50 π 50 sin 2 π 50 cos 2 π 50 2 π 50 , ・・・, sin 10 π cos 10 π 10 π が得られます.これらの点をプロットして線分で結ぶと, 3次元のグラフができます. もちろん,グラフは2次元のディスプレイに表示されるのですから, 実際には3次元ではありません.

式 switch_expression の値が文字列 'sinc' である場合の実行文は,


[X Y] = meshgrid(-8:0.5:8);
R = sqrt(X.^2 + Y.^2) + eps;
Z = sin(R) ./ R;
mesh(X, Y, Z)

です.最後の実行文は,テキストでは mesh(Z) となっていますが, 上のように修正してください. その方がわかりやすいですし,式とグラフが一致するので.

1行目で使われている meshgrid 関数は, 引数にベクトルをとり,それを異なった方向(縦および横)に複製して並べ, 2つの行列を作ります. 関数 ex1701sub では,生成された2つの行列を, それぞれ変数 X と Y に代入しています. ベクトル -8 -7.5 8 を入力引数として,

X = -8 -7.5 8 -8 -7.5 8 -8 -7.5 8
Y = -8 -8 -8 -7.5 -7.5 -7.5 8 8 8

となります.この meshgrid 関数については 5.4.1 節で学習します.

2行目と3行目の実行文では, 1行目の meshgrid 関数で生成された X と Y の値を使って, テキスト p.34 の (1.44) 式から Z の値を計算しています. 2行目の実行文で, 最後についている + eps を無視(これは machine epsilon と呼ばれ,扱うことのできる最小の正の数です)すると,

R = ( -8 ) 2 ( -7.5 ) 2 8 2 ( -8 ) 2 ( -7.5 ) 2 8 2 ( -8 ) 2 ( -7.5 ) 2 8 2 + ( -8 ) 2 ( -8 ) 2 ( -8 ) 2 ( -7.5 ) 2 ( -7.5 ) 2 ( -7.5 ) 2 8 2 8 2 8 2
Z = sin R R

となります. 行列に対してドット(.)つきの演算が実行されていることに注意してください. 2行目の実行文では, 行列 X および Y の各要素を2乗し,それから2つの行列を加えて, 最後に各要素の平方根をとります. この計算により得られた結果は行列であり, 変数 R に代入されています. 3行目の実行文では.行列 R の各成分に対して sin の値を計算し, 得られた値を対応する R の成分で割り算しています.

実行文の4行目では mesh 関数が使われています. この関数は,式 switch_expression の値が文字列 'helix' の場合に用いた plot3 関数と同様に, 3次元のグラフの描画に用いられます. (X, Y, Z) で与えられる点をメッシュ状につないで表示します. この関数については,meshgrid 関数とともに,5.4.1 節で解説されています.

式 switch_expression の値が文字列 'helix' である場合と, 'sinc' である場合に実行される動作を説明してきました. 例題 ex1701 から関数(サブルーチン) ex1701sub が呼び出されるとき, もうひとつの場合は文字列 'exit' です. したがって,ex1701sub の中で case exit という記述すればいいのですが, ここでは異なる方法を使うことにします. それが,関数 ex1701sub の定義(テキスト p.33)12行目と13行目の,


otherwise
  close

です. これは文字通り「その他(otherwise)」の場合に実行する動作を記述しています. どの case にも該当しない入力引数で関数 ex1701sub が呼び出されるときには, この動作が実行されます. ここでの実行文は close だけですね. これは関数です.入力引数なしでこの関数が呼ばれたときには, 現在の Figure ウィンドウが閉じられます.

それでは以下のように,入力される文字列が 'helix' である場合と, 最後の「その他」の場合だけ書いてみましょう. 書き終えたら,ex1701sub.m というファイル名で保存します. 保存先は,後で作成する例題 ex1701 のプログラムと同じディレクトリにします. 授業用ディレクトリに保存しておけばよいでしょう.

ex1701sub

テキストでは関数定義の最後がひとつの end で終わっています. これは switch 文の終了を意味します. これに対して上図では,最後に end が2つあります. 最初の end は switch 文の終わりを意味し, 次の end は関数定義の終わりを意味します. 関数定義のフォーマットでエディターを起動すると, 一組の fuction と end があらかじめ書かれています. 定義で switch 文を使えば,switch 文の終わりを示すために, もうひとつ end が書かれることになります. 関数定義の終わりを意味する end は常に必要ではなく, この関数 ex1701sub ではテキストのように省略可能です.

入力される文字列が 'helix' である場合と, 最後の「その他」の場合だけ書いたら, この関数 ex1701sub をキーインして実行してみましょう. コマンドウィンドウで,


ex1701sub('helix')

とキーインしてください. これは,文字列 'helix' を入力として, ex1701sub 関数を呼び出しています. 次のようなグラフが描かれたウィンドウが出現すれば, 'helix' の case について, switch 文は正しく動作しています.

helix

続いて,文字列 'exit' を入力として, ex1701sub 関数を呼び出してください. 入力が「その他」の文字列である場合に, 開かれていたグラフのウィンドウが閉じられることを確認します. コマンドウィンドウで,


ex1701sub('exit')

とキーインしてください. グラフのウィンドウが閉じられたでしょうか. 閉じられれば,ex1701sub は全体として正しく動作しています.

部分的に書いた関数 ex1701sub が正しく動作していることを確認したら, 関数の残りの部分を書いて完成させましょう. 関数への入力が文字列 'sinc' である case を switch 文に追加します. 関数を上書き保存して,動作をテストしてください. 入力が 'helix' および 'exit' である場合はすでにテストしましたから, 入力が 'sinc' である場合を新たにテストすればいいですね.

関数 ex1701sub の作成が終わったら, 例題 ex1701 のプログラムを始めます. コマンドウィンドウで直接にキーインするのではなく, エディターを開いてプログラムを書きます. エディターを開くには,MATLAB の [ホーム] タブにある, [新規スクリプト] をマウスで左クリックします.

例題 ex1701 で作ろうとしているのは, 「ボタンをマウスでクリックすることによって, 描画したい関数のグラフを選択するプログラム」(p.31)です. このプログラムを実行すると, テキスト p.33 に図示(ex1701-a)されている Figure ウィンドウが表示されます. ここで,[Helix] というボタンを押すとテキスト p.34 の図(ex1701-b), [Sinc] というボタンを押すとテキスト p.35 の図(ex1701-c)に示されているグラフが描画されます. [Exit] というボタンを押すと Figure ウィンドウが閉じられます.

例題 ex1701 のプログラムは,

という3つの部分から構成されています. これらについて順に説明します.

プログラムの1行目は Figure ウィンドウの描画です. コマンドウィンドウに


figure('Position', [150 150 500 400])

とキーインしてみてください. プログラムを理解するために, その一部をキーインして実行するのはよい方法です. このようにキーインすると, 下図のように何も描かれていない Figure ウィンドウが表示されます.

Figure window

ここで使われている figure 関数は Figure ウィンドウを作成する関数です. 実は,単に figure とキーインすれば, Figure ウィンドウが現れます. そうすると, 例題 ex1701 でこの関数に与えている引数は何を意味しているのでしょうか. これは,作成する Figure ウィンドウのサイズと,表示する位置を指定しているのです.

Figure ウィンドウには, 色,名前,表示位置とサイズなどの プロパティ があります. これらプロパティにはあらかじめデフォルト値が設定されています. 単に figure とキーインすると, これらデフォルト値を使って Figure ウィンドウが作成されます. あるプロパティ(たとえば,ウィンドウの位置とサイズ)に関して特定の値を使いたければ, figure 関数を使うときに,そのプロパティの名前と値を引数として与えます. 例題 ex1701 での


figure('Position', [150 150 500 400])

という関数呼び出しは, Position プロパティの値を [150 150 500 400] として, Figure ウィンドウを作成することを要求しています. プロパティ名は引用符(')で囲みます. プロパティは複数指定できます.一般的な書式は


figure('名前1', 値1, '名前2', 値2, '名前3', 値3, ...)

となります. プロパティの「値」は単一の数値であるとは限りません. 数値を集めたベクトルであったり,文字列であったりします.

Figure ウィンドウの Position プロパティの値は, 4つの数値をひとまとめにしたベクトルです. 例題 ex1701 では [150 150 500 400] ですね. この4つの数値は [left bottom width height] という意味です. 単位はピクセルです. すなわち,ウィンドウの左下の座標値が (150, 150), 幅が 500 ピクセル,高さが 400 ピクセルということです. 座標の原点は (1, 1) で, これはコンピュータのスクリーンの左下です. 原点は (0, 0) でないとおかしいような気がしますが, 1ピクセルは小さな正方形であることに注意してください. コンピュータのスクリーンは小さな正方形が敷き詰められています. それぞれに色を表示することで, スクリーン全体としての画像が出力されます. いちばん左下のピクセルの位置が (1, 1) となります.

コンピュータのスクリーンのサイズを知りたければ,


get(0, 'ScreenSize')

とキーインしてください. 下図のようにスクリーンのサイズが表示されます. これは私のノート PC のスクリーンです.

ScreenSize

表示される4つの数値は Figure ウィンドウの Position プロパティの値と同じ意味です. 幅が1440ピクセル,高さが900ピクセルであることがわかります. 左下の座標値は (1, 1) です.

MATLAB のバージョンが R2014b 以降なら, groot とキーインすると, 下図のようにスクリーンサイズの情報が得られます.

groot function

以上の説明で,例題 ex1701 のプログラムの1行目,


h0 = figure('Position', [150 150 500 400])

の右辺の意味は理解できたと思います. それでは,これを左辺の変数 h0 に代入しているのは, どういう意味なのでしょうか? テキストの p.34 に, 「ex1701 の1行目ではウィンドウを1つ生成します. このウィンドウを指し示す情報,すなわち, ハンドル を h0 に格納しています」 という説明があります.ちょっと抽象的で難しいですね. 以前のバージョンの MATLAB では, 例題 ex1701 の1行目をキーインして実行すると, 下図(MATLAB R2013b での出力です)のようにハンドル番号と呼ばれる数値が出力されました. 変数 h0 の情報を whos 関数(1.3.1 節で説明しました)で表示してやると, 数値を格納した変数と同じであることがわかります. (よくわからなければ,適当な変数に数値を代入して比べてみてください.) Figure ウィンドウには識別のための数値(ハンドル番号)が与えられ, この数値がハンドルだったのです.

figure handle MATLAB R2013b

MATLAB R2014b 以降のバージョンでは,この仕様が変更され, Figure ウィンドウのハンドルは オブジェクトハンドルと呼ばれるものになりました. 例題 ex1701 の1行目を MATLAB R2016a でキーインして実行すると, 次のような出力がされます.

figure handle MATLAB R2016a

Figure ウィンドウのプロパティが表示されるのですね. 変数 h0 の情報を whos 関数で表示してやると, 数値を格納した変数とは異なることがわかります. こうなると,ハンドルは単に数値ではなく, Figure ウィンドウを指し示す「情報」という説明(テキスト p.34)がしっくりきます.

プログラムの2行目から6行目では, 1行目ですでに描画した Figure ウィンドウの中に, グラフを表示するための座標軸を描きます. 2行目の axes 関数は, すでに描かれている Figure ウィンドウの中に座標軸を描く関数です. 1行目で用いた figure 関数と同様に, プロパティの値を指定して呼び出します. プロパティの値が指定されなかったときにはデフォルト値が用いられます. プロパティの値の指定方法は figure 関数と同じです.一般的な書式は.


axes('名前1', 値1, '名前2', 値2, '名前3', 値3, ...)

となります.例題 ex1701 では, Parent と Units と Position という3つのプロパティについて値が指定されています. プログラムの2行目から5行目のそれぞれで最後についている3つのピリオド(...)は, ひとつの実行文を複数行に分けることを意味するのでした(1.2.2節). 3つのピリオドの前に半角スペースをひとつ空けることをお勧めします.

Parent プロパティには,座標軸を描く Figure ウィンドウのハンドルを指定します. この例題の1行目で描いた Figure ウィンドウのハンドルは, 変数 h0 に格納されていたのでした. したがって,ここでの Parent プロパティの値は h0 となっています. 座標軸を描くための Figure ウィンドウを, 座標軸に対する 親ウィンドウ と呼びます. 一般に,配置したいオブジェクト(この例題では座標軸)の上位にあるオブジェクトを 親オブジェクト, 親オブジェクトの中に配置されるオブジェクト(この例題では座標軸)を 子オブジェクト と呼びます. これについては第5章で学びます.

Unit プロパティには, 座標軸の位置やサイズを指定するときの測定単位を指定します. ここでは単位をピクセル(pixels)としています.

Position プロパティは, Figure ウィンドウの中に描く座標軸の位置と大きさを指定します. このプロパティは, プログラム1行目での Figure ウィンドウの描画でも指定しました. 描く座標軸の位置と大きさは [60 100 400 280] となっています. この4つの数値は [left bottom width height] という意味です. 原点 (1, 1) は親ウィンドウの左下です. 親ウィンドウの左下を (1, 1) として, 描く座標軸の左下の位置は (60, 100), 幅が 400 ピクセル,高さが 280 ピクセルということです.

Figure ウィンドウのハンドルを変数 h0 に格納したのと同様に, 座標軸のハンドルを変数 h1 に格納します(プログラムの2行目).

ここまでの説明が理解できたら, エディターで例題 ex1701 のプログラムを6行目まで書き, ex1701.m という名前で保存してください. この例題のように長いプログラムでは, 全体を一度に書いて実行してはいけません. 部分的に書き,テストを繰り返しながら完成させます

6行目まで書いた例題 ex1701 のプログラム ex1701.m を実行しましょう. コマンドウィンドウで,ex1701 とキーインします. すると,下のような Figure ウィンドウが表示されます. 1行目だけを実行した場合と比べてみてください. 座標軸が追加されていますね.

Figure window plus axes

プログラムの7行目から27行目では,Figure ウィンドウの内部, 座標軸の下に,3つのボタンを描画しています. 7行目から13行目が [Helix], 14行目から20行目が [sinc], 21行目から27行目が [Exit] というボタンの描画です.

長いプログラムは部分的に書き, テストを繰り返しながら完成させるという原則に従って, [Helix] というボタンだけを描画してみましょう. ボタンは,Figure ウィンドウの中に描かれるという点では, すでに描いた座標軸と同じです. これらを比較すると, プロパティの指定方法が類似していることに気がつくでしょう.

座標軸を描くための関数が axes だったのに対して, ボタンを描くための関数は uicontrol です. テキスト p.35 で説明されているように, これはボタンなどの GUI オブジェクトを作るための組み込み関数です. この例題で作るプッシュボタンの他に,ラジオボタン, トグルボタン,チェックボックス,リストボックス, などを作ることができます. MATLAB では,これらを User Interface Control オブジェクトと呼びます. 関数名 uicontrol はここから来ています. GUI オブジェクトについては第7章で学習します.

作成する GUI オブジェクトの種類は Style プロパティで指定します. プッシュボタンを作るには, このプロパティの値を 'Pushbutton' とします(プログラムの9行目).

ボタンの中に表示する文字列は String プロパティで指定します. ここでは 'Helix' です(プログラムの11行目).

Parent プロパティの値が h0 であるのは, 座標軸を描画したときと同じです. 座標軸もボタンも,同じ親ウィンドウの中に描画します. 親ウィンドウのハンドルが変数 h0 に格納されています.

Position プロパティは,Figure ウィンドウの中での, [Helix] ボタンの位置とサイズを指定しています. これは [60 20 100 40] となっています. すでに描いた座標軸の左下の位置は (60, 100) でしたね. [Helix] ボタンの左下は (60, 20) ですから, 左端の位置は座標軸とそろえていますね. ボタンの幅は 100,高さは 40 です. ボタンの幅と高さは,後で作る2つのボタンでも同じです.

テキストでは Callback プロパティの値が [ ] で囲まれていますが, これは不要です.


'Callback', 'ex1701sub(''helix'')' ...

としてください. テキスト(改訂版)の第1刷では関数名に誤植がありました. 関数名が sub1701 となっているテキストを持っている人は, ex1701sub に訂正してください. すぐ後で説明するように, helix を囲んでいるのは二重引用符(")ではなく, 2つの単一引用符('')です.注意してください.

Callback プロパティは, ユーザーが User Interface Control オブジェクトを操作したときに実行される関数を指定します. ここでは,ボタンが押されたときに,文字列 'helix' を引数として, 前もって自作しておいた関数 ex1701sub を実行します. 関数 ex1701sub の動作をテストしていたとき,コマンドウィンドウで


ex1701sub('helix')

とキーインすると, テキスト p.34 のグラフが表示されたことを思い出してください. Callback プロパティの値として設定されているのは, このコマンドを単一引用符で囲んだ文字列ですね(プログラムの12行目). ただし,'helix' の部分が ''helix'' なっています. これは,引用符がひとつだと, 'ex1701sub(' までが Callback プロパティの値だと解釈されてしまうためです. 単一引用符で囲まれた文字列の中で, さらに単一引用符を使いたいときには, このように単一引用符を2つ重ねます

他のプロパティの値がボタンの描画においてただちに反映されるのに対して, このプロパティで指定された関数は,ボタンの描画の段階では実行されません. こうした関数は,電話番号だけ伝えておいて, 後で相手から電話をかけてもらうという行為に似ているため, コールバック関数(callback function) と呼ばれます. MATLAB に ex1701sub という関数を伝えておくと, ボタンが押されたときに MATLAB が ex1701sub を呼び出してくれるということです.

ここまでの説明が理解できたら, エディターを開いて,例題 ex1701 のプログラムを13行目まで書き, ex1701.m を上書き保存してください. コマンドウィンドウで,ex1701 とキーインして実行すると, 下図のように,[Helix] というボタンが追加されたことがわかります.

Figure window plus axes plus [Helix] button

ここで [Helix] ボタンを押すと, 文字列 'helix' を引数として関数 sub1701 が呼び出され, 下図のグラフが描画されます.

Helix graph

[Helix] ボタンのハンドルは変数 h2 に格納されます.

プログラムの残りの部分は, [Sinc] ボタンと [Exit] ボタンの描画ですね. これは [Helix] ボタンの描画とほとんど同じですから, この部分のプログラムは容易に理解できるでしょう. まずは [Sinc] ボタンを描画するところまでプログラムを書いて実行してください. 動作に問題がなければ,[Exit] ボタンの描画を行い, 例題 ex1701 のプログラムを完成させてください.

1.8 MATLAB は難しい!?:関数のまとめ

1.8.1 22世紀まであと何日?:datenum

MATLAB には非常に多くの組み込み関数があります. この入門テキストだけでもかなりの数の関数を学習します. どのような関数が用意されているのか, 経験を積んで少しずつ知っていきましょう.

キーワードを使って関数を検索するには, lookfor 関数を使います. これについては第4章で学習します.

コマンドウィンドウで help とキーインすると, 「ヘルプトピック」が表示されます. テキストではこれを「関数のグループ」と書いていますが, 関数だけでなく,MATLAB で使うあらゆる文(たとえば, for ループ)についての解説にアクセスできます.

テキストでの記述と異なり,現在の MATLAB では, メニューの「ヘルプ」の中に「プロダクトのヘルプ」はありません. かわりに「ドキュメンテーション」にアクセスしてください. 下図は MATLAB R2017a でのドキュメンテーションへのアクセスです. コマンドウィンドウで doc とキーインすることによっても, ドキュメンテーションにアクセスできます.

access to documantaion

1.8.2 関数のまとめ

この章で学習した書式,演算記号,定数,関数については, ここに挙げられているリストを見て, 使い方を説明できるかチェックしましょう.

関数の使い方を調べるには help 関数を使います(p.41). 学習済みの関数について,たとえば


help('plot')

というように,使い方を調べてみるとよいでしょう.

学習の初期段階では,ヘルプを読んでも意味のわからないところばかりです. しかし, ヘルプを読むスキルは学習が進むにつれて重要になってきますから, 意味のわかるところを探して読み, 説明の書き方になじんでいってください