拟合参数间强相关性是否导致圆形Moffat函数拟合失效？

问题描述 #

尝试用lmfit将圆形Moffat函数拟合至通量剖面，拟合效果极差；拟合报告显示参数beta与alpha的相关系数高达0.999，且二者误差极大，疑问是否为参数相关性导致拟合失败。

代码实现 #

from lmfit import Model
import matplotlib.pyplot as plt

def PSF_fit(r,beta,alpha):
    R=((beta-1)/(3.14*alpha*alpha))*((1+ (r / alpha)**2)**(-beta))
    return R

moffat = Model(PSF_fit)

print(f'parameter names: {moffat.param_names}')
print(f'independent variables: {moffat.independent_vars}')

params = moffat.make_params(beta=2.85,alpha=1.5)

result = moffat.fit(Intensity, params,r=r)

print(result.fit_report())

plt.rcParams.update({'font.size': 14})
plt.figure(figsize = (10,7))
plt.plot(r,Intensity, 'bo')

plt.plot(r, result.init_fit, '--', label='initial fit')
plt.plot(r, result.best_fit, '-', label='best fit')
plt.xlabel('Radius(pixels)')
plt.ylabel('Saturation corrected profile')
plt.legend()
plt.title('Circular Moffat Fit')
plt.show()

拟合报告 #

parameter names: ['beta', 'alpha']
independent variables: ['r']
[[Model]]
Model(PSF_fit)
[[Fit Statistics]]
# fitting method   = leastsq
# function evals   = 248
# data points      = 28
# variables        = 2
chi-square         = 127948.959
reduced chi-square = 4921.11382
Akaike info crit   = 239.961102
Bayesian info crit = 242.625511
[[Variables]]
beta:   18.2260915 +/- 1845366.32 (10124860.41%) (init = 2.85)
alpha:  9.61025099 +/- 553166.734 (5756007.15%) (init = 1.429321)
[[Correlations]] (unreported correlations are < 0.100)
C(beta, alpha) = 0.999

拟合结果 #

拟合曲线与实测数据偏差极大，最佳拟合完全无法匹配通量剖面的分布特征。

问题根源 #

参数beta与alpha的高相关性（0.999）正是拟合失败的核心原因：

1. 这两个参数在当前模型形式下存在强共线性——调整其中一个参数时，按比例调整另一个参数几乎不会改变拟合误差，导致拟合算法无法收敛到唯一的最优解，最终参数估计值偏离合理范围，误差被无限放大。
2. 现有Moffat函数的参数化方式放大了这种相关性：当beta较大时，(r/alpha)^2项的影响会被指数放大，此时alpha的缩放和beta的增减可以产生几乎一致的曲线形状，尤其当数据仅覆盖通量剖面的中心区域（缺乏尾部数据）时，beta无法被有效约束。

解决方案 #

1. 给参数添加合理约束 #

通过设置参数的上下限，强制拟合算法在物理合理的范围内搜索最优值：

params = moffat.make_params(beta=2.85, alpha=1.5)
# 根据天文PSF的典型值设置范围
params['beta'].set(min=2, max=20)
params['alpha'].set(min=0.5, max=5)

2. 重新参数化模型 #

改用**半高宽（FWHM）**代替alpha作为参数，因为FWHM的物理意义更明确，且与beta的相关性更低。FWHM与alpha的转换关系为：FWHM = 2*alpha*sqrt(2^(1/beta)-1)，重新定义模型：

import numpy as np
from lmfit import Model

def PSF_fit(r, beta, fwhm):
    alpha = fwhm / (2 * np.sqrt(2**(1/beta) - 1))
    norm = (beta - 1) / (np.pi * alpha**2)
    return norm * (1 + (r/alpha)**2)**(-beta)

moffat = Model(PSF_fit)
params = moffat.make_params(beta=2.85, fwhm=3.0) # 根据数据估算FWHM初始值

3. 优化数据覆盖范围 #

确保数据包含通量剖面的尾部区域（即较大半径r的部分），因为Moffat函数的尾部下降速率由beta主导，只有足够的尾部数据才能有效约束beta参数，降低与alpha的相关性。

4. 优化初始参数 #

根据数据先估算半高宽，再反推alpha的初始值；或者参考同类型PSF的典型参数值设置初始值，帮助拟合算法更快收敛到合理区域。

内容的提问来源于stack exchange，提问作者Arghya Chakraborty