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

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

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

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

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

このQ&Aに関連する最新のQ&A

A 回答 (3件)

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


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

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

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以下になることを仮定しています。
この50と80の差はyのふらつきを考慮したものです。
また、データは xについてソートされているものとしています。

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

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

長くなってすみません。
    • good
    • 0
この回答へのお礼

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

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

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


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

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

私の読んでいる本では、EMアルゴリズムという
アルゴリズム(とGibbs Samplerというアルゴリズム)
を使っているようです。
    • good
    • 0
この回答へのお礼

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

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

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



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

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

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

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

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Qガウスフィッティングについて

時間 t に対するデータ列 yi (i=1, 2, …) に
最小二乗法を用いて
ガウス関数 y(t) = a + b exp( -(t-c)^2 / 2σ^2 ) を
当てはめたいのです。
これまでは、エクセルなどのソフトを用いて
未知数 a, b, c, σを求めていたのですが、
これをソフトなどを用いずに理論式で解くにはどうすれば良いでしょうか?

ひとつのやり方として、
ガウス関数の右辺の a を左辺に移行して、両辺の対数を取り
対数データ列 ln(yi - a) との誤差を求めて最小二乗法…
という手を考えたのですが、
これだと、誤差の総和を求める際、Σの中にパラメータの a が
ln(yi - a) という形で残ってしまい、
誤差の偏微分で得られるパラメータに関する方程式が、
非線形な連立方程式となって解くことができません。
(説明が分かり難くて申し訳ありません)

その他にもいろいろ考えたのですが、
どうにもうまくいきません。
とにかく、パラメータの a が邪魔で、
これを b, c, σを用いずに先に求める、
という方向で考えているのですが、
何か良い方法がありましたら、よろしくお願いします。

時間 t に対するデータ列 yi (i=1, 2, …) に
最小二乗法を用いて
ガウス関数 y(t) = a + b exp( -(t-c)^2 / 2σ^2 ) を
当てはめたいのです。
これまでは、エクセルなどのソフトを用いて
未知数 a, b, c, σを求めていたのですが、
これをソフトなどを用いずに理論式で解くにはどうすれば良いでしょうか?

ひとつのやり方として、
ガウス関数の右辺の a を左辺に移行して、両辺の対数を取り
対数データ列 ln(yi - a) との誤差を求めて最小二乗法…
という手を考えたのですが、
これだと、誤差の総和を求める際、Σの...続きを読む

Aベストアンサー

 最小化したいのは
E=Σ(yi-y(t))^2
ですね。パラメータをベクトルで
p=(a,b,c,σ)
と書くことにします。

 どうやっても非線形で、反復計算が必要です。「理論式」という観点で言い換えるなら、「最小二乗解pに収束する漸化式」を作るところまでしか出来ない、ってことです。
 その漸化式は、たとえば「第n回目の反復で得られたパラメータの推定値p[n]を使って、(∂E/∂p)(p[n])を係数とするpの一次式でEを近似し、これを使ってp[n]よりも最小二乗解に近いp[n+1]を計算する」というものですから、「適切な初項を与えれば最小二乗解に収束する漸化式」でしかありません。(他にもアプローチはありますが、初項が必要なのは同じ。)
 ここで「適切な初項」というのは、「初項で決まるガウス曲線がもしデータとかけ離れていたら、漸化式が収束しない」ってことです。では初項が「適切」かどうかを予め判定する手段はというと、これまた一般にはありませんで、計算してみて収束していくかどうかを見るしかない。そういう事情です。

 なので、実地に計算するには、漸化式の初項(つまりpの最初の推測値)を決める方法が必要。これはイロイロ工夫するしかありません。

 以下は「適切な初項を決める手段」の話なので、現段階ではご興味ないのかも知れませんが:

 「パラメータaをなんとか推定したうえで、対数を取って線形最小二乗法でb,c,σを出す」というのは一つの方法です。ただしこのままだと
F=Σ(ln(yi)-ln(y(t)))^2
を最小化することになり、Eの最小化とはかなり違います。むしろ
G=Σ(yi (ln(yi)-ln(y(t)))^2
でEをよく近似できます。( ln(1+x)≒x を使っています。)
 もちろん、これは飽くまで初項を推定する方法であって、それを使って「Eを最小化する漸化式」で反復計算をすることになります。

 別の手段として、もしデータyiによくフィットする曲線y(t)が存在する(データがほぼ綺麗にガウスカーブになっている)なら、データの最大値付近に放物線をフィッティングするなどしてcの推測値を決め、さらに直線
y=(最大値-最小値)/2
とデータの折れ線グラフとの交点の間隔2w(半値幅)からσを推定して、線形最小二乗法でa,bを決める、などなど。

 最小化したいのは
E=Σ(yi-y(t))^2
ですね。パラメータをベクトルで
p=(a,b,c,σ)
と書くことにします。

 どうやっても非線形で、反復計算が必要です。「理論式」という観点で言い換えるなら、「最小二乗解pに収束する漸化式」を作るところまでしか出来ない、ってことです。
 その漸化式は、たとえば「第n回目の反復で得られたパラメータの推定値p[n]を使って、(∂E/∂p)(p[n])を係数とするpの一次式でEを近似し、これを使ってp[n]よりも最小二乗解に近いp[n+1]を計算する」というものですから、「適切な初項を与...続きを読む

Q任意の形の波形を関数でフィッティングするためにはどうしたら良いですか?

任意の形の波形を関数でフィッティングするためにはどうしたら良いですか?

ある実験より得られた波形があります。

その波形が多項式によってフィッティング可能であると仮定してフィッティングを行いたいのですが、
どうすれば良いのでしょうか?

例えば、その波形が10乗までの関数で表されると考えたとしても
y=A + Bx^1 + Cx^2....Kx^10


とフィッティングパラメータが11個もあるわけですよね?
そうすると、フィッティングをかけたとしても、一通りだけでなく幾通りも解が出てきてしまうのではないでしょうか?

昔、電子工学の授業でこういう任意の波形のフィッティング式の作り方に関して
習ったことがあるのですが、完全に忘れてしまいました。
どなたかご存じでしたら、参考書などを教えて頂けますでしょうか?

Aベストアンサー

実用的には,先の二人の回答者が考えてみえるパラメトリック回帰ではなく,ノンパラメトリック回帰になるんでしょうね.
例えば,RBFとかKrigingとか.
設計最適化における応答曲面近似においては,一般的に利用されています.
一方,化学の世界の機器分析では,ニューラル・ネットワークのような学習的手法が利用されます.
イオンクロマト分析のスペクトル解析などに使われています.これは,ピークの出る位置が先験的に分かっていて,繰り返し出てくる場合などに応用できます.

その先には,EMアルゴリズムとかベイジアン・ネットワークとかがあるようですが,
私は詳しくないので,他の方に説明はお任せします.

参考書:洋書ですが,David Ruppert,M.P.Wand,R.J.Carroll(2003):Semiparametric Regression,Cambridge University Press (ペーパーバックで8000円くらい)

Q曲線と曲線の交点を通る曲線の求め方(曲線群)

皆様、こんにちは。

円A:f(x,y)と円B:g(x,y)の交点を通る円の方程式は全て
kf(x,y)+lg(x,y)=0の形で表せると習ったのですが、

これの応用で
円A:f(x,y)と円B:g(x,y)の交点を通る三次曲線は全て
(ax+b)f(x,y)+(cx+d)g(x,y)=0・・・・(1)
の形で表せるのでしょうか?

もし2円の交点を通る3次曲線が全て(1)で表せるのでしたら
その証明方法なども教えてください。

よろしくお願いします。

Aベストアンサー

#1さんがすでにご指摘の通り、
三次曲線をすべて (1) で表すことはできません。
そして、kf(x,y)+lg(x,y)=0 も、
二次曲線すべてではなく、円だけを表しているにすぎません。
でも、vigo24 さんの着眼点はなかなかおもしろいと思います。
これをヒントにして、次のように考えてみました。

f1(x,y) = x^2 + y^2 + l1 x + m1 y + n1
f2(x,y) = x^2 + y^2 + l2 x + m2 y + n2
とおきます。そして、f1(x,y) = 0 と f2(x,y) = 0 が
円であり、2点で交わるものとします。
(円でない場合もありますし、2点で交わらない場合もあります。
詳細は http://okwave.jp/qa3076718.html をご覧下さい。)

f3(x,y) = f1(x,y) - f2(x,y) = (l1-l2)x + (m1-m2)y + (n1-n2) とおきます。
すると、f1(x,y) = 0 かつ f2(x,y) = 0 ならば f3(x,y) = 0 が成り立ちます。
これを利用して、次のような方程式を考えてみます。
f4(x,y) = a f1(x,y) + b f2(x,y) + (cx+dy) f3(x,y) = 0

まず、f4(x,y) = 0 が二次以下の曲線になることはすぐにわかります。
そして、f1(x,y) = 0 と f2(x,y) = 0 の交点を (x1,y1) , (x2,y2) とすると、
f1(x1,y1) = 0 かつ f2(x1,y1) = 0 ですから、f3(x1,y1) = 0 となり、
f4(x1,y1) = a×0 + b×0 + (c x1 + d y1)×0 = 0 となります。
ですから、f4(x1,y1) = 0 , f4(x2,y2) = 0 となり、
曲線 f4(x,y) = 0 は f1(x,y) = 0 と f2(x,y) = 0 の交点を通ります。

さて、問題となるのが、f4(x,y) = 0 が2点を通る二次以下の曲線を
すべて網羅しているかどうか、ということです。
このことについて厳密な証明はしていませんが、
二次以下の曲線は ax^2 + bxy + cy^2 + dx + ey + f = 0 で表せます。
つまり、パラメータは6個(0でない係数で割れば5個)必要です。
そして、曲線 f4(x,y) = 0 のパラメータは、a,b,c,d と交点2個で、
やはり6個(0でない係数で割れば5個)あることになります。
ですから、f4(x,y) = 0 の形で、
2点を通る二次以下の曲線をすべて表せているものと思われます。

ここからいろいろ発展させることができると思います。
例えば、円と円の交点ではなく、二次曲線と二次曲線の交点で、
交点が4個あるとすれば、k1 f1(x,y) + k2 f2(x,y) = 0 の形で
交点を通る二次曲線を表せるかもしれません。
vigo24 さんのもともとの質問である、円と円の交点を通る三次曲線も、
(a1 x + a2 y + a3) f1(x,y) + (b1 x + b2 y + b3) f2(x,y)
+ (c1 x^2 + c2 xy + c3 y^2) f3(x,y) = 0
の形で表すことが可能ではないかと思われます。
(上記の形ではパラメータが1つ多すぎるので、
不要なパラメータが含まれていると思います。)

#1さんがすでにご指摘の通り、
三次曲線をすべて (1) で表すことはできません。
そして、kf(x,y)+lg(x,y)=0 も、
二次曲線すべてではなく、円だけを表しているにすぎません。
でも、vigo24 さんの着眼点はなかなかおもしろいと思います。
これをヒントにして、次のように考えてみました。

f1(x,y) = x^2 + y^2 + l1 x + m1 y + n1
f2(x,y) = x^2 + y^2 + l2 x + m2 y + n2
とおきます。そして、f1(x,y) = 0 と f2(x,y) = 0 が
円であり、2点で交わるものとします。
(円でない場合もありますし、2...続きを読む

Qデータ列の曲線によるフィッティング

データ列(xi,yi)(i=1,...,n)を関数
y=α(x+γ)^2 exp(-β(x+γ)) (α>0,β>0,γ>0)
でフィッティングしたいです。対数をとって普通に最小2乗法で解こうとしたら得られる連立方程式が線形でなくて解けませんでした。
どうしたらいいのでしょうか?

Aベストアンサー

完全な意味での最小二乗法ではないけれど、無理矢理線形の問題に帰着して簡単に計算するstomachman我流のやりかたがあります。(超平面法と勝手に呼んでいます。)この方法は、「データがモデルに良く合う」という時に限り「最小二乗解とほぼ一致する」という性質があります。

モデルをデータと区別して、
f(t) = α((t+γ)^2) exp(-β(t+γ))
と書きましょう。
両辺をtで微分すると
f'(t) = 2α(t+γ) exp(-β(t+γ))-βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2α(t+γ)^2) exp(-β(t+γ))-(t+γ)βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2f(t)-(t+γ)βf(t)
(t+γ)f'(t) = (2-βγ)f(t)-βtf(t)
という微分方程式が得られます。
 この両辺をtで積分すると、(t=x1~xiまで)
∫(t+γ)f'(t)dt = (2-βγ)∫f(t)dt-β∫tf(t)dt
左辺を部分積分して
(xi+γ)f(xi)-(x1+γ)f(x1)-∫f(t)dt = (2-βγ)∫f(t)dt-β∫tf(t)dt
整理すると
xif(xi)-x1f(x1) = (3-βγ)∫f(t)dt-β∫tf(t)dt-γ(f(xi)-f(x1))
です。
ここで、f(xi)の代わりにyiを使っても、「データとモデルが良く合う」のなら大差は出ない筈。だから
∫f(t)dtの代わりにyをx1~xiまで数値積分したものJi、
∫tf(t)dtの代わりにxyをx1~xiまで数値積分したものKi、を用いることにして、残差をεiと書くことにし、
xiyi-x1y1 = (3-βγ)Ji-βKi-γ(yi-y1)+εi (i=1,2,....,n)

改めて
Yi = xiyi-x1y1
Li = yi - y1
A = (3-βγ)
B = -β
C = -γ
とおくと、
Yi = A Ji + B Ki + C Li + εi (i=1,2,....,n)
という線形問題になったわけで、これは簡単に解けます。
 ところで、A,B,Cはβとγだけで表される一方、αが出てきませんね。
つまり
・A,B,Cからβ,γを求めると1通りに決まりません。これは気にしないで、たとえばAを無視してβ,γを計算します。ついでに、こうして得たβ,γを使って(3-βγ)を計算してみます。もしこれがAとまるで合わないようなら、「データがモデルに良く合う」という条件が満たされていない。yiに含まれるノイズが大きすぎるんです。
・αを求めるには、
f(t) = α((t+γ)^2) exp(-β(t+γ))
を再び使います。γとβは既に分かっているから、簡単な線形最小二乗法の問題ですね。

より高精度の最小二乗解を求めるには、超平面法で得た解を初期値にして、非線形最小二乗法(ガウス・ニュートン法など)を使って改良すると良いです。(教科書としてはいつも、中川・小柳「最小二乗法による実験データ解析」東京大学出版会をお薦めしています。)

完全な意味での最小二乗法ではないけれど、無理矢理線形の問題に帰着して簡単に計算するstomachman我流のやりかたがあります。(超平面法と勝手に呼んでいます。)この方法は、「データがモデルに良く合う」という時に限り「最小二乗解とほぼ一致する」という性質があります。

モデルをデータと区別して、
f(t) = α((t+γ)^2) exp(-β(t+γ))
と書きましょう。
両辺をtで微分すると
f'(t) = 2α(t+γ) exp(-β(t+γ))-βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2α(t+γ)^2) exp(-β(t+γ))-(t+γ)βα((t+γ)^2) exp(-β(t+...続きを読む

Q柿2個、りんご4個、みかん6個の中から6個を取り出す方法は何通りあるか

柿2個、りんご4個、みかん6個の中から6個を取り出す方法は何通りあるか?ただし、取り出されない果実があっても良い。

この問題が分かりません。

Aベストアンサー

質問者様がこの問題が分からないように私にもこの問題が何を問うものなのかが分かりません。
『6個を取り出す方法は何通りあるか?』だけでは、回答出来ないですね。
おそらく他に何かの条件が有るのだと思いますので、それを記載してほしいです・・・。


人気Q&Aランキング

おすすめ情報