セットアップ日記

_ Gaussianでfit

たくさんのデータをgaussianでfitする必要があって、どうやってやろうか考えた。グラフソフトを使って、手でするのはデータが多すぎて現実的ではない。Rはプログラムを組めるので、可能ではあるだろうが、gaussianのfitをどうするかまだ知らないので、別の方法を試みることにした。調べてみると、GSLというライブラリがあり、科学計算に用いるいろいろな道具が揃っているらしい。さらに、これがrubyからたたけるらしい。というわけで、rubyでgaussian fitをしてみた。まずは、

aptitude install libgsl-ruby

でインストール。Histogramだとfit_gaussianというmethodがあるようだが、これだとxyデータをそのまま扱えない。いろいろと探していたら、MultiFitでできることが分かった。こんな感じだ。

require "gsl"
x=[0,1,2,3,4,5,6]
y=[0,1,2,3,2,1,0]
gx=GSL::Vector[*x]
gy=GSL::Vector[*y]
guess=[gy.min,gy.max,gx.mean,2]
y0,int,x0,var=GSL::MultiFit::FdfSolver.fit(gx,gy,"gaussian",guess)[0].to_a

いろいろと分からないことが多いが、配列を返してくるので、先頭にあるパラメータを取り出して、これがGSL::Vectorなので、rubyの通常の配列に戻して、結果を得ている。いろいろなデータを試してみると、明らかにピークがあるデータはうまくいくのだが、ピークがあるか無いかわからないようなデータだと、変な値になってしまう。ピークがなさそうなら、強度が0に近い値になって欲しいのだが。あとはこの辺りをなんとかしないと。