MATLABの基本操作からデータ処理・グラフ作成まで
2017.05.26 MATLAB TA
MATLABとは？
MATLABは、科学・技術計算に特化したプログラミング言語・環境です。世界中の多くの大学や企業で、最先端の研究・開発にMATLABが使用されています。
東工大では、計算機室等の学内の設備で利用可能なほか、皆さんの個人所有PCにもインストールが可能なTAHライセンスが締結されています。個人で購入するととても高価なソフトウェアですので、この機会にインストールして使ってみましょう。
アカウントの取得
アカウントの取得には、m.titech.ac.jpドメインのメールアドレスが必要です。アカウント取得の手順はGSICのウェブサイトに詳しく記載されています。
インストール・アクティベーション
インストーラをダウンロードします。MathWorksアカウントへのログインが必要です。
各種OS向けのインストールガイドはこちらをご覧ください：
インストールが完了したら、「MATLABのアクティベーション」にチェックを入れて完了します。アクティベーション画面が起動するので、画面の指示に従ってライセンスをアクティベートしてください。
この際、「アカウントの取得」で使用したアクティベーションキーを使用します。
画面の基本操作
MATLABを起動すると、「ファイルエクスプローラ」、「ワークスペース」、「コマンドウィンドウ」の3つが現れます。
「ファイルエクスプローラ」
Windowsのエクスプローラと同じ感覚でファイルを操作できます。
「ワークスペース」
定義した変数のリストが表示されます。変数名をダブルクリックすると「変数エディタ」が開き、変数を直接編集することもできます。
「コマンドウィンドウ」
MATLABコマンドを入力して実行することができます。MATLABのコマンドを実行する方法は、
- 「コマンドウィンドウ」から1つずつ実行する方法
- 「スクリプト」と呼ばれる、複数のコマンドをつなげて1つのファイルに記述したものを実行する方法
の、大きく2つあります（「スクリプト」については後述）。
まず関数電卓として使ってみる
MATLABでは、通常の数式と同じようにして、変数に数値を代入できます。つぎの式を「コマンドウィンドウ」に入力して、Enterを押してください。
a = 1
a に1を代入した結果が表示されましたか？うまくいかない場合は、日本語入力がONになっていないか、スペースや記号が全角になっていないかどうかを確認してください。入力は英数入力モードで半角文字を使います。
コマンドの最後にセミコロン「;」をつけると、結果を表示せずに実行させることができます。
b = 2.3;
ここまでで定義した a と、b に2をかけたものを足し合わせて、新しい変数 c に代入してみましょう。掛け算の記号にはアスタリスク「*」を、割り算にはスラッシュ「/」を使用します。
c = a + 2*b;
変数の値は、変数名をそのまま入力して Enter を押すことで表示できます。
c
MATLABの特徴は、ベクトル・行列計算が簡単に行えることです。まずはベクトルを作ってみましょう。数値をカンマ区切りで列挙して角かっこで囲むことで、行ベクトルができます。つぎの2つの書き方のうちどちらを書いても同じです（カンマの必要はありません）。
v = [1, 2, 3]
v = [1 2 3]
列ベクトルを作る際には、セミコロン「;」を区切り文字として使います。
w = [1; 2; 3]
この行ベクトルと列ベクトルの表記を組み合わせることで、行列を作ることができます。
A = [1 2; 3 4; 5 6]
A = [1 2;
3 4;
5 6]
行列やベクトルの演算も、スカラー変数と同様に可能です。
[1 2; 3 4] * [5 6; 7 8]
MATLABでは、範囲を指定したベクトルの生成をよく使用します。「:」を使って、「初項:刻み幅:末項」と書きます。
x = 1:10
y = 0:0.1:0.5
零行列、単位行列などの特殊な行列は、コマンドを使って生成することができます。
I = eye(3)
O = zeros(3)
Q = zeros(4,2)
行列やベクトルの転置には、シングルクォート「'」を使います。厳密には「'」は共役転置で、単なる転置を行う際にはドットをつけた「.'」を使います。
A % もとの行列
B = A' % 共役転置
B = A.' % 転置
ベクトルや行列の結合も簡単にできます。ブロック行列を作る感覚で以下のように記述します。
C = [B; v; 0 1 2]
D = [A w]
また、行列の要素を取り出すには 変数名(行, 列) のようにします。
D(1,2)
行と列の指定に範囲ベクトルを使うと、部分行列を取り出すことができます。全部の行や列を指定したいときは、数字を書かずに「:」だけを書きます。また、最後の行や列まで取り出したいときは、 end という特別な記号を使います。これらは行数や列数がたくさんあるデータを扱うときに便利です。
D(1:2, 2:3)
D(:,2:3)
D(1:2, 2:end)
MATLABでは、円周率 πや複素数を簡単に扱えます。円周率には pi を、虚数単位には 1i または 1j を使います。
pi
1+2i
ただし、コンピュータの数値計算の原理上、若干の誤差が生じます。たとえば以下の例を計算すると答えは0ですが、非常に小さな虚数が結果に含まれてしまっています。(気になる人は「浮動小数点演算 誤差」で検索してみましょう)
1+exp(pi*1i)
ここで説明した基本事項はMATLABの機能のほんの一部ですが、これらを組み合わせることで基礎的なデータ処理が可能です。まずは手を動かして慣れることが重要ですので、説明したことを活用していろいろ計算をしてみてください。数学でよく使う関数のMATLABでの利用方法については、こちらのページを参照してください：
スクリプトの作成
ここまではコマンドウィンドウを使って対話的に計算を実行してきました。スクリプトを使うと、一連の処理の流れをファイルにまとめておくことができます。これによって、同じ処理を複数のデータに対して何度も繰り返す際に、打ち込む手間が省け、ミスの防止にもつながります。
スクリプトを作成するには、メニューの「新規作成」→「スクリプト」を選ぶか、「ファイルエクスプローラ」で右クリックして「新規ファイル」→「スクリプト」をクリックします。
スクリプトエディタが開くので、実行したいコマンドを順番に書いていきます。今回は以下の画像のようなベクトル同士の演算を行ってみます。スクリプトを書く際は適宜空行を入れ、見やすくなるように心がけましょう。
ここで、clear というコマンドが1行目にありますが、これはワークスペースの変数（=現在使われている変数）をすべて消去するコマンドです。スクリプトの冒頭に書くと、前回の実行結果の混入を防ぐことができるので、 clear を書く癖をつけましょう。
出来上がったらファイルを保存します。今回は「myscript.m」というファイル名で保存します。MATLABのスクリプトには「.m」という拡張子を使います。
スクリプトを実行するには、画面上部にある実行ボタンを押すか、コマンドウィンドウにスクリプトの名前を打ち込んで実行します。エディタに入力できる状態で「F5」キーを押しても実行できます。
データの取り込み・入力・保存
MATLABには、さまざまな種類のファイルのデータを扱うコマンドが用意されています。以下に一例を示します。
MATLAB形式のファイルの読み込みと書き込みを行うコマンドです。「.mat」という拡張子を使います。MATLABで処理した後のデータは.matファイルに保存しておくと、データの再利用ができて便利です。
CSV形式のファイルを読み書きします。CSVとは、つぎのようなカンマ区切りのテキストで表現されるデータ形式のことです。
type data.csv
A = csvread('data.csv')
Excelファイルのデータを読み書きします。「.xls」ファイルと「.xlsx」ファイルの両方に対応しています。
A = xlsread('data.xlsx')
A = xlsread('data.xls')
また、変数エディタを使用するとExcelのような感覚で変数の値を入力することができます。「ワークスペース」にリスト表示されている変数の名前をダブルクリックすると、変数エディタが開きます。
変数エディタはExcelからの貼り付けにも対応しています。Ctrl+C, Ctrl+V 等のショートカットや、右クリックであらわれるメニューを使ってコピーと貼り付けを行います。
csvやxlsから取り込んだデータは行列として使用することができます。実践的な例を見ていきましょう。
例題1. 桜の開花日の推定
桜は2月1日からの平均気温の合計が400度を超えた日、もしくは最高気温の合計が600度を超えた日に開花するという、よく知られた経験則があります。
この例題では、気象庁が公開している気温のデータを使って、2017年の開花日を確かめてみます。なお、ここでは400度、600度を使用するルールをそれぞれ「400度則」、「600度則」と呼ぶことにします。
T = csvread('sakura.csv');
このデータは、各行が1日を表しており、1列目に平均気温、2列目に最高気温が格納されています。まずは400度則を調べるスクリプトを書いてみます。
day_count = 0;
sum_temp = 0;
for k=1:size(T,1) % size(T,1)はTの行数(日数)。日数分繰り返す。
sum_temp = sum_temp + T(k,1); % 平均気温の和の値を更新
if sum_temp > 400
day_count = k; % k日目が答え
break % for ～ end のループから抜ける
end
end
day_count
さて、実行結果を見ると2月1日を1日目として52日目ですから、3月24日が結果になります。ちなみに実際の開花日は3月21日でした。
ここで、 for ～ end は繰り返し同じ処理を行わせるための構文です。ここでは、データの日数分だけ for ～ end の内容を実行させるようになっています。
また、 if は条件分岐のための構文です。ここでは、平均気温の和 sum_temp が400より大きくなったかどうかで処理を分岐させています。上記のスクリプトでの処理は、 sum_temp > 400 が真のとき if ～ end で囲まれた部分の処理を実行し、それ以外の場合は素通りします。
600度則も試してみましょう。まずは下にある答えを見ずに上記のスクリプトを書き換えてみてください。今度は平均気温ではなく最高気温の和であることに注意してください。
day_count = 0;
sum_temp = 0;
for k=1:size(T,1) % size(T,1)はTの行数
sum_temp = sum_temp + T(k,2); % 平均気温の和の値を更新
if sum_temp > 600
day_count = k; % k日目が答え
break % for ～ end のループから抜ける
end
end
day_count
この場合の結果は3月20日になりました。いずれの経験則でも実際の開花日に近い日付が推定ができていることがわかります。
高度な使い方を身につけると、以下のようにループを使わずに目的の動作を実現できます。これはMATLABコードの高速化において重要なテクニックですので、使い慣れてきたら挑戦してみてください。
day_count = find(cumsum(T(:,1))>400,1) % 平均気温の累積和が400を超える最初のインデックスを計算
day_count = find(cumsum(T(:,2))>600,1) % 最高気温の累積和が600を超える最初のインデックスを計算
グラフの作成
レポートや論文を作成する際に、実験やシミュレーションで得られたデータをそのまま貼ってしまうと、読みづらく結果が伝わらない文書が出来上がってしまいます。したがって、グラフを活用してデータを可視化する必要があります。
ここでは、基本的なデータの準備とプロットの仕方、見た目の微調整に便利な機能を紹介します。
例題2. 桜の開花日推定結果の可視化
400度則と600度則それぞれの結果を表すプロット(グラフ)を作成して、先ほどの桜の開花日推定の結果を可視化します。今回は平均気温の累積和を日付に対してプロットして、その上に400度の線を重ねて描いてみます。まずは累積和のデータを準備しましょう。
エラーが出ましたね。どこが間違っているかわかりますか？ここでエラーが出ているのは SUM_AVG(k) = ... の行ですね。エラーメッセージを読んでみると、「添字インデックス」の値が間違っているようです。
SUM_AVG(k-1) の部分に k=1 を代入してみると、ベクトルの0番目の値を要求するコードになってしまっていることがわかります。MATLABでは、行列やベクトルの最初の要素の番号は1ですので、0番目の値は存在しません。したがって、このようなエラーが発生するのです。
今回は、あらかじめ1行目のデータを入れておいて、ループ箇所を2行目からにすることで修正します。プログラムを書く際には、こうしたエラーを解決するためのアイディアも必要になります。
clear
T = csvread('sakura.csv');
SUM_AVG = zeros(size(T,1), 1);
SUM_MAX = zeros(size(T,1), 1);
SUM_AVG(1) = T(1,1); % 2月1日(初日)の値を格納しておく
SUM_MAX(1) = T(2,1);
for k=2:size(T,1) % 2行目からスタート
SUM_AVG(k) = SUM_AVG(k-1) + T(k,1);
SUM_MAX(k) = SUM_MAX(k-1) + T(k,2);
end
ここでデータを見ておきましょう。気温の累積和がきちんと計算できていますね。
SUM_AVG
それではこのデータを早速折れ線グラフにプロットしてみましょう。プロットするにはx軸のデータとy軸のデータが必要になります。今回はx軸のデータとして2月1日を1とした日付番号、y軸のデータとして平均気温の累積和を使用します。見やすいようにグリッドを表示するように grid on コマンドを使用しています。
figure;
plot(1:size(T,1), SUM_AVG);
grid on
プロットのウィンドウを閉じるには、おなじみの閉じるボタンをクリックするか、コマンドウィンドウで close all を実行します。
つぎに、400度の線を引きます。 plot コマンドを繰り返し実行するとそれまでのプロットが消えてしまいます。そこで、2回目の plot コマンドの前に hold on を実行して、プロットを保持するようにします。
figure;
plot(1:size(T,1), SUM_AVG);
grid on; hold on;
plot(1:size(T,1), 400*ones(size(T,1), 1));
これで、気温の累積和のグラフの上に線を引くことができました。
つぎは凡例と軸のラベルを出してみましょう。データをプロットする際は、どの線が何の値を示しているのか凡例をつけ、各軸が何を表しているのかラベルをつけることで、わかりやすくすることが重要です。
figure;
plot(1:size(T,1), SUM_AVG);
grid on; hold on;
plot(1:size(T,1), 400*ones(size(T,1), 1));
legend('平均気温の累積和', '開花条件')
xlabel('日数 [日]'); ylabel('気温 [℃]')
凡例にLaTeX形式の数式を書くこともできます。少し発展的な使い方になるので詳細は割愛しますが、以下の書き方をまねて試してみてください。
figure;
plot(1:size(T,1), SUM_AVG);
grid on; hold on;
plot(1:size(T,1), 400*ones(size(T,1), 1));
legend({'$$\sum_{i=1}^k T_i$$', '開花条件'}, 'interpreter', 'latex')
xlabel('日数 [日]'); ylabel('気温 [℃]')
ここまでで作ったプロットをそのままレポートに貼ってしまうと、おそらく減点を食らいます。以下にレポート向けのプロット作成のチェックリストを示します。
- プロット内の文字は印刷後に読める大きさになっているか？
- 線は細すぎないか？
- 凡例や軸ラベルは適切か？
- 凡例がデータにかぶっていないか？
これらの点に気を付けて最後にプロットの仕上げを行います。仕上げに便利なのが、「プロパティエディター」です。
線の太さ、点線や破線等の形式、フォントの大きさや種類などを変更すると、以下のように最初よりも見やすいプロットができあがります。
ほかにも、つぎのようなコマンドを使うと様々なプロットを作成することができます。以下はほんの一例です。
- subplot : 1つのFigure内に複数のプロットを入れる
- loglog : 対数グラフを作成する
- semilogx, semilogy : 片対数グラフを作成する
- polarplot : 極座標プロットを作成する
その他の機能や詳しい使い方については、こちらのページをご覧ください。
最小二乗法
基礎的でかつ重要なデータ処理として、測定データのモデルへのフィッティングが挙げられます。ここでは最小二乗法を使って誤差の自乗和を最小化するようなモデルのパラメータを求め、測定データをモデルにフィッティングします。
最もシンプルな例として、単純な1次関数へのフィッティングを行ってみます。使用するデータは 
に平均0分散0.1の正規分布に従う乱数を足し合わせたものです。
に平均0分散0.1の正規分布に従う乱数を足し合わせたものです。clear
load linear.mat
figure; plot(x,y,'*')
1次関数 
の場合は誤差の自乗和が最小となるパラメータ 
はつぎのように求められます。
の場合は誤差の自乗和が最小となるパラメータ はつぎのように求められます。
これをMATLABで実装するとつぎのようになります。結果を見ると乱数を加える前の傾き・切片をうまく推定できていることがわかります。
n = length(x);
den = n*sum(x.^2)-sum(x)^2;
a = (n*sum(x.*y)-sum(x)*sum(y))/den
b = (sum(x.^2)*sum(y)-sum(x.*y)*sum(x))/den
figure; plot(x,y,'*',x,a*x+b,'r-')
MATLABに Curve Fitting Toolbox や Optimization Toolbox をインストールすると、直線以外にも曲線や曲面に対するフィッティングが簡単に（全く頭を使わずに）できるようになりますが、データ処理方法の原理を理解した上で使うようにしましょう。
MATLAB Academyでの無料「MATLAB入門」コース
この資料では、MATLABの基礎的な使用方法のうちのほんの一部しか取り上げていません。もう少し詳しい使い方を知りたい場合は、MATLAB Academyの「MATLAB入門」コースを一通りこなしてみるとよいでしょう。
受講にはMathWorksアカウントが必要です。GSICのウェブサイトに記載されている方法にしたがって、アカウントを作成しましょう。
困ったとき・質問があるときは
ドキュメンテーションを読む
うまくいかないとき、いつでも人に聞ける状況ではないと思います。そんな時に役立つのがドキュメンテーションです。使い方のわからないコマンドがあるときは、help コマンドや doc コマンドを使って確認しましょう。help コマンドを使うと、コマンドウィンドウに簡単な使用方法が出てきます。
help eig
これだけでは詳しい使い方はわからないので、「eig のリファレンス ページ」をクリックするか、 doc eig コマンドを実行してドキュメンテーションを読むようにしましょう。ドキュメンテーションはWeb上にも用意されているので、検索エンジンから調べることも可能です。
この資料を書いたMATLAB TAも、普段ドキュメンテーションを読んで確認しながらコードを書いています。
MATLAB TAに聞く
ドキュメンテーションを参照してもよくわからない場合は、ぜひMATLAB TAをご活用ください。オフィスアワーをほぼ毎日開いておりますので、お気軽にお越しください。オフィスアワーの時間帯や場所については、こちらをご覧ください。また、メールやTwitterによるお問い合わせも大歓迎です。
メールアドレス: sim_edu@citl.titech.ac.jp
Twitter: @matlab_titech
MATLAB Answersで質問する
MATLAB開発元のMathWorks公式の質問コミュニティサイトです。MATLABユーザ同士で問題解決の助け合いができます。日本語でも質問ができるので、ぜひ活用してみてください。
(このドキュメントは、MATLABのライブエディターを使用して作成されました。)