MATLAB基本操作復習
先日のExcelの例題
x
y
0.278
0.398
0.522
0.743
0.911
0.937
1.022
1.398
1.521
1.843
（１） ｘとｙを単純にプロットする
2.4739
3.6955
4.1761
5.499
7.124
7.4567
8.2219
12.481
13.69
18.079
← 「散布図」のグラフを描いてみる
（２）近似直線をあてはめる
● これと同じ作業をMATLABで
やってみる。
MATLABのホームタブ、変数メ
ニュー、新規変数で表示される
シートに、Excelから10個の点の
データをコピーペースとして変数名
を付けるだけ。
プロットタブでtestdataの散布図
• 変数testdataをマウスで選び、プロットタブで
scatter(散布図)を描くだけ。コマンドでやるな
らscatterコマンドを使う。
Figureウィンドウのツールメニュー
• データの統計
• 基本的な近似
• 以上、簡単でした。TeXの大レポートのグラフは、
MATLABで描いても良いです。
• グラフの細かい調整はplottoolsで行います。
• ヘルプの「プロットツール」を見て下さい。「グラフ
の注釈」も見て下さい。グラフ資料の作り方の参
考になります。
>> testdata
testdata =
0.2780
0.3980
0.5220
0.7430
0.9110
0.9370
1.0220
1.3980
1.5210
1.8430
2.4739
3.6955
4.1761
5.4990
7.1240
7.4567
8.2219
12.4810
13.6900
18.0790
>> x=testdata(:,1)
>> y=testdata(:,2)
x=
y=
0.2780
0.3980
0.5220
0.7430
0.9110
0.9370
1.0220
1.3980
1.5210
1.8430
2.4739
3.6955
4.1761
5.4990
7.1240
7.4567
8.2219
12.4810
13.6900
18.0790
コマンドで行う場合
観測データプロットと近似式あてはめ
>> y=[2.4739 3.6995 4.1761 5.599 7.124 7.4567
8.2219 12.481 13.69 18.079]
>> x=[0.278 0.398 0.522 0.743 0.911 0.937
1.022 1.398 1.521 1.843]
>> plottools
プロットツールを起動
散布図 oは、アルファベット小文字オー
>> plot(x,y,'o')
1次式（直線）近似の当てはめ結果
>> p=polyfit(x,y,1)
１次の多項式p y=9.7340x-1.0183
p = 9.7340 -1.0183
>> plot(x, polyval(p,x), '-', x,y,'o')
グラフ描画
＊ polyval(p,x)は多項式pの定義域xに対する値域のベクトル
つまり、近似直線上の値
plotコマンド
•
•
•
•
２次元線形プロット関数
plot(y)
yの配列番号を横軸、yを縦軸に
plot(x,y) xを横軸、yを縦軸に
plot(x,y，仕様指定)
仕様指定で線種、マーカ、色を指定
• plot(x,y，仕様指定, ‘LineWidth’，線幅ポイント数)
線幅も指定
仕様指定
• 線種、マーカ、色指定の3組をシングルクォー
ト「’」で囲む 。 ’-*r’ は、赤い色で実線でアスタ
リスクを結ぶ。
• 線種： 「-」実線、「:」点線、「-.」鎖線、「--」破線
• マーカ： 「.」「o」「x」「+」「*」「s」四角形「d」ダイ
ヤ「^」「v」「>」「<」上下左右矢印「p」五角形「h」
六角形
• 色： 「y」黄「m」マゼンタ「c」シアン「r」赤「g」緑
「b」青「w」白「k」黒
ひとつのfigureウィンドに複数のグラフ
• plot(x1,y1,x2,y2,・・・)
(xn,yn)の組み合わせで複数のグラフを描く
• plot(x1,y1,仕様1,x2,y2,仕様2，・・・,
‘LineWidth’，線幅ポイント数)
• hold onコマンドを使う方法
hold on
plot(・・・)
plot(・・・)
２次元線形以外のプロット
•
•
•
•
•
loglog(x,y)
semilogx(x,y)
semilogy(x,y)
polar(theta,r)
bar（x,y)
両対数プロット
x軸対数
y軸対数
極座標
棒グラフ
３次元プロットは、メッシュプロットや表面プロット
グラフの装飾や調整
•
•
•
•
•
title(‘文字列’)
xlabel(‘文字列’)
ylabel(‘文字列’)
ｔext(x, y, ‘文字列’)
axis([xmin, xmax, ymin, ymax])
＊FigureウィンドウのGUIツールでもできる。
GUIとCUI
• 散布図を描いて近似直線を引く、はマウスを
使ってもできるし、コマンドでもできることが分
かった。
• マウスでできるなら、わざわざコマンドを打っ
てやることはない（コマンドは知っていないと
打てないから）、だろうか？
コンピュータリテラシ 第14回
MATLAB（2a）
コンピュータでの数値計算：数値の表現と扱い
今日のテーマ：
コンピュータによる数値情報の扱い
コンピュータの演算処理速度
• コンピュータの演算性能の目安
← １秒間に実行できる浮動小数点演算（実数
どうしの計算）の回数
単位 FLOPS
（floating-point operations per second）
コンピュータが扱う数値
• コンピュータが扱う数値（数値とは「計算の対
象になるもの」）は２種類
整数型数値 ： 小数点がつかない
実数型数値 ： 小数点がつく
• 整数と実数といっても数学で習ったのとは異
なる
整数は実数に含まれるものではない
数値もコンピュータの中では有限桁の２進数
で表される → 扱える範囲が有限
整数型数値
• 正の整数と負の整数を表すため、2の補数表
現を用いると便利
1bit
０
31bits
００････････････・・・・・・・・・・・・・・・・・０
正数については、その数の大きさ
負数については、その数の２の補数
その整数の符号
符号は、１のとき負の数
2の補数は、その数を2進数で表し、すべてのビッ
トの0と1を反転させて１を足した数
２の補数表現
• 3bitsの場合
０１１ ＝
３
０１０ ＝
２
００１ ＝
１
０００ ＝
０
１１１ ＝
- １
１１０ ＝
- ２
１０１ ＝
- ３
１００ ＝
- ４
正の数と負の数の足し算、負の数と負の数の足し算が、（負数
を２の補数表現すれば）通常の2進数の足し算としてできる！
4bitsで表せる整数の範囲
０１１１ 正の整数の最大値
１１１１
０１１０
１１１０
０１０１
１１０１
０１００
１１００
００１１
１０１１
００１０
１０１０
０００１
１００１
００００
１０００ 負の整数の最大値
正の数と負の数の足し算、負の数と負の数の足し算が、（負数
を２の補数表現すれば）通常の2進数の足し算としてできる！
正整数の最大値に１を足すと負整数の最大値になる。。。
２の補数表現
• 4bitsの場合 - 8~7の範囲の整数
• 8bitsの場合 -128~127の範囲の整数
• コンピュータで扱う整数は、32bitsや64bits
表せる範囲の整数を越えて計算しようとする→
桁あふれ（オーバーフロー）という現象
整数は、コンピュータ内で数えたり、（文字コード
のように）区別したりするのに使う。
計算には向かない。
実数型数値
• 仮数部と指数部を分けて表す浮動小数点形
式で表現
→ 大きな数も小さな数も表現して計算できる
• 仮数部の有効桁数は決まっているので、近
似値しか表せない
• （例）10進数の浮動小数点形式で、仮数部4
桁（=有効数字）、指数部2桁が使えるなら、次
の範囲の数を表せる
- 9.999×1099 ~ +9.999×1099
IEEE754規格の倍精度（64bits）実数
1bit
０
11bits
52bits
００・・・・・・０
００････････････・・・・・・・・・・・・・・・・・０
仮数部 ： 2進52桁は、10進数の15桁くらいの有効数字
指数部（符号含む）
仮数部の符号
符号は、１のとき負の数
Mathematicaはがんばる
• Mathmaticaでは、厳密整数や浮動小数点数の
桁数に制約がない。
• それはMathematicaが異常、すごくがんばってい
る。
• 普通は、コンピュータでの計算結果は疑わしいと
思って下さい^^;
• それが証拠に、Excelだって、3.3-3.2もちゃんと計
算できません ← セルの書式設定で、数値の
有効桁数いくつ以上から計算結果がおかしくな
るか観察してみましょう。
誤差と精度に対する考え方
• 計算や計測などで「本当の値とのずれ」のこ
とを誤差と呼ぶ
• 絶対誤差: 誤差の絶対的な量がどれくらいか
• 相対誤差: ある値について、誤差がどれくらい
含まれているか
たとえば「絶対誤差が1度の温度計」があったものとする。
室温(20度とか)を計測 → 相対誤差 5% → ぼろい
熔融金属(1000度とか)を計測 → 相対誤差 0.1% → 優秀
有効数字
• 精度: 「どれくらい正確か」を表す言葉だが曖昧
• 有効数字(有効桁数): 「数字が何桁正しいか」
先の例で
23±1度 → 有効数字1桁
1042±1度 → 有効数字3桁
• 倍精度実数型の場合、有効数字は10進数16桁
くらい。優秀だと思いますか?
「2進の」「浮動小数点」表現での数値計算
→ さまざまな気をつけるべき状況
実数計算の精度
• 有効桁数が決まっていることで、実数の計算では精
度が落ちることがある
（例）今日のポイント！
・ 10で割ったり、0.1を掛けるのは、（人にとっては容
易な計算だが、）0.1は2進表現では循環小数になる
ためコンピュータは精度を落とす計算
・ 実数どうしの計算の精度で、桁数がすごく大きな数
どうしや、すごく大きな数と小さな数を計算するときは
（すべての桁の数が扱えず）桁落ちするのを観察しよ
う
今日のポイント！
数値計算するとき
• 計算には、（整数でなく）倍精度実数を使って
いる
• 計算問題は、人が計算する方法と、計算機に
解かせる方法とでは、かなり違うことがある
(コンピュータなら無理やり力ずくで計算するこ
とも可能で頼もしい、一方で人の予想外とこ
ろで容易に計算精度を落とすことがある)
• コンピュータの計算には、誤差を回避する工
夫が必要なことがある
コンピュータの計算結果は、すべて近似計算と思え。
本日のMATLAB実習
(1) MATLABで計算誤差について観察・体験しよ
う。コンピュータでの計算結果が誤差まみれ
（になりかね）な（い）ことをぜひ体験しておきま
しょう。
(2) 円周率の値の近似値を複数の方法で求めて
みる計算などを例題に用いています。
(3) そのために、MATLABで簡単なプログラムを動
かします。プログラムの中身に興味のある人は
それも見て下さい。
＊ Mathematica、R、Excelも、プログラミングの機能を備えています。
準備１
(1) ドキュメントのMATLABフォルダ内にrenshuフォルダを
作っていない人は、先に作りましょう。
(2)クラスのページのリンクから圧縮フォルダmatfunc.zipをダ
ウンロードし中身のファイルを確認しましょう（.mファイルが５
つ）。
(3)５つのファイルをMATLAB/renshuフォルダにコピーしま
しょう。
準備２
●MATLAB起動
●カレントディレクトリをrenshuフォルダに移る。
（先ほど保存した.mファイル5つがあることを
確認）
本日の教科書
MATLAB/Scilabで理解する数値計算、櫻井鉄也著、
東京大学出版会、ISBN978-4-13-062450-3
• １章～２章を使用。主に「第2章 有限桁の数
値」のところの内容がコンピュータの数値計
算の「基本のき」を取り扱っています。
• MATLABはプロプライエタリソフトですが、
scilabはMATLABと互換性の高いフリーソフト
です（構文が違うところも多そうですが、比較
しながら使えるようです）。
MATLABプログラミング
～分岐の基本構文～
• 条件分岐の基本構文
if 式
命令
end
「式」は条件式で、式が成り立って
いるとき「命令」実行する。ifを打っ
てからendを打つまでは、改行しな
がら命令を書いてよい。
>> k=5
>> if k>0
disp(k)
end
if 式
命令
elseif 式
命令
else
命令
end
式が成り立たないときは、elseifで
次の条件式で判定する。それを１
回以上繰り返した後、最後までど
の条件式も成り立たなかった場合
を扱うにはelseを用い（省略も可）、
if文はendで終わる。
MATLABプログラミング
～繰り返しの基本構文～
• 繰り返しの基本構文
for ループ変数 = ループ変数の範囲
命令
end
ループ変数の値を「範囲」で変化さ
せながら「命令」繰り返し実行する。
ループ変数の変化の範囲は、ベク
トルで与える。
>> for k=1:10
disp(k)
end
>> v=5:-1:1
>> v
>> for k=v
disp(k)
end
while 条件
命令
end
条件が成り立つ間、「命令」を繰り返
し実行する。「命令」を実行したら、
また条件を調べる、というようにして。
>> a=10;
>> while a>=1;
disp(a);
a=a/2;
end
MATLABプログラミング
～ ファンクションMファイル（入出力あり） ～
• 関数（ファンクションMファイル）の定義
>> edit func.m
エディタでfuncの関数定義を書き保存する。func.mのあるディ
レクトリでは関数funcを（読み込まなくても）利用できる。キー
ワードfuncionを使って定義する関数には、入力を渡して実行
し、出力を受け取ることができる。
● func.mの中身を打ち込もう
以下のように使う。
>> [C, d] = func(3, sqrt(2))
function[A,b]=func(x,y)
A=eye(2,2)*x;
b=[y;y^2];
プログラミング
～ スクリプトMファイル（入出力なし） ～
• 関数（スクリプトMファイル）の定義
>> edit enshuritu.m
エディタでMATLABインタプリタに打つものをmファイルに書い
ておき、実行することができる。実行時に入力は渡せない。
>> enshuritu
n=[100 1000 10000 100000 1000000];
m=length(n);
for k=1:m
z=rand(n(k),2);
r=z(:,1).^2+z(:,2).^2;
p(k)=sum(r<1)/n(k)*4;
end;
disp(p);
乱数によるπの値のシミュレーション
を例題に。
• 1×1のマス内に点(x,y)を0～1の一様乱数で
沢山打つものとする
• 打った点が1/4円内に入る確率や回数は、外
側の面積１の正方形と1/4円の面積π/4の比
に従う
• 打った点が1/4円内に入る確率 ～ 「1/4円の
面積÷全体の面積」 = π/4
• これを利用してπを求めてみよう
y
１
０
１
x
enshuritu.m
• 乱数を使って円周率の近似値を求める。試行
回数をn= 100, 1000, 10000・・・と増やすと真
の円周率の値に近づくのが分かる。
n=10（10回乱数を発生）の場合だと、こんな計算をやっている
>> z=rand(10,2)
１０個分の点の座標（x,y）,0=<x,y<1のデータを乱数で生成
>> z(:,1)
>> z(:,1).^2
「.」（ピリオド）は行列の要素ごとに演算を行うのでした
>> r=z(:,1).^2+z(:,2).^2;
１０個の点それぞれの原点からの距離を計算し、
>> r
>> p=sum(r<1)/10*4;
>> p
その距離が１より小さいものを合計する。その合計は、nが大
きくなると半径１の1/4円の面積（π/4）と、外側の一辺が１の正
方形の面積（１）の比に近づくことを利用してπの値を計算する
他のプログラミング言語でやると
public class Test {
public static void main(String[] args) {
int n = Integer.parseInt(args[0]);
int count = 0;
for(int i = 0; i < n; ++i) {
double x = Math.random(), y = Math.random();
if(x*x + y*y <= 1.0) ++count;
}
System.out.printf("%.20g¥n", 4.0 * count / n);
}
}
乱数発生（サイコロを振る）回
数が増やしてみた。乱数の値
は毎回違うからこんなもの。C
MATLABのプログラムは簡潔に書ける
から好き、という人が多いらしい
> java Test 1000
3.1400000000000000000
> java Test 10000
3.1308000000000000000
> java Test 100000
3.1414400000000000000
> java Test 1000000
3.1408160000000000000
コンピュータは力任せに、
そもそも近似計算をする道具
• コンピュータはいくら計算しても疲れないので、
連続量をサンプリングする精度を細かくして
いったり（細かくしても誤差はあるので、このと
き真の値との間に生じる誤差を量子化誤差と
いう）、シミュレーションの試行回数を多くする
ような計算に用いることが多い。
そもそもの近似計算：量子化
• 本来連続的なものを「個数」に置き換える（サ
ンプリング、標本化に伴う）ことの誤差を『量
子化誤差』という
• 多くの場合、量子化誤差はアナログ量をディ
ジタル量に変換するところで現れる
• たとえば「0.1kg単位で表示されるデジタル体
重計」は50gの量子化誤差
量子化誤差の例
• 「半径1mの円の面積=π（m2）」のはずだから、
• 「2m四方の紙を1cm単位に分けて、円内に
入っているものの数を数える」ことで面積が求
まるはず。
• 上下/左右対称だから1/4円についてやって4
倍する
← 実際この方法でπの値を求めてみると、
１ｍ
Ｎ
１ｍ
・
・
・
２
１
1マス
１ ２
・・・
Ｎ
1/n×1/n （ｃｍ２）
enshurituR.m
n=[100 1000 10000];
m=length(n);
for k=1:m
count=0;
for x=0:n(k)-1
for y=0:n(k)-1
if x*x+y*y<=n(k)*n(k)
count=count+1;
end;
end;
end;
p(k)=4*count/(n(k)*n(k));
end;
disp(p);
乱数を使うenshuritu.m
は毎回結果が違うが、こ
ちらは毎回同じ結果
個数を増やして近似（→
精度を上げる）すればす
るほど、真の値に近づく
ただし精度を上げすぎてもあまり
意味がないこともある
テーマ:
サイコロを振る
シミュレーションと大数の法則
そもそもの近似計算：
・乱数を使ってのシミュレーションに基づく数値は、乱数を
使っているからもともと毎回違う誤差がある（この誤差に特に
名前はついていない）
・ただ、「サイコロを振る回数Nを増やすと正確な値に近づく
（大数の法則）」
・誤差は 1/√N に比例（中心極限定理）
・実際にNを増やしながら試してみよう（100倍にすると誤差が
1/10になるか）
• このような、「サイコロを振って値を求める」方
法もさまざまなところで使われている。『モン
テカルロ法』
大数の法則（サイコロシミュレーション）
• saikoro.mで、大数の法則の様子を見てみよう
• サイコロの目が1/6の確率で出るというような
確率も、試行回数を増やすとゆらぎながらそ
の値に近づく。
• このスクリプトMファイルの出典。
「信号処理」「画像処理」のためのMATLAB入
門、高井信勝著、工学社、ISBN978-4-875933C3055
二進数による実数の表現
を見てみよう！
２進の浮動小数点による数の表現
>> f=[1 1];
>> f .* 0.5 .^(1:length(f))
>> sum(f .* 0.5 .^(1:length(f)))
0.11(2)= (1/2)+(1/2)^2
>> f=[0 1 0 1 0 1];
>> f .* 0.5 .^(1:length(f))
>> sum(f.*0.5.^(1:length(f)))
0.010101(2)
コンピュータ内部での実数は、２進の指数表現。
1.1(2) × 2-1
仮数部が、「１．-----」（小数点より左が必ず１になるよ
うな形）に正規化した後に、数の大きさを指数部で表
している。
1.0101(2) × 2-2
IEEE754規格の倍精度（64bits）実数
1bit
０
11bits
52bits
００・・・・・・０
００････････････・・・・・・・・・・・・・・・・・０
仮数部 ： 2進52桁は、10進数の15桁くらいの有効数字
指数部（符号含む）
仮数部の符号
符号は、１のとき負の数
次の計算を打ち込んでやってみよう
何が起きたか、考えよう
浮動小数点による数の表現の限界による誤差
オーバーフロー、アンダーフロー
>> realmax
最大の数
>> realmin
最小の数
realmaxやrealminの周辺で起きる
>> (1e150)^2
>> (1e200)^2
無限大
>> 1/Inf
ゼロ
オーバフローとアンダーフローを観察
>> (1e-150)^2
>> (1e-200)^2
ゼロ
>> 1/(1e-200)^2
無限大
実数の丸め誤差による計算結果の狂い
>> x=1；
<= x=1からはじめ
て、 1+xが１と等し
くない間xを半分に
していく。
>> while (1+x ~= 1)
x=x/2;
end
>> x
x=
1.1102e-016
<= While文が停止した時の
ｘ の値？？？
実数の丸め誤差：0.1の２進表記
• 10進の0.1は、2進表
記ではどのように表せ
るか、自分で手計算し
てみよう。
• コンピュータでの実際
のdouble型実数0.1が
どのように表されてい
るか、bitex.mで調べよ
う。
>> bitex(0.1)
ans =
1001100110011001100110011001100
110011001100110011010
「bitex(数値)」は、その数値が倍精度実数で表
現されたときの仮数部の52ビットぶんを見
せてくれます。
0.1の二進表現に係わる
計算の狂い
>> x=0.1;
>> while x<1
x=0.1からはじめて、x が 1 より小さい間
0.1を足す計算のはずが・・・
x = x + 0.1; disp(x);
disp(bitex(x));
end
●上のwhile文を実行したら、どんなまずい事が起るか。
また、なぜそんな事になるのか簡単に説明せよ。
fl.m
• 関数ｆｌは、数値ｘを10進の有効桁数５桁に丸め
る関数である。 使用例：
>> pi
>> format long
>> pi
>> fl(pi)
←有効桁５桁
• 倍精度実数計算ではだいたい１０進数１６桁の
有効桁数であるが、何桁であっても有限桁数
が決まっていればそれによる誤差が生じること
を関数flを使って観察しよう。
有効桁数に伴う情報落ち
>> x1=121730; x2=3.1481;
>> x1+x2
ans =
1.217331481000000e+005
有効桁数１６桁の正しい答え
>> fl(x1+x2)
ans =
121730
赤： 誤差の大きな結果
有効桁数が５桁しかないと、結果は
こうなる。x2の値は無いのと同じ。
>> a=8.1651; b=8.1653;
有効桁限界あたりの怪しい計算
例： ２数の平均
>> (a+b)/2
ans =
8.16520000000000
正しい答え
同じ計算を有効桁数５桁
ですると、誤差が大きい
>> fl(fl(a+b)/2)
ans =
有効桁数５桁でも工夫に
よって正しい答え
>> fl(a+fl(fl(b-a)/2))
ans =
平均値であり得ない値！
8.165200000000001
8.164999999999999
赤： 誤差の大きな結果。なぜ誤差
が大きいのか？
緑：計算の工夫によって正しい答え
に近い値。どういう工夫をしたのか？
中間結果のオーバフロー
>> a=2.0e200;
>> b=1.0e200;
>> c=a^2 + b^2
c=
Inf
cは非常に大きな値のため
オーバーフロー
C= 𝑎2 ＋𝑏 2
を計算することはで
きないのか？
© Copyright 2024 ExpyDoc
