如何在三维空间中进行带固定点的多项式拟合

9
我有一组x、y、z坐标在三维空间中的点,还有另一个名为“charge”的变量,它表示在特定的x、y、z坐标上沉积的电荷量。我希望对这些数据进行加权(按探测器中沉积的电荷量进行加权,更多的电荷量对应更高的权重),使其经过给定点(即顶点)。
现在,当我尝试在2D中实现时,我尝试了各种方法(将顶点移到原点并对所有其他点进行相同的转换并强制拟合通过原点,使顶点具有非常高的权重),但是没有一个好像Jaime在这里提供的答案好:如何使用固定点拟合多项式 它使用拉格朗日乘数法,在本科高级多变量课程中我有一点模糊的了解,但除此之外并不太清楚,并且似乎这段代码的转换并不像仅添加z坐标那么容易。(请注意,即使代码没有考虑沉积的电荷量,它仍然给出了最佳结果)。我想知道是否有同样算法的3D版本。我也联系了答案的作者,但没有收到回复。
这里有关于我的数据和我在2D中试图做什么的更多信息:如何为拟合加权散点图上的点? 以下是我使用一种方式将顶点强制设为原点,然后设置fit_intercept = False对数据进行拟合的代码。我目前正在尝试这种方法来处理2D数据,因为我不确定Lagrange乘数法是否存在3D版本,但是在3D中有线性回归的方法可以实现,例如:在3D中拟合直线
import numpy as np
import sklearn.linear_model

def plot_best_fit(image_array, vertexX, vertexY):
    weights = np.array(image_array)
    x = np.where(weights>0)[1]
    y = np.where(weights>0)[0]
    size = len(image_array) * len(image_array[0])
    y = np.zeros((len(image_array), len(image_array[0])))
    for i in range(len(np.where(weights>0)[0])):
        y[np.where(weights>0)[0][i]][np.where(weights>0)[1][i]] = np.where(weights>0)[0][i]
    y = y.reshape(size)
    x = np.array(range(len(image_array)) * len(image_array[0]))
    weights = weights.reshape((size))
    for i in range(len(x)):
        x[i] -= vertexX
        y[i] -= vertexY
    model = sklearn.linear_model.LinearRegression(fit_intercept=False)
    model.fit(x.reshape((-1, 1)),y,sample_weight=weights)
    line_x = np.linspace(0, 512, 100).reshape((-1,1))
    pred = model.predict(line_x)
    m, b = np.polyfit(np.linspace(0, 512, 100), np.array(pred), 1)
    angle = math.atan(m) * 180/math.pi
    return line_x, pred, angle, b, m
image_array是一个numpy数组,vertexXvertexY分别表示顶点的x和y坐标。以下是我的数据:https://uploadfiles.io/bbhxo。我无法创建一个玩具数据,因为没有简单的方法复制这个数据,它是通过Geant4模拟中微子与氩核相互作用产生的。我不想消除数据的复杂性。而且,这个特定事件恰好是我的代码不能处理的事件,我不确定是否可以生成一份数据,使得我的代码不能处理它。

理想情况下两者都要(在我的情况下,只有一个点需要通过)。但正如我在问题中所说的,如果我能够得到一个没有权重的三维拉格朗日乘数(因为它在二维时效果最好),那就足够了。 - Always Learning Forever
不过,对于我的数据来说,那并不起作用。 - Always Learning Forever
是的,我的意思是从数据集中的每个数据点中减去约束点的坐标。为什么这对你的数据无效? - Andrew Guy
这是一个非常好的问题,我不确定。如果你想的话,你可以自己查看数据。顶点的坐标是:129 255。目前只是在二维平面内,一旦我找到一个合适的方法,我会尝试在三维空间中实现。 - Always Learning Forever
我可以在这里包含我的代码,但我不认为我能模拟我的数据。我的代码对大多数数据都有效,但对某些样本无效。这是其中一个无法工作的样本。 - Always Learning Forever
显示剩余10条评论
1个回答

6

这更像是一种使用基本优化的手工制作解决方案。它很直接。只需测量一个点到所拟合线的距离,并使用基本的optimize.leastsq最小化加权距离即可。下面是代码:

# -*- coding: utf-8 -*
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
from scipy import optimize
import numpy as np

def rnd( a ):
    return  a * ( 1 - 2 * np.random.random() ) 

def affine_line( s, theta, phi, x0, y0, z0 ):
    a = np.sin( theta) * np.cos( phi )
    b = np.sin( theta) * np.sin( phi )
    c = np.cos( theta )
    return np.array( [ s * a + x0, s * b + y0, s * c + z0 ] )

def point_to_line_distance( x , y, z , theta, phi, x0, y0, z0 ):
    xx = x - x0
    yy = y - y0
    zz = z - z0
    a = np.sin( theta) * np.cos( phi )
    b = np.sin( theta) * np.sin( phi )
    c = np.cos( theta )
    r = np.array( [ xx, yy, zz ] )
    t = np.array( [ a, b, c ] )
    return np.linalg.norm( r - np.dot( r, t) * t )

def residuals( parameters, fixpoint, data, weights=None ):
    theta, phi = parameters
    x0, y0, z0 = fixpoint
    if weights is None:
        w = np.ones( len( data ) )
    else:
        w = np.array( weights )
    diff = np.array( [ point_to_line_distance( x , y, z , theta, phi , *fixpoint ) for x, y, z in data ] )
    diff = diff * w
    return diff

### some test data
fixpoint = [ 1, 2 , -.3 ]
trueline = np.array( [ affine_line( s, .7, 1.7, *fixpoint ) for s in np.linspace( -1, 2, 50 ) ] )
rndData = np.array( [ np.array( [ a + rnd( .6), b + rnd( .35 ), c + rnd( .45 ) ] ) for a, b, c in trueline ] )
zData = [ 20 * point_to_line_distance( x , y, z , .7, 1.7, *fixpoint ) for x, y, z in rndData ]

### unweighted
bestFitValuesUW, ier= optimize.leastsq(residuals, [ 0, 0],args=( fixpoint, rndData ) )
print bestFitValuesUW
uwLine = np.array( [ affine_line( s, bestFitValuesUW[0], bestFitValuesUW[1], *fixpoint ) for s in np.linspace( -2, 2, 50 ) ] )

### weighted ( chose inverse distance as weight....would be charge in OP's case )
bestFitValuesW, ier= optimize.leastsq(residuals, [ 0, 0],args=( fixpoint, rndData, [ 1./s for s in zData ] ) )
print bestFitValuesW
wLine = np.array( [ affine_line( s, bestFitValuesW[0], bestFitValuesW[1], *fixpoint ) for s in np.linspace( -2, 2, 50 ) ] )

### plotting
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1, projection='3d' )
ax.plot( *np.transpose(trueline ) ) 
ax.scatter( *fixpoint, color='k' )
ax.scatter( rndData[::,0], rndData[::,1], rndData[::,2] , c=zData, cmap=cm.jet )

ax.plot( *np.transpose( uwLine ) ) 
ax.plot( *np.transpose( wLine ) ) 

ax.set_xlim( [ 0, 2.5 ] )
ax.set_ylim( [ 1, 3.5 ] )
ax.set_zlim( [ -1.25, 1.25 ] )

plt.show()

返回的是
>> [-0.68236386 -1.3057938 ]
>> [-0.70928735 -1.4617517 ]

黑色点表示修正点。蓝色表示原始线。橙色和绿色分别表示未加权和加权拟合。数据根据与线的距离着色。

结果


快速问题:返回的两个值是什么?[-0.68236386 -1.3057938],它们是theta和phi吗? - Always Learning Forever
还有,“original line”是什么意思? - Always Learning Forever
@AlwaysLearningForever 是的,返回的参数是 thetaphi 的最佳拟合值。由于这是通用数据,我通过提供给定行并添加随机数来生成它。这是“原始线”和拟合,因此应该接近这条线。 - mikuszefski
@AlwaysLearningForever 在这种情况下,我/我们有优势知道适合的样子。这可能有助于稍微调整一下。例如,可以看到,根据随机点的分布方式,加权会产生更多或更少的影响。 - mikuszefski

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接