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

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

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

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

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

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

No.1 回答者： alfeim 回答日時： 2001/08/31 22:15

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

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

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

この回答へのお礼

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

お礼日時：2001/09/01 01:17

No.2 回答者： paatje 回答日時： 2001/08/31 23:32

話としては、混合モデル(Mixture model)の推定

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

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

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

この回答へのお礼

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

お礼日時：2001/09/01 01:20

No.3 ベストアンサー 回答者： hazenoki 回答日時： 2001/09/01 00:06

パルスを探すアルゴリズムの具体的な文献はしりませんが、

たとえば、ｘが小さい方から順にｙの値を見ていき最大値を探し
ｙが一定値以下になったら次の最大値を探し始めるようなのでは
どうでしょうか。

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

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)
などがあります。ふつうは最小二乗法かな？ データ処理や統計処理の教科書に
載っていると思います。たぶん。

長くなってすみません。

この回答へのお礼

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

お礼日時：2001/09/01 01:33