最小二乗法？

i 個の測定点 (x[i],y[i]) を，最小二乗法などを用いて下記の式にフィッティングさせようと考えています。Visual Basic で作成した測定プログラムの中で使用したいのですが，具体的にどのようなアルゴリズムでフィッティングを行えばいいのか分かりません。

Y = A * sin(X - C)^2 + B

実測する x[i] の範囲は狭く，例えば -15°～ +15°まで 0.2°毎の計 151 プロット，といった感じです。そして定数 A，B，C の内，最も高い精度で求めたい定数は C です。測定の段階で x の範囲を狭めているのは，正確な C （通常 1°未満）を求めるためです。

この測定は x[i] にはほとんど誤差が含まれませんが y[i] には誤差があります。y[i] の含まれている誤差は試料によってまちまちなので，一概には言えません。目視ではほとんど誤差が分からない綺麗なカーブの場合，逆に目視で辛うじて下に凸の曲線が分かる程度の場合，どちらもあり得ます。

考え方だけでも構いませんので，どうかご教授下さい。よろしくお願いいたします。

No.11 ベストアンサー 回答者： stomachman 回答日時： 2001/12/13 05:36

No.9のコメントについて、概要だけです。

●自由度については、モデルにパラメータが３つあるんですから、m=3が正解でしょう。
なお、重みの式から見ると、ワンセットの計測においては、どのサンプルのノイズも標準偏差が同じであると仮定して良いようですね。

●Cが含む誤差の評価方法について。
モデル方程式を、得られたA,B,Cの周りで線形近似したものを考えます。j回目の計測セットでサンプル点k (k=1,2,....,151; x=0.2(k-1)-15)で得た測定値と、モデルとの重み付き残差の縦ベクトルをε[j]とします。そして、
∂f/∂A = (sin(0.2(k-1)-15-C))^2
∂f/∂B = 1
∂f/∂C = -2A[j](sin(0.2(k-1)-15-C))(cos(0.2(k-1)-15-C))
(k=1,2,.....,151)
というヤコビアン行列J[j]を作ります。計測のセットjごとにA[j],B[j]の値は異なることにご注意下さい。すると、
ε[j]= w[j]J[j] (ΔA,ΔB,ΔC)'
（ここに'は転置、w[j]はdiag[√W(j),√W(j),.....]という対角行列で、計測のセットjごとに異なる。）
と書けます。
　さらに、j=1,2,...,N についてε[j]を縦につないだベクトルをεとし、同じくJ[j]を縦につないだ行列をJとし、wはw[j]を右下に繋いだ対角行列diag[√W(1),√W(1),....., √W(2),√W(2),....., ...... ,√W(N),√W(N),.....]　（全部で151N個の対角要素が並びます。）とすれば、j回の計測を全部まとめて
ε= wJ (ΔA,ΔB,ΔC)'
と表現できるでしょう。そしてεの要素はどれも同じ正規分布に従うと仮定できます。（このために重みwを掛けたんです。）
　この線形問題において、εからΔCへの「誤差の伝播」を考えれば良い。これは教科書に書いてあると思いますよ。ΔCが正規分布に従うと仮定して、（平均は0の筈ですから）その標準偏差を求める問題です。
　ΔCの分布に対してどんな危険率（つまり信頼性限界）を設定するかは結果の利用目的によりますが、たとえばC±2(ΔCの標準偏差)、という風に表すことが出来ます。

この回答へのお礼

お礼が遅くなり，大変申し訳ございませんでした。
自分でも本を読みつつ，stomachman さんのご回答について少しずつ理解を進めておりました。
私の中で，まだ完全に問題を解決したわけではありませんが，長い間スレッドを開けておくのも失礼かと思い，とりあえずこのスレッドは一旦締めさせて頂きます。もう少し勉強をすすめ，分からない部分がまとまりましたら，またお伺いすることがあるかもしれませんが，その時もまたどうかよろしくお願いいたします。

お礼日時：2001/12/20 23:39

No.10 回答者： stomachman 回答日時： 2001/12/13 04:16

No.9のコメントまだ読んでません。

その前にまず、stomachman我流の線形解法（超平面法）で解くとどうなるかをご紹介しておくことにします。

f(x) = A(sin(x-C))^2 + B
ですから、
f(x) =(A/2+B)-(A/2)cos(2x-2C)
すなわち
f(x) = (A/2+B)-(A/2)cos(2C)cos(2x)-(A/2)sin(2C)sin(2x)
と書き換えることができます。ここで
P=(A/2+B)
Q=-(A/2)cos(2C)
R=-(A/2)sin(2C)
とおくと、
f(x) = P + Q cos(2x) + R sin(2x)
ですから、パラメータP,Q,Rを線形最小二乗法で決めることができます。

で、P,Q,RからA,B,Cを計算するのは簡単です。
A=2√(Q^2+R^2)
B=P-A/2
そしてCは二通りの答が出てきます。
C=Arcsin(-2R/A)/2, C=Arccos(-2Q/A)/2

両者を比べてみてください。余程ノイズが大きくない限り、非常に良く一致することと思います。

No.9 回答者： stomachman 回答日時： 2001/12/11 03:47

No.8のコメントを拝見しました。

CはA,Bと断じて無関係。ということですね。

だとすれば、こういう方法が考えられます。

実験をM回繰り返したとして、それぞれについて、Cを定数とする最小二乗法の正規方程式を作ります。つまり、M組の線形最小二乗法の問題が作られる訳です。問題jの残差二乗和をS(C,j) (j=1,2,....,M)と書くことにしましょう。

これら全ての問題に対して、同じCを与えて、S(C,j) (j=1,2,....,M)を計算します。
実験ごとに誤差の程度が違いますから、それを勘案した重みW(j) (j=1,2,....,M)を用意して、
T(C) = ΣS(C,j) W(j)　　(j=1,2,....,M)
を計算します。
このT(C)を最小にするようなCを、例えば四分法で求めるのです。

重みWをどう与えるべきか、については、これはもう既に教科書で勉強されたことを応用すれば解決できそうに思います。

この回答への補足

いつもお世話になり，大変ありがとうございます。
今まで教えていただいたことにつきましては，大筋で理解したつもりです。
しかし以下の３点につきまして，依然疑問が残っております。
平行して私も勉強いたしますが，何かヒントでもいただければ幸いです。

※３における標準偏差および重みの計算方法は正しいでしょうか？
特に，自由度 n-m の与え方が分からず m=2 としておりますが，これで良いのでしょうか。

※４にて C の信頼限界を算出する具体的な式が分かりません。

※５にて，重みを考慮した標準偏差を計算したいのですが（信頼限界を算出するため），
具体的な計算の式が分かりません。


--------

■ クロスニコル測定におけるデータ解析方法

測定は M 回行う。M 回中 j 回目の測定結果は x[i,j]，y[i] に格納する。

解析を始める前に，まず最初に以下のようなサブルーチンを作成する。

｜最小残差和を計算するサブルーチン（Calc_S と命名）
｜１．引数として C，および j を受け取る。
｜２．この C を用い，残差二乗和が最小となる A，B を，線形最小二乗法で求める（※１）。
｜３．残差二乗和 S(C,j) を返す。

このサブルーチンを作った上で，次のように場合わけを行って解析を行う。

｜（１）測定を 1 回のみ行ったときのアルゴリズム
｜１．初期値として，Cm を x[i] の最小値，Cp を x[i] の最大値，Cc = (Cm + Cp) / 2 とおく。
｜２．C を計算する都度 Calc_S をコールしながら，S(C) が最小となる C を求める（※２）。
｜３．標準偏差を計算し（※３），カイ二乗検定によって C の信頼限界を求める（※４）。

｜（２）測定を M 回行ったときのアルゴリズム
｜１．Calc_S をコールし，j 回目の測定の残差二乗和 S(C,j) を求める。
｜２．標準偏差を基に，j 回目の測定の重み W(j) を求める（※３）。
｜３．T(C)=ΣS(C,j)W(j) とおき，T(C) が最小となる C を求める（※２）。
｜４．重みを考慮した標準偏差を計算し（※５），C の信頼限界を求める（※４）。

※１　線形最小二乗法の計算方法
∂Σε[i]^2 / ∂A = 0 ∧ ∂Σε[i]^2 / ∂B = 0 を解く（連立方程式）。
A = …， B = … として求まる。

※２　四分法のアルゴリズム
１．Cm，Cc，Cpの，計3点について，Calc_S をコールして S(C) を求める。
２．(Cc+Cm)/2，(Cc+Cm)/2 の計 2 点について，Calc_S をコールして S(C) を求める。
３．5 点のうち，もっとも S が小さい時の C を新たな Cc，その両端を新たな Cm，Cp とする。
４．C が所望の精度で求まるまで，２から繰り返す。

※３　標準偏差σ(j) および重み W(j) の計算方法
標準偏差 σ(j) = Sqrt (S(C,j) / (n - 2))
重み W(j) = 1 / σ(j)^2

※４　C の信頼限界を求めるアルゴリズム
カイ二乗検定（計算方法不明）

※５　重みを考慮した標準偏差の計算方法
（計算方法不明）

補足日時：2001/12/12 21:53

No.8 回答者： stomachman 回答日時： 2001/12/09 04:10

siegmund先生のご指摘の通り、ドリフトの問題を解決できれば一番すっきりしますね。

どうにもならない場合、むしろ例えば「Aがある一定範囲に入った場合のCだけを採用し、他の実験は捨てる」ぐらいの方が確かなのかもしれません。

また、No.6に書き間違いがあるので訂正です。

●「四分法」の説明は以下のように訂正します。
(1) Cm=-1, Cc=0, Cp=1 とします。S(Cm), S(Cc), S(Cp)を（線形最小二乗法で）計算します。すると、S(Cm), S(Cc), S(Cp)のうちで一番大きいやつがもしS(Cc)だったら、これはどうも状況がおかしいので、この方法は失敗です。（Cの範囲を大幅に見誤ったとか、モデルが全然データに合っていないと思われます。）まともな状況なら、S(Cm)かS(Cp)が最大になります。
(2) さらにS((Cc+Cm)/2)、S((Cc+Cp)/2)を計算します。
(3) そして、Cm, Cc, Cpを変更します。５個の値Cm, (Cc+Cm)/2,Cc, (Cc+Cp)/2, Cpのなかから引きつづく３つを選んで新しいCm, Cc, Cpの値とし、そのときS(Cm), S(Cc), S(Cp)（これらはもう既に計算してありますね）のうちのS(Cc)が最小になるように、引きつづく３つを選ぶのです。
そして(2)に戻る。

●もうひとつ、細かいことですが
＞　C=-1～1を2^N分割　たとえば32768の等間隔に
は
＞　C=-1～1を2^N分割　たとえば65536個の等間隔に
の誤りでした。

この回答への補足

大変申し訳ございません。うっかり siegmund さんの名前を呼び捨てにしてしまいました。コピー＆ペーストしている内に，どこかで消えてしまったようです。本当に，大変失礼いたしました。

補足日時：2001/12/11 01:13

この回答へのお礼

ご回答ありがとうございます。大変申し訳ございませんが，stomachman さんへのお礼と，siegmund へのお礼とを兼ねさせていただきます。また，返事が遅くなり大変失礼いたしました。

> Aが経時的にドリフトしてしまう、という事情があったのですね。Bもおそらく動くのでしょう。

A，B，C につきまして，それぞれの相関関係は全くありません。強いて言うならば，A と B との間については，多少の相関がある可能性はあります。ただし C だけは完全に独立です。

今回取り扱っている式は物質の旋光度測定のモデルでして，

「光源 → 偏光子 → 旋光性を持つ試料 → 検光子 → ホトマル」

という光学系における，検光子の回転角度（x[i]）と検出する光の強度（y[i]）との関係を示しています。各定数については，A は入射光の強度に比例する値（入射光の安定性によりドリフトを持つ），B はホトマルの暗電流などに起因する一定のノイズ，C は試料による旋光度です。よって C は，A や B とは完全に独立に振る舞うものと考えられます。この測定において，求めたい値は C だけであり，A や B については求めても決して使われない値です。

厳密に言いますと，今回の測定は旋光度測定の原理を応用した別の測定でして，光源にはパルスレーザーを用いております。そして，この出力不安定が解析を困難にしております。単パルスごとの強度のバラツキについては，一定時間積算をとれば解決する問題だと思っておりますが，長い周期での出力強度のドリフトについては，目下の所解決策が見あたりません。光学系の途中にビームスプリッターを設置して，A の値を別途モニターできればよいのですが，すぐに光学系を変えられるという訳ではありません。

結局の所，A がドリフトしてしまうために，もはやサインカーブにフィッティングしても意味がないかと思い，例えば y = ax^2 + bx + c とかにフィッティングしてもいいのかな？　などと思ってしまったりします。しかし |x| < 15 の範囲で二次関数では，いかにもモデルから外れすぎているとも思い，正直なところ，そこら辺の数学的（？）な判定については非常に悩むところがあります。

> 必要なのはCを振ることだけです。

振るのが C だけで良いとなると，さらに stomachman さんのおっしゃる四分法のアルゴリズムを考えると，四分法によって C を求めるのが最も簡便そうですね。結局の所，以下のようなアルゴリズムを組もうと考えているのですが，これで問題ありませんよね？

１．まず cos(C)=1，sin(C)=C とおいて線形最小二乗法で解く。（大雑把に A，B，C が求まる。）
２．四分法によって，正確な C を求める。
３．線形最小二乗法で A と B を求める。
４．所望の精度で C が求まるまで，再び２に戻る。
５．残差二乗和を元にカイ二乗検定を行い，C の信頼限界を求める。

> 四分法を高速でやるためには、

高速化のアルゴリズムを読んでいて，なるほどと思いました。但し今回は，測定範囲と測定間隔とを自由に設定できるようにプログラムを組んでしまったので，これは後々の参考とさせていただきます。

> １回のfittingに於いて求められたCの値に含まれる誤差の評価方法についてでした。

stomachman さんが過去の回答でおっしゃっていた「最小二乗法による実験データ解析」を入手し，早速読んでおりました。一回のフィッティングにおける C の信頼限界を求めるには，最適な A，B，C を求めた後，その時の残差二乗和を元にカイ二乗検定を行えばいいんですよね。「最小二乗法による実験データ解析」の§3.5 を参考に致しました。

お礼日時：2001/12/10 03:32

No.7 回答者： siegmund 回答日時： 2001/12/09 00:56

siegmund です．

> 線形最小二乗法 → 解析的に解が求まる方法
> 非線形最小二乗法 → 値を変えながら繰り返し計算する方法
stomachman さんの書かれているとおりで，
正規方程式が，決めるべき未知係数(A,B,C，など)について線型の場合が
線型最小自乗法です．
この場合は連立１次方程式を解けばいいのですから，解析的に解が求められます．
線型でない場合は解析的な解が求められる場合もありますが，
たいていはうまくいきません．
それで，いろいろ数値的方法が必要になるわけです．
A,B,C，などの値の見当をつけてそのまわりで連立方程式を線型化して解くのが
ニュートン法です．

stomachman さん：
> 必要なのはCを振ることだけです。
> Cを決めれば、その状態で残差二乗和S(C)を最小にするA,Bは線形最小二乗法
>（つまり正規方程式が連立１次方程式になって、これを解くこと）
> で簡単に求められるからです。
おっしゃる通りです．
私のやり方はちょっと頭が悪かったです．

A のドリフトはちょっと深刻なように思えます．
A がドリフトするように見えると言うことは，
時間が経過すると同じ X の値に対しても Y の値が推定される誤差以上に
系統的に変化してしまうんでしょうか？
そうすると，Y は X のみで決まるのではなくて，
まだ他に何かに対する(例えば，温度とか)依存性も持っていることになります．
今までの議論は Y が X のみで決まることが前提でしたから，
前提が怪しくなってしまいます．
最小自乗法はランダムな誤差に対する処理方法ですから，
系統的にずれるような話はまた別ですね(stomachman さんも書かれています)．
数学的検討よりは，測定系を検討してドリフトの原因を探る方が
本質のように思われます．
たぶん，トライされてなかなか解決がつかないんだろうとは思いますが．．．

No.6 回答者： stomachman 回答日時： 2001/12/08 20:31

Aが経時的にドリフトしてしまう、という事情があったのですね。

Bもおそらく動くのでしょう。だとすれば、仰るとおり、実験毎にA,B,Cを求めるのが良さそうです。

一般に、非線形最小二乗法というのは、正規方程式（残差二乗和が極値をとる、ということを表す方程式）がA,B,Cの一次式にならない場合に適用されます。そういう場合には、正規方程式のA,B,CをA+ΔA, B+ΔC,C+ΔCで置き換え、A,B,Cとして適当な「推定値」を（定数として）入れ、ΔA, ΔB, ΔCの２次以上の項を全て無視してしまう。そうするとΔA, ΔB, ΔCを変数とする連立１次方程式が得られます。これを解いて、A+ΔA, B+ΔC,C+ΔCを新しい「推定値」とする。これを繰り返して「推定値」を改良していく訳です。

ご質問の場合には、Cだけについてこの処方が必要になる。ですから、No.5 siegmund先生の
＞あるいは，A0，B0，C0 の周囲で A,B,C を適当なステップで振って
＞残差の自乗和を直接計算する手もあります．
＞いわば，３次元の表にするわけです．
＞今の範囲に残差の自乗和の極小値があれば，その周囲でステップを細かくしてまた振る．
＞A,B,C の値を 10個ずつトライしますと，
＞１つのステップについて 1000 通りの場合の残差の自乗和を計算するわけですが，
＞今のパソコンならあっという間でしょう．

については、必要なのはCを振ることだけです。Cを決めれば、その状態で残差二乗和S(C)を最小にするA,Bは線形最小二乗法（つまり正規方程式が連立１次方程式になって、これを解くこと）で簡単に求められるからです。
Cを有効数字４桁で決めたければ少なくとも数万通りについて調べることになります。これを端から順に試すのでは効率が悪いですね。
小数点以下１桁目について探す
小数点以下2桁目について探す
　：
という風に絞り込んでいくと速くなります。

これを徹底して行うのが以下の方法です。

(1) Cm=-1, Cc=0, Cp=1 とします。S(Cm), S(Cc), S(Cp)を（線形最小二乗法で）計算します。すると、S(Cm), S(Cc), S(Cp)のうちで一番大きいやつがもしS(Cc)だったら、これはどうも状況がおかしいので、この方法は失敗です。（Cの範囲を大幅に見誤ったとか、モデルが全然データに合っていないと思われます。）まともな状況なら、S(Cm)かS(Cp)が最大になります。
(2) さらにS((Cc+Cm)/2)、S((Cc+Cp)/2)を計算します。そして、Cm, Cc, Cpを変更します。５個の値Cm, (Cc+Cm)/2,Cc, (Cc+Cp)/2, Cpのなかから引き?く３つを選んで新しいCm, Cc, Cpの値とし、そのときS(Cm), S(Cc), S(Cp)（これらはもう既に計算してありますね）のうちのS(Cc)が最小になるように、引き?く３つを選ぶのです。
そして(1)に戻る。

　これで、（２進法で）１桁ずつCを絞り込んでいくことができます。stomachmanはこの方式（「四分法」と呼んでいます）を極めて多量のデータ処理を行うような、色々な目的に実際に利用しています。

　ついでに、四分法を高速でやるためには、sin(x[i]+C)をテーブルルックアップで計算できるようにしておくと良いです。つまりC=-1～1を2^N分割　たとえば32768の等間隔に分割します。-1, -1+1/32768, -1+2/32768, ..... , 1 という32769個の数値の中から、S(C)を最小にするCを探そうという訳です。一方、x[i]のサンプリング間隔も1/32768の整数倍にしておきます。そして、θ=-16～16の範囲についてθを1/32768刻みで変えて(sin(θ))^2を計算した表を作っておけば、毎回sinの計算をしなくても済む訳です。

＞例えば No 2 での補足で，私は数回の測定から信頼限界を求めたと書きましたが，
＞この時は数回測定した結果から標準偏差，平均誤差をもとめて，Student の t 分布
＞を用いて，○○%信頼限界といった値を計算しました。stomachman さんの
＞おっしゃる方法が具体的にはどういう方法なのか，よく分かりません。

No.4でstomachmanが説明したのは、１回のfittingに於いて求められたCの値に含まれる誤差の評価方法についてでした。yにランダムなノイズがあるために、ノイズが違えばCも違ってくる。そのCのばらつきは（僅かなので）ノイズの標準偏差と比例する標準偏差をもつ正規分布で近似できます。詳しくは教科書をご参照下さい。

さて、ご質問で言う「信頼限界」は、ドリフトを起こす系において繰り返し求めたCの値の分布の話だったんですね。これはまた別の話です。
信頼性限界を云々する以前に、Cの値がどう変化するのか。たとえば、AとCの相関はないか、などを調べて、Cの値が本当にランダムなのかどうかをチェックする所から始めなくてはいけません。（AやBとCに緩やかな相関が認められることはかなりありそうです。）
　その上で、一体何を測りたいのかを良く考え直して、たとえばA,Bを使ってCを補正するなど、A,BとCが無相関になるようにしてやりますと、初めてCの信頼性限界について議論が可能になります。

No.5 回答者： siegmund 回答日時： 2001/12/06 14:09

siegmund です．

出馬依頼だけでは何なので
(いや，出馬依頼は反則らしいんで，決してそういうつもりでは...)，
私もちょっとだけ．

C が 1゜未満というのでしたら，
stomachman さんの sin の展開で C の１次まででもかなり見当がつきそうです．
これでしたら，普通の線型最小二乗法ですし，
パラメーターの数の心配はいりません．
で，A,B,C の値の見当がついたら(A0，B0，C0 とします)，
A = A0 + δA
B = B0 + δB
C = C0 + δC
として，δA，δB，δC について線型化して最小二乗法で
δA，δB，δC を求める．
で，修正された A,B,C を A0,B0,C0 として同じことを繰り返し，
δC が所用の精度(10^(-4) 位？)になるまでやります．
非常に運が悪いなどでなければ，数回の繰り返しでかなりの精度が得られます．

あるいは，A0，B0，C0 の周囲で A,B,C を適当なステップで振って
残差の自乗和を直接計算する手もあります．
いわば，３次元の表にするわけです．
今の範囲に残差の自乗和の極小値があれば，その周囲でステップを細かくしてまた振る．
A,B,C の値を 10個ずつトライしますと，
１つのステップについて 1000 通りの場合の残差の自乗和を計算するわけですが，
今のパソコンならあっという間でしょう．
力業でエレガントではありませんが，こういう手もあります．

この回答へのお礼

ご回答をありがとうございました。siegmund さんのご回答と stomachman さんのご回答を総合し，以下のような手順で，今回の質問の内容を解決できるのではないかと考えています。

１．まず cos(C)=1，sin(C)=C とおいて線形最小二乗法で解く。
２．１の結果を用い，線型最小二乗法でδA，δB，δC を求める。
３．２を繰り返して C の制度を高める。
４．∂S/∂C を求めて C の信頼限界を求める。

この手順を踏む上で，まずは１と２で用いる線形最小二乗法について，本を探して自分でも勉強したいと思います。

> A,B,C を適当なステップで振って残差の自乗和を直接計算する手もあります．

この方法はわかりやすいですね。このアルゴリズムならすぐにでもプログラミングできそうです。
この度は，大変ありがとうございました。

お礼日時：2001/12/07 03:06

No.4 回答者： stomachman 回答日時： 2001/12/06 03:58

あやや。

stomachman読み違えたんですね。
y[i] = A (sin(x[i] - C))^2 + B+ε[i]
から、 |C| < 1°を求める。なお、|x[i]|＜15°
こういう問題でしたか。

非線形最小二乗法でCを追いつめていくのは勿論、真っ当な方法です。つまり、Cをちょっと変えては線形最小二乗法でA,Bを決め、残差二乗和S(C)=Σ(ε[i] ^2)が小さくなるようにCを修正していく。

でも、それじゃ回答として詰まらないので、まっとうでないやり方を考えましょう。
|C| が小さいのが特徴ですね。だから
sin(x[i] - C)= sin(x[i]) cosC - cos(x[i]) sinC
cosC=1-(C^2)/2+(C^4)/4! -…
sinC = C-(C^3)/3! + (C^5)/5! -…
を利用する手はどうでしょう。
sin(x[i] - C)≒sin(x[i]) - Ccos(x[i])
としてしまうんです。すると、
y[i] ≒ A (sin(x[i]) - Ccos(x[i]))^2 + B+ε[i]
= A ((sin(x[i]))^2 + A(C^2)(cos(x[i]))^2 -2ACsin(x[i])Ccos(x[i])+ B+ε[i]
ですから、
a=A,
b=A(C^2),
c=-2AC,
d=B
とおけば、線形最小二乗法の問題に帰着します。
a,b,cからAとCを決める訳で、得られるパラメータの個数が１個余計ですけど、そこは適当にごまかすことにしましょう。
　もちろん、こうやって求めたCを利用して、さらに非線形最小二乗法で追いつめていくというのもアリです。

> 測定を数回繰り返して C の平均値と信頼限界を求めよう
同じ対象を何度も測定する、という意味でしょうか。だったら測定値を全部ぶち込んで一度に解いた方が良いですよ。ベクトルx[i]の異なる成分に同じ値が何度も現れたって一向に構いません。測定毎にノイズ（yの測定誤差）の標準偏差が異なるのであれば、その逆数を重みにした重み付き最小二乗法を行うのが適切かと思います。

> フィッティングにおける相関係数は，どのようにして求めればよいのでしょうか。

んーと。Cの誤差を評価するのが目的なら、Cを微小に動かした時に残差二乗和がどう変化するか、つまり∂S/∂C を求めて、残差二乗和の極小値から、（たとえばε[i] が正規分布だと仮定した場合の）信頼性限界にぶつかるにはCがどれだけ動けるかを調べればよいのです。
　相関係数というのは、この場合モデルA (sin(x[i] - C))^2 + Bとy[i]との相関ですね。普通に考えれば良いと思いますよ。

この回答への補足

ご回答ありがとうございます。最小二乗法については y = ax + b の a と b を求める以外には用いたことがなく，線形／非線形と接頭語のつく最小二乗法には全く馴染みがありません。stomachman さんのご回答の文面から判断すると，

線形最小二乗法 → 解析的に解が求まる方法
非線形最小二乗法 → 値を変えながら繰り返し計算する方法

という漠然とした認識をもったのですが，正しいでしょうか？

> 何度も測定する、…中略…だったら測定値を全部ぶち込んで一度に解いた方が良いですよ。

この件についてですが，実は y = A (sin(x - C))^2 + B 式における A の値に問題がありまして，測定系上の都合により，長時間測定中していると A の値はゆっくりと約±20%くらいの変動を持つのです。この A の変動による影響を最小限に抑えるため，x の測定順序を 0, 0.2, -0.2, 0.4, -0.4, … という具合に，ほぼ左右対称になるように測定したりしています。この様な状況下でも，全部の測定結果を同時に解析した方が良いのでしょうか？　正直言いまして，私は数学的な処理について全く自信がありません。

> 残差二乗和の極小値から、信頼性限界にぶつかるにはCがどれだけ動けるかを調べればよい

感覚的には雰囲気が分かりますが，具体的にどういった計算をすればいいのか見えてきません。例えば No 2 での補足で，私は数回の測定から信頼限界を求めたと書きましたが，この時は数回測定した結果から標準偏差，平均誤差をもとめて，Student の t 分布を用いて，○○%信頼限界といった値を計算しました。stomachman さんのおっしゃる方法が具体的にはどういう方法なのか，よく分かりません。

どちらにせよ，私はもう少し勉強しないといけないことがよく分かりましたので，少々勉強をして出直してきます。この度は，ご丁寧にありがとうございました。

補足日時：2001/12/07 03:03

No.3 回答者： kony0 回答日時： 2001/12/06 01:49

stomachmanさんの表記を借りて、

S(A,B,C)=Σ(ε[i] ^2)という、R^3(A,B,C)→R^1(S)の関数を考えて
minS(A,B,C)という（非線形）最適化の問題を解いてしまえばいいのではないでしょうか？
つまり、最適化問題を解く数値解法としてニュートン法、準ニュートン法などを使うもよし、
∂S/∂A=0,∂S/∂B=0,∂S/∂C=0の連立方程式を解くのもよし、
と思うのですが。（目的関数はろくに考えてないですけど凸関数っぽいし）

＃もともと最小二乗法の回帰係数もこの考えで導出されてますよね？というだけなのです。まともに考えず直感勝負で。（汗）

この回答へのお礼

ご回答ありがとうございます。大変申し訳ございませんが，私の勉強不足が原因で，おっしゃっている内容をフォローしきれておりません。とりあえず，「ニュートン法」および「準ニュートン法」というキーワードで，図書館で文献をあたってみます。

お礼日時：2001/12/07 03:03

No.2 回答者： stomachman 回答日時： 2001/12/05 21:41

呼ばれて飛び出てジャジャジャジャン。

stomachmanです。

モデルの式がちょっと曖昧ですね。
(a) y[i] = A (sin(x[i] - C))^2 + B+ε[i] 　(x[i]= (-15+ 0.2(i-1))度、i=1,2,....,151）
ということでしょうか、いやそれとも
(b) y[i] = A sin((x[i] - C)^2) + B+ε[i] 　(x[i]= (-15+ 0.2(i-1))度、i=1,2,....,151）
なのかなあ？
ま、それはどうでもいいや。

x[i]の測定の範囲を±15度に絞っているということは、|C|＜15度、などと予想されていると思われます。そしてCの精度として１度程度の誤差を許すのであれば、これはもう簡単です。
たとえばC=-15,-14,......,15度のそれぞれについて、
S(C)=Σ(ε[i] ^2)　　　(i=1,2,.....,151)
を作ってみて、S(C)が一番小さくなるCを見つければよいのです。Cが決めてあれば、これは線形最小二乗法の問題ですから、実に簡単。
　そうやって見つけた、S(C)が最小になるCをminCと書くことにします。さらに
<minC-1,S(minC-1)>, <minC,S(minC)>, <minC+1,S(minC+1)>の３点を通る放物線の頂点を探せば、一層高精度でCの最適値が決まりますね。

　さらに、こういったデータが何組も次々と与えられるというのでしたら、予めsinの表を0.2度刻みで作っておけば、計算の手間が随分省けます。

この回答への補足

ご回答，どうもありがとうございます。

> モデルの式がちょっと曖昧ですね。

確かに曖昧でした。申し訳ございません。フィッティングさせたい関数は (a) の方です。

> そしてCの精度として１度程度の誤差を許すのであれば、

これも，私の言葉が足りなくてご迷惑をお掛けしました。C の範囲は通常 |C| < 1°です。最小二乗法によって少なくとも C を小数点以下４桁くらいまで求め，測定を数回繰り返して C の平均値と信頼限界を求めようと考えております。測定範囲を -15°～ +15°としているのは，ノイズが多い場合でもちゃんと下に凸のカーブを得るためです。測定範囲を狭めてしまうと，ノイズが多かった場合にカーブが見えてきません。

補足ついでに，もう一つ質問を追加してもよろしいでしょうか？　このフィッティングにおける相関係数は，どのようにして求めればよいのでしょうか。残差平方の平均値から求めるのでしょうか？

どうかよろしくお願いいたします。

ＰＳ。
stomachman さんの過去の回答履歴を見て，かなり感動しました。こちらの方もかなり参考になりそうです。

補足日時：2001/12/06 01:31