我想利用下面的代码来拟合一些2D数据到椭圆,这是我从this帖子中得到的。
import numpy as np
import matplotlib.pyplot as plt
alpha = 5
beta = 3
N = 500
DIM = 2
np.random.seed(2)
# Generate random points on the unit circle by sampling uniform angles
theta = np.random.uniform(0, 2*np.pi, (N,1))
eps_noise = 0.2 * np.random.normal(size=[N,1])
circle = np.hstack([np.cos(theta), np.sin(theta)])
# Stretch and rotate circle to an ellipse with random linear tranformation
B = np.random.randint(-3, 3, (DIM, DIM))
noisy_ellipse = circle.dot(B) + eps_noise
# Extract x coords and y coords of the ellipse as column vectors
X = noisy_ellipse[:,0:1]
Y = noisy_ellipse[:,1:]
# Formulate and solve the least squares problem ||Mx - b ||^2
M = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(M, b)[0].squeeze()
# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))
# Plot the noisy data
plt.scatter(X, Y, label='Data Points')
# Plot the original ellipse from which the data was generated
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse')
# Plot the least squares ellipse
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
但是,我需要使用完整的通用方程:Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 1,我从here得到,包括椭圆的倾斜和平移。
在上面的代码中,M矩阵:
M = np.hstack([X**2, X * Y, Y**2, X, Y])
假设我拟合的是Ax^2 + Bxy + Cy^2 + Dx + Ey = 1。
**问题:**我如何在F系数中“添加”,或者使用这种方法是不可能的?
1条答案
按热度按时间b4lqfgs41#
事实上,这是可能的,在优化的一点知识。我们知道F的平凡解肯定是
F=1
,其他参数等于0,但我们不希望这样,可能存在除平凡解之外的另一个解。我们来做个理论分析。在标准最小二乘法中,我们有
| 标准最小二乘法|
| - -----|
|
|
但是有几种方法可以扩展这种方法,以便它可以满足每个问题的需求。在你的例子中,我们知道的是,平凡的解决方案重置了参数
A, B, C, D and E
,也就是说,最简单的解决方案是认为这些对于米不能是Zero
因此,一个解决方案是添加一个正则化器,它添加到我们的目标函数中,如果我们的参数距离零,则添加一个越来越小的误差。其中一种方法是:以获得参数之和的倒数(参数之和越高,正则化器的惩罚就越小),以
Regularizer = alpha/sum(|x|)
的方式。sum(|x|)
是线性代数中的L1范数。增加alpha
以增加正则化函数的强度所以我们的LSTSQ是
| 我们的修正最小二乘|
| - -----|
|
|
但是要小心,因为
x
的范数不能为零,我们将在代码中处理这个问题。那么,让我们来编码
第一步
从scipy导入
minimize
第二步
创建自己的minimize函数
(the调整后的变量总是在其他参数之前,这就是为什么x在前面)
第三步
创建参数并应用于minimize函数
第四步
执行相同的plot代码,现在使用
x[5]
(F param)。结果
椭圆由-0.337x^2 +-0.143xy +-0.543y^2 +-0.0325x +-0.0331y +5.33 = 1给出
Scipy最小化
基本上,“最小化”函数通过迭代来工作,其中它计算成本函数具有最低值的位置,因此它在每次迭代时调整参数,直到它在给定容差和允许的最大迭代次数的情况下找到最佳参数。
查看更多文档:https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html