scipy curve_fit拟合多变量数组非线性回归报错求解

问题背景 #

• 需求：使用scipy.optimize.curve_fit拟合自定义多变量函数的9个待求系数，所有输入自变量均为形状(1000,)的一维numpy数组，真实观测值ddphi11形状同为(1000,)
• 前置验证：手动逐点代入初始系数计算的预测值与真实值匹配度符合预期，函数逻辑、初始系数、原始数据无基础错误

问题复现 #

自定义目标函数 #

import numpy as np 
import matplotlib.pyplot as plt
import tensorflow as tf
import math
    
# 函数自变量：dphi1,dphi2,phi1,phi2,M；待拟合系数：c1~c9
def ddphi1(dphi1,dphi2,phi1,phi2,M,c1,c2,c3,c4,c5,c6,c7,c8,c9):
    return (-(c1*np.sin(phi1-phi2)*np.cos(phi1-phi2)*dphi1**2)-(c2*np.sin(phi1-phi2)*dphi2**2)+(c3*np.cos(phi1-phi2)*np.sin(phi2))+(c4*np.cos(phi1-phi2)*(dphi2-dphi1))-(c5*np.sin(phi1))+c6*M-c7*dphi1)/(c8-(c9*np.cos(phi1-phi2)*np.cos(phi1-phi2))) 

手动验证代码 #

初始系数取值：

p = [0.5625, 0.375, 27.590625000000003, 0.09375, 55.18125, 62.5, 0.425, 1, 0.5625]

逐点计算预测值并绘制对比曲线：

n = len(time1)
y = np.empty(n)

for i in range(n):
    y[i] = ddphi1(dphi11[i],dphi22[i],phi11[i],phi22[i],M[i],p[0],p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8])

plt.plot(time1, ddphi11)
plt.plot(time1, y, 'r')

手动计算结果可视化符合预期，证明基础逻辑无问题。

报错的自动拟合代码 #

from scipy.optimize import curve_fit

g = [0.56, 0.37, 27.63, 0.094, 55.18, 62.5, 0.625, 1, 0.56]
c,cov =curve_fit(ddphi1,(dphi1,dphi2,phi1,phi2,M),ddphi11,g)
print(c)

触发报错信息：

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (5,) + inhomogeneous part.

报错定位在curve_fit内部将xdata转换为浮点数组的步骤，即传入的自变量集合无法被解析为形状规则的numpy数组。

报错根因 #

curve_fit对输入的xdata格式有明确要求：

• 单变量拟合时，xdata为形状(n_samples,)的一维数组即可
• 多变量拟合时，xdata必须为形状(n_features, n_samples)的规则二维数组，直接将多个一维数组打包为元组传入时，numpy无法保证能将其解析为形状统一的二维数组，就会触发形状不匹配错误。
  同时原目标函数将5个自变量设为独立位置参数，和curve_fit传入二维数组时的参数传递逻辑不匹配，也会导致运行异常。

修复方案 #

• 使用np.vstack将5个长度一致的一维自变量按行堆叠，得到形状为(5, 1000)的规则二维数组作为xdata传入
• 调整目标函数的参数接收逻辑：第一个形参接收完整的二维自变量数组，在函数内部按行拆分出5个独立自变量后再参与计算
• （可选）通过bounds参数给待拟合系数设置合理取值范围，避免拟合过程中分母为0、参数溢出等计算异常

修复后可运行代码 #

from scipy.optimize import curve_fit
import numpy as np

# 调整目标函数参数结构
def ddphi1(x, c1,c2,c3,c4,c5,c6,c7,c8,c9):
    # 拆分二维输入数组为独立自变量
    dphi1, dphi2, phi1, phi2, M = x
    return (-(c1*np.sin(phi1-phi2)*np.cos(phi1-phi2)*dphi1**2)
            -(c2*np.sin(phi1-phi2)*dphi2**2)
            +(c3*np.cos(phi1-phi2)*np.sin(phi2))
            +(c4*np.cos(phi1-phi2)*(dphi2-dphi1))
            -(c5*np.sin(phi1))
            +c6*M - c7*dphi1) / (c8 - c9*np.cos(phi1-phi2)**2)

# 初始系数
g = [0.56, 0.37, 27.63, 0.094, 55.18, 62.5, 0.625, 1, 0.56]
# 堆叠自变量为规则二维数组
x_data = np.vstack([dphi1, dphi2, phi1, phi2, M])
# 执行拟合
c, cov = curve_fit(ddphi1, x_data, ddphi11, p0=g)
print("拟合得到的系数：", c)

内容的提问来源于stack exchange，提问作者Carlos Raymundo