12

私はいくつかのパラメーター推定を行おうとしていますが、約30変数にわたって予測方程式の二乗誤差を最小化するパラメーター推定を選択したいと思います。方程式が線形の場合、30個の偏導関数を計算し、それらをすべてゼロに設定して、線形方程式ソルバーを使用します。しかし残念ながら、方程式は非線形であり、その導関数も非線形です。

方程式が単一の変数を超えている場合は、ニュートン法（ニュートンラプソンとも呼ばれます）を使用します。Webには、単一変数の関数に対してニュートン法を実装するための例とコードが豊富にあります。

約30個の変数があるとすると、ニュートン法を使用してこの問題の数値解をプログラムするにはどうすればよいですか？私は閉じた形の方程式を持っており、一次および二次導関数を計算できますが、そこからどのように進めるかはよくわかりません。私はウェブ上で多数の処理を見つけましたが、それらはすぐに重い行列表記になります。ウィキペディアで適度に役立つものを見つけましたが、それをコードに変換するのに問題があります。

分解が心配なのは、行列代数と行列反転です。一次方程式ソルバーを使用して行列を反転することはできますが、正しい行と列を取得したり、転置エラーを回避したりすることなどが心配です。

具体的には：

  • 変数をそれらの値にマッピングするテーブルを操作したいと思います。そのようなテーブルを引数として与えられた場合に二乗誤差を返すようなテーブルの関数を書くことができます。また、任意の変数に関する偏導関数を返す関数を作成することもできます。

  • 表の値について妥当な開始見積もりがあるので、収束について心配する必要はありません。

  • 推定値（各変数の値のテーブル）、関数、および偏微分関数のテーブルを使用して新しい推定値を生成するループを作成する方法がわかりません。

最後は私が助けて欲しいものです。直接の助けや良い情報源へのポインタは大歓迎です。

編集：私は閉じた形で一次および二次導関数を持っているので、それらを利用して、シンプレックス検索のようなよりゆっくりと収束する方法を避けたいと思います。

4

7 に答える 7

4

NumericalRecipesリンクが最も役に立ちました。誤差推定を象徴的に微分して30個の偏導関数を生成し、ニュートン法を使用してそれらをすべてゼロに設定しました。コードのハイライトは次のとおりです。

__doc.findzero = [[function(functions, partials, point, [epsilon, steps]) returns table, boolean
Where
  point       is a table mapping variable names to real numbers 
              (a point in N-dimensional space)
  functions   is a list of functions, each of which takes a table like
              point as an argument
  partials    is a list of tables; partials[i].x is the partial derivative
              of functions[i] with respect to 'x'
  epilson     is a number that says how close to zero we're trying to get
  steps       is max number of steps to take (defaults to infinity)
  result      is a table like 'point', boolean that says 'converged'
]]

-- See Numerical Recipes in C, Section 9.6 [http://www.nrbook.com/a/bookcpdf.php]




function findzero(functions, partials, point, epsilon, steps)
  epsilon = epsilon or 1.0e-6
  steps = steps or 1/0
  assert(#functions > 0)
  assert(table.numpairs(partials[1]) == #functions,
         'number of functions not equal to number of variables')
  local equations = { }
  repeat
    if Linf(functions, point) <= epsilon then
      return point, true
    end
    for i = 1, #functions do
      local F = functions[i](point)
      local zero = F
      for x, partial in pairs(partials[i]) do
        zero = zero + lineq.var(x) * partial(point)
      end
      equations[i] = lineq.eqn(zero, 0)
    end
    local delta = table.map(lineq.tonumber, lineq.solve(equations, {}).answers)
    point = table.map(function(v, x) return v + delta[x] end, point)
    steps = steps - 1
  until steps <= 0
  return point, false
end


function Linf(functions, point)
  -- distance using L-infinity norm
  assert(#functions > 0)
  local max = 0
  for i = 1, #functions do
    local z = functions[i](point)
    max = math.max(max, math.abs(z))
  end
  return max
end
于 2009-01-08T20:42:38.480 に答える
3

ニュートン法を使用する代わりに、Nelder-Mead のシンプレックス法が理想的にこの問題に適しており、C の Numerical Recpies で参照されています。

ロブ

于 2008-12-24T20:47:51.517 に答える
3

関数最小化アルゴリズムを求めています。ローカルとグローバルの 2 つの主要なクラスがあります。問題は最小二乗であるため、ローカル最小化アルゴリズムとグローバル最小化アルゴリズムの両方が同じ一意のソリューションに収束する必要があります。ローカル最小化はグローバルよりもはるかに効率的であるため、それを選択します。

多くの局所最小化アルゴリズムがありますが、最小二乗問題に特に適しているのはLevenberg-Marquardtです。そのようなソルバーを手元に持っていない場合 (たとえば、MINPACK から)、おそらくニュートンの方法で問題を解決できます。

x <- x - (hessian x)^-1 * grad x

ここで、線形ソルバーを使用してベクトルを掛けた逆行列を計算します。

于 2010-03-16T06:36:58.570 に答える
3

Numerical Recipes in C Web ページで必要なものを見つけることができるかもしれません。オンラインで利用できる無料版があり ます。 ここ(PDF) は、C で実装された Newton-Raphson 法を含む章です。また、 Netlib (LINPack など)で利用できるものを確認することもできます。

于 2008-12-24T20:34:44.117 に答える
1

十分な解決策があると思うかもしれませんが、私にとって、これについて考える最も簡単な方法は、最初に1変数の場合でそれを理解し、次にそれを行列の場合に拡張することです。

1変数の場合、一次導関数を二次導関数で割ると、次の試行ポイントまでの（負の）ステップサイズが得られます（例：-V / A）。

N変数の場合、1次導関数はベクトルで、2次導関数は行列（ヘッセ行列）です。導関数ベクトルに2次導関数の逆数を掛けると、結果は次の試行ポイントへの負のステップベクトルになります（例：-V *（1 / A））。

二階導ヘッセ行列を取得できると思います。それを反転するためのルーチンが必要になります。さまざまな線形代数パッケージにはこれらがたくさんあり、非常に高速です。

（このアイデアに精通していない読者のために、2つの変数がxとyであり、表面がv（x、y）であると仮定します。その場合、1次導関数はベクトルです。

 V = [ dv/dx, dv/dy ]

二階導関数は行列です：

 A = [dV/dx]
    [dV/dy]

また：

 A = [ d(dv/dx)/dx, d(dv/dy)/dx]
    [ d(dv/dx)/dy, d(dv/dy)/dy]

また：

 A = [d^2v/dx^2, d^2v/dydx]
    [d^2v/dxdy, d^2v/dy^2]

これは対称的です。）

表面が放物線状（一定の2階導関数）である場合、1ステップで答えが得られます。一方、2階導関数が非常に一定でない場合は、振動が発生する可能性があります。各ステップを半分（または一部）にカットすると、安定するはずです。

N == 1の場合、1変数の場合と同じことを行うことがわかります。

幸運を。

追加：コードが必要でした：

 double X[N];
// Set X to initial estimate
while(!done){
    double V[N];    // 1st derivative "velocity" vector
    double A[N*N];  // 2nd derivative "acceleration" matrix
    double A1[N*N]; // inverse of A
    double S[N];    // step vector
    CalculateFirstDerivative(V, X);
    CalculateSecondDerivative(A, X);
    // A1 = 1/A
    GetMatrixInverse(A, A1);
    // S = V*(1/A)
    VectorTimesMatrix(V, A1, S);
    // if S is small enough, stop
    // X -= S
    VectorMinusVector(X, S, X);
}
于 2009-01-14T15:46:19.410 に答える
1

すでに偏導関数があるので、一般的な勾配降下アプローチはどうでしょうか?

于 2009-01-08T15:55:38.530 に答える
0

私の意見では、Particle Swarm 法などの確率的オプティマイザを使用します。

于 2009-01-14T18:17:05.377 に答える