ガウシアンフィッティングのアルゴリズム

質問

c言語等を用いて実験データの解析を行おうと思っています。

データの形式は（x・・・時間軸、y・・・値）の列になっていて、
これをグラフにすると、パルス状の波形が連続する形になっています。

求めたい情報は、各パルス波形の山に対応する時刻の羅列
なのですが、現在ではガウス曲線近似機能のついたグラフソフトで、
一個一個手作業で求めています。

これを、解析プログラムを作って自動化しようと思っているのですが、
ガウス曲線（鋭い立ち上がりの）に近い部分を自動的に検出する方法、
また、フィッティングを数値計算的に行う方法がわかりません。

このようなアルゴリズムを考える上で参考になるようなHPや文献を
ご存知の方がいらっしゃれば、教えていただきたいと思っています。

回答No.3

hazenoki · ベストアンサー

パルスを探すアルゴリズムの具体的な文献はしりませんが、
たとえば、ｘが小さい方から順にｙの値を見ていき最大値を探し
ｙが一定値以下になったら次の最大値を探し始めるようなのでは
どうでしょうか。

言葉で説明すると難しいのでコードで示すと以下のようになると
思います。（適当に書いたのでバグがあるかも。）

int i;
int flag_peak=0; /* 0: 谷 1: 山 */
float y_max=0;
float x_peak=0;
float x[1000], y[1000]; /* data */
/* ここで x と y の配列にデータをセットする。 */
for(i=0; i<1000; i++){
　if( y[i]<50 ){
　　if( flag_peak ==1 ) printf("peak %d %d\n", x[i], y[i]);
　　flag_peak=0;
　　y_max=0;
　}
　if( y[i]<80 ) continue;
　if( y_max<y[i] ) {
　　y_max = y[i];
　　x_peak = x[i];
　}
　flag_peak=1;
}

ここで、パルスは上向きで最大値が常に80以上になる、
パルス間の谷は常に50以下になることを仮定しています。
この５０と８０の差はｙのふらつきを考慮したものです。
また、データは xについてソートされているものとしています。

このパルスの位置は最大値の場所と定義していますが、
フィッティングする場合は、この値を フィッティングの初期値として
使うことになると思います。

フィッティングは、最小二乗法や 最尤推定法 (maximum likelihood method)
などがあります。ふつうは最小二乗法かな？ データ処理や統計処理の教科書に
載っていると思います。たぶん。

長くなってすみません。

お礼

プログラムまで考えていただいて本当にありがとうございます！
パルスのフィッテングを行った後、次のパルスへと移行するときの
考え方の参考になりそうです。

回答No.2

paatje

話としては、混合モデル(Mixture model)の推定
というやつだと思いますが、
私も現在勉強中なので、アルゴリズムに関しては、
詳しく教えて差し上げられません。

「混合モデル」や「Mixture model」で
検索されてみてはいかがでしょうか？

私の読んでいる本では、EMアルゴリズムという
アルゴリズム（とGibbs Samplerというアルゴリズム）
を使っているようです。

お礼

alfeimさんの教えて下さった文献とともにしらべてみます。
混合モデルという言葉ははじめて聞きました。
有難う御座います。

回答No.1

alfeim

#無学なため、ガウス曲線なるものを知りませんが・・・

鋭い立ち上がり位置を検出、というのであれば、サンプルデータを数値微分してやれば一発では？

んで、アルゴリズム関連の書籍としては
・奥村晴彦『C言語による最新アルゴリズム事典』技術評論社，1991年，ISBN4-87408-414-1，2400円
がオススメ・・・というか定番ですね。
なお、ソースコードはVectorでダウンロード可能です(http://www.vector.co.jp/soft/data/prog/se002453.html)。

お礼

早速の解答、有難う御座います。
微分は試行錯誤では行ってみたのですが、
自力では上手くいきませんでした。