20

X、Y、およびZを解きたい3 つの未知数を持つ 4 つの非線形方程式があります。方程式は次の形式です。

F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ

...ここでa、bとは、4 つの方程式の のc各値に依存する定数です。F

これを解決する最善の方法は何ですか？

4

2 に答える 2

52

これには 2 つの方法があります。

  1. 非線形ソルバーを使用する
  2. 問題を線形化し、最小二乗法で解く

設定

したがって、あなたの質問を理解しているように、4 つの異なるポイントで F、a、b、および c を知っており、モデル パラメーター X、Y、および Z を反転したいと考えています。3 つの未知数と 4 つの観測データ ポイントがあるので、問題は過剰決定です。したがって、最小二乗法で解いていきます。

この場合、逆の用語を使用する方が一般的であるため、式をひっくり返しましょう。それ以外の：

F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ

かきましょう：

F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)

F、X、Y、およびZ4 つの異なる点 (例: )を知っている場所F_0, F_1, ... F_i。

方程式自体ではなく、変数の名前を変更しているだけです。（これは何よりも私の考えやすさのためです。）

線形解

実際、この方程式を線形化することは可能です。a^2、b^2、a b cos(c)、については簡単に解くことができますa b sin(c)。これを少し簡単にするために、もう一度ラベルを付け直しましょう。

d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)

これで、方程式ははるかに簡単になりました: F_i = d + e X_i + f Y_i + g Z_i. d、e、f、の最小二乗線形反転を行うのは簡単gです。a次に、 、b、およびを取得できますc。

a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)

さて、これを行列形式で書きましょう。の 4 つの観察結果を翻訳します (これから書くコードは任意の数の観察結果を取得しますが、現時点では具体的なものにしておきましょう):

F_i = d + e X_i + f Y_i + g Z_i

の中へ：

|F_0|   |1, X_0, Y_0, Z_0|   |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2|   |1, X_2, Y_2, Z_2|   |f|
|F_3|   |1, X_3, Y_3, Z_3|   |g|

または: F = G * m(私は地球物理学者なのでG、"Green's Functions" とm"Model Parameters" に使用dします。通常、 の代わりに "data" にも使用しFます。)

Python では、これは次のように変換されます。

def invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

非線形解

scipy.optimize@Joeが提案したように、を使用してこれを解決することもできます。で最もアクセスしやすい関数はscipy.optimize、scipy.optimize.curve_fitデフォルトでレーベンバーグ・マルカート法を使用する関数です。

Levenberg-Marquardt は「山登り」アルゴリズムです (この場合は下り坂になりますが、この用語はとにかく使用されます)。ある意味では、モデル パラメーターの最初の推測 (デフォルトでは ではすべて 1 scipy.optimize) を行いobserved - predicted、パラメーター空間の の勾配を下に向かって下ります。

警告:適切な非線形反転法を選択し、最初の推測を行い、メソッドのパラメーターを調整することは、まさに「闇の芸術」です。やることでしか身につかないし、うまくいかないこともたくさんあります。レーベンバーグ・マルカート法は、パラメーター空間がかなり滑らかな場合に適した一般的な方法です (これは滑らかである必要があります)。他の状況で優れている他の多くの方法があります (シミュレーテッド アニーリングのようなより一般的な方法に加えて、遺伝的アルゴリズム、ニューラル ネットワークなどを含む)。ここではその部分について掘り下げるつもりはありません。

一部の最適化ツールキットが修正scipy.optimizeしようとして、処理しようとしない一般的な落とし穴が 1 つあります。モデル パラメータのマグニチュードが異なる場合 (例: a=1, b=1000, c=1e-8)、マグニチュードが類似するように再スケーリングする必要があります。そうscipy.optimizeしないと、 の「ヒル クライミング」アルゴリズム (LM など) は局所勾配の推定値を正確に計算できず、非常に不正確な結果が得られます。今のところa、 、b、および のc大きさが比較的似ていると仮定しています。また、基本的にすべての非線形メソッドでは、最初の推測を行う必要があり、その推測に敏感であることに注意してください。デフォルトはかなり正確p0な.curve_fita, b, c = 1, 1, 1a, b, c = 3, 2, 1

警告が邪魔にcurve_fitならないように、関数、観測が行われたポイントのセット (単一のndim x npoints配列として)、および観測された値が渡されることを期待します。

したがって、関数を次のように書くと:

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

に渡す前に、わずかに異なる引数を受け入れるようにラップする必要がありますcurve_fit。

手短に：

def nonlinear_invert(f, x, y, z):
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

スタンドアロン 2 つの方法の例:

完全な実装を提供するために、次の例を示します

  1. 関数を評価するためにランダムに分散された点を生成し、
  2. それらのポイントで関数を評価します (設定されたモデル パラメーターを使用)。
  3. 結果にノイズを追加し、
  4. 次に、上記の線形方法と非線形方法の両方を使用して、モデル パラメーターを反転します。
import numpy as np
import scipy.optimize as opt

def main():
    nobservations = 4
    a, b, c = 3.0, 2.0, 1.0
    f, x, y, z = generate_data(nobservations, a, b, c)

    print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
    print linear_invert(f, x, y, z)

    print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
    print nonlinear_invert(f, x, y, z)

def generate_data(nobservations, a, b, c, noise_level=0.01):
    x, y, z = np.random.random((3, nobservations))
    noise = noise_level * np.random.normal(0, noise_level, nobservations)
    f = func(x, y, z, a, b, c) + noise
    return f, x, y, z

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

def linear_invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

def nonlinear_invert(f, x, y, z):
    # "curve_fit" expects the function to take a slightly different form...
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

main()
于 2013-10-23T14:35:21.663 に答える