２変数関数の近似曲線

締切済

質問者：jirotaku
質問日時：2007/10/27 02:13
回答数：6件

計測データの２変数（X、Y）から、一つの値Z(X,Y)が決まる関数の近似曲線を作りたいのですが、方法（ツール）を教えて頂けないでしょうか？有料ソフトでも構いません。
「これを使えば出来るのでは？」ではなく、実際にした事があり「これを使えば出来た」という情報が嬉しいです。

通報する

この質問への回答は締め切られました。

質問の本文を隠す

回答 (6件)

最新から表示
回答順に表示

No.6

回答者： stomachman
回答日時：2007/10/29 23:41

ANo.5 のstomachmanです。

　B,つまりfittingをやりたいということですね。
　測定を行った点を(Xj, Yj)とし、そこでの測定値をMjとします。(j=1,2,…, J)

　さて、モデルとして、例えば多項式
Z(X,Y) = a[1]+a[2]X+a[3]Y+a[4]XY
を考えているとしましょう。X,Yに具体的に測定点の座標(Xk, Yk)を代入した式
Z(Xj, Yj) = a[1]+a[2]Xj+a[3]Yj+a[4]XjYj
は、未知数a[1]～a[4]を含んでいます。一方、1, Xj, Yj, XjYj はいずれも定数です。だから、右辺は変数a[1]～a[4]の一次式になっています。

ここで、変数の個数よりも測定点の個数Jの方が多いことが必要です。

　測定値とのずれεjをj=1,2,…, J について
εj = Z(Xj, Yj) - Mj
と定義します。そして、それらの二乗和
E = Σ(εj ^2) 　　（Σはj=1,2,…,Jについての和。以下同様)
を最小にするようなa[1]～a[4]を決定することを考えます。そのようなa[1]～a[4]は
∂E/∂a[n] = 0　　(n=1,2,…,4)
という連立方程式を解くことによって得られます。
さて、
∂E/∂a[n]
=(∂/∂a[n]) Σ(εj ^2)
=　Σ(∂/∂a[n]) (εj ^2)
=　2Σεj (∂εj /∂a[n])
= 2Σ((Z(Xj, Yj) - Mj)(∂Z(Xj, Yj) /∂a[n]))
であり、たとえば
∂Z(Xj, Yj) /∂a[4] = (XjYj)
です。
従って、
2Σ(Z(Xj, Yj) - Mj)=0
2Σ(Z(Xj, Yj) - Mj)Xj=0
2Σ(Z(Xj, Yj) - Mj)Yj=0
2Σ(Z(Xj, Yj) - Mj)(XjYj)=0
という式ができ、これをさらに展開して整理すると、
a[1]Σ1+a[2]ΣXj+a[3]ΣYj+a[4]ΣXjYj=ΣMj
a[1]ΣXj+a[2]Σ(Xj^2)+a[3]Σ(XjYj)+a[4]Σ((Xj^2)Yj)=Σ(XjMj)
a[1]ΣYj+a[2]Σ(XjYj)+a[3]Σ(Yj^2)+a[4]Σ(Xj(Yj^2))=Σ(YjMj)
a[1]Σ(XjYj)+a[2]Σ((Xj^2)Yj)+a[3]Σ(Xj(Yj^2))+a[4]Σ((Xj^2)(Yj^2))=Σ(XjYjMj)
という連立方程式になります。
これを行列で表せば、ANo.1のように逆行列を使ってa[1]～a[4]が計算できる訳です。

　最初から行列を使うと、もうちょっとスマートに表せます：
J行4列の行列P[j,n]　(j=1,2,…,J; n=1,2,3,4)を
P[j,1]=1
P[j,2]=Xj
P[j,3]=Yj
P[j,4]=XjYj
とし、縦ベクトルm　を
m[j] = Mj
とすると、上記の連立方程式は
(P' P)a = P'm
と書けます。（ここに ' は転置行列を表します。）なので、
a = (P' P)* (P'm)
です。（ただし　()* は逆行列を表します。）
Excelの関数で言えば、
(P' P)* P'm
は
=MMULT(MINVERSE(MMULT(TRANSPOSE("P"), "P")), MMULT(TRANSPOSE("P"),"m"))
ですね。

　このサイトで「最小二乗法 MINVERSE」というキーワードで検索すれば、他にも参考になる情報が見つかるでしょう。

- 5
- 件

通報する

No.5

回答者： stomachman
回答日時：2007/10/29 05:44

「近似」と仰る意味によって話が違います。

A: とびとびの(X,Y)における「計測データ」しかないので、計測しなかった(X,Y)における値をソレナリに推定したい。
B: 理論的にZ(X,Y)の関数の形が分かっているが、未知の係数を幾つか含んでいる。「計測データ」に含まれる誤差の影響を除いて、この関数の係数を決めたい。

どちらでしょうか。Aの場合は補間法(内挿法, interpolation)を使い、Bの場合はフィッティング(当てはめ, fitting)を行います。

いずれにしても、曲面がどんな種類の関数であるか（これをモデルと言う）を決めておいて、そのパラメータ（係数）を計算するというやり方をします。

A: 補間法で、（X,Y）が格子点になっている場合によく使われるのは、双n次多項式補間（ラグランジュ補間）です。たとえば双3次補間では、Xについて3次式、Yについても3次式であるような多項式
Z(X,Y) = (a[3,3] X^3 + a[2,3] X^2 + a[1,3] X + a[0,3])Y^3 + (a[3,2] X^3 + a[2,2] X^2 + a[1,2] X + a[0,2]) Y^2 + (a[3,1] X^3 + a[2,1] X^2 + a[1,1] X + a[0,1])Y +(a[3,0] X^3 + a[2,0] X^2 + a[1,0] X + a[0,0])
を使って、4個の格子点（計測点）を四隅とする矩形の内部での値を計算します。この式の16個の係数a[3,3], a[3,2], …は、前記の矩形を囲む4×4個の格子点での測定値を使って計算します。すなわち、曲面Zがこれら16個の格子点において測定値と一致するように係数を決めるのです。（上記の式は未知の係数a[0,0], a[0,1], …について１次式になっていますから、連立1次方程式の問題であり、簡単です。）格子点4個で囲まれる矩形ひとつずつについて、16個の係数を計算することになります。
　最も簡単なのは双1次補間で、格子点4個で囲まれる矩形の内部を、四隅の4個の格子点での計測値だけで決まる係数で表します。すなわち、
Z(X,Y) = (a[1,1] X + a[0,1])Y +(a[1,0] X + a[0,0])
です。（実用上、これで十分であるような場合も多々あります。）
　なお、測定点が格子になっていない場合には、スプライン補間が使えます。（幾つもバリエーションがありますが。）

B: フィッティングにおいては、全ての計測値を使ってひとつの式Z(X,Y)の係数を決めます。そうして決めた曲面Z(X,Y)は、各測定点(X,Y)に於ける測定値と必ずしも一致しません。このずれは、測定値の方に誤差がある、と考えるのです。
　フィッティングにはANo.1にある最小二乗法を使うのが普通です。特にモデルが多項式
Z(X,Y) = a[0,0]+a[1,0]X+a[0,1]Y+a[1,1]XY+a[2,0](X^2)+a[0,2](Y^2)+…
であるとき、未知の係数a[0,0], a[0,1], …を決める訳ですが、この式は（変数の数や次数が幾つだろうと関係なく）a[0,0], a[0,1], …について１次式になっていますから、線形最小二乗法で簡単に計算できます。
　なお、各測定点における測定誤差のばらつきが分かっている場合には、重み付き線形最小二乗法を使うのが適切です。

　という訳で、どういうツールを使うか以前に、何をやりたいかをはっきりさせる必要があります。