2016-11-04 10 views
2

いくつかのデータに私の表面方程式を適合させたいと思います。私はすでにscipy.optimize.leastsqを試しましたが、境界を指定することができないため、使用できない結果が出ます。私はまた、scipy.optimize.least_squaresをしようとしたが、それは私にエラーを与える:3次元表面のためにPythonに収まる最小四角形

ValueError: too many values to unpack 

私の式は、上記の式のようになるようにA、B、Cが発見されなければならない

f(x,y,z)=(x-A+y-B)/2+sqrt(((x-A-y+B)/2)^2+C*z^2) 

パラメータ以下の点は、X、Y、Zのために使用される場合にゼロに近い、できるだけ:

[ 
    [-0.071, -0.85, 0.401], 
    [-0.138, -1.111, 0.494], 
    [-0.317, -0.317, -0.317], 
    [-0.351, -2.048, 0.848] 
    ] 

境界は> 0、B> 0、C> 1

であろう

どのように私はそのようなフィット感を得る必要がありますか?それを行うためのPythonの最も良いツールは何ですか?私は3Dサーフェスをどのようにフィットさせるかの例を検索しましたが、関数フィッティングを含むほとんどの例は、ラインまたはフラットサーフェスのフィッティングに関するものです。

答えて

5

私はこの問題を、scipyの一般的なoptimize.minimizeメソッドとscipyのoptimize.least_squaresメソッドでどのように解決できるかについて、より一般的な例を提供するためにこの答えを編集しました。


まず問題を設定できます:

import numpy as np 
import scipy.optimize 

# =============================================== 
# SETUP: define common compoments of the problem 


def our_function(coeff, data): 
    """ 
    The function we care to optimize. 

    Args: 
     coeff (np.ndarray): are the parameters that we care to optimize. 
     data (np.ndarray): the input data 
    """ 
    A, B, C = coeff 
    x, y, z = data.T 
    return (x - A + y - B)/2 + np.sqrt(((x - A - y + B)/2) ** 2 + C * z ** 2) 


# Define some training data 
data = np.array([ 
    [-0.071, -0.85, 0.401], 
    [-0.138, -1.111, 0.494], 
    [-0.317, -0.317, -0.317], 
    [-0.351, -2.048, 0.848] 
]) 
# Define training target 
# This is what we want the target function to be equal to 
target = 0 

# Make an initial guess as to the parameters 
# either a constant or random guess is typically fine 
num_coeff = 3 
coeff_0 = np.ones(num_coeff) 
# coeff_0 = np.random.rand(num_coeff) 

これは厳密には最小二乗法ではありませんが、どのようにこのようなものについては? この解決法は、問題にハンマーハンマーを投げるようなものです。おそらくSVDソルバーを使って最小二乗法を使ってより効率的に解を得る方法がありますが、答えを探しているのであれば、scipy.optimize.minimizeがあなたの答えを見つけるでしょう。

# =============================================== 
# FORMULATION #1: a general minimization problem 

# Here the bounds and error are all specified within the general objective function 
def general_objective(coeff, data, target): 
    """ 
    General function that simply returns a value to be minimized. 
    The coeff will be modified to minimize whatever the output of this function 
    may be. 
    """ 
    # Constraints to keep coeff above 0 
    if np.any(coeff < 0): 
     # If any constraint is violated return infinity 
     return np.inf 
    # The function we care about 
    prediction = our_function(coeff, data) 
    # (optional) L2 regularization to keep coeff small 
    # (optional) reg_amount = 0.0 
    # (optional) reg = reg_amount * np.sqrt((coeff ** 2).sum()) 
    losses = (prediction - target) ** 2 
    # (optional) losses += reg 
    # Return the average squared error 
    loss = losses.sum() 
    return loss 


general_result = scipy.optimize.minimize(general_objective, coeff_0, 
             method='Nelder-Mead', 
             args=(data, target)) 
# Test what the squared error of the returned result is 
coeff = general_result.x 
general_output = our_function(coeff, data) 
print('====================') 
print('general_result =\n%s' % (general_result,)) 
print('---------------------') 
print('general_output = %r' % (general_output,)) 
print('====================') 

出力は次のようになります。

==================== 
general_result = 
final_simplex: (array([[ 2.45700466e-01, 7.93719271e-09, 1.71257109e+00], 
     [ 2.45692680e-01, 3.31991619e-08, 1.71255150e+00], 
     [ 2.45726858e-01, 6.52636219e-08, 1.71263360e+00], 
     [ 2.45713989e-01, 8.06971686e-08, 1.71260234e+00]]), array([ 0.00012404, 0.00012404, 0.00012404, 0.00012404])) 
      fun: 0.00012404137498459109 
     message: 'Optimization terminated successfully.' 
      nfev: 431 
      nit: 240 
     status: 0 
     success: True 
      x: array([ 2.45700466e-01, 7.93719271e-09, 1.71257109e+00]) 
--------------------- 
general_output = array([ 0.00527974, -0.00561568, -0.00719941, 0.00357748]) 
==================== 

は、私はあなたがする必要があるすべては、実際の最小二乗にこれを適応する文書で発見された残差を計算する関数を指定します。

# =============================================== 
# FORMULATION #2: a special least squares problem 

# Here all that is needeed is a function that computes the vector of residuals 
# the optimization function takes care of the rest 
def least_squares_residuals(coeff, data, target): 
    """ 
    Function that returns the vector of residuals between the predicted values 
    and the target value. Here we want each predicted value to be close to zero 
    """ 
    A, B, C = coeff 
    x, y, z = data.T 
    prediction = our_function(coeff, data) 
    vector_of_residuals = (prediction - target) 
    return vector_of_residuals 


# Here the bounds are specified in the optimization call 
bound_gt = np.full(shape=num_coeff, fill_value=0, dtype=np.float) 
bound_lt = np.full(shape=num_coeff, fill_value=np.inf, dtype=np.float) 
bounds = (bound_gt, bound_lt) 

lst_sqrs_result = scipy.optimize.least_squares(least_squares_residuals, coeff_0, 
               args=(data, target), bounds=bounds) 
# Test what the squared error of the returned result is 
coeff = lst_sqrs_result.x 
lst_sqrs_output = our_function(coeff, data) 
print('====================') 
print('lst_sqrs_result =\n%s' % (lst_sqrs_result,)) 
print('---------------------') 
print('lst_sqrs_output = %r' % (lst_sqrs_output,)) 
print('====================') 

出力はここにある:あなたが書いた

==================== 
lst_sqrs_result = 
active_mask: array([ 0, -1, 0]) 
     cost: 6.197329866927735e-05 
     fun: array([ 0.00518416, -0.00564099, -0.00710112, 0.00385024]) 
     grad: array([ -4.61826888e-09, 3.70771396e-03, 1.26659198e-09]) 
     jac: array([[-0.72611025, -0.27388975, 0.13653112], 
     [-0.74479565, -0.25520435, 0.1644325 ], 
     [-0.35777232, -0.64222767, 0.11601263], 
     [-0.77338046, -0.22661953, 0.27104366]]) 
    message: '`gtol` termination condition is satisfied.' 
     nfev: 13 
     njev: 13 
    optimality: 4.6182688779976278e-09 
     status: 1 
    success: True 
      x: array([ 2.46392438e-01, 5.39025298e-17, 1.71555150e+00]) 
--------------------- 
lst_sqrs_output = array([ 0.00518416, -0.00564099, -0.00710112, 0.00385024]) 
==================== 
+0

コードは私に正しい結果を与えます!しかし、私はあなたが境界条件をどのように課したのか本当に理解していないのですか? – UN4

+0

私はそれをするのを忘れた。 ¯\\ _(ツ)_ /¯。最初の推測では正の値を使用していたので、ちょうどうまくいったと思うし、L2の正則化は係数が小さくなるよう奨励したので、良い答えが得られたときに停止した。メソッドは非常に一般的ですが、if np.any(coeff <0):return np.infを含めることができ、制約を課すことになります。 – Erotemic

+0

ああ、そうだよ。ありがとう! – UN4

関連する問題