如何在 3D 中用固定点进行多项式拟合

Question

我有一组 3D 中的 x、y、z 点 space 和另一个名为 charge 的变量，它表示沉积在特定 x、y、z 坐标中的电荷量。我想对该数据进行加权（根据检测器中沉积的电荷量加权，这恰好对应于更多电荷的更高权重），使其通过给定点，即顶点。

现在，当我为 2D 执行此操作时，我尝试了各种方法（将顶点带到原点并对所有其他点进行相同的转换并强制拟合通过原点，给顶点真的很重）但是 none 和 Jaime 在这里给出的答案一样好：How to do a polynomial fit with fixed points

它使用拉格朗日乘子的方法，我从本科高级多变量课程中隐约熟悉它，但其他的不多，而且看起来代码的转换并不像刚才那样简单添加 z 坐标。（请注意，即使代码没有考虑存入的电荷量，它仍然给了我最好的结果）。我想知道是否有相同算法的 3D 版本。我也在 Gmail 中联系了答案的作者，但没有收到他的回音。

这里有一些关于我的数据和我在 2D 中尝试做的事情的更多信息：

这是我执行此操作的代码，我强制顶点位于原点，然后适合数据设置 fit_intercept=False。我目前正在为 2D 数据采用这种方法，因为我不确定拉格朗日乘数是否有 3D 版本，但是有线性回归方法可以在 3D 中执行此操作，例如，此处：Fitting a line in 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 数组，vertexX 和 vertexY 分别是顶点的 x 和 y 坐标。这是我的数据：https://uploadfiles.io/bbhxo。我无法创建玩具数据，因为没有复制此数据的简单方法，它是由 Geant4 模拟中微子与氩核相互作用产生的。我不想摆脱数据的复杂性。而这个特定事件恰好是我的代码不起作用的事件，我不确定我是否可以专门生成一个数据，这样我的代码就不起作用了。

Answer 1

这更像是一个使用基本优化的手工制作解决方案。这是直截了当的。一个只是测量点到要拟合的线的距离，并使用基本 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()

哪个returns

>> [-0.68236386 -1.3057938 ]
>> [-0.70928735 -1.4617517 ]

固定点显示为黑色。蓝色的原始线。未加权和加权拟合分别为橙色和绿色。数据根据到线的距离着色。

如何在 3D 中用固定点进行多项式拟合

How to do a polynomial fit with fixed points in 3D

python

curve-fitting

data-fitting