高斯拟合约化卡方值远小于1的拟合技术问题咨询

开展拟合优度评估过程中发现指标异常：已知*约化卡方（原文笔误写为reduced R²）*用于表征拟合优度，理想取值应在1附近，但本次拟合得到的约化卡方值仅为0.033，现有3个待解决的技术问题：

• 当前拟合模型中纳入实验数据不确定度的处理方式是否正确？
• 如何调整缩放规则使约化卡方值趋近于1？查阅资料得知需选用更合理的缩放方式，若同时处理多个拟合模型，应当对每个模型独立缩放，还是所有模型采用统一缩放标准？
• 除约化卡方、决定系数（Coefficient of determination）外，还有哪些统计参数可用于多个候选拟合模型间的最优模型判定？

import numpy as np
from lmfit import Minimizer, Parameters, fit_report
import scipy.constants as sc

x=np.array([-0.55140294, -0.53524022, -0.51157048, -0.46087502, -0.36867469,
       -0.36463724, -0.33394184, -0.31759441, -0.24301006, -0.19108247,
       -0.17992331, -0.15978266, -0.12565861, -0.06353306, -0.03416749,
        0.        ,  0.05381462,  0.05741924,  0.12256307,  0.13431629,
        0.15070565,  0.24975708,  0.30714032,  0.36017666,  0.36062918,
        0.44028426,  0.46225328,  0.49792282,  0.51374084,  0.54654304,
        0.5681194 ])
Data=np.array([3.10022439e-07, 5.28936553e-07, 1.08389150e-06, 4.51866520e-06,
       4.03761662e-05, 4.51393005e-05, 8.40351313e-05, 1.11528073e-04,
       3.84782855e-04, 7.40352481e-04, 8.34990761e-04, 1.01873368e-03,
       1.31826299e-03, 1.85815838e-03, 2.02222472e-03, 2.09575003e-03,
       1.98385694e-03, 1.91161400e-03, 1.40295968e-03, 1.28660846e-03,
       1.12519421e-03, 3.50500171e-04, 1.43286865e-04, 5.16806731e-05,
       4.98973395e-05, 8.13812022e-06, 4.47760445e-06, 1.66754498e-06,
       1.07843394e-06, 3.95451986e-07, 1.97595319e-07])
uncertainities=np.array([3.45420178e-06, 3.69177339e-06, 4.15460196e-06, 5.79355068e-06,
       1.25738539e-05, 1.30260341e-05, 1.69780223e-05, 1.94687535e-05,
       3.41877881e-05, 4.68776601e-05, 4.97061736e-05, 5.47758901e-05,
       6.29249004e-05, 7.43568853e-05, 7.73802945e-05, 7.85165928e-05,
       7.49519891e-05, 7.44931618e-05, 6.24757861e-05, 5.97435073e-05,
       5.57673389e-05, 3.15179837e-05, 2.03474708e-05, 1.29556740e-05,
       1.29044986e-05, 6.53164560e-06, 5.52338978e-06, 4.36260354e-06,
       3.99640819e-06, 3.45133957e-06, 3.21225716e-06])

def fGaussian(x,sig,A,xc,y0):
    y=A/((2*sc.pi)**0.5*sig)*np.exp(-(x-xc)**2/(2*sig**2))+y0
    return (y)

def fcn2min_fGauss(params, x, data, uncertainities):
    sig = params['sig']
    A= params['A']
    xc= params['xc']
    y0= params['y0']
    model = fGaussian(x,sig,A,xc,y0)
    return (model - data)/uncertainities

# 初始化高斯拟合参数
params_fGauss= Parameters()
params_fGauss.add('sig', value=10**-3, min=0)
params_fGauss.add('A', value=10**-4, min=0)
params_fGauss.add('xc', value=0, min=0)
params_fGauss.add('y0', value=0, min=0)
print ('******************** Gauss Parameters **********************')
params_fGauss.pretty_print()
print ('************************************************************')

print('fitters initiation')
fcn_args=(x, Data,uncertainities)
fitter_fGauss = Minimizer(fcn2min_fGauss, params_fGauss, fcn_args=fcn_args,nan_policy='propagate')

# 采用信赖域反射法执行最小二乘最小化
result_least_squares_Gauss = fitter_fGauss.minimize(method='least_squares',jac='3-point',loss='soft_l1')
print('fGauss fit report ---------------------------------')
print(fit_report(result_least_squares_Gauss))

[[Fit Statistics]]
    # fitting method   = least_squares
    # function evals   = 37
    # data points      = 31
    # variables        = 4
    chi-square         = 0.91407466
    reduced chi-square = 0.03385462
    Akaike info crit   = -101.238737
    Bayesian info crit = -95.5027884
[[Variables]]
    sig:  0.13202573 +/- 2.6268e-04 (0.20%) (init = 0.001)
    A:    6.9992e-04 +/- 1.6989e-06 (0.24%) (init = 0.0001)
    xc:   0.00130681 +/- 3.3337e-04 (25.51%) (init = 0)
    y0:   1.4545e-24 +/- 2.4519e-07 (16857606207555559424.00%) (init = 0)
[[Correlations]] (unreported correlations are < 0.100)
    C(sig, y0) = -0.314
    C(A, y0)   = -0.149
    C(sig, xc) =  0.124
    C(sig, A)  =  0.115

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