1.
tan(Y)=(a(X/b)/((X/b)^2-1)の式に
データ配列を最小二乗法によりカーブフィットさせたいのですが、線形代数での解き方をご教授お願いします。

2.
(マトリックスを使う場合は、どのように展開すれば
良いのでしょうか。)

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

A 回答 (3件)

ふふふ。

皆さん、非線形と決めてかかってますね。
じつは、データがたくさんあって、しかもモデルと良く合うことが分かっている場合には、実用上問題ない精度で線形計算ができます。
元のモデルを変形すると
tan(Y)X^2=(abX)+tan(Y)b^2
ですよね。ここで、
A= b^2, B= ab
yi = (tan(Yi) Xi^2), xi = tan(Yi) (i=1,2,....,N)
とおけば
yi= A xi+ B + 誤差
というモデルをフィッティングすることで、A, Bが得られます。A,Bからa,bを求めるのはできますよね。
 どうしても本来の問題をきちんと解くと言う場合には、やはり非線形最小二乗法が必要ですが、上記の方法でa,bの非常に良い近似値を出しておけば、あとは「ガウス-ニュートン法」、いやさ「ニュートン法」のようなシンプルかつ高速な方法で収束させることができます。

この方法は、縦軸を(tan(Y) X^2)、横軸をtan(Y)にとったグラフを描けばデータが直線上に乗り、この空間での最小二乗法をやる、という仕組みなのです。非常に応用が広いんですよ。私的に開発したもので「超平面法」と読んでいます。
    • good
    • 0

 ▼返答遅くなりました。

さあ回答しようとPCの前に座ったものの、はっと気がつくともう朝という状態が二晩続きました。P(^.^)

 さて、前回の質問や円の方程式もそうですが、展開するとaやbというパラメータについて線形結合になっています。このような場合は、代数的に解いたり、直接法(Gauss-Jordan法とか)で求めることができます。ところが、今回はtanY=abx/(x^2-b^2)ですからパラメータについて非線形です。これは直接法で解けないので、反復法で求める必要があります。反復法というのは、解に適当な初期値を設定した後、ある基準によって徐々に解を改善して、誤差が許容範囲内に収まれば終了するという方法です。SOR法、最急降下法、共役勾配法などがあります。
 a、b、評価関数(偏差平方和のような)からなる3次元空間があるとします。評価関数を高さ方向にとると、aとbの値によってうねうねとうねる曲面をイメージできるでしょう。山あり谷ありの曲面ですが、ある場所に評価関数最小の点があります。その点のaとbが求めたい最尤推定値です。しかし、この曲面上には霧が立ち込めており、どこに評価関数最小の点があるのか最初は見えません。そこで適当な場所に立ち、その場所の評価関数の値や勾配を調べて、次のステップはどれくらいの距離、どの方向へ進もうかと考えながら谷底(評価関数最小の点)へと降りていく、というのが反復法のイメージです。その「次のステップはどれくらいの距離、どの方向へ進もうか」という増分の求め方に反復法ごとの特徴があります。
 逆Hesse法では、極値(最尤推定値といえる)近傍では評価関数が2次関数で十分近似できることを利用しています。Hesse行列(評価関数の2階偏微分)の1/2をα、評価関数の1階偏微分の1/2をβとして、Σαδ=βを解いて増分δを求めます。しかし、「評価関数が2次関数で近似できる」という仮定があまり良くない場合には収束が悪くなります。そこで、反復初期は最急降下法を、極値近傍では逆Hesse法を利用するようにしたアルゴリズムが
Levenberg-Marquardt法(単にMarquardt法ともいう)です。最急降下法では、増分をδ=定数xβで求めますが、この定数を求める際にHesse行列の対角項を利用するところに特徴があり、非線形最小2乗問題の解法の標準といえます。
 大雑把というとこんな感じなのですが、Numerical Recipes in C(技術評論社)に詳細な解説がありますので参照してみてください。Levenberg-Marquardt法のソースコード(C言語)も載っていますので便利ですよ。わかりにくいところがありましたら、また回答いたしますので、この質問は閉じなくてもいいですよ。それでは、健闘を祈ります。

この回答への補足

お手数をおかけしております。

非常に、参考になります。
分かったような気にさせてくれますね。(^_^;)
ご教授頂いた 概要をガイドラインとして、もうすこし勉強してから再度質問させて頂きます。

補足日時:2000/09/18 00:08
    • good
    • 0

 こんばんは、TCMです。


 非線形のフィッティングですね。これは繰り返し操作が伴うので、質問#8014のように簡単には解けません。キャラクタ環境だけで、これの説明をするのはかなりやっかいです。ということでまず確認したいのですが、omaiさんはこのフィッティングについて何らかの情報はお持ちでしょうか?
 例えば「逆Hesse法で計算することは知っているのだが、Hesse行列の作り方がわからない」とか「Levenberg-Marquardt法で解きたい」とか。2でマトリクスという文言が出ていますので、何らかの資料を手元に持っておられるとお見受けするのですが。ちょっとポイントを絞ってみるというのはいかがでしょうか?

この回答への補足

TCMさん、お世話になっております。

資料は、webです。
泥縄ですが、コメントにありました"逆Hesse法"とか
"Levenberg Marquardt"をwebで検索して調べてます。
しかしみな難解で 直には理解できそうにありません(~_~;)

以前、円の方程式 (x-a)^2 + (y-b)^2 = c^2 を
A=[x^2+y^2]
B=[2x,2y,1]
C=[a,b,c^2-a^2-b^2]
とおいて
C=inv(B'*B)*B'*A
で近似した事がある程度なんで、同様に解けることを期待してたんです。

よろしければ、まず"逆Hesse法"と
"Levenberg Marquardt法"の大雑把な特徴を教えて下さい。

補足日時:2000/09/13 12:06
    • good
    • 0

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

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

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

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

Aベストアンサー

絶対値があるので、x<a1 と a1≦x<a2 と a2≦x の3通りの場合分け
が必要です。0<b1<b2ですから、与式の両辺に b1b2 をかけておいて
 b2|(x-a1)|>b1|(x-a2)| と変形してからやるといいです。
考えとしては絶対値の外し方[x<0のときlxl=-x,0≦xのときlxl=x]を使い
ます。
1.x<a1 のとき・・・x-a1もx-a2も負になるからマイナスをつけてはずす
   -b2(x-a1)>-b1(x-a2) →両辺に-1をかけてb2(x-a1)<b1(x-a2)
   これを解いて、 x<(a1b2-a2b1)/(b2-b1) ・・・(1)
   ここで a1 と (a1b2-a2b1)/(b2-b1) の大小関係を調べると
   両方に(b2-b1)をかけた式で a1(b2-b1)-(a1b2-a2b1)=-a1b1+a2b1
   =b1(-a1+a2)>0 となるので a1>(a1b2-a2b1)/(b2-b1) となります
   したがって、ここでの解は(1)の解でよいことになります。
2.a1≦x<a2 のとき・・・x-a1は正、x-a2は負だから
   b2(x-a1)>-b1(x-a2)
   これを解いて、x>(a1b2+a2b1)/(b1+b2)
   ここで、1.のときと同様にして (a1b2+a2b1)/(b1+b2) とa1,a2
   との大小関係を考えると、省略しますが、
     a1<(a1b2+a2b1)/(b1+b2)<a2 となり、
   ここでの解は (a1b2+a2b1)/(b1+b2)<x<a2・・・(2)
3.a2≦x のとき・・・x-a1もx-a2も正だから
   b2(x-a1)>b1(x-a2)
   これを解いて x>(a1b2-a2b1)/(b2-b1)
   同様に a2 と (a1b2-a2b1)/(b2-b1) の大小関係を調べると、また
   省略しますが a2>(a1b2-a2b1)/(b2-b1) となり
   ここでの解は a2≦x・・・(3)

以上、(1)~(3)が解となります。
各場合について、数直線をかいて考えるといいでしょう。

絶対値があるので、x<a1 と a1≦x<a2 と a2≦x の3通りの場合分け
が必要です。0<b1<b2ですから、与式の両辺に b1b2 をかけておいて
 b2|(x-a1)|>b1|(x-a2)| と変形してからやるといいです。
考えとしては絶対値の外し方[x<0のときlxl=-x,0≦xのときlxl=x]を使い
ます。
1.x<a1 のとき・・・x-a1もx-a2も負になるからマイナスをつけてはずす
   -b2(x-a1)>-b1(x-a2) →両辺に-1をかけてb2(x-a1)<b1(x-a2)
   これを解いて、 x<(a1b2-a2b1)/(b2-b1) ・・・(1)
   ここで a1 と (...
続きを読む

Q2直線 x/a+y/b=1, x/a+y/b=2(a>0, b>0)の

2直線 x/a+y/b=1, x/a+y/b=2(a>0, b>0)の間の距離を求めよ。

という問題の解説に、

2直線は平行だから、第一の直線上の点(1、0)を通る。よって、ここからbx+ay=2abまでの距離を求める

と、ありました。

なぜ(1,0)を通るのですか?

Aベストアンサー

誤記なんてレベルでは済まないですよ。
 A.2直線は平行である。
 B.第一の直線が点(1,0)を通る。or B'.第一の直線上のどこかの点を第二の直線が通る。
 C.AがB(またはB')の根拠になっている。
このうち正しいのはAだけです。

第一の直線は点(a,0)を通る。
また、2直線は平行だから、点(a,0)から第二の直線までの距離を求めればよい。
とでも書くのなら良いのですが、論理が滅茶苦茶ですね。

Q数学の線形代数の問題なのですが、n×nの2つのマトリックスA,Bがあり

数学の線形代数の問題なのですが、n×nの2つのマトリックスA,Bがあります。AとBの積の行列式はAの行列式とBの行列式の積となるようです。
すなわち、det(AB)=det(A)det(B) です。これは任意のn(1以上の整数)で成り立つのでしょうか。
テキストを見たのですが、省略されているようです。n=2の場合は、計算が簡単なので確かめられますが、高次だったらどうなるでしょうか。
よろしくお願いします。

Aベストアンサー

A=[a_i](n)(aiは列ベクトル)
B=[b_i,j](n,n)(b_i,jは実数)とすると
det(AB)=det(Σa_i・b_i,1,Σa_i・b_i,2,・・・,Σa_i・b_i,k,・・・,Σa_i・b_i,n)
これを
det(a_1,a_2,・・・,a_k+b,・・・,a_n)=det(a_1,a_2,・・・,a_k,・・・,a_n)+det(a_1,a_2,・・・,b,・・・,a_n)
(a_i、bともに列ベクトル)を用いて和の形に展開して、
det(a_1,a_2,・・・,α・a_k,・・・,a_n)=α・det(a_1,a_2,・・・,a_k,・・・,a_n)
を用いて行列Bの要素を行列式から外に出すと、
det(AB)=det(a_1,a_2,・・・,a_k,・・・,a_n)b_1,1・b_2,2・・・b_n,n+・・・
+det(a_σ(1),a_σ(2),・・・,a_σ(k),・・・,a_σ(n))b_σ(1),1・b_σ(2),2・・・b_σ(k),k・・・b_σ(n),n
+・・・
(ここでdet(a_σ(1),a_σ(2),・・・,a_σ(k),・・・,a_σ(n)=sgn(σ)det(a_1,a_2,・・・,a_k,・・・,a_n)を用いる)
=det(a_1,a_2,・・・,a_k,・・・,a_n)Σsgn(σ)_σ(1),1・b_σ(2),2・・・b_σ(k),k・・・b_σ(n),n
(行列式の定義がdet(B)=Σsgn(σ)_σ(1),1・b_σ(2),2・・・b_σ(k),k・・・b_σ(n),n なので)
det(AB)=detA・detB
となります。

A=[a_i](n)(aiは列ベクトル)
B=[b_i,j](n,n)(b_i,jは実数)とすると
det(AB)=det(Σa_i・b_i,1,Σa_i・b_i,2,・・・,Σa_i・b_i,k,・・・,Σa_i・b_i,n)
これを
det(a_1,a_2,・・・,a_k+b,・・・,a_n)=det(a_1,a_2,・・・,a_k,・・・,a_n)+det(a_1,a_2,・・・,b,・・・,a_n)
(a_i、bともに列ベクトル)を用いて和の形に展開して、
det(a_1,a_2,・・・,α・a_k,・・・,a_n)=α・det(a_1,a_2,・・・,a_k,・・・,a_n)
を用いて行列Bの要素を行列式から外に出すと、
det(AB)=det(a_1,a_2,・・・,a_k,・・・,a_n)b_1,1・b_2,2...続きを読む

Q∫〈x〉=b^x /( a^x+b^x )   (1<a<b) のときlim_

関数∫(x) = b^x /( a^x+b^x)  (1<a<b)のときの 
lim_(x→∞)∫(x)、lim_(x→-∞)∫(x)求めよ。

関数の記号の入力の仕方が解らなくて間違えているかもしれません。ごめんなさい。よろしくお願いします。手書きですが、問題画像付けました。

Aベストアンサー

x→∞のとき
分母と分子をb^xで割る。
|a/b| < 1 だから(a/b)^xが0に収束することを利用

x→-∞のとき
分母と分子をa^xで割る
(b/a)^xが0に収束することを利用。

こういう問題はより大きくなる方で割ればいい。

Qy=A(x-B)/(x-C)についてのカーブフィッティング

24個のデータ(x,y)について、y=A(x-B)/(x-C)という関数を用いて近似曲線を描いたときのA,B,Cの値を知りたいのですが、どうすればよいのでしょうか?

Aベストアンサー

よく知られたn次式の最小二乗法のように一般的な方法でさっと出来るものではない場合と理解します。理論的な最小二乗法の式が構築出来る場合もあるとは思いますが、泥臭くやる方が結局得な場合もあるので下記に示します。EXCEL前提です。

私はこういう場合、先ず実測データ(meas) [x, y] を B, C 列にプロットした後、A, B, Cをパラメータとするために例えばA2, A4, A6に適当な3ヶの数字を入れます。A1, A3, A5に"A" "B" "C"と書いておきます。このバラメータ値を使って各x値に対する理論値(theory)をD列に計算します。E列に実験と理論の差の二乗和 (C2-D2)^2 , (C3-D3)^2, ... を計算します。その合計値 [=sum(E2:E100)] をA8に入れ、A7に"SUM"とでも書きます。さらに実験と理論の曲線を同じグラフ上に書きます。
パラメータ A, B, C を変更しますと、SUMは増減します。グラフも見ながら、差が減少する方向にA, B, Cを繰り返し変更していきます。
3ヶもパラメータがあるので場合によっては大変ですが、その関数形ですと、さほど悩まずに収束させることが出来るだろうと思います。EXCELですとSUMの有効数字を10桁でも可能ですから、相当の精度でパラメータを決めることが出来ます。

よく知られたn次式の最小二乗法のように一般的な方法でさっと出来るものではない場合と理解します。理論的な最小二乗法の式が構築出来る場合もあるとは思いますが、泥臭くやる方が結局得な場合もあるので下記に示します。EXCEL前提です。

私はこういう場合、先ず実測データ(meas) [x, y] を B, C 列にプロットした後、A, B, Cをパラメータとするために例えばA2, A4, A6に適当な3ヶの数字を入れます。A1, A3, A5に"A" "B" "C"と書いておきます。このバラメータ値を使って各x値に対する理論値(theory)をD列に計...続きを読む


人気Q&Aランキング

おすすめ情報