カーブフィッティング (回帰分析) の練習
線形最小二乗フィッティング: R, python, excel
非線形最小二乗フィッティング: mathematica
目次
1
最小二乗法
2
2
線形フィッティング (線形回帰) とは
2.1 線形フィティングの要件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3
線形最小二乗法の理論の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.2
3
直線偏光測定のモデル
5
4
線形フィッティング – 統計解析言語 R 編
6
5
4.1
4.2
4.3
R の作業ディレクトリ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
R の使い方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
データ入力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
6
8
4.4
4.5
4.6
フィッティング . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
9
10
グラフの作成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ひとまとめにする . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
線形フィッティング – 科学計算パッケージを使う
10
5.1
5.2
5.3
使用するシステム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4
5.5
モデル構築とフィッティング . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
14
5.6
スクリプトの作成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
データ入力
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
データ読み込み . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
グラフ表示
11
12
12
6
線形フィッティング – Excel 編
16
7
非線形フィッティング – mathematica 編
17
17
7.1
データ入力
7.2
7.3
モデルの入力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
フィッティングの実行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
18
19
A Mathematica の基礎
A.1 起動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
21
A.2 式の入力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 単純な定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
21
A.4 関数定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.5 リスト . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
22
測定可能な 2 つの量の間の関係が理論式であらわされていて，その式が未知のパラメータを含ん
でいる場合，いろいろな条件で得た複数の測定データから未知パラメータの値を求めることができ
る． これをカーブフィッティング (curve fitting), あるいは回帰分析 (regression analysis) という．
理論式はモデル式ともいう．
測定誤差のため，普通，モデル式を測定データに完全に一致させることはできないので，できる
だけ多くの測定データを使って，パラメータの値を推定するしかない． 推定計算の指針は，デー
タとモデルの「不一致の程度」が最小になるようなパラメータ値を求めるということである． 「不
一致の程度」をあらわす量としては，残差 (測定データとモデルの差) の二乗和が使われることが
多い．測定誤差がランダムに生じるものであるならば，残差の二乗和を最小にするようなパラメー
タ値が最も確からしいパラメータ値であるということが言えるからである (ガウス・マルコフの定
理)．この考え方を最小二乗法という．
1
最小二乗法
モデル式が従属変数 y を独立変数 x の関数としてあらわすようなものであるとする．未知パラ
メータ a を含めると，y は x と a の 2 変数関数として，
y = f (x, a)
(1)
のように書ける．m 組の測定値を x1 , y1 , x2 , y2 , ..., xm , ym とあらわす．i 番目の測定値の残差は yi
と f (xi , a) の差である．これを ei であらわすと，
y1
y2
ym
= f (x1 , a) + e1
= f (x2 , a) + e2
..
.
= f (xm , a) + em
(2)
残差の二乗和を J とする． すなわち，
J = e1 2 + e2 2 + ... + em 2
(3)
ei をひとまとめにして m 次元列ベクトル e にすると，J はベクトル e のノルムの 2 乗であり，e
と e の内積でもある．
e1
e2
e=
J = ∥e∥2 = e · e
(4)
.. ,
.
em
2
a を変えれば J は変化するが，J ≧ 0 であるから J の値には下限がある．未知パラメータの推定
とは，a の関数である J の最小点を求めることである．
パラメータが複数あっても同じである．n 個のパラメータ a1 , a2 , ..., an のがあるなら，これらを
まとめて一つの n 次元列ベクトル a とみなせばよい．
a1
a2
a=
..
.
(5)
an
f (x, a) は f (x, a) に変えればよい．
f (x, a) がどういう式であらわされるかに関係ない方法は非線形フィッティング (nonlinear fitting)
という． f (x, a) が a1 , a2 , ... の 1 次式であらわされる場合に特化した方法は線形フィッティング，あ
るいは線形回帰という．特に，f (x, a) が x についても 1 次式になっている場合は，f (x, a) = a1 x+a2
という形になってしまうが，この場合は直線回帰 (linear regression) という．
直線回帰は計算量が少ないので手計算でもなんとかなるし，統計機能を備えた電卓や表計算ソフ
トのグラフ作成機能を使えば手軽におこなえる．一般の線形フィッティングでは，面倒な行列計算
が必要なので，行列計算の機能を備えたコンピュータープログラムを使用する．非線形フィッティ
ングには，複雑なアルゴリズムによる大量の計算をおこなうコンピュータープログラムを使う．
後で見るように，線形フィッティングは１次方程式を解くことに帰着するので，正しい答えが得
られるどうかは単純な条件で決まる．非線形フィッティングは試行錯誤によって J の最小点を探し
出すので，正しい答えが得られたかどうかを判断することが難しい場合がある．また，問題によっ
て最適なアルゴリズムが異なるかもしれない．そのため，未知パラメータになんらかの置き換えを
おこなうことで，線形な式にすることができるかどうかを検討するのが得策である．
線形フィッティング (線形回帰) とは
2
2.1
線形フィティングの要件
フィッティングをおこなうにはモデル式を用意する必要がある．線形モデルの場合は，
f (x) = a1 g1 (x) + a2 g2 (x) + ... + an gn (x)
(6)
のように，x の関数 g1 (x), g2 (x), ..., gn (x) を係数とする未知パラメータ a1 , a2 , ..., an の線型結合
の形である．f (x1 ) = a1 g1 (x1 ) + a2 g2 (x1 ) + ... + an gn (x1 ) などを式 (2) に代入すると，
y1
= a1 g1 (x1 ) + a2 g2 (x1 ) + ... + an gn (x1 ) + e1
y2
= a1 g1 (x2 ) + a2 g2 (x2 ) + ... + an gn (x2 ) + e2
..
.
ym
(7)
= a1 g1 (xm ) + a2 g2 (xm ) + ... + an gn (xm ) + em
g1 (x1 ) などはモデル式の n 個ある項に m 個ある測定点の値を代入して得られる mn 個の値である．
式 (7) は a1 , a2 , ... を未知数とする線形連立方程式とほぼ同じである．そこで，ベクトルと行列を
使って簡潔に表現しよう．
3
y1 , ..., ym と e1 , ..., em を m 次元列ベクトル，a1 , ..., an を n 次元列ベクトル，g1 (x1 ), ..., gn (xm )
を m 行 n 列の行列にまとめる．
y1
y2
y=
.. ,
.
e1
e2
e=
.. ,
.
ym
g1 (x1 )
g1 (x2 )
G=
..
.
em
g2 (x1 )
g2 (x2 )
..
.
g1 (xm ) g2 (xm )
···
···
..
.
gn (x1 )
gn (x2 )
..
.
(8)
· · · gn (xm )
すると，式 (7) は次のようにまとめられる．
y = Ga + e
(9)
G は計画行列 (design matrix) と呼ばれる．その i 行 j 列の成分は Gij = gj (xi ) である．
残差の二乗和 (J = e · e) を最小にするという意味でこの方程式 (9) を「解く」ことができる．こ
れが線形最小二乗法 (linear least-squares method) であり，この解を最小二乗解という．ただし，
未知パラメータを決めるには，未知パラメータの数よりも多くの測定値が必要であるから，a の次
元 n よりも式の数 m の方が多くなければならない．いいかえれば，G は縦長の行列になっている
ことが必要である．
2.2
線形最小二乗法の理論の概要
a の微小変化 da に対する J の変化 dJ を計算しよう．式 (9) から e = y − Ga であるから，e の
変化は de = −Gda となる．また，J = e · e であるから，dJ = 2de · e．よって，
(
)
dJ = 2(G da) · (G a − y) = 2 da · G∗ (Ga − y)
となる．ただし，∗ で随伴行列をあらわす．J の極小点では a のどの向きへの微小変化に対しても
J の変化は 0 であるから，G∗ (Ga − y) はゼロでなければならない．よって，
G∗ Ga = G∗ y
(10)
が得られる 1 ．
G∗ G は n × n の半正定値対称行列である 2 ．(10) が一意的に解けるには，G∗ G は正則でなけれ
ばならないので，G はフルランク (階数が n に等しい) でなければならない．これは G の各列が
独立ということであるが，式 (8) を見ればわかるように，それには，少なくともモデル中の関数
g1 (x), ..., gn (x) が独立でなければならない．これはモデル式が冗長ではないということである． ま
た，少なくとも n 個の独立な測定点が必要である．だから，同じ測定値をコピーして増やすのは
だめ．
n が小さい場合は，式 (10) を直接解いて a を求めることができる．直線回帰はそういう場合に
あたる． n が大きい場合にも正確に解けるコンピュータープログラムでは，行列を巧妙に処理し
て解くようになっている 3 .
√
(10) は次のように考えても得られる．a が自由に動くと，Ga は n 次元の超平面 s を形成する． J = ∥y − Ga∥ だ
から，s 上で y との距離が最小になる場所が Ga である．したがって，Ga − y は s に垂直なベクトルになるが，s は G の n
本の縦ベクトルが張るのであるから，Ga − y はこれらの縦ベクトルとの内積がすべて 0 である．よって，G∗ (Ga − y) = 0.
2 内積の性質から (Ga) · (Ga) ≧ 0 であるが，内積を行列積で書くと，a∗ (G∗ G)a ≧ 0 となる．これが G∗ G の半正定
値性．対称性の方は，(G∗ G)∗ = G∗ (G∗ )∗ = G∗ G
3 QR 分解や特異値分解が利用される．QR 分解では，G の像空間の正規直交基底を構成して，式 (10) を解ける形に
1式
4
3
直線偏光測定のモデル
偏光フィルターは，入射光の偏光の向きに対して回転させることによって光の透過率が変化する
ものである．理想的な偏光フィルターでは，回転角度 x に対する透過率 t の関係は次の式であらわ
される．
t = P cos2 (x − b) + P ′
(11)
b は偏光フィルターへの入射光の「最大偏光方向」の角度である．係数 P と P ′ は入射光の偏光状
態で決まる 4 ．実際には偏光フィルターには偏光によらない反射や吸収があるし，直接測定するの
は偏光フィルターを通過した光の強さを光センサーで何らかの量に変換したものであるが，これら
の過程が線形性を有し，入射光強度が一定である場合は，測定値 y は
y = a cos2 (x − b) + c
(12)
となる．パラメータの a と c は入射光強度やセンサーの感度に依存するし，c はセンサーのゼロ点
のずれを含むかもしれない．しかし，b にはそのような依存性は無いと考えられる．
式 (12) は a と c について線形であるが，b は非線形な関数 cos の中に入っているので，b について
は非線形である． このままでは非線形フィッティングをおこなう必要があるので， 変形を試みる．
倍角公式と加法定理を使って式 (12) を変形すると，次のようになる．
y=
a cos 2b
a sin 2b
a
cos 2x +
sin 2x + + c
2
2
2
よって，次のように，パラメータを A, B, C に置き換えれば線形になる．
y = A cos 2x + B sin 2x + C,
a
a
a
A = cos 2b, B = sin 2b, C = + c
2
2
2
√
1
a
a = 2 A2 + B 2 , b = arctan(A, B), c = C −
2
2
(13)
(14)
(15)
ここで，arctan(A, B) は平面のベクトル (A, B) の方位角である．プログラミング言語や数値計算パッ
ケージの多くはこの関数が計算できる 5 ． 普通の arctan を使う場合は，A > 0 なら arctan(B/A),
それ以外で B < 0 なら −(π/2) − arctan(A/B), それ以外なら (π/2) − arctan(A/B) とする．
する．
特異値分解は，行列 A に対して，適当な直交行列 U,V と対角行列 S によって A = USV−1 となるようにすることで
ある．S の対角成分は A によって一意的に決まる．任意の行列 A に対して「擬似逆行列」A+ が一意的に決まることも知
られている．A の特異値分解によって，A+ = VS+ U−1 となる． ここで，対角行列 S の擬似逆行列 S+ は，S のゼロで
ない成分をすべて逆数にしてから転置するだけで得られる．擬似逆行列によって最小二乗解は a = G+ y で求められるが，
これが一意解であるのは G がフルランクの場合である．
なお，行列 A に対して，BAB = B, ABA = A, (AB)∗ = AB, (BA)∗ = BA を満たす行列 B を A の擬似逆行列とい
い，A+ と書く．
4 形式的に，強さ S の光を，強さ S の直線偏光成分と，それに独立な強さ S の成分とに分けて考えることができ，式
1
2
(11) の係数 P と P ′ は，P = S1 /S ，P ′ = S2 /(2S) = (1 − P )/2 とあらわせる．
5 C 言語など多くの言語にある atan2 という関数は atan2(B,A) のように使用する． Y 座標が先，X 座標が後である．
Excel や 計算サイト http://keisan.casio.jp の atan2 や，Mathematica の ArcTan は X 座標が先．
5
線形フィッティング – 統計解析言語 R 編
4
R 言語は統計解析に特化したシステムなので，データ・フィッティングは比較的楽にできる．
4.1
R の作業ディレクトリ
Windows ではスタートメニューやタイルで R を起動する．Linux ではメニューやデスクトップ・
アイコンで起動するか，端末エミュレーターの中で R と打って起動する．R を終了するには，R
の中で quit() と打つ．このとき，作業内容を保存するかどうかを答える必要がある．
ファイルの読み書きのために「作業ディレクトリ」が重要である．端末エミュレーターで R を
起動する場合は R を起動したときの端末の現在のディレクトリが R の作業ディレクトリになる．
これは R を起動する前に pwd コマンドで見ることができる．また，cd コマンドによって変更す
ることもできる．しかし，R を GUI で起動すると，作業ディレクトリがわからないかもしれない．
R では次のコマンドで作業ディレクトリを確認することができる．
getwd()
これは ”getwd” という名前の関数を引数無しで呼び出す式であり，この式の評価結果が表示される．
getwd の wd は”working directory”から来ている．また，dir() で作業ディレクトリにあるファイル
とサブディレクトリのリストが得られる．サブディレクトリのみのリストは list.dirs(recursive=0)
で得られる．作業ディレクトリを変更する必要がある場合は setdir 関数を使う．たとえば，現在
の作業ディレクトリの中にある xxxx という名前のサブディレクトリに移るには，setwd("xxxx")
とすればよいし，現在のディレクトリの親ディレクトリに移るには，setwd("..") とすればよい．
これらは現在の作業ディレクトリからの相対パスである．ユーザーのホームディレクトリに移るに
は，setwd("~") とする．Windows では，以上のような操作を GUI でもおこなうことができる．
これ以降は，
「データ入力」のセクションと「ひとまとめにする」のセクション以外は読むのを省
略することもできる．
4.2
R の使い方
R では簡単な数値計算がすぐにできる．四則演算の演算子は +, -, *, / である．べき乗は ^ でも
** でもよい．式をくくる括弧にはすべて ( ) を使う．関数呼び出しはすべての引数を括弧で囲ん
で関数名の後ろに付ける．
(1+2)/2
sin(0.1)
tan(0.1)
複数の数値からなるデータは c という関数を使って，たとえば，
v <- c(1,3,5,4)
v
v+10
v*10
6
のように書くことができる．これは４つの数値からなるデータに v という名前を付け，v の値や各
要素に数値を各要素に足したりかけたりした結果を表示している．このようなデータは「ベクト
ル」と呼んでいる．たくさんの要素からなるベクトルを簡単に作る方法もいろいろ用意されてい
る．重要なのは等差数列で，コロン”:”を使って，たとえば
0:40
と書くと，0 から 40 までの公差 1 の等差数列 (0, 1, ..., 40) ができる．次のように，このベクトルに
0.5 をかけると，(0, 0.5, ..., 20) になり，さらに 10 を引くと，(−10, −9.5, ..., 10) になる．
x <- 0:40*0.5-10
x
同じことは，seq という関数を使ってもできる．
x <- seq(-10,10,0.5)
x
seq はたぶん sequence から来ている 6 ．
数値データをプロットする関数は plot である．y = x2 のグラフを描いてみる．
x <- seq(-10,10)
y <- x^2
plot(x,y)
21 点しか計算していないので，表示もそれだけである．
「任意の数値に対してその２乗を計算する」といった計算方法自体は「関数」である．関数を定
義してそれに名前を付け，後でその関数を使って計算をおこなうことができる．
f <- function(x){x^2}
f(2)
f(10)
f(2)+f(10)
f(0:10)
関数定義のところは，
function(仮引数){仮引数を含んだ式}
という形式である．中括弧の中には関数の定義が入る．簡単な計算では結果をあらわす式だけでよ
いが，もっと複雑な処理の記述もできる．
関数定義を使うと，plot 関数によって曲線も描くことができるが，左端と右端を指定する必要
がある．
plot(f, -10, 10)
これにデータ点を重ねるには，points 関数を使う．
x <- seq(-10,10,2)
points(x, f(x))
6 関数の説明を見るには，?seq
のようにすればよい．Unix では，カーソルキーで表示をスクロールし Q で元に戻る．
7
4.3
データ入力
測定データは任意のテキストエディタを使って入力して適当なファイル名で保存する．保存する
場所は R の作業ディレクトリにする．次のは練習用データである．# はコメント行を示している．
# レーザーの偏光 4 月 1 日
# 回転角 (度) : 透過光強度 (ボルト)
0 0.809
30
60
90
0.403
0.039
0.077
120
150
180
0.477
0.836
0.794
210
240
0.400
0.040
270
300
330
0.081
0.481
0.843
360
0.782
このファイルを R で読み込む．ファイル名を xxyy.dat とすれば，
data = read.table("xxyy.dat")
読み込んだデータに data という名前を付けたが，日本語のほうがよければ，生データ=read.table("xxyy.dat")
のようにしてもよい．測定データのプロットは次のようにすればすぐにできる．
plot(data)
データ・プロットをよく見れば，誤った測定点が見つかる可能性がある．
data
で data を表示すると，1 列目に回転角，２列目に透過光強度が入っていることがわかる．data か
ら各列のデータを取り出すには，次のように２重の大括弧を使う．
x = data[[1]]
y = data[[2]]
回転角が独立変数なので，それに x という名前を付け， 透過光強度には y という名前を付けた．
4.4
フィッティング
あらためて，データプロットを作っておく．今度は横軸とグラフ全体に名前を付けよう．
plot(x,y,xlab="角度 (度)",main="テストデータ")
8
角度の単位を度にしたので，モデル式 (13) は，
(
)
(
)
A cos (π/90)x + B sin (π/90)x + C
という形になる．線形フィッティングには lm 関数を使う．
r <- lm( y ~ cos(x*pi/90)+sin(x*pi/90) )
lm に指定するのは，従属変数 y と~ とモデル式から未知パラメータを取ってしまったものである．
名前 r を使って計算結果を見ることができる．
lm 関数が計算したさまざまな結果全体を見るには，summary(r) とすればよいが，未知パラメー
タの値だけ取り出すには次のようにする．
p = r$coefficients
p
p の各要素に対応しているモデル式の項が表示され，定数項 C がはじめに来て，次に cos の係数
A，sin の係数 B の順になっていることがわかるので，
C <- p[[1]]
A <- p[[2]]
B <- p[[3]]
式 (15) を使って a, b, c の計算もやっておく．R では atan2 関数が使える．
a <- 2*sqrt(A^2+B^2)
b <- 0.5*atan2(B,A)*(180/pi)
c <- C-a/2
a;b;c
4.5
グラフの作成
理論曲線を描くために，フィッティングの結果を入れて関数を定義する．(12) 式に戻ろう．
f <- function(x) a*cos((x-b)*pi/180)^2+c
角度は度からラジアンに直した．準備ができたので，データ点と理論曲線を描く．
plot(f,0,360,xlab="角度 (度)",main="テストデータのフィッティング")
points(x,y)
grid()
始めの plot では少し飾りを付けた．また，grid() で目盛線を入れている．
グラフは PNG か JPEG ファイルに保存するのが簡単である．保存してから印刷したりレポー
トに入れたりすればよい．
savePlot("xxxxyyyy.png","png")
9
4.6
ひとまとめにする
データの読み込みからグラフの保存までを１つの関数に入れ，analyze という名前をつけよう．
この関数には，データファイルの名前を引数として与えることにする．下の関数定義を PDF から
R へコピー＆ペーストしても動くかもしれない 7 ．
analyze <- function(datafile) {
data <- read.table(datafile)
# データ読み込み
x <- data[[1]]; y <- data[[2]]
# x 列， y 列
r = lm(y ~ cos(x*pi/90)+sin(x*pi/90)) # フィッティング
p <- r$coefficients
# 未知係数の値
C <- p[[1]]; A <- p[[2]]; B <- p[[3]]
a <- 2*sqrt(A^2+B^2)
b <- 0.5*atan2(B,A)*(180/pi)
c <- C-a/2
print(c(a=a,b=b,c=c))
# グラフを描く
h <- max(y)*1.1-min(y)*0.1
# 結果をまとめて表示
# 上に 10% 余裕
l <- min(y)*1.1-max(y)*0.1
# 下に 10% 余裕
plot(function(x) a*cos((x-b)*pi/180)^2+c, 0, 360,
xlab=" 角度 (度)", ylab=" 強度 (V)", ylim=c(l,h),
main=paste(datafile, ", b=", b) )
points(x, y)
grid()
savePlot(paste(datafile,"png", sep="."), type="png")
}
これで，次のようにするだけで，パラメータ値を計算して表示し，グラフを画像ファイルに保存
してくれる．
analyze("xxyy.dat")
R を終了するときに，”Save workspace image? [y/n/c]: ” と訊かれたら，y で答える．すると，次に
同じディレクトリで R を起動すると，前におこなった定義がそのまま使える．この情報は ”.RData”
というファイルに保存されている．何もかも終わったら，このファイルを消去すること．
5
線形フィッティング – 科学計算パッケージを使う
線形フィッティングは数値計算法の定番であるので，他のさまざまな数値計算法と同様，Excel
でも Fortran でも何でも使っておこなうことができる．計算は関数を１つか２つ呼び出すだけの
ことであるが，データ入力と結果の可視化に手間がかかるようだと困る．データを並べる順番を確
かめたり結果の妥当性を確認したりしながらの試行錯誤も必要になる．
7 これが入っているファイルを入手することができた場合は，R
10
で一度だけ source(”ファイル名”) のようにすればよい．
5.1
使用するシステム
パソコンで動く数値解析用ソフトとして，Matlab が普及し，他のものでも， Matlab を手本に
して作られているのが多い．それらは Fortran で培われた数値計算サブルーチンを誰でも簡単に使
えるようにするために，対話実行環境と高度な表示機能を装備したものと考えればよい．Matlab
と同じ使い方ができるフリーソフトとしては，Scilab や Octave がある．汎用プログラミング言語
でも，適当なライブラリーによって Matlab と同様の機能が使えることが多いが，使い方を早く
マスターするには，インタープリターがよい．ここでは ’Python’ 言語を使ってフィッティングを
やってみよう．
Python 言語の文法 8 は理解していなくてもよい． 数値処理とグラフィックスのために ’pylab’
9
というライブラリーを使う．実行環境として ’spyder’ を使うことにする． spyder の中でデータ
入力も計算も表示もプログラミングもできる 10 ．spyder を起動するには，たとえば端末に次のよ
うにタイプする．
spyder
あるいは，メニューから起動できるかもしれない．
spyder のウィンドウは大きく３つの枠に分かれているが，そのうち Console に式やコマンドを
打ち込めば，その場で計算の実行ができる．1+2 とか 2*3 とか 12/-4 とか sin(0.1) とか tan(pi/4)
などとやってみれば，ただの電卓としての使い方がわかるが，注意すべき点が少しある．
• 1, 2, 3 は 1.0, 2.0, 3.0 とは異なる． 前者は「整数型」後者は「実数型」である． Python 2
では，整数型どうしの割り算は切り捨てによって整数型になってしまう． これは Python 3
では異なる． 面倒なので，実数計算の際には 1, 2, 3 は 1., 2., 3. のように，小数点を付けて
書くのがよい．
• べき乗には ^ ではなく ** を使う． ^ は違う意味に使われる．
• 式やコマンドは左端に空白を入れないで書く．Python の文法では，行の始まりを右に下げ
ることでブロックの始まりを示し，左に戻すことでブロックの終わりを示す．これはプログ
ラム中で構造を表現するのに使われるので，勝手に空白をいれてはいけない．
Variableble explorer では変数の値を見ることができる． Console に
a=sin(0.1)
と打って，Variable explorer を見ると，a の値が表示されているだろう． さらに，
del(a)
と打てば，a が消去されるので，次に a と打っても “NameError” になる． del は delete のこと．
上向き矢印キー ⟨ ↑ ⟩ で前のコマンドを再現し，⟨enter⟩ で再実行することができる． これは「ヒ
ストリー機能」である．
IPython console という別の種類のコンソールを使うこともできる．それには，Console ウィン
ドウの下側のタブで IPython console を選び，さらに，すぐ上のグレー領域を右クリックし，Open
an IPython console をクリックする．普通のコンソールではグラフは別ウィンドウに表示される
が，IPython console ではグラフがウィンドウ内に表示される．
8 http://www.python.jp/doc/release/
9 http://docs.scipy.org/doc/
10 spyder
を準備するには，Linux ではパッケージシステムを使って spyder を導入する．Windows と Mac では，
Anaconda を web で探して導入する． 簡単．
11
データ入力
5.2
Editor は python 言語のスクリプトやデータファイルを書く場所である． ここにデータを入れ
て保存しよう． 何か書いてあるのは消してしまい，以下のようにデータを書き込む．これは練習
用データである．
# レーザーの偏光 4 月 1 日
# 回転角 (度) : 透過光強度 (ボルト)
0
30
60
0.809
0.403
0.039
90
120
0.077
0.477
150
180
210
0.836
0.794
0.400
240
270
0.040
0.081
300
330
360
0.481
0.843
0.782
# を付けて注意書きを入れることができる．File メニューから，Save as... で保存する．このと
き，どのホルダーに保存するかをよく考えて選び，必要なら保存場所を新しく作る．ファイル名に
は英語アルファベットでも数字でもひらがなでも漢字でもなんでも使える． Windows のプログ
ラムで開く可能性がある場合は ’.txt’ という拡張子を付けるが，そうでないなら，拡張子はいら
ない．
なお，測定点は規則的に並んでいなくても問題ない．等間隔になっていなくてもよいし，順番に
並んでいる必要もない． 同じ測定点で複数回測定してデータを追加してもかまわない．
5.3
データ読み込み
ここから先はスクリプト作成のところまで手動操作が続くが，スクリプトができてしまえば，そ
れを使って簡単にデータ処理をおこなうことができ，面倒な手動操作もスクリプトの内容も忘れる
ことができる． そこで，5.6 節のスクリプトの作成のところまで約２ページ分を飛ばしてもよい．
データを保存したら，Console でデータをファイルから読み込む．Console は「作業ディレクト
リ」を情報として持っていて，作業ディレクトリにあるファイルは ファイル名だけで指定するこ
とができる．それを見るには，Console で %pwd と打つ．
作業ディレクトリをファイルのあるディレクトリに合わせるには，ファイルを表示している Editor
ウィンドウのタブ (左上の出っ張り) を右クリックし，Set console working directory をクリックす
ればよい．
データファイルを読んでみよう． Console で下のように打つ．コマンド，関数，ファイル名など
は，入力の途中で ⟨tab⟩ キーを押せば残りが自動的に入るし，関数は括弧まで入れた時点で，その
説明が Object inspector に表示される． このような入力補助機能を利用して，記憶力を節約す
るのがよい．
12
data = loadtxt( "data01a.txt" )
読み込んだデータには data という変数名が付いたので，Variable explorer にその内容の一部が表
示されるが，全部を見るには，IPython console で data と打てばよい．data は 2 成分の行が 13 行
集まってできている．data の特定の要素を行番号と列番号で指定することができる．行と列はい
ずれも 0 番目から始まることになっているので，たとえば，data[0,1] は最初の行の右側の要素で
あり，data[12,0] は最後の行の左側の要素である．
行全体や列全体を表現することもできる．data[0]，data[1]，... data[12] は各行をあらわす．
data[0,:]，data[1,:]，... data[12,:] も各行をあらわす．data[:,0]，data[:,1] は各列をあらわす．左の
列は独立変数なので x とし，右の列は従属変数なので y とする．
x = data[:,0]
y = data[:,1]
これは，行と列を入れ替える「プロパティ」T を使って
x,y = data.T
としても同じである．
このデータをグラフにしよう．グラフを見てデータの間違いを発見することもある．
plot(x, y, ’o’)
ここで，’o’ はグラフの形式を指定するもので，小文字の o は丸でマークすることを示している．
データ点を線で結びたいなら，’o-’ を使うとよい．
モデル構築とフィッティング
5.4
(13) 式がモデル式であるので，それにもとづいて計画行列 G を作る．まず．
g1 = cos(2*x*pi/180)
g2 = sin(2*x*pi/180)
で G の左端の列とまんなかの列ができる．三角関数では角度はラジアンにしなければならないの
で，π/180 をかけている．3 番目の項は定数項なので，計画行列の 3 番目の列はすべて 1 にする．
すべて 1 のデータを作る方法はいろいろあるが，ここでは，既存のデータ x から次のようにして
作る．
g3 = x-x+1
x-x で x の中身をすべて 0 にしたものが作られ，それに 1 を加えることで，すべての要素が 1 に
なる．データ数は x と同じである 11 ．g1, g2, g3 を列ベクトルとみなして，それらを横につない
で G を作る．
G = column_stack((g1,g2,g3))
この column_stack 関数は 2 重の括弧を使うことに注意．
線形フィッティングの計算をまとめておこなってくれる便利な関数 ’lstsq()’ があるので，これを
使おう．最小二乗とは’least square’ なので，こういう名前になっている．
11 g3
= ones like(x) でも同じ．
13
A, B, C = lstsq( G, y )[0]
[0] が付いているのは，lstsq() が返す計算結果から必要なものだけを取り出すためである． A, B,
C と打てば，3 つのパラメータの推定値がわかる．
元の非線形モデルのパラメータは次のようにして得られる．
a = 2 * sqrt( A ** 2 + B ** 2 )
b = math.atan2( B, A ) * 90 / pi
c = C - a/2
atan2() に math. を付けなければならないのは，特殊なものとして別扱いになっているためである．
5.5
グラフ表示
パラメータ A, B, C が決まったので，モデルのグラフを表示することができる． 測定データは
0 から 360 まであったので，この間で 200 点ほど計算しよう．
XX = linspace(0, 360, 200)
YY = A*cos(pi/90*XX) + B*sin(pi/90*XX) + C*ones_like(XX)
XX には 0 と 360 の間の 200 個の値が等間隔に入る． YY は XX に対応する関数値が入る 12 ．こ
れを測定データと一緒にプロットする．
plot(XX, YY, x,y,’ro’,grid())
グラフの上にボタンがあれば，どれかを押して画像をセーブすることができる．あるいは，グラフ
を右クリックすると画像セーブができるかもしれない．
5.6
スクリプトの作成
同じことを繰り返さなくて済むように，データファイルを読み込んでフィッティングをおこなう
操作を固めて，
「スクリプト」を作る．
File メニューで New file... を選ぶ． Editor ウィンドウが新たに開くので，次の通り書き込ん
で， ’polarization.py’ という名前を付けて保存する．
12 GG = column stack((cos(pi/90*XX), sin(pi/90*XX), XX-XX+1)); YY = inner(GG,[A,B,C]) としてもよい．
inner は内積を計算する関数．
14
各行の左側のスペースの数に注意する必要がある．全角文字が紛れ込んでエラーになるので，
PDF からのコピー＆ペーストはだめ 13 ．
# -*- coding: utf-8 -*u"""
偏光測定用線形フィッティング・モジュール
モジュール読み込み: import polarization
使用法は fit クラスを見ること
"""
import pylab as P
class fit:
u""" 使用例:
r1=polarization.fit(’datafile1’)
--ファイルから読み込んで処理
r1.save()
--datafile1.png に結果を保存
"""
def __init__(me, filename):
x, y = P.loadtxt(filename).T
xr = P.pi/90*x #ラジアン
G = P.column_stack(( P.cos(xr), P.sin(xr), xr-xr+1 ))
results = P.lstsq(G, y)
A, B, C = results[0]
a = 2*P.sqrt(A**2 + B**2)
b = P.math.atan2( B, A ) * 90 / P.pi
c = C - a/2
me.xy, me.ABC, me.abc = (x, y), (A, B, C), (a, b, c)
me.name, me.results = filename, results
me.plot()
def plot(me, npoints=180, **kwopt):
(x, y), (A, B, C) = me.xy, me.ABC
X = P.linspace( x.min(), x.max(), npoints )
Xr = (P.pi/90) * X
#in dadian
Y = A*P.cos(Xr) + B*P.sin(Xr) + C*P.ones_like(Xr)
P.plot( x, y, ’ko’, X, Y, ’k-’, **kwopt )
P.grid( True )
P.title( me.__str__() )
P.show()
me.area, me.fig = P.gca(), P.gcf()
13 polarization.py
をダウンロードすることができるかどうかは担当者に尋ねること．
15
def __str__(me):
abc = ’a={:.5g}, b={:.5g}. c={:.5g}’.format(*me.abc)
return me.name + ’, ’ + abc
def save(me):
me.fig.savefig(me.name + ’.png’)
これは fit というクラスを定義しているモジュールであり，手動操作をなぞってはいるが，少し
余分な要素がある． たとえば，pylab という外部のモジュールを使用することを import 文で指定
している． また，数学関数やグラフを作成する関数には ’P.’ という修飾子が付いている．複数の
データの同時処理ができるような考慮もしてある．パラメータ値とファイル名がグラフの中に入る．
polarization モジュールを使うのは大変簡単である．spyder の Console で以下のように操作
する．
import polarization
f = polarization.fit("data01a.txt")
f.save()
f.save() でグラフがファイルに保存される． ファイル名はデータファイルの名前から作られる．
data_0505-1415, a=0.87823, b=-17.615. c=-0.00039712
0.8
0.6
0.4
0.2
0.0
0
6
50
100
150
200
250
300
350
400
線形フィッティング – Excel 編
Excel のような表計算プログラムでは，LINEST という関数で線形フィッティングをおこなうこ
とができる． 表を埋めたりグラフを表示するのに手間がかかるが，結果を着実に得ることができる．
次の図のように，測定値として，x と y の値を縦向きに並べる． ここでは三角関数を使うの
で，x はラジアンに直している． その右にモデル行列を入れる． ここでは「パラメータ係数」と
称している．定数項はモデル行列ではすべて 1 であるが，表計算の場合はそれを入れる必要は無
い．
「フィットパラメータ」の２行下の３つのセルをセレクトして，そこに LINEST 関数を入れる．
16
=LINEST(C3:C15,D3:E15,1)
LINEST 関数には，y の範囲とモデル行列の範囲を指定する． 定数項の計算もさせるために，3 番
目のパラメータを１にする． この式を配列式として確定するため，control キーと shift キーを押
しながら enter する．パラメータの推定値が表示されるが，定数項以外の順序がなぜか逆になる．
推定曲線は推定値を使って別に計算して表示しなければならない． なぜなら，測定点だけで計
算しても滑らかな線にはならないからである． 図では，表のこの部分は上の方の一部だけが見え
ている．グラフは「散布図」として作成する． ２系列のデータを一つのグラフに入れる方法を調
べるには，とにかくいろいろ試してみればよい．
F
G
H
I
J
K
0.5
C
C
0.45
)
0.4
0.35
0.3
0.25
(
A
B
C
D
E
1
2 degree x
y
cos
sin
0.4165
1
3
30 0.524 0.2264
0.5 0.866
4
60 1.047 0.0318
-0.5 0.866
5
90 1.571 0.0432
-1 1E-016
6
120 2.094 0.2223
-0.5 -0.866
7
150
2.618
0.4315
0.5
-0.866
8
180 3.142 0.4371
1 -2E-016
9
210 3.665 0.2548
0.5 0.866
10
240 4.189 0.0408
-0.5 0.866
11
270 4.712 0.0215
-1 4E-016
12
13
300 5.236 0.2149
-0.5 -0.866
330 5.76 0.4266
0.5 -0.866
14
360 6.283 0.4128
1 -5E-016
15
16
17
B
A
C
18
-0.11 0.1982 0.2294
19
20
y
21 degree x
0.4275
22
6 0.105 0.401
23
24
12 0.209 0.3669
18 0.314 0.3268
25
24 0.419 0.2824
26
30 0.524 0.2358
27
0.2
0.15
0.1
0.05
0
0
50
100
150
200
250
300
350
400
degree
非線形フィッティング – mathematica 編
7
数学ソフトの Mathematica を使って非線形フィッティングをおこなってみよう． どういう計算
をしているかがわかるようにする． 答えは一意的に求まる．
まずは Mathematica を起動して，ノートブックを表示する． 起動法や基本的な操作法は附録に
書いたので，適宜参照してほしい．
7.1
データ入力
データを次のように入力する．ノートブックで，dat= とタイプし，続けて，挿入メニューの表・
行列 の 新規作成で，ダイアログを出し，
「列数」を 2 にして OK ボタンを押す． 測定値を下のよ
うに入力する． 左側は回転角度，右側は透過光強度である． 行を追加するには，入力済みデータ
の一番下の行で，⟨control⟩-⟨enter⟩ を打つ．
17
0
0.809
30
60
0.403
0.039
90
120
150
0.077
0.477
0.836
dat= 180
210
0.794 ;
0.400
240
270
300
0.040
0.081
0.481
330
360
0.843
0.782
空白データが無いようにする． データがすべて入ったら，表の右側にセミコロンを付けて，⟨shift⟩-
⟨enter⟩ を打つ． 間違えたらデータを修正してから ⟨shift⟩-⟨enter⟩ を打てばよい．dat とすると，内
容を確認することができる． 回転角度と透過光強度からなる複数のリストが全体としてひとつの
リストにまとめられている．
Mathematica では，何か動作をさせるには ⟨shift⟩-⟨enter⟩ を打つ必要があることを忘れないこと．
次に，今入力したデータをグラフ表示する．
ListPlot[dat]
0.8
0.6
0.4
0.2
50
100
150
200
250
300
350
念のため，これから使う名前を洗浄する． 何か変なことが起こるとき，名前の洗浄をしなかっ
たのが原因ということは多い．
Clear[x,a,b,c]
7.2
モデルの入力
モデル関数を定義する．
f[ x , a , b , c ] := a Cos[(Pi/180)(x-b)]∧ 2 + c;
モデルに適当なパラメータを入れてみて，グラフでデータと比較する．
18
Show[Plot[f[x,0.8,160,0], {x,0,360}], ListPlot[dat]]
0.8
0.6
0.4
0.2
50
100
150
200
250
300
350
dat の中の角度と強度を別々のリストにする．
{X, Y}=Transpose[data];
残差の二乗和を関数にする．
J[a , b , c ] := Total[(Y-f[X,a,b,c])∧ 2]
7.3
フィッティングの実行
次のようにして，J(a, b, c) の最小点を求める． a, b, c の初期値として，上で推定した値を使う．
r = FindMinimum[J[a,b,c], {a,0.8}, {b,160}, {c,0}]
r から結果を取り出す．
s = r[[2]]
グラフを見る．
Show[Plot[f[x,a,b,c]/.s, {x,0,360}], ListPlot[dat], Frame->True]
0.8
0.6
0.4
0.2
0.0
0
50
100
150
200
250
300
350
もう一度結果を見る．
s
ファイルメニューの「保存」でノートブックを保存する． ただしこのファイルは Mathematica
が無いと役に立たない．グラフを右クリックし，
「形式を選択してグラフィックスを保存...」でグ
19
ラフを保存することができる． ファイルの種類はワープロ用なら PNG がよい．用途によっては，
PDF や EPS がよいかもしれない．
20
A
Mathematica の基礎
以下は Mathematica を使う場合に最低限必要な情報である．
A.1
起動
Mathematica を起動する方法はいろいろである． メニューやアイコンから起動できればそれ
でよい． 「端末」や「コマンドプロンプト」で Mathematica とタイプして起動することもできる
かもしれない．
起動の途中でソフトウェアライセンス条項に対する同意を求める画面が出た場合は，それを受け
入れる必要がある．
起動したら，白紙のノートブックが出れば準備完了であるが，Mathematica のバージョンによっ
ては，メニューかボタンを使ってノートブックを開く必要があるかもしれない．
メニューなどを日本語にするには，おそらく Edit メニューの Preference で設定することがで
きる． この場合，Mathematica を起動し直す必要がある．
A.2
式の入力
ノートブックが無事に開いたとしよう． 式を打ち込んだら ⟨shift⟩-⟨enter⟩ で計算ができる． 数
値の四則計算をしてみてほしい． かけ算は 3 4 のようにしてもよいし，3*4 のようにしてもよい．
割り算は 10/5 のようにする．整数どうしの割り算では結果が分数になる場合があるが，小数点
を付けた数の割り算では近似値が求められる． たとえば，10/12 は
5
6
となるが，10.0/12.0 は
0.833333 のようになる． 近似値を求めるには，式の右側に //N を付けることもできる． たと
えば，10/12//N とすればよい．
Mathematica に組み込まれている関数の名前はすべて大文字で始まることになっている． たと
えば，指数関数は Exp である． 引数は中括弧 [ , ] で囲む． よって，exp(x) は Exp[x] と
あらわす．正弦関数は Sin である． 角度の単はラジアンである． 円周率 π は Pi とあらわす．
よって，たとえば，角度 15 度の正弦を求めるには，Sin[15 Pi/180] とする． 近似値でよい場
合は，Sin[15 Pi/180]//N，あるいは Sin[15.0 Pi/180] のようにすればよい．
A.3
単純な定義
式に名前を付けると，その後は式のかわりに名前を使うことができる． たとえば，x としても，
41
x になるだけだが，x = 41/263; としたあとは，x とすれば， 236
になるし，x//N とすれば，
0.155894 と計算してくれる．一度使った名前を別の式に付けることもできる． この場合，前の式
は忘れられる．x が不要になったら，Clear[x] として記憶を消すのがよい．
名前は英語アルファベットやギリシャ文字を任意に組み合わせて勝手に作ることができる． 大
文字と小文字の区別がある． Mathematica にはじめからある名前は，大文字からはじまることに
なっている．
21
A.4
関数定義
f (x) = x2 − 2x − 1 というような関数を作るには次のようにする (この場合も ⟨shift⟩-⟨enter⟩ を
打つ)．
f[x_ ] := x^2-2x-1;
x_ という表現は「パターン」といい，
「ここに何かが来る． 何が来るかわからないが，それを
仮に x と呼ぶ」という意味である． 右辺は，仮に付けられた名前 x を使って，どのような計算
をするのかを指定する．
これは関数の定義である． 計算の実行は「関数呼び出し」といい，関数名と実際の計算対象を
指定する． たとえは，今定義した関数 f を使って，f[Pi] とすると，−1 − 2π + π 2 となるし，
√
√
f[2^(1/2)] とすると，1 − 2 2 となる (2^(1/2) は 21/2 , すなわち 2)．
A.5
リスト
複数のものをまとめてひとつにしたものを「リスト」という．ひとつにまとめたいものを中括弧に
入れてあらわす． たとえば，{1, 1/2, Pi} は，整数 1, 有理数 21 ，および無理数 π からなるリス
トである． リストはまとめて計算することができる． たとえば，{1, 1/2, Pi} 2 は，{2, 1, 2π}
となる．
a={1, 1/2, Pi}; としてから，a とすると，リストが表示される．2 番目の要素だけを見るに
は，a[[2]] とする． 2 重の大括弧が必要である．
複数のリストからなるリストを作ることもできる．
b={{4/3, 7/3}, {Pi, 2 Pi}, {2^(1/3), 2^(1/2)}};
22
© Copyright 2024 ExpyDoc
