python - curve_fit 包含 numpy 数组的 2D 函数 --> 形状 (3,3,9) 和 (3,1) 未对齐 : 9 (dim 2) ! = 3 (dim 0)

Question

这是我的第一个问题。我希望我能为你们提供足够的信息。

我在将 curve_fit 与使用 numpy 数组的函数一起使用时遇到了问题。使用具有固定参数的函数效果很好。

任务:

我想找到两个向量。我所知道的是向量如何围绕彼此和坐标系旋转，以及一个向量与坐标系 y 轴之间的最终角度。

问题:

my_func 是我想使用的函数，但只有在将其放入 curve_fit 时才会抛出错误。

my_func2 应该是有效问题的分割。据我所知，主要的功能是相同的，除了 my_func2 的公式越来越不复杂。

这是我写的代码:

import numpy as np
from scipy.optimize import curve_fit

'''The function I want to optimize. 2D-xdata with 6 free parameters.'''
def my_func(X, hx, hy, hz, lx, ly, lz):

    # These are the two independently measured angles (2D-xdata).
    phi, alpha = X

    # y-axis of the coordinate system.
    yaxis = np.array([[0],
                      [1],
                      [0]])

    # First wanted vector h (first three parameters).
    h = np.array([[hx],
                  [hy],
                  [hz]])

    # Second wanted vector l (last three parameters).
    l = np.array([[lx],
                  [ly],
                  [lz]])

    # Projection matrix.
    Pxy = np.array([[1, 0, 0],
                    [0, 1, 0],
                    [0, 0, 0]])

    # Needed to generate the rotation matrix around the unknown vector h.
    h_norm = h / np.linalg.norm(h)
    n1, n2, n3 = h_norm[0][0], h_norm[1][0], h_norm[2][0]

    # Rotation matrix for rotation around the vector h by angle alpha.
    R_h = np.array([[n1 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n1 * n2 * (1 - np.cos(alpha)) - n3 * np.sin(alpha), n1 * n3 * (1 - np.cos(alpha)) + n2 * np.sin(alpha)],
                    [n1 * n2 * (1 - np.cos(alpha)) + n3 * np.sin(alpha), n2 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n2 * n3 * (1 - np.cos(alpha)) - n1 * np.sin(alpha)],
                    [n1 * n3 * (1 - np.cos(alpha)) - n2 * np.sin(alpha), n2 * n3 * (1 - np.cos(alpha)) + n1 * np.sin(alpha), n3 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha)]])

    # Rotate the vector l around the vector h by angle alpha.
    l_rot = np.dot(R_h, l)

    # Rotation matrix for rotation around x-axis by angle phi.
    R_x = np.array([[1, 0, 0],
                    [0, np.cos(phi), -np.sin(phi)],
                    [0, np.sin(phi), np.cos(phi)]])

    # Rotate the vector l_rot around the x-axis by angle phi.
    l_final = np.dot(R_x, l_rot)

    # Project the vector l_final into the xy-plane.
    l_final_xy = np.dot(Pxy, l_final)

    # Get the angle between the projected vector l_final_xy and the y-axis.
    angle = np.arccos(np.vdot(l_final_xy, yaxis) / (np.linalg.norm(l_final_xy)))

    # Return angle in degree.
    return angle * 180 / np.pi


'''A simplified version of the function above with less complex formulas.'''
def my_func2(X, a1, a2, a3, b1, b2, b3):

    # Represents phi and alpha of my_func.
    x1, x2 = X

    # Represents the first wanted vector of my_func.
    va = np.array([[a1],
                   [a2],
                   [a3]])

    # Represents the second wanted vector of my_func.
    vb = np.array([[b1],
                   [b2],
                   [b3]])

    # Represents the rotation matrix of my_func. It depends on the x-data and the parameters.
    M1 = np.array([[x1 * a1, x2 * b1, 0],
                   [0, x1 * a2, x2 * b2],
                   [x2 * b3, 0, x1 * a3]])

    # Some simplified math with the wanted vectors and the generated matrix.
    v_new = np.vdot(np.dot(M1, va), vb)

    return v_new

下面是一些测试。

# Some x-data: phi and alpha.
xdata = [[0, 0, 0, 30, 30, 30, 60, 60, 60],
         [0, 90, 180, 0, 90, 180, 0, 90, 180]]

# Some y-data.
ydata = [10, 11, 12, 13, 14, 15, 16, 17, 18]

# Test if my_func works as expected.
print(my_func([np.pi / 4, np.pi / 4], 1, 0, 0, 1, 1, 1))

此行打印出正确的 135.0。我还测试了其他值，结果看起来总是正确的。

print(curve_fit(my_func2, xdata, ydata)[0])

此行打印 [-0.88635298 2.75337506 0.66050304 0.13882423 0.01404608 0.02166652]。所以简化问题的拟合是可行的。

print(curve_fit(my_func, xdata, ydata)[0])

此行抛出以下错误:

l_rot = np.dot(R_h, l) ValueError:形状 (3,3,9) 和 (3,1) 未对齐:9 (dim 2) != 3 (dim 0)

所以最后一个问题是:为什么我只在使用 curve_fit 时遇到尺寸问题，我该如何规避？

Answer 1

感谢Mr.T我注意到我对曲线拟合的理解是错误的。我认为它会遍历传递的 x 和 y 数据并将每组值输入到函数中。

实际上曲线拟合进入孔 x 数据，它是一个列表/数组。该函数必须能够处理该列表/数组本身，如 Mr.T在他的第二条评论中解释道。

为了解决我的问题，我刚刚在我的函数中添加了一个 for 循环，该循环遍历 x 数据列表/数组。现在返回值列表而不是单个值。每个 x 数据集一个值。

我不确定这是否是最好的解决方案，但工作程序如下。

import numpy as np
from scipy.optimize import curve_fit

'''The function I want to optimize. 2D-xdata with 6 free parameters.'''
def my_func(X, hx, hy, hz, lx, ly, lz):

    angles = []

    for index in range(len(X)):

        # These are the two independently measured angles (2D-xdata).
        phi, alpha = X[index]

        # y-axis of the coordinate system.
        yaxis = np.array([[0],
                          [1],
                          [0]])

        # First wanted vector h (first three parameters).
        h = np.array([[hx],
                      [hy],
                      [hz]])

        # Second wanted vector l (last three parameters).
        l = np.array([[lx],
                      [ly],
                      [lz]])

        # Projection matrix.
        Pxy = np.array([[1, 0, 0],
                        [0, 1, 0],
                        [0, 0, 0]])

        # Needed to generate the rotation matrix around the unknown vector h.
        h_norm = h / np.linalg.norm(h)
        n1, n2, n3 = h_norm[0][0], h_norm[1][0], h_norm[2][0]

        # Rotation matrix for rotation around the vector h by angle alpha.
        R_h = np.array([[n1 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n1 * n2 * (1 - np.cos(alpha)) - n3 * np.sin(alpha), n1 * n3 * (1 - np.cos(alpha)) + n2 * np.sin(alpha)],
                        [n1 * n2 * (1 - np.cos(alpha)) + n3 * np.sin(alpha), n2 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n2 * n3 * (1 - np.cos(alpha)) - n1 * np.sin(alpha)],
                        [n1 * n3 * (1 - np.cos(alpha)) - n2 * np.sin(alpha), n2 * n3 * (1 - np.cos(alpha)) + n1 * np.sin(alpha), n3 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha)]])

        # Rotate the vector l around the vector h by angle alpha.
        l_rot = np.dot(R_h, l)

        # Rotation matrix for rotation around x-axis by angle phi.
        R_x = np.array([[1, 0, 0],
                        [0, np.cos(phi), -np.sin(phi)],
                        [0, np.sin(phi), np.cos(phi)]])

        # Rotate the vector l_rot around the x-axis by angle phi.
        l_final = np.dot(R_x, l_rot)

        # Project the vector l_final into the xy-plane.
        l_final_xy = np.dot(Pxy, l_final)

        # Get the angle between the projected vector l_final_xy and the y-axis.
        angle = np.arccos(np.vdot(l_final_xy, yaxis) / (np.linalg.norm(l_final_xy)))

        angles.append(angle * 180 / np.pi)

    # Return angle in degree.
    return angles

# Some x-data: phi and alpha.
xdata = [[0, 0],
         [0, 90],
         [0, 180],
         [30, 0],
         [30, 90],
         [30, 180],
         [60, 0],
         [60, 90],
         [60, 180]]

# Some y-data.
ydata = [10, 11, 12, 13, 14, 15, 16, 17, 18]

print(curve_fit(my_func, xdata, ydata)[0])