0

我正在尝试学习如何使用 python // 优化适合更高维度的numpy数据scipy。成功地使用 scipy 来适应表格的一行后y = f(x)，我尝试扩展逻辑以适应表格的椭圆z = f(x, y)；两者都如下所示。我希望这种方法可以推广到适合更高维度（即球体）的形状。

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

## LINE EXAMPLE
npoints = 30
x = np.linspace(-5, 5, npoints)
m, b = 3, -4 # slope, y-intercept
initial_parameter_guess = (2.5, -1)
y = m * x + b # exact line
noise = np.random.uniform(-1, 1, size=x.size)
yn = y + noise # line with noise

def get_residuals(prms, x, y):
    """ """
    return (y - (prms[0] * x + prms[1]))**2

def f_error(prms, x, y):
    """ """
    resid = get_residuals(prms, x, y)
    return np.sum(resid)

result = minimize(f_error, x0=initial_parameter_guess, args=(x, yn))
# print(result)
yf = result.x[0] * x + result.x[1]

fig, ax = plt.subplots()
ax.scatter(x, yn, color='b', marker='.')
ax.plot(x, yf, color='r', alpha=0.5)
ax.grid(color='k', alpha=0.3, linestyle=':')
plt.show()
plt.close(fig)

线拟合示例

将此逻辑应用于椭圆的情况，

## ELLIPSE EXAMPLE
npoints = 75
theta = np.random.uniform(0, 2*np.pi, size=npoints)
a, b = (3, 5)
initial_parameter_guess = (2.5, 6)
xnoise = np.random.uniform(-1, 1, size=theta.size)
ynoise = np.random.uniform(-1, 1, size=theta.size)
x = a**2 * np.cos(theta)
xn = x + xnoise
y = b**2 * np.sin(theta)
yn = y + ynoise

def get_residuals(prms, x, y):
    """ """
    return 1 - ((x/prms[0])**2 + (y/prms[1])**2)

def f_error(prms, x, y):
    """ """
    resid = get_residuals(prms, x, y)
    return np.sum(resid)

result = minimize(f_error, x0=initial_parameter_guess, args=(xn, yn))
# print(result)

minimize算法viascipy没有找到最优参数；以下输出显示为print(result)：

      fun: -4243.611573066522
 hess_inv: array([[41.44400248, 39.68101343],
       [39.68101343, 37.99343048]])
      jac: array([-1496.81719971,  2166.68896484])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 174
      nit: 1
     njev: 42
   status: 2
  success: False
        x: array([-1.51640333,  2.93879038])

我已经看到了这个问题的另一种解决方案，它使用最小二乘的矩阵公式，但是这个例子对我来说有点难以理解。我见过的几乎所有方法都基于发布链接中的方法。我更喜欢使用minimize，除非由于我目前不知道的原因，线性代数方法更好。

无论如何，我上面的方法可以以一种可行的方式进行调整/调整，并且可以推广到更高的维度吗？

4

1 回答 1

3

代码有两个问题。不是最小化残差之和，而是最小化残差平方和。其次，椭圆方程应定义为x=a*np.cos(theta)和y=b*np.sin(theta) 

npoints = 75
theta = np.random.uniform(0, 2*np.pi, size=npoints)
a, b = (3, 5)
initial_parameter_guess = (2.5, 6)
xnoise = np.random.uniform(-1, 1, size=theta.size)
ynoise = np.random.uniform(-1, 1, size=theta.size)
x = a * np.cos(theta)
xn = x + xnoise
y = b * np.sin(theta)
yn = y + ynoise
def get_residuals(prms, x, y):
    """ """
    return 1 - ((x/prms[0])**2 + (y/prms[1])**2)

def f_error(prms, x, y):
    """ """
    resid = get_residuals(prms, x, y)
    return np.sum(np.square(resid))
result = minimize(f_error, x0=initial_parameter_guess,args=(xn, yn))

result
  fun: 5.85099318913613
 hess_inv: array([[ 0.07025572, -0.02902779],
       [-0.02902779,  0.12040811]])
      jac: array([-5.96046448e-08,  1.19209290e-07])
  message: 'Optimization terminated successfully.'
     nfev: 48
      nit: 10
     njev: 12
   status: 0
  success: True
        x: array([3.35248219, 5.13728323])
于 2019-12-14T16:32:06.467 回答