カーブフィッティング ( 回帰分析 ) の練習

カーブフィッティング ( ^回帰分析 ) ^の練習

線形最小二乗フィッティング

: R, python, excel

非線形最小二乗フィッティング

: mathematica

目次

1

^{最小二乗法}

2

2

線形フィッティング

(

線形回帰

) 3

2.1

線形フィティングの要件

. . . 3

2.2

線形最小二乗法の理論の概要

. . . 4

2.3

定数項

. . . 5

3

^{直線偏光測定のモデル}

5 3.1

線形モデルへの変換

. . . 6

3.2

特別な場合（離散フーリエ変換）

. . . 6

4

^{線形フィッティング}

–

^{統計解析言語}

R

^編

7 4.1 R

^{の作業ディレクトリ}

. . . 7

4.2 R

^の使い方

. . . 7

4.3

データ入力

. . . 9

4.4

フィッティング

. . . 10

4.5

グラフの作成

. . . 11

4.6

ひとまとめにする

. . . 11

5

^{線形フィッティング}

–

^{スクリプト言語}

python

^編

12 5.1

基本的操作法

. . . 12

5.2

^{スクリプト作成}

. . . 13

5.3

^{ライブラリーの使用}

. . . 14

5.4

^{計画行列作成}

. . . 14

5.5

^{データ入力と読み込み}

. . . 15

5.6

フィッティング

. . . 16

5.7

グラフの完成

. . . 17

5.8

解析をおこなう関数の作成

. . . 18

6

^{線形フィッティング}

– Excel

^編

19 7

^{非線形フィッティング}

– mathematica

^編

20 7.1

データ入力

. . . 20

7.2

モデルの入力

. . . 21

7.3

フィッティングの実行

. . . 22

付録

A Mathematica

^の基礎

23 A.1

起動

. . . 23

A.2

式の入力

. . . 23

A.3

^{単純な定義}

. . . 23

A.4

^関数定義

. . . 23

A.5

^リスト

. . . 24

測定可能な

2

つの量の間の関係が理論式であらわされていて，その式が未知のパラメータを含んでいる場 合，いろいろな条件で得た複数の測定データから未知パラメータの値を求めることができる． これをカーブ フィッティング

(curve fitting),

あるいは回帰分析

(regression analysis)

という．理論式はモデル式ともいう．

多くの場合，測定誤差のため，モデル式を測定データに完全に一致させることはできないので，できるだけ 多くの測定データを使ってパラメータの値を推定する． これはデータとモデルの「不一致の程度」が最小に なるようなパラメータの値を求めるということである． 「不一致の程度」としては，残差

(

^{測定データとモデル} の差

)

^{の二乗和が使われる}

^*1

．測定誤差がランダムに生じるものであるならば，残差の二乗和を最小にするよ うなパラメータ値が最も確からしいパラメータ値であるということが言えるからである

(

^{ガウス・マルコフの} 定理

)

．この考え方を最小二乗法という．

1 ^{最小二乗法}

モデル式が従属変数

y

を独立変数

x

の関数としてあらわすようなものであるとする．未知パラメータ

a

を含 めると，次のように，

y

は

x

と

a

の

2

変数関数として書ける．

y= f(x,a) (1)

n

組の測定値を

(x₁,y1), (x₂,y2), ..., (x_n,yn)

とあらわす．

i

番目の測定値の残差を

ei

であらわすと，

y1 = f(x₁,a) +e1

y2 = f(x₂,a) +e2

...

y_n = f(x_n,a) +e_n

(2)

残差の二乗和を

J

は，

J=e12+e22+...+en2 (3)

*1

空間の

2

点間の距離を座標の差の二乗和から求めるのと同じである．

e1, ...,en

をひとまとめにして

n

次元列ベクトル

e

にする．

J

はこのベクトルの大きさの

2

^{乗であり，}

e

と

e

の 内積でもある．

e=

















 e1

e2

... en



















, J=kek²=e·e (4)

パラメータ

a

を変えれば，

f(x_i,a)

の値が変わり，

J

も変化する．ただし，

J=0

であるから

J

の値には下限 がある．未知パラメータの推定とは，

a

の関数である

J

の最小点を求めることである．

パラメータが複数あっても同じである．

m

個のパラメータ

a1,a2, ...,am

があるなら，これらをまとめて一つ の

m

次元列ベクトル

a

とみなし，

a=

















 a1

a2

... a_m



















(5)

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) =a1x+a2

という形になってしまう が，この場合は直線回帰

(linear regression)

^という．

直線回帰は計算量が少ないので手計算でもなんとかなるし，統計機能を備えた電卓や表計算ソフトのグラフ 作成機能を使えば手軽におこなえる．一般の線形フィッティングでは，連立方程式を解かなければならないの で，その機能を備えたコンピュータープログラムを使用することが多い．非線形フィッティングには，複雑な アルゴリズムによって大量の計算をおこなうコンピュータープログラムを使う．

後で見るように，線形フィッティングは連立１次方程式を解くことに帰着するので，答えが得られるどうか は単純な条件で決まる．非線形フィッティングは試行錯誤によって

J

の最小点を探し出すので，正しい答えが 得られたかどうかを判断することが難しい場合がある．また，問題によって最適なアルゴリズムが異なるかも しれないし，使用するアルゴリズムによって答えが異なる．そのため，未知パラメータになんらかの置き換え をおこなうことで，線形な式にすることができるかどうかを検討するのが得策である．

2 ^{線形フィッティング} ( ^線形回帰 )

2.1 線形フィティングの要件

m

個の未知パラメータに対して線形のモデル式は，

f(x) =a1g1(x) +a2g2(x) +...+amgm(x) (6)

のように書ける．これは

m

個の関数

g1(x),g2(x), ...,gm(x)

を係数とする

m

個の未知パラメータ

a1,a2, ...,am

の線型結合である．

f(x1) =a1g1(x1) +a2g2(x1) +...+a_mg_m(x1)

などを式

(2)

^{に代入すると，}

y1 =a1g1(x1) +a2g2(x1) +...+amgm(x1) +e1

y2 =a1g1(x2) +a2g2(x2) +...+amgm(x2) +e2

...

yn =a1g1(x_n) +a2g2(x_n) +...+amgm(x_n) +en

(7)

g1(x₁)

などはモデル式の

m

個の関数にそれぞれ

n

個の測定点の値を代入して得られる

m×n

個の値である．

式

(7)

は

a1,a2, ...

を未知数とする線形連立方程式とほぼ同じであるから，ベクトルと行列を使って簡潔に表

現しよう．

y1, ...,yn

と

e1, ...,en

を

n

次元列ベクトル，

a1, ...,am

を

m

次元列ベクトル，

g1(x₁), ...,gm(x_n)

を

n

行

m

列の行列にする．

y=













 y1

... yn















, e=













 e1

... en















, a=













 a1

... am















, G=















g1(x1) · · · gm(x1)

... ... ...

g1(x_n) · · · gm(x_n)















(8)

すると，式

(7)

は次のようにまとめられる．

y=Ga+e (9)

G

は計画行列

(design matrix)

^{と呼ばれる．その}

i

行

j

列の成分は

G_{i j} =g_j(x_i)

である．

残差の二乗和

(J=e·e)

を最小にするという意味でこの方程式

(9)

を「解く」ことができる．これが線形最

小二乗法

(linear least-squares method)

で，その解を最小二乗解という．ただし，未知パラメータを決めるに

は，未知パラメータの数

m

よりも測定値の数

n

が多いことが必要である．よって，

G

は縦長の行列である．

2.2 線形最小二乗法の理論の概要

a

の微小変化

da

に対する

J

の変化

dJ

を計算しよう．式

(9)

^から

e = y−Ga

であるから，

e

の変化は

de=−Gda

となる．また，

J=e·e

であるから，

dJ= 2de·e

．よって，

dJ=−2(Gda)·(y−G a) =−2da·

G^∗(y−Ga)

となる．ただし，

^∗

で随伴行列をあらわす．

J

の極小点では

a

のどの向きへの微小変化に対しても

J

の変化は

0

であるから，

G^∗(y−Ga)

はゼロでなければならない

^*2

．よって，

G^∗Ga=G^∗y (10)

G^∗G

は

n×n

の半正定値対称行列である

^*3

．

(10)

が一意的に解けるには，

G^∗G

は正則でなければならない ので，

G

はフルランク

(

^階数が

m

に等しい，単射

)

でなければならない．これは

G

の各列が独立ということで あるが，式

(8)

を見ればわかるように，それには，少なくともモデル中の関数

g1(x), ...,gm(x)

が独立でなけれ ばならない．これはモデル式が冗長ではないということである． また，少なくとも

m

個の独立な測定点が必 要である．だから，同じ測定値をコピーして増やすのはだめ．

*2G^∗(y−Ga) = 0

の左辺は

m

次元ベクトルで，その各成分，G の

i

列目を

g_i

として，

g_i·(y−Ga)

となる．右辺は

0

ベクトルだか ら，この内積はすべて

0

である．これは，G の全縦ベクトルと

y−Ga

が垂直であることを意味する．よって，線形変換

G

の像で ある

m

次元超平面に点

y

から下ろした垂線の足が

Ga

である．

*3

内積の性質から

a·(G^∗G)a= (Ga)·(Ga)=0

である

(

半正定値

)

．また，(G

^∗G)^∗=G^∗(G^∗)^∗=G^∗G

である

(

対称

)

．

パラメータ数

m

が小さい場合は，式

(10)

^{を直接解いて}

a

を求めることができる．直線回帰はそういう場合 にあたる．

m

が大きい場合にも使われるコンピュータープログラムでは，計画行列を巧妙に処理して解くよ うになっている

^*4.

2.3 ^定数項

パラメータの

1

つが定数項になっている場合は，モデル式は

h(x,a,b) = f(x,a) +b

のように書ける．ここで，定数項はパラメータ

b

であり，他のパラメータは

a

，

f(x,a)

は

a

をパラメータとす る

x

の関数である．

J=Pn

i=1

y_i−f(x_i,a)−b2

を最小にする

b

は容易に求められ，

b=1 n

n

X

i=1

(y_i−f(x_i,a))

となる．そこで，平均値

(y1 +...+yn)/n

を

y

とし，平均値

f(x1,a) +...+ f(xn,a)

/n

を

f(x,a)

とし，

Yi=yi−y,F(xi,a) =h(xi,a)−f(x,a)

と置くと，

yi= f(x,a,b) +ei

は

Yi=F(xi,a) +ei

となって，未知パラメータが

1

つ減る．線形モデルでは，

f(x,a)

は

f(x,a) =a1g1(x) +· · ·+amgm(x)

のように計算することができる．よって，まず

a

を求め，その結果を使って

b

の値を計算すればよい．

3 ^{直線偏光測定のモデル}

偏光フィルターは，入射光の偏光の向きに対して回転させることによって光の透過率が変化するものであ る．理想的な偏光フィルターでは，回転角度

x

に対する透過率

t

の関係は次の式であらわされる．

t=Pcos²(x−b) +P⁰ (11)

b

は偏光フィルターへの入射光の「最大偏光方向」の角度である．係数

P

と

P⁰

は入射光の偏光状態で決ま る

^*5

．実際には偏光フィルターには偏光によらない反射や吸収があるし，直接測定するのは偏光フィルターを 通過した光の強さを光センサーで何らかの量に変換したものであるが，これらの過程が線形性を有し，入射光 強度が一定である場合は，測定値

y

は

y=acos²(x−b) +c (12)

となる．パラメータの

a

と

c

は入射光強度やセンサーの感度に依存するし，

c

はセンサーのゼロ点のずれを含 むかもしれない．しかし，

b

にはそのような依存性は無いと考えられる．

*4QR

分解や特異値分解が利用される．

QR

分解はグラム・シュミットの直交化のアイデアを行列の分解に適用するもので，直交行 列

Q

と上三角行列

R

によって

G=QR

と分解され，a

=R⁻¹Q^∗y

となる．

特異値分解では，行列

G

は２つの直交行列

U,V

と対角行列

S

によって

G=USV^∗

と分解される．ただし

S

は

G^∗G

のすべて の固有値の平方根からなる対角行列である．S のゼロでない成分をすべて逆数に変えたものを

S⁺

とし，G

⁺ =VS⁺U^∗

とする．

G⁺

は

G

の「擬似逆行列」という．G が縦長でフルランクの場合，式

(10)

の解は

a=G⁺y．

*5

形式的に，強さ

S

の光を，強さ

S1

の直線偏光成分と，それに独立な強さ

S2

の成分とに分けて考えることができ，式

(11)

の係数

P

と

P⁰

は，P

=S₁/S

，P

⁰=S₂/(2S) = (1−P)/2

とあらわせる．

3.1 ^{線形モデルへの変換}

式

(12)

は

a

と

c

について線形であるが，

b

は非線形な関数

cos

の中に入っているので，

b

については非線形 である． そこで，パラメータの変換によって線形の式にすることができるかどうかを検討する．倍角公式と 加法定理を使って式

(12)

を変形すると，

y=acos2b

2 cos2x+asin2b

2 sin2x+a 2 +c

よって，次のように，パラメータを

A, B,C

に置き換えれば線形になる．

y=Acos2x+Bsin2x+C, (13)

A= a

2cos2b, B=a

2sin2b, C=a

2 +c (14)

a= 2p

A²+B², b=1

2arctan(A,B), c=C−a

2 (15)

ここで，

arctan(A,B)

は平面のベクトル

(A,B)

の方位角である． プログラミング言語や数値計算パッケージ

の多くはこの関数が計算できる

^*6

．普通の

arctan

^{を使う場合は，}

A>0

なら

arctan(B/A),

^{それ以外で}

B<0

なら

−(π/2)−arctan(A/B),

それ以外なら

(π/2)−arctan(A/B)

とする．

3.2 特別な場合（離散フーリエ変換）

式

(13)

^は

y

を周期

π

の

2

つの成分と定数成分に分解した形になっているから，

A, B,C

は「周波数成分」で ある．特に，独立変数

x

の値が正確に全周期にわたって等間隔に取られた場合には，次の離散フーリエ変換に よってこれらを求めることができる．

A= 2 N

m

X

i=1

yicos(2xi), B= 2 N

m

X

i=1

yisin(2xi), C= 1 N

m

X

i=1

yi (16)

xi,yi

はデータ点，

m

はデータ点数である．

xi

は

2Nπ(i+h)/m

の形

(N

は偏光フィルターの回転数，

h

は任意 の定数

)

に書けなければならない．偏光フィルターの回転数は

1/2

の整数倍ならよい．離散フーリエ変換が直 交変換であることから，この計算が最小二乗法になっていることが言える

^*7

．また，線形最小二乗法の式

(10)

にモデルを代入して計算すれば，上式が直接得られる．

一般の線形フィッティングはコンピュータープログラムで簡単にできるので，わざわざ

x

の取り方に制約を 付けて離散フーリエ変換で計算するメリットはあまり無い．

*6C

言語など多くの言語にある

atan2

という関数は

atan2(B,A)

のように使用する．

Y

座標が先，

X

座標が後である．

Excel

や 計算サイト

http://keisan.casio.jp

の

atan2

や，Mathematica の

ArcTan

は

X

座標が先．

*7

「パーセバルの等式」が成り立つ．

4 ^{線形フィッティング} – ^{統計解析言語} R ^編

R

言語は統計解析に特化したシステムなので，データ・フィッティングは比較的楽にできる．

4.1 R の作業ディレクトリ

Windows

ではスタートメニューやタイルで

R

を起動する．

Linux

ではメニューやデスクトップ・アイコン

で起動するか，端末エミュレーターの中で

R

と打って起動する．

R

を終了するには，

R

の中で

quit()

^と打 つ．このとき，作業内容を保存するかどうかを答える必要がある．

ファイルの読み書きのために「作業ディレクトリ」が重要である．端末エミュレーターで

R

を起動する場合 は

R

を起動したときの端末の現在のディレクトリが

R

の作業ディレクトリになる．これは

R

^{を起動する前に}

pwd

コマンドで見ることができる．また，

cd

コマンドによって変更することもできる．しかし，

R

^を

GUI

^で 起動すると，作業ディレクトリがわからないかもしれない．

R

では次のコマンドで作業ディレクトリを確認す ることができる．

1 getwd()

これは

”getwd”

という名前の関数を引数無しで呼び出す式であり，この式の評価結果が表示される．

getwd

の

wd

は

”working directory”

から来ている．また，

dir()

で作業ディレクトリにあるファイルとサブディレ

クトリのリストが得られる．サブディレクトリのみのリストは

list.dirs(recursive=0)

で得られる．作業 ディレクトリを変更する必要がある場合は

setdir

関数を使う．たとえば，現在の作業ディレクトリの中にあ る

xxxx

という名前のサブディレクトリに移るには，

setwd("xxxx")

とすればよいし，現在のディレクトリ の親ディレクトリに移るには，

setwd("..")

とすればよい．これらは現在の作業ディレクトリからの相対パ スである．ユーザーのホームディレクトリに移るには，

setwd("~")

^とする．

Windows

^{では，以上のような} 操作を

GUI

でもおこなうことができる．

4.2 R ^の使い方

R

では簡単な数値計算がすぐにできる

^*8

．四則演算の演算子は

+,-,*,/

^{である．べき乗は}

^

^でも

**

^でも よい．式をくくる括弧にはすべて

( )

を使う．関数呼び出しはすべての引数を括弧で囲んで関数名の後ろに 付ける．

1 (1+2)/2

2 sin(0.1)

3 tan(0.1)

複数の数値からなるデータは

c

という関数を使って作る．たとえば，

1 v <- c(1,3,5,4)

*8

これ以降は，4.3 の「データ入力」と

4.6

の「ひとまとめにする」のセクション以外は読むのを省略することもできるが，推奨はし

ない．

2 v

3 v+10

4 v*10

のようにすると，４つの数値からなるデータを作って

v

^{という名前を付け，}

v

の各要素に数値を足したりかけ たりした結果を表示する．このようなデータは「ベクトル」と呼んでいる．たくさんの要素からなるベクトル を簡単に作る方法もいろいろ用意されている．重要なのは等差数列である．コロン

”:”

を使って，たとえば

1 0:40

と書くと，

0

から

40

までの公差

1

の等差数列

(0,1, ...,40)

ができる．次のように，このベクトルに

0.5

をかけ ると，

(0,0.5, ...,20)

になり，さらに

10

^{を引くと，}

(−10,−9.5, ...,10)

になる．

1 x <- 0:40*0.5-10

2 x

同じことは，

seq

という関数を使ってもできる．

1 x <- seq(-10,10,0.5)

2 x

seq

はたぶん

sequence

から来ている

^*9

．

数値データをプロットする関数は

plot

^である．

y=x²

のグラフを描いてみる．

1 x <- seq(-10,10)

2 y <- x^2

3 plot(x,y)

21

点しか計算していないので，表示もそれだけである．

「任意の数値に対してその２乗を計算する」といった計算方法自体は「関数」である．関数を定義してそれ に名前を付け，後でその関数を使って計算をおこなうことができる．

1 f <- function(x){x^2}

2 f(2)

3 f(10)

4 f(2)+f(10)

5 f(0:10)

関数定義のところは，

*9

関数の説明を見るには，

?seq

のようにすればよい．

Unix

では，カーソルキーで表示をスクロールし

Q

で元に戻る．

1 function(

^仮引数

){

^{仮引数を含んだ式}

}

という形式である．中括弧の中には関数の定義が入る．簡単な計算では結果をあらわす式だけでよいが，もっ と複雑な処理の記述もできる．

関数定義を使うと，

plot

関数によって曲線も描くことができるが，左端と右端を指定する必要がある．

1 plot(f, -10, 10)

これにデータ点を重ねるには，

points

^{関数を使う．}

1 x <- seq(-10,10,2)

2 points(x, f(x))

4.3 ^{データ入力}

測定データは任意のテキストエディタを使って入力して適当なファイル名で保存する．保存する場所は

R

の作業ディレクトリにする．次のは練習用データである．

#

はコメント行を示している．

1 #

^{レーザーの偏光}

4

^月

1

^日

2 #

^回転角

(

^度

) :

^{透過光強度}

(

^ボルト

)

3 0 0.809

4 30 0.403

5 60 0.039

6 90 0.077

7 120 0.477

8 150 0.836

9 180 0.794

10 210 0.400

11 240 0.040

12 270 0.081

13 300 0.481

14 330 0.843

15 360 0.782

このファイルを

R

で読み込む．ファイル名を

xxyy.dat

とすれば，

1 data = read.table("xxyy.dat")

読み込んだデータに

data

という名前を付けたが， 日本語のほうがよければ， 生データ

=read.table("xxyy.dat")

のようにしてもよい．測定データのプロットは次のようにすればすぐにできる．

1 plot(data)

データ・プロットをよく見れば，誤った測定点が見つかる可能性がある．

1 data

で

data

^{を表示すると，}

1

列目に回転角，２列目に透過光強度が入っていることがわかる．

data

^{から各列の} データを取り出すには，次のように２重の大括弧を使う．

1 x = data[[1]]

2 y = data[[2]]

回転角が独立変数なので，それに

x

という名前を付け， 透過光強度には

y

という名前を付けた．

4.4 ^{フィッティング}

あらためて，データプロットを作っておく．今度は横軸とグラフ全体に名前を付けよう．

1 plot(x,y,xlab="

^角度

(

^度

)",main="

^{テストデータ}

")

角度の単位を度にしたので，モデル式

(13)

は，

Acos

(π/90)x

+Bsin

(π/90)x +C

という形になる．線形フィッティングには

lm

関数を使う．

1 r <- lm( y ~ cos(x*pi/90)+sin(x*pi/90) )

lm

に指定するのは，従属変数

y

と

~

とモデル式から未知パラメータを取ってしまったものである．名前

r

^を 使って計算結果を見ることができる．

lm

関数が計算したさまざまな結果全体を見るには，

summary(r)

とすればよいが，未知パラメータの値だ け取り出すには次のようにする．

1 p = r$coefficients

2 p

p

の各要素に対応しているモデル式の項が表示され，定数項

C

がはじめに来て，次に

cos

の係数

A

，

sin

の係 数

B

の順になっていることがわかるので，

1 C <- p[[1]]

2 A <- p[[2]]

3 B <- p[[3]]

式

(15)

を使って

a,b,c

の計算もやっておく．

R

では

atan2

関数が使える．

1 a <- 2*sqrt(A^2+B^2)

2 batan2(B,A)*(180/pi)

3 c <- C-a/2

4 a;b;c

4.5 グラフの作成

理論曲線を描くために，フィッティングの結果を入れて関数を定義する．

(12)

式に戻ろう．

1 f <- function(x) a*cos((x-b)*pi/180)^2+c

角度は度からラジアンに直した．準備ができたので，データ点と理論曲線を描く．

1 plot(f,0,360,xlab="

^角度

(

^度

)",main="

テストデータのフィッティング

")

2 points(x,y)

3 grid()

始めの

plot

では少し飾りを付けた．また，

grid()

で目盛線を入れている．

グラフは

PNG

か

JPEG

ファイルに保存するのが簡単である．保存してから印刷したりレポートに入れた りすればよい．

1 savePlot("xxxxyyyy.png","png")

4.6 ^{ひとまとめにする}

データの読み込みからグラフの保存までを１つの関数に入れ，

analyze

という名前をつけよう．この関数に は，データファイルの名前を引数として与えることにする．下の関数定義を

PDF

^から

R

^{へコピー＆ペース} トしても動くかもしれない

^*10

．

1 analyze <- function(datafile) {

2 data <- read.table(datafile) #

^{データ読み込み}

3 x <- data[[1]]; y <- data[[2]] # x

^列，

y

^列

4 r = lm(y ~ cos(x*pi/90)+sin(x*pi/90)) #

^{フィッティング}

5 p <- r$coefficients # $

^{で未知係数の値を得る}

6 C <- p[[1]]; A <- p[[2]]; B <- p[[3]]

7 a <- 2*sqrt(A^2+B^2)

8 b <- 0.5*atan2(B,A)*(180/pi)

9 c <- C-a/2

*10

これが入っているファイルを入手することができた場合は，R で一度だけ

source(”

ファイル名”) のようにすればよい．

10 print(c(a=a,b=b,c=c)) #

^{結果をまとめて表示}

11 #

^{グラフを描く}

12 h <- max(y)*1.1-min(y)*0.1 #

^上に

10%

^余裕

13 l <- min(y)*1.1-max(y)*0.1 #

^下に

10%

^余裕

14 plot(function(x) a*cos((x-b)*pi/180)^2+c, 0, 360,

15 xlab="

^角度

(

^度

)", ylab="

^強度

(V)", ylim=c(l,h),

16 main=paste(datafile, ", b=", b) )

17 points(x, y)

18 grid()

19 savePlot(paste(datafile,"png", sep="."), type="png")

20 }

これで，次のようにするだけで，パラメータ値を計算して表示し，グラフを画像ファイルに保存してくれる．

1 analyze("xxyy.dat")

R

を終了するときに，

”Save workspace image? [y/n/c]: ”

と訊かれたら，

y

で答える．すると，次に同じ ディレクトリで

R

を起動すると，前におこなった定義がそのまま使える．この情報は

”.RData”

という隠し ファイルに保存されている．何もかも終わったら，このファイルを消去する．

5 ^{線形フィッティング} – ^{スクリプト言語} python ^編

データ解析をプログラミング言語でおこなうにはいろいろな選択肢がある．豊富な科学計算機能と素早いテ ストができることと正確な情報が得やすいことから

Python

を使うことにしよう．

Python

のデータ処理ライ

ブラリの

1

つである

pandas^*11

と

pylab^*12

を使って線形フィッティングをやってみる。

Python

の文法の理

解はいらない

^*13

。

Python

を実行する環境はいろいろあるが，

spyder

というユーザー・インタフェース を使 うことにする

^{*14 *15}

。

5.1 基本的操作法

spyder

を起動するには，たとえば端末に次のようにタイプする．

1 spyder

あるいは，メニューから

spyder

を起動できるかもしれない．

*11http://pandas.pydata.org/

*12http://docs.scipy.org/doc/

*13http://www.python.jp/doc/release/

*14spyder

を準備するには，Linux ではパッケージシステムを使って

spyder

を導入する．Windows と

Mac

では，Anaconda を

web

で探して導入する． 簡単．

*15spyder

を使わず，

15

ページのデータファイルの形式を見，スクリプトファイル

henkou.py

をダウンロードし，

19

ページの脚注

の最後の３行を見れば，解析だけはできる．

spyder

のウィンドウは大きく３つの枠に分かれているが，そのうち

IPython console

^{に式やコマンドを} 打ち込めば，その場で計算の実行ができる．

1+2

^とか

2*3

^とか

12/-4

などとやってみれば，ただの電卓とし ての使い方がわかるが，注意すべき点が少しある．

• 1, 2, 3

^は

1.0, 2.0, 3.0

とは異なる． 前者は「整数型」後者は「実数型」である．

Python 2

^では，整 数型どうしの割り算は小数部が切り捨てられてしまう．実数の結果が必要な場合は，整数でも

2.0

とか

2.

のように小数点を付けて書く．

•

べき乗には

^

^ではなく

**

^を使う．

^

は違う意味に使われる．

2×10⁹

のようなものは

2*10**9

^のよ うに書くよりも，

2e9

のように書く方がよい．

•

行の左の空白（インデント）の量は文法上意味があり，

if

文のような複合文のブロックの始まりで空白 を増やし終わりで元に戻す．

• log(2.718)

のように，関数呼び出しはいつも小括弧を付けて書く．

Variable explorer

で は 変 数 の 値 を 見 る こ と が で き る．

IPython console

^に

a=b=0.1

^{と 打 っ て，}

Variable explorer

^{を見ると，}

a

^と

b

の値が表示されているのがわかる．また，

del a, b

^{と打てば，}

a

と

b

^{が消去される．}

上向き矢印キー

h

↑

i

で前のコマンドを再現し，修正してから

henteri

キーで再実行することができる． こ れは「ヒストリー機能」である．

5.2 ^{スクリプト作成}

Editor

ウィンドウにはプログラムなどを書いて保存することができる．まず，

File

^{メニューの}

New file...

で新しい

Editor

ウィンドウを開く．すでに開いている

Editor

ウィンドウは閉じてしまってよい．

Spyder

では，

Editor

ウィンドウに書いたプログラムのどの部分でも任意に実行してテストすることがで

きる．たとえば，

hF9i

を押すと，あらかじめ選択した部分が

IPython console

で実行される．もう１つの 方法は，

hCtrlihEnteri

，または

Run current cell

である．こちらはマークによって区切られた範囲を実行 するものである．便利なので練習しておこう．

Editor

ウィンドウ中に次のように書く

(

初めから書いてある ものは置いておく． 左端から書く． 文字列は半角のシングルクォート

’

かダブルクォート

”

で囲む

)

．

1 #%%

2 a='

^たら

'

3 print a

4 #%%

5 a=a+a

6 print a

#%%

^の行は，

Editor

ウィンドウを「セル」という部分に区切るための記号である．今書いたものは

2

つのセ ルに分かれている．境界は横線でわかる．上側のセルの中のどこかにカーソルを置いておき，

Run

メニューで

Run current cell

を実行する．するとそのセルの中身が

IPython console

で実行される．

Run current cell

は

hCtrlihEnteri

を押してもできるし，ツールバーにもこれをおこなうボタンがある．

上側のセルと下側のセルをいろいろな順番で実行してみてほしい．試し終わったら，この意味のないプログラ ムは区切り記号と一緒に消してしまう．

5.3 ^{ライブラリーの使用}

Python

は汎用プログラミング言語なので，原子的なもの以外の機能はプログラムモジュール

(

ライブラ

リー

)

として提供され，必要なモジュールを指定して使用する．数値計算や図表示のさまざまな機能を使用す

るには，

pylab

モジュールが適している．また，

pandas

モジュールを併用するともう少し便利になる．これら

のモジュールを使用するため，

Editor

ウィンドウに次のコードを書き，

Run current cell

^{で実行し，エラー} にならないことを確認しておく

^*16

．

1 from pylab import *

2 from pandas import *

3 #%%

ここでファイル名を付けて保存しておく．ファイル名は

henkou.py

のように適当でよいが，拡張子として

.py

を付けておく．何か書き足したときには，忘れないで保存をおこなう．

なお，コマンド，関数，ファイル名などは，入力の途中で

htabi

キーを押せば残りが自動的に入る．また，

関数名を入力した直後に

hcontroli

キーと

I

^{を押すと説明が}

Object inspector

に表示される．関数呼び出 しを括弧まで書くだけで自動的に説明があらわれるようになっている場合もある．このような入力補助機能を 利用して，記憶力を節約するのがよい．

5.4 ^{計画行列作成}

線形モデルの式

(13)

^{を使って式}

(8)

^{の計画行列}

G

を構成するには，独立変数

x

の値リスト

x1,x2, ...,x_n

が 必要である．

python

では，単純なリストはたとえば

[1, 2, 3]

^{のように表現する}

^*17

^{．しかし，このリスト} 構造は一般的すぎて扱いが面倒なので，実数のベクトルや行列を効率よく扱うことができるリッチな機能を備 えた

array

や

DataField

というタイプのデータ構造を使う

^*18

．

array([1,2,3])

などとすれば変換することがで きる．そこで，

[1,2,3]*2

^の結果と

array([1,2,3])*2

の結果の違いを見てほしい．用途の違いがわかる．

x

のリストを与えれば計画行列を作ってくれる関数を定義しておこう

^*19

．

Editor

ウィンドウに次のコード を入力する

(

先ほどの区切り記号の下に続けて

)

．

def ...

^{の行は左端から書く．}

2

行目以降には自動的に空 白が入る．これらの空白を変更してはいけない．

4

行目で改行した後

5

行目にできる空白は消しておく．括弧 の種類に注意．関数名

make_g

は

make

とアンダースコアと

g

である．

*16import *

はそれぞれのモジュールで定義されたすべてのものを読み込んでしまい，同じ名前のものがたまたま複数あると取り返

しがつかないので，一般的には推奨されない．ここではそのようなことが起きる可能性を無視し，手軽さを優先している．

*17

これは整数

1, 2, 3

を並べてひとまとめにしたものであり，

a=[1,2,3]

のようにして全体にひとつの名前を付けることができる．

その中のひとつずつの要素は，0 から始まる要素番号を使って，

a[0], a[1]

のように書き表す．

a=[[1,2],[3,4,5]]

のような リストのリストは要素を

a[0][0], a[1][2]

のようにあらわす．

*18array

は

numpy

モジュールによって，DataField は

pandas

モジュールによって定義されている．

numpy

モジュールは

pylab

モジュールの一部である．

*19

この程度の計算は実質的に

1

行のプログラムですんでしまうが，それでもわざわざ新しい関数を定義するのは，関数名を見ればど

んなことをするのかわかるからである

(子供の名前を見れば少なくとも親の思いはわかる．)．

1 def make_g(xx):

2 u"

測定点リストから計画行列を作る関数：定数項省略

"

3 rr = array(xx) * pi/90

4 return DataFrame( [ [cos(r),sin(r)] for r in rr ] )

5 #%%

カーソルをこの定義の中に入れて，

Run current cell

で実行する．この実行は

make_g

^{という名前の関数の} 定義

(

すなわちどのように計算すればいいかということの記憶

)

^を

python

^{が行っただけである}

^*20

^．

こ の 関 数 の 動 作 を テ ス ト す る た め，

IPython console

で， 次 の よ う に， 測 定 点 の リ ス ト と し て

[0,15,30,45]

^を

make_g

に与えて実行してみる．

1 make_g([0,15,30,45])

これで

4

^行

×2

列の配列ができる．これは次のような計算をしたのである．各行は各測定点に対応している．

















cos(2×0°) sin(2×0°) cos(2×15°) sin(2×15°) cos(2×30°) sin(2×30°) cos(2×45°) sin(2×45°)

















5.5 データ入力と読み込み

測定データはテキストファイルに入れておくのがよい．テキストファイルは作るのも変更するのも容易だか らである．

File

メニューの

New file

で新しい

Editor

ウィンドウを開き，ここに測定データを入れて保存し よう．何か書いてあるのを全部消してしまってから，次のように書き込む

(

これは例である

)

．はじめの行には

x

^と

y

を列名として書く．２行目から偏光フィルターの角度と透過光強度のデータを空白で区切って書き込ん でいく．各行の右側に

#

で注釈を入れてもよい．空行や注釈のみの行は許されないかもしれない．

1 x y #

^{スクロース溶液}

1, 5/16

2 0 0.809

3 30 0.403

4 60 0.039

5 90 0.077

6 120 0.477

7 150 0.836

8 180 0.794

9 210 0.400

*20

１行目のキーワード

def

は

‘define’(定義)

ということで，関数定義の開始をあらわす．この行には関数名と仮パラメータを書 く．関数のパラメータにどのようなデータを想定するのかを明示するような文法は無い．２行目の文字列は

‘docstring’

といい，

IPython console

で

make_g ?

と打ったりした時に表示される説明文である．３行目はこの関数に与えられたリストの全成分

に

π/90

を掛けている．三角関数には円周角

(

ラジアン

)

で角度を与えなければならない．４行目は上の結果を使って行列のすべ

ての行を作り全体をまとめて結果としている．array や

DataFrame

はリストや行列データの特性を数値計算向きにするものであ

る．未知パラメータは

A,B,C

と

3

つあるので，計画行列は本来

3

列になるが，C は定数項であるため別扱いになる．

10 240 0.040

11 270 0.081

12 300 0.481

13 330 0.843

14 360 0.782

これを

File

メニューの，

Save as...

で保存する．ファイル名には英語アルファベットと数字が使えるが，ひ らがなや漢字を使うこともできる．説明の都合上，今はファイル名を

data10

としておこう．もし拡張子を付 けるなら，テキストファイルという意味で

.txt

を付けるとよい．

Desktop

に保存しよう．

ファイルに保存されているデータを計算に使うには，それを

python

に読み込んでやる必要があるが，ま ず，

IPython console

で

pwd

と打って作業ディレクトリを表示し，これがデータファイルの保存場所と異 なる場合は，作業ディレクトリをファイルの保存場所に変更する

^*21

．

データファイルを読込むために，

IPython console

^{で次のように打つ．}

1 file = 'data10'

2 data = read_table(file, sep='\s+', comment='#')

ファイル名はいろいろ変わるので，代わりに

“file”

のような変数を使うのがよい．

sep='\s+'

^は，

read_table

関数に，データが１つ以上の空白によって区切られていることを教えるために必要なパラメータである．ま

た，

comment='#'

は注釈文字を指定するパラメータである．読み込んだデータには

“data”

という変数名が

付いたので，

data

と打てば中身を確認することができる．

ｘ

-y

グラフを表示するには次のようにする．

1 data.plot(x='x', y='y', kind='scatter')

'scatter'

^とは

Excel

などの「散布図」と同じ意味である．測定の間違いがわかる場合があるから，グラフ

を表示してよく見ることは重要である．

5.6 フィッティング

線形フィッティングを

pandas

モジュールの

ols

関数でおこなう．

IPython console

で試してみよう．

1 ols(x=make_g(data.x), y=data.y)

ols

には，

x=

で計画行列を与え，

y=

で 測定値リストを与える．

data.x

^は

data

^中の

x

列を意味する．

data.y

は

y

列．

ols

の結果は多くの数値からなるので，それらを一旦適当な変数名でまとめて保存するのがよい．そこから パラメータの推定値を取り出して

A, B, C

に入れ，それらの標準誤差を

dA, dB, dC

に入れる．

*21Editor

ウィンドウのタブ

(左上の出っ張り)

の右クリックメニューで

Set console working directory

を選ぶと，作業ディ

レクトリがファイルの保存ディレクトリに変更される．ファイル保存場所と作業ディレクトリが一致している場合は，ファイル名

のみでファイル指定ができる．pwd は

print working directory

のこと.

1 fit = ols(x=make_g(data.x), y=data.y)

2 A, B, C = fit.beta

3 dA,dB,dC = fit.std_err

ここで，

A, dA, B, dB, C, dC

と打てば，

3

つのパラメータの推定値と標準誤差がわかる．

式

(15)

から，元の非線形モデルのパラメータは次のようにして得られる

(

角度はラジアンから度に変換

)

^．

1 a = 2 * sqrt( A ** 2 + B ** 2 )

2 b = math.atan2( B, A ) * 90 / pi

3 c = C - a/2

atan2()

には

math.

を付けなければならない．これらのパラメータの値の標準誤差は

A, B, C

の標準誤差か

ら求められるが，ここでは

b

が大事なので

b

の標準誤差だけを求める．誤差の伝播則で計算すると

∆b= s ∂b

∂A∆A

!2

+ ∂b

∂B∆B

!2

= 2 a²

q

(B∆A)²+ (A∆B)²

となる

^*22

．よって，

b

の標準誤差は次のように計算すればよい．

1 db = sqrt( (B*dA)**2 + (A*dB)**2 ) * (2/a**2) * (180/pi)

5.7 グラフの完成

決定したパラメータを使ってモデルのグラフを表示しよう．それには，

fit.predict

という関数

(

オブジェク ト

fit

のメソッドともいう

)

を使って，測定点以外での値を計算する．細かい間隔で計算しないとなめらかな 曲線にならないので，

1°

ごとに計算する

^*23

．

1 xfine = arange(data.x.min(), data.x.max())

2 plot(xfine, fit.predict(x=make_g(xfine)))

グラフのタイトルを入れることができるが，そこに計算結果を入れておく．

1 title('{}: a={:.4}, b={:.4}(sigma={:.2}), c={:.4}\n'

2 .format(file, a, b, db, c))

.format

の前の文字列の中の中括弧の部分が

file, a, b, db,

および

c

の値に置き換えられて表示される．中括 弧の中には，

{:.4}

のように，表示精度を入れることができる．

*22

式

(14)

または式

(15)

から，∂b/∂A

=−2B/a²

，∂b/∂A

= 2A/a²

．∆

b,∆A,∆B

は各パラメータの標準誤差．ただし，これらの式 では角度の単位はラジアンである．

*23

測定データの独立変数の最大値と最小値を

data.x.min()

と

data.x.max()

で取得し，arange 関数によって，最小値と最大値

の間で間隔

1

の数値列を作り，この数値列を

make_g

関数に与えて細かい間隔での計画行列を組み立てる．それを

fit.predict

に与えると，すでに

fit

の中に収められている情報を使って，関数値リストが計算される．plot 関数には横軸データとして

xfine

も与える必要がある．

グラフに画像保存のボタンがあれば，それでファイルに保存することができるが，プログラムで保存するに は，

gcf().savefig(file+'.png')

^{のようにやればよい．}

file+'.png'

は，データファイルの名前に拡張

子として

‘.png’

を追加したものを意味する．

gcf()

はいま表示している図の情報．

5.8 解析をおこなう関数の作成

データを読み込んでグラフを作るまでの手続きを

analyse

という名前の関数にまとめよう．この関数はデー タファイル名を与えて呼び出すようにする．

Editor

^{ウィンドウの}

make_g

の定義に続いて，次のコードを入 力する．

1 def analyse(file):

2 u"file

で与えられたデータを解析し結果を図ファイルに出力

"

3 data = read_table(file, sep='\s+', comment='#')

4 fit = ols(x=make_g(data.x), y=data.y) #

^推定計算

5 A, B, C = fit.beta #

^推定値

6 dA,dB,dC = fit.std_err #

^{推定値標準誤差}

7 a = 2 * sqrt( A ** 2 + B ** 2 ) #

^振幅

8 b = math.atan2( B, A ) * 90 / pi #

^偏光方向

9 c = C - a/2 #

^最低値

10 db=sqrt( (B*dA)**2 + (A*dB)**2 )*(2/a**2)*(180/pi) #b

^誤差

11 data.plot(x='x', y='y', kind='scatter') #

^{測定点プロット}

12 xfine = arange(data.x.min(),data.x.max()) #1°

^{ごとに分割}

13 plot(xfine, fit.predict(x=make_g(xfine)) ) #

^{曲線プロット}

14 title('{}: a={:.4}, b={:.4}(sigma={:.2}), c={:.4}\n'

15 .format(file, a, b, db, c))

16 gcf().savefig(file+'.png') #

^{画像ファイルへ保存}

17 return {'a':a, 'b':b, 'db':db, 'c':c}

18 #%%

return

文はキーと値の形式で結果を返している．

Run current cell

^{をしてから，}

IPython console

^で

1 analyse('data10')

と打って呼び出す．結果が

data10.png

というファイルに保存されることを確認し，うまくいったら，ここま でのプログラムを

1

つのセルにまとめるため，区切り記号をすべて消すか

#%

に変え，このプログラムを保存 する．以後は，

IPython console

^で

from henkou import analyse

^{とするだけで，}

analyse

^{関数が使える} ようになる．

なお，以上の計算では単位の表示を考慮していないが，報告書では当然単位の明示が必要である．

Linux

や

Mac OS X

では，このプログラムを独立した完全なコマンドにすることもできる．それには，プログラムの

1

行目に

#! /usr/bin/python

と書き

(これはshebang

という)．最後のところに次のコードを追加して保存する．

import sys

if __name__ == "__main__":

for f in sys.argv[1:]:

analyse(f)

最後に，このファイルに実行権を付ける

(たとえばファイルのプロパティーで)．端末を開いて，

./henkou.py data10 data11 data12

のように打てば，複数のデータファイルの解析が一度にできる．データは任意のテキストエディタで入力することができる．

6 ^{線形フィッティング} – Excel ^編

Excel

のような表計算プログラムでは，

LINEST

という関数で線形フィッティングをおこなうことができ

る． 表を埋めたりグラフを表示するのに手間がかかるが，結果を着実に得ることができる．

次の図のように，測定値として，

x

と

y

の値を縦向きに並べる． ここでは三角関数を使うので，

x

はラジア ンに直している． その右にモデル行列を入れる． ここでは「パラメータ係数」と称している．定数項はモデ ル行列ではすべて

1

であるが，表計算の場合はそれを入れる必要は無い． 「フィットパラメータ」の２行下の ３つのセルをセレクトして，そこに

LINEST

関数を入れる．

=LINEST(C3:C15,D3:E15,1)

LINEST

関数には，

y

の範囲とモデル行列の範囲を指定する． 定数項の計算もさせるために，

3

番目のパラ

メータを１にする． この式を配列式として確定するため，

control

キーと

shift

キーを押しながら

enter

する．

パラメータの推定値が表示されるが，定数項以外の順序がなぜか逆になる．

推定曲線は推定値を使って別に計算して表示しなければならない． なぜなら，測定点だけで計算しても滑

らかな線にはならないからである． 図では，表のこの部分は上の方の一部だけが見えている．グラフは「散布

図」として作成する． ２系列のデータを一つのグラフに入れる方法を調べるには，とにかくいろいろ試して

みればよい．

測定値 パラメータ係数

degree x y cos sin

0.4165 1

30 0.524 0.2264 0.5 0.866 60 1.047 0.0318 -0.5 0.866 90 1.571 0.0432 -1 1E-016 120 2.094 0.2223 -0.5 -0.866 150 2.618 0.4315 0.5 -0.866 180 3.142 0.4371 1 -2E-016 210 3.665 0.2548 0.5 0.866 240 4.189 0.0408 -0.5 0.866 270 4.712 0.0215 -1 4E-016 300 5.236 0.2149 -0.5 -0.866 330 5.76 0.4266 0.5 -0.866 360 6.283 0.4128 1 -5E-016

フィット・パラメータ

B A C

-0.11 0.1982 0.2294 degree x y

0.4275 6 0.105 0.401 12 0.209 0.3669 18 0.314 0.3268 24 0.419 0.2824 30 0.524 0.2358

0 50 100 150 200 250 300 350 400

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

0.5 列 C

列 C

角度（degree）

透過率(任意単位)

A B C D E F G H I J K

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

7 ^{非線形フィッティング} – mathematica ^編

数学ソフトの

Mathematica

を使って非線形フィッティングをおこなってみよう． どういう計算をしてい るかがわかるようにする． 答えは一意的に求まる．

まずは

Mathematica

を起動して，ノートブックを表示する． 起動法や基本的な操作法は附録に書いたの

で，適宜参照してほしい．

7.1 データ入力

データを次のように入力する．ノートブックで，

dat=

とタイプし，続けて，挿入メニューの表・行列 の 新 規作成で，ダイアログを出し， 「列数」を

2

にして

OK

ボタンを押す． 測定値を下のように入力する． 左側は 回転角度，右側は透過光強度である． 行を追加するには，入力済みデータの一番下の行で，

hcontroli-henteri

を打つ．

　

dat=

　

0 0.809 30 0.403 60 0.039 90 0.077 120 0.477 150 0.836 180 0.794 210 0.400 240 0.040 270 0.081 300 0.481 330 0.843 360 0.782

;

空白データが無いようにする． データがすべて入ったら，表の右側にセミコロンを付けて，

hshifti-henteri

を 打つ． 間違えたらデータを修正してから

hshifti-henteri

を打てばよい．

dat

とすると，内容を確認すること ができる． 回転角度と透過光強度からなる複数のリストが全体としてひとつのリストにまとめられている．

Mathematica

では，何か動作をさせるには

hshifti-henteri

を打つ必要があることを忘れないこと．

次に，今入力したデータをグラフ表示する．

ListPlot[dat]

50 100 150 200 250 300 350

0.2 0.4 0.6 0.8

念のため，これから使う名前を洗浄する． 何か変なことが起こるとき，名前の洗浄をしなかったのが原因 ということは多い．

Clear[x,a,b,c]

7.2 モデルの入力

モデル関数を定義する．

f[ x_, a_, b_, c_ ] := a Cos[(Pi/180)(x-b)]^∧2 + c;

モデルに適当なパラメータを入れてみて，グラフでデータと比較する．

Show[Plot[f[x,0.8,160,0], {x,0,360}], ListPlot[dat]]

50 100 150 200 250 300 350

0.2 0.4 0.6 0.8

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 50 100 150 200 250 300 350

0.0 0.2 0.4 0.6 0.8

もう一度結果を見る．

s

ファイルメニューの「保存」でノートブックを保存する． ただしこのファイルは

Mathematica

^{が無いと役}

に立たない．グラフを右クリックし， 「形式を選択してグラフィックスを保存

...

」でグラフを保存することがで

きる． ファイルの種類はワープロ用なら

PNG

がよい．用途によっては，

PDF

や

EPS

がよいかもしれない．