[ 岡山大学 | 理学部 | 地球科学科 | 地球および惑星大気科学研究室 ]

2019年度 大気科学演習1

gnuplot

フィッティング

データを任意の関数でフィティングする.

マウナロアで測定されたCO2濃度を図にしてみる

gnuplot> reset
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2

データの各列がそれぞれ何であるのかは,ファイルの中を見たら書いてあるので,各自で確認してください. 欠損値が -999.99 になっているので,

gnuplot> set datafile missing '-999.99'
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2

データを関数でフィッティングするときは,(1)関数を定義,(2)初期値の設定,(3)フィッティング,の順番でおこなう.(2)は省略することもできるが,初期値がダサいとフィッティングに失敗することがある.

とりあえず,やってみる

gnuplot> f(x)=a*x+b
gnuplot> fit f(x) '/mnt/data/co2/insitu.txt' using 0:8 via a,b
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, f(x)

1行目で関数を定義,2行目でフィッティング,3行目は結果の表示,である. フィッティングをすると,なんかいろいろなことが表示される(ちゃんと読んだら,フィッティングして決まったパラメタの値とか,その当てはまり具合とか,そういうことが書いてあるとわかる). 一次関数でフィッティングしたら,こういう結果になるだろうという結果になっている.

関数を変えて試してみる.

gnuplot> g(x)=a+b*exp(x/c)
gnuplot> a=330
gnuplot> b=80
gnuplot> c=600
gnuplot> fit g(x) '/mnt/data/co2/insitu.txt' using 0:8 via a,b,c
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, g(x)

けっこういい感じに見える.

CO2は季節変化しているので,それもフィッティングしてみる.

gnuplot> h(x)=a+b*exp(x/c)+d*sin(e*x+f)
gnuplot> d=3
gnuplot> e=0.5
gnuplot> f=1
gnuplot> fit h(x) '/mnt/data/co2/insitu.txt' using 0:8 via a,b,c,d,e,f
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, h(x)

これはあんまりよくないように見える. 季節変化が明らかに過小評価されている. おそらく,与えた初期値がいまいちだったので,最適解を見つけることに失敗したのだと予想される. こういうときは,初期値を調整して再度試してみる,というのが常套手段である. 答えに近い初期値を与えれば,解はちゃんと見つかる. とはいえ,答えがわからない場合も多い(というか,答えを知っていたらフィッティングなどする必要がない). 答えに近い初期値を与えればよいというのは理屈では正しくとも役に立たないことが多かったりする.

ここでは,ちょっと違う方法で答えを探してみる. まず,フィッティングを全域でなく狭い範囲に限定しておこなう

gnuplot> d=3
gnuplot> e=0.5
gnuplot> f=1
gnuplot> fit [0:100] h(x) '/mnt/data/co2/insitu.txt' using 0:8 via a,b,c,d,e,f
gnuplot> set xrange [0:100]
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, h(x) 

fit の後ろに範囲を書くと,その範囲に限定してパラメータのフィッティングがおこなわれる. ここでは [0:100] の範囲でフィッティングをおこなっている. 結果を見ると,季節変化をうまく抽出できているように見える. ちなみに,全体に広げるとずれている

gnuplot> set xrange [0:600]
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, h(x) 

次は,[0:100] の区間でフィッティングした結果を初期値にして,もうちょっと広い区間でフィッティングをおこなう. フィッティングをおこなうと,それぞれのパラメタ(ここでの例だと a, b, c, d, e, f)にフィッティングで求められた最適解が代入される. したがって,そのまま次のフィッティングをおこなえば,先のフィッティングで得られた結果が初期値となる.

gnuplot> fit [0:200] h(x) '/mnt/data/co2/insitu.txt' using 0:8 via a,b,c,d,e,f
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, h(x)

式のグラフがいまいちなのを修正する

gnuplot> set samples 2000
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, h(x)

最後は範囲の制限を解除してフィッティングする

gnuplot> fit h(x) '/mnt/data/co2/insitu.txt' using 0:8 via a,b,c,d,e,f
gnuplot> plot '/mnt/data/co2/insitu.txt' using 0:8 with points pointtype 7 pointsize 0.2, h(x)

フィッティングの結果はfitを実行するたびに表示されるが,それと同じものが fit.log という名前のファイルに書き込まれている. gnuplotを終了して見てみる.

gnuplot> exit
$ cat fit.log

季節変化の振幅(d)は 2.9 (ppm)で,その周期は 2 pi / e = 11.99 (month) となっている. 周期が1年になっていることは,そうであるべきものが正しく求まっているということで,とてもよいことである. 繰り返しになるが,計算機が出した結果を鵜呑みにしてはいけない. 常識やその他の知識に照らし合わせて,おかしなことが起こっていないかどうか,確認できることは確認しなければならない.

補足

何が目的であるのかに依存するが,周期は1年であるとわかっているのであれば,パラメタ e を求める必要はなく,最初から関数を

h(x)=a+b*exp(x/c)+d*sin(2.0*pi*(x/12.0)+f)

で与えるのが正解,という場合もあるだろう.




Last Updated: 2019/11/10, Since: 2019/11/10.
This page is generated by Makefile.rd2html.