Python Grid Finder三参数寻优结果异常，寻求问题排查帮助

我尝试编写一个Grid Finder，基于x3data和y3data两个列表寻找a、b、c三个参数的最优值，函数定义在代码顶部：

delx=np.zeros(1)
dely=np.zeros(1)  
    
def func(x, a, b, c):
   #y = 100*a*x + 100*b
   
   #y = y - 360.*np.floor(y/360.)    # subtracting 360 times the nearest integer of y/360 from y
   
   delx = np.cos(a*x + b) - c 
   dely = np.sin(a*x + b) 
   rad = np.sqrt(delx*delx+dely*dely)
   y=np.arccos(delx/rad)
   ysign1 = np.sign(dely)
   ysign2 = (1. - ysign1)/2.
   y = ysign2 * 2. * np.pi + (ysign1 *y)
   y = y*180./np.pi      # in degrees 
   return y


ndat=len(x3data)
ngrid=50

chimin=1.e10
    
apar_best = 0.0
bpar_best = 0.0
cpar_best = 0.0

for i in range(ngrid):
    for j in range(ngrid):
        for k in range(ngrid):
            apar=10.*i/float(ngrid-1)
            bpar=1.*j/float(ngrid-1)
            cpar=1.*k/float(ngrid-1)
            chi=0.
            for l in range(ndat):
                diff=(y3data[l]-func(x3data[l],apar,bpar, cpar))*np.pi/180.
                chi=chi+(np.sin(diff))**2
                if(chi<chimin):
                    print(i, j, k, apar, bpar, cpar, np.sqrt(chi/float(ndat)))
                    chimin=chi
                    apar_best=apar
                    bpar_best=bpar
                    cpar_best = cpar
    
print("best fit results:")                                          

print(apar_best, bpar_best, cpar_best, np.sqrt(chimin/float(ndat)))           
print("a =" , apar_best)
print("b =" , bpar_best)
print("c =" , cpar_best)

运行得到的最优拟合结果如下：
0.0 0.8367346938775511 0.673469387755102 6.264895969973369e-05
a = 0.0
b = 0.8367346938775511
c = 0.673469387755102

我始终得到相同的a、b、c值，且结果不正确，不确定问题出在哪里，希望得到帮助。

Grid Finder逐行说明 #

• 1. 将x3data数组的长度声明为变量ndat，网格大小声明为ngrid。
• 2. 将chimin（最小卡方值）设置为一个极大的数值。
• 3. 声明三个变量apar_best、bpar_best、cpar_best并初始化为0.0。
• 4. 开始外层for循环，i从0到ngrid-1遍历。
• 5. 第一个for循环内嵌套第二个for循环，j从0到ngrid-1遍历。
• 6. 第二个for循环内嵌套第三个for循环，k从0到ngrid-1遍历。
• 7. 针对每组i、j、k值，计算apar=10*i/(ngrid-1)、bpar=j/(ngrid-1)、cpar=k/(ngrid-1)。
• 8. 将chi值初始化为0。
• 9. 开始内层for循环，l从0到ndat-1遍历。
• 10. 在此循环内，计算y3data与func输出值的差值并存储到diff变量中。
• 11. 通过计算diff正弦值的平方累加得到chi值。
• 12. 如果当前chi值小于chimin，则打印相关参数及拟合指标，并更新chimin和最优参数值。

问题排查与修复建议 #

1. 卡方更新时机错误
   当前代码在内层遍历单个数据点时就判断chi < chimin并更新最优参数，此时chi仅累加了部分数据点的贡献，并非完整的总卡方值。这会导致第一个数据点的chi就可能替换掉初始的极大chimin，后续参数即使更优也无法被正确识别。
   修复：将if(chi<chimin)的判断和更新代码移到内层循环（遍历所有数据点）之后：

   chi=0.
   for l in range(ndat):
       diff=(y3data[l]-func(x3data[l],apar,bpar, cpar))*np.pi/180.
       chi=chi+(np.sin(diff))**2
   # 所有数据点计算完成后再判断更新
   if(chi<chimin):
       print(i, j, k, apar, bpar, cpar, np.sqrt(chi/float(ndat)))
       chimin=chi
       apar_best=apar
       bpar_best=bpar
       cpar_best = cpar
   

2. 参数网格范围可能不合理

   • apar的范围是0_{10，bpar和cpar的范围是0}1，如果真实参数不在这些区间内，网格寻优必然找不到正确结果。需要根据业务场景或数据特性调整参数范围，比如扩大/缩小apar的范围，或者调整bpar、cpar的取值区间。
   • 网格密度ngrid=50，若参数范围过大，步长会太粗，可能错过最优解；若范围合适，步长可以保证精度，可根据实际情况调整。

3. 全局变量与局部变量混淆
   代码开头定义了全局的delx=np.zeros(1)和dely=np.zeros(1)，但函数func内部又重新定义了同名的局部变量。虽然此处未引发错误，但容易导致变量作用域混淆，建议删除开头的全局变量定义，避免意外问题。

4. 卡方统计量合理性验证
   当前用(np.sin(diff))**2作为卡方项，diff是角度差转弧度。需要确认这个统计量是否符合你的数据分布和拟合目标，常规卡方统计量通常是(观测值-预测值)^2/方差，如果你的数据是角度类周期性数据，当前定义可能有合理性，但需结合业务场景验证。

内容的提问来源于stack exchange，提问作者Rory Fitzpatrick