物理学特別実験

最小2乗法を用いたモデル関数の最適化(フィッティング)の練習

手順

以下では情報端末室のPCを利用することを想定しています。自分のPCを用いる場合, "z: ドライブ"を各自のPCの"C:ドライブ"等に読み替えて下さい。また，gnuplotをgnuplot homepageなどからダウンロードしてPCにインストールして下さい。Mathematicaは個人で使う場合には購入すル必要があります。大学の端末室で利用しましょう。
  1. z:ドライブに fitfile という名前のフォルダを作り, その中に次のファイルをダウンロードする(windowsではマウスを右クリックし, 「リンク先を名前を付けて保存」をクリック).
    最小2乗法復習ファイル (ファイル名 lsfit.pdf)
    Mathematica用ファイル (ファイル名 lsfitting.nb)
    dataファイル (ファイル名 12.dat)
    dataファイル (ファイル名 22.dat)
    dataファイル (ファイル名 32.dat)
    dataファイル (ファイル名 42.dat)
  2. 「最小2乗法復習ファイル」を読んで, 多項式のフィッティングについて復習しよう. 要点は□で囲ってある部分に出ている.

gnuplotによる最適化

gnuplotはデータフィッティングの機能ももつグラフ作成ソフトです. 様々な関数を利用することも簡単に出来る.
  1. gnuplotを起動する.
  2. fitting dataの入っているディレクトリを指定する
    gnuplot> cd 'z:\fitfile'
  3. フィットするデータをグラフに表示
    gnuplot> plot '12.dat'
  4. フィットする関数を入力する
    gnuplot> f(x)= a*x+b
    (a, b はfitting 係数, gnuplotが備える関数はタイトルバーの[Functions]に掲載されている. xの2乗はx**2と書く)
  5. 関数f(x)をデータにフィットする
    gnuplot> fit f(x) '12.dat' via a,b
    via の後にはfittingで決める係数を記入する. 何回かの繰り返しの後係数が誤差とともに表示される.
  6. データとフィットされた関数を図示
    gruplot> plot f(x), '12.dat'
うまく計算できないときのヒント
  1. 計算が中断して,
    「 undefined value during function evaluation　」
    などと表示される場合がある. フィッティング係数(a, b, ...等)の初期値が良くないことが原因となっていることがある. フィッティング係数の値は,
    gnuplot> show all
    により確認できる. 出来るだけよい初期値を与えるためには, データ点とフィッティング関数をよく見比べながら, 係数の初期値を決めてやる必要がある. 例えば,
    gnuplot> a=4; b=5;
    等のように, 具体的に数値を代入してやり
    gnuplot> plot f(x), 'データファイル'
    として, データと関数を同時に描かせながら, 初期値を出来るだけよい値にしてやるとうまくいく.
  2. グラフの表示範囲が適切でない
    グラフ上で右クリックしたままマウスを移動することにより、表示部分を変更ができる。しかし、うまくやらないと画面からグラフが見えなくなってしまう。 そのような場合には、
    set xrange [0:10] ( 0＜x＜10の例)
    set autoscale (範囲を自動で決める)
    コマンドを使うとよい。
  3. 軸のラベル
    set xlabel 't (s)'
    など。
  4. その他

Mathematicaプログラムによる最適化

  1. 上で用いた文字変数 data, x, yはMathematicaではListと呼ばれており, いくつかの数字をまとめたもので, これは行列やベクトルと同様な扱いができる. 「data」と入力して[Shift]+[Enter] keyでその中身が表示される. 例えば {{1,2},{3,4},....}のように, 各成分は{ }で一つにまとめられ, それが更に{ }でまとめられる. 各要素には 例えば data[[2,1]]とやることでアクセスできる. これは2番目の中括弧の1番目のデータということで, 上の例では3ということになる.
  2. list1=Table[・・・・・, {n,1,ndim}] では, 1からndimまでのndim個のデータ・・・・・からなるListを作り, それにlist1と名前を付けている.
  3. Sum[・・・・, {i,1,ndata}] では, iの値を1からndataまで変えながら,・・・・を加えている.
  4. 「.」(ピリオド)は内積と似ている. a={{1,2},{3,4}}, b={c,d}のとき,
    a.b={c+2d, 3c+4d}, b.a={c+3d,2c+4d}
    となる.
Mathematicaの練習
1. 変数cに適当な数値を代入, それを表示しなさい.
2. 変数dに3変数のListを代入し, それを表示しなさい.
3. 変数eに3×3のListを代入し, それを表示しなさい.
4. 変数fにeの逆行列を代入し, それを表示しなさい.
5. eとfの積を表示しなさい.
6. 1から100までの和を求め, gに代入しなさい.
7. sin x (-π< x < π)のグラフを描きなさい.
8. x2+y2 (-1< x < 1, -1< y < 1)のグラフを描きなさい.

解答例 (大文字小文字に注意)
1. c=5; Print[c]
2. d={1, 2, 3};Print[d]
3. e={{1, 2, 3}, {4, 5, 6}, {7, 8, 10}}; Print[e]; Print[TableForm[e]]
4. f=Inverse[e]; Print[f]
5. Print[e.f]
6. g=Sum[i, {i,1,100}] または, g=0; Do[g=g+i,{i,1,100}]; Print[g]
7. Plot[Sin[x], {x, -Pi, Pi}]
8. Plot3D[x^2 + y^2, {x, -1, 1}, {y, -1, 1}]

参考書
スティーブン ウルフラム Mathematica --A System for Doing Mathematics by Computer-- アジソンウエスレイ


課題

(a) 22.dat, (b) 32.dat, (c) 42.dat に入っているデータについて適切な関数をフィットし, その関数, 求めた係数, その誤差を示すとともに, データ点, 関数を図にまとめよ. ただし,
32.dat は横軸が時間を表していると仮定して, 減衰時間を誤差とともに求めよ. (減衰時間は, 縦軸の値が1/eの大きさになるまで経過する時間とせよ)
42.dat は横軸が周波数を表していると仮定し, 中心周波数, 半値全幅(ピークとバックグラウンドの中間を取る位置での幅)を誤差とともに求めよ.


注意 誤差を求める場合には次の誤差伝播の法則を利用しよう.
誤差の伝播 x, y の分散がそれぞれσx2, σy2で, z がx, y の関数 z=f(x, y) と表される場合に, z の分散は
σz2=(∂ f/∂ x)2σx2+(∂ f/∂ y)2σy2
と表される.

tips