3d 表面的最小二乘拟合 python
Least squares fit in python for 3d surface
我想用我的曲面方程拟合一些数据。我已经尝试过 scipy.optimize.leastsq 但由于我无法指定界限，所以它给了我一个无法使用的结果。我也试过 scipy.optimize.least_squares 但它给了我一个错误：
ValueError: too many values to unpack
我的等式是：
f(x,y,z)=(x-A+y-B)/2+sqrt(((x-A-y+B)/2)^2+C*z^2)
应找到参数 A、B、C，以便当以下点用于 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]
]
边界将是 A > 0，B > 0，C > 1
我应该如何获得这样的身材？ python 中最好的工具是什么？我搜索了有关如何拟合 3d 表面的示例，但大多数涉及函数拟合的示例都是关于直线或平面拟合的。
我编辑了这个答案，以提供一个更通用的示例，说明如何使用 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])
====================
我想用我的曲面方程拟合一些数据。我已经尝试过 scipy.optimize.leastsq 但由于我无法指定界限，所以它给了我一个无法使用的结果。我也试过 scipy.optimize.least_squares 但它给了我一个错误：
ValueError: too many values to unpack
我的等式是：
f(x,y,z)=(x-A+y-B)/2+sqrt(((x-A-y+B)/2)^2+C*z^2)
应找到参数 A、B、C，以便当以下点用于 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]
]
边界将是 A > 0，B > 0，C > 1
我应该如何获得这样的身材？ python 中最好的工具是什么？我搜索了有关如何拟合 3d 表面的示例，但大多数涉及函数拟合的示例都是关于直线或平面拟合的。
我编辑了这个答案，以提供一个更通用的示例，说明如何使用 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])
====================