scipy.optimize.curve_fit 引发运行时错误

Question

这是我第一次使用curve_fit，我还没有找到与我的问题相匹配的示例。我的问题是，我是否正确使用了 curve_fit 数据格式？如果是，那么我的问题是数学问题，我必须纠正它（我还没有发现任何错误......）。 这是我的代码：

import numpy as np
import math as m
import scipy.optimize as sci

def f4(X,Tx,Ty,Tz,i,j,k):
    res=[]
#X looks like [x1,y1,z1,x2,y2,z2....]
    for n in range(np.shape(X)[0]/3):
        xr=X[n*3]+Tx+m.cos(j)*m.cos(k)*X[n*3]-m.cos(j)*m.sin(k)*X[n*3+1]+m.sin(j)*X[n*3+2]
        yr=X[n*3+1]+Ty+(m.sin(i)*m.sin(j)*m.cos(k)+m.cos(i)*m.sin(k))*X[n*3]+(-m.sin(i)*m.sin(j)*m.sin(k)+m.cos(i)*m.cos(k))*X[n*3+1]-m.sin(i)*m.cos(j)*X[n*3+2]
        zr=X[n*3+2]+Tz+(-m.cos(i)*m.sin(j)*m.cos(k)+m.sin(i)*m.sin(k))*X[n*3]+(m.cos(i)*m.sin(j)*m.sin(k)+m.sin(i)*m.cos(k))*X[n*3+1]+m.cos(i)*m.cos(j)*X[n*3+2]
        res.append(xr)
        res.append(yr)
        res.append(zr)
    return res

xdata2=[998.362,5000.052,99.4,997.862,5000.052,99.4,998.362,5000.052,98.9,997.862,5000.052,98.9]
ydata2=[999.46555,4999.801,99.4,999.112,5000.15455,99.4,999.46555,4999.801,98.9,999.112,5000.15455,98.9]
p0=[1,0,0,0,0,0.8]

popt,pcov=sci.curve_fit(f4,xdata2,ydata2,p0)

它引发：RuntimeError：未找到最佳参数：函数调用次数已达到 maxfev = 1400。

根据我的计算，p0 非常接近解。 我的代码还可以吗，还是应该继续寻找数学错误？

如果有人好奇，我试图找到 3 个平移和 3 个旋转，我必须将其应用于一组点 Xp 以获得一组点 Xr。

感谢任何建议

编辑 1/04：

我尝试了 unutbu 的方式来改变我的功能：

def f2(X,Tx,Ty,Tz,i,j,k):
    res=[]
    for n in range(np.shape(X)[0]):
        xr=X[n][0]+Tx+m.cos(j)*m.cos(k)*X[n][0]-m.cos(j)*m.sin(k)*X[n][1]+m.sin(j)*X[n][2]
    yr=X[n][1]+Ty+(m.sin(i)*m.sin(j)*m.cos(k)+m.cos(i)*m.sin(k))*X[n][0]+(-m.sin(i)*m.sin(j)*m.sin(k)+m.cos(i)*m.cos(k))*X[n][1]-m.sin(i)*m.cos(j)*X[n][2]
    zr=X[n][2]+Tz+(-m.cos(i)*m.sin(j)*m.cos(k)+m.sin(i)*m.sin(k))*X[n][0]+(m.cos(i)*m.sin(j)*m.sin(k)+m.sin(i)*m.cos(k))*X[n][1]+m.cos(i)*m.cos(j)*X[n][2]
    aux=[xr,yr,zr]
    res.append(aux)
res=np.array(res)
return res

我在数组中添加了两个点：

xdata3=np.array([[998.362,5000.052,99.4],[997.862,5000.052,99.4],[998.362,5000.052,98.9],[997.862,5000.052,98.9],[999.112,4999.801,98.9],[999.112,4999.801,99.4]])
ydata3=np.array([[999.46555,4999.801,99.4],[999.112,5000.15455,99.4],[999.46555,4999.801,98.9],[999.112,5000.15455,98.9],[1000.112,4999.801,98.9],[1000.112,4999.801,99.4]])

我试用了这个功能：

test=f2(xdata3,1,0.00001,0.00001,0.00001,0.00001,0.8)

In [17]:test
Out[17]:
array([[-1891.88925911,  9199.80185261,   198.87092002],
   [-1892.73761247,  9199.44317456,   198.87091992],
   [-1891.88926411,  9199.80185761,   197.87092002],
   [-1892.73761747,  9199.44317956,   197.87091992],
   [-1890.4366777 ,  9199.91400129,   197.87091663],
   [-1890.4366727 ,  9199.91399629,   198.87091663]])

也改变了 p0：

p0=[1,0.00001,0.00001,0.00001,0.00001,0.8]

然后我尝试了curve_fit：

test=f2(xdata3,1,0.00001,0.00001,0.00001,0.00001,0.8)

我有一个不同的错误：

error: Result from function call is not a proper array of floats. 

不确定发生了什么，因为我的测试变量看起来不错。

Answer 1

我认为问题仅在于您选择的模型函数。 数据根本不适合该模型。

如果您想使用平移和旋转来拟合数据，也许您需要 正在寻找的更像是f4（下），其中

• 找到对象的质心Xbar。 X-Xbar 将对象移回更靠近原点的位置。
• 围绕原点旋转X-Xbar（使用np.dot(R, X-Xbar)），然后
• 由T + Xbar 翻译。

（请注意，与原始 f4，X 也不翻译。）

---

import numpy as np
import scipy.optimize as optimize
from math import cos, sin

def f4(X,Tx,Ty,Tz,i,j,k):
    T = np.array([Tx, Ty, Tz])[:, np.newaxis]
    X = X.reshape(-1,3).T
    Xbar = X.mean(axis=1)[:, np.newaxis]
    X = X - Xbar
    R = np.array([
        (cos(j)*cos(k), -cos(j)*sin(k), sin(j)),

        (sin(i)*sin(j)*cos(k)+cos(i)*sin(k), 
         -sin(i)*sin(j)*sin(k)+cos(i)*cos(k), 
         -sin(i)*cos(j)),

        (-cos(i)*sin(j)*cos(k)+sin(i)*sin(k), 
         cos(i)*sin(j)*sin(k)+sin(i)*cos(k), 
         cos(i)*cos(j))])
    result = np.dot(R, X) + T + Xbar
    return result.T.ravel()

p0=[1,0,0,0,0,0.8]

xdata2 = np.array([998.362,5000.052,99.4,997.862,5000.052,99.4,998.362,5000.052,98.9,997.862,5000.052,98.9])
ydata2 = [999.46555,4999.801,99.4,999.112,5000.15455,99.4,999.46555,4999.801,98.9,999.112,5000.15455,98.9]

popt, pcov = optimize.curve_fit(f4, xdata2, ydata2, p0)
print(popt)

产量

[  1.17677500e+00  -7.42250000e-02  -4.02305065e-37  -1.28557241e-20
  -1.61334110e-20  -7.85398164e-01]

Answer 2

对于其中一些拟合例程，您应该不从 0.0 开始初始参数。可能想尝试将它们设置为小但非零（1e-6）。