我有一些要拟合的数据,然后执行卡方检验以获得拟合优度。显然,我要应用的拟合度不能很好地拟合数据(这本身不是问题,我不一定期望如此),但是返回的值scipy.stats.chisquare会建议几乎完美的配合,这显然是错误的。

到目前为止,我所做的是定义一个描述我要应用的拟合的函数(正弦拟合),然后使用scipy.optimize.curve_fit从popt获取拟合参数,然后在先前定义的函数以生成拟合。

然后,我将测量的数据和拟合的数据放入scipy.stats.chisquare中,以尝试拟合,但返回的p值为1.0,这是不正确的。我的假设是,在scipy.stats.chisquare中使用scipy.optimize.curve_fit生成的值存在一些问题,但是如果是这种情况,我不明白为什么这是一个问题或如何解决。

我将测量数据存储在两个列表中,以下分别称为“时间”和“费率”

import numpy as np
import math
%matplotlib inline
import matplotlib.pyplot as plt
from statistics import stdev
import scipy


time =[309.6666666666667, 326.3333333333333, 334.6666666666667, 399.9166666666667, 416.5833333333333, 433.25, 449.91666666666663, 466.58333333333337, 483.25, 499.91666666666663,]

rate = [0.298168, 0.29317, 0.306496, 0.249861, 0.241532, 0.241532, 0.206552, 0.249861, 0.253193, 0.239867]

def oscillation(t,A,C):
    return(A*np.cos((2*np.pi*(t-x0))/(t0))+C)
t0 = 365.25
A = 0.35/2
x0 = 152.5
C = 0.475

popt, pcov = curve_fit(oscillation, time, rate, p0=[A,C])


rate_fit = []

for t in time:
    r = oscillation(t, popt[0],popt[1])
    rate_fit.append(r)

print(scipy.stats.chisquare(rate, f_exp=rate_fit))

plt.plot(time,rate, '.')
plt.plot(time,rate_fit,'--')


上面的输出是一个拟合,在绘制时看起来确实最适合数据,但显然不是完美拟合,因此另一个输出的p值为0.99999999999458533,这显然是错误的

最佳答案

您仅适合两个参数AC,因此强制相位和周期。
如果您还适合阶段和期间,那么您会更适合:

python - scipy.stats.chisquare没有提供输入数据所期望的结果-LMLPHP

同样在这种情况下,我的p值为1.0。

固定x0t0时p值为1.0的原因是,您的结果是可以使用x0t0的那些值进行的最佳拟合。强迫这些值,很可能会导致整体更差。为了比较,我免费使用x0t0

A = -3.45840427e-02
C = 2.65142203e-01
x0 = 1.88838771e+02
t0 = 2.61112538e+02


将其与t0 = 365.25x0 = 152.5进行比较。

当然,您有(物理)原因要修复,例如t0到一年,但是在这种情况下,您不必担心情节看起来很糟;您的p值仍会考虑到这一点。

但是,更可能的原因是您还忘记了ddof中的scipy.stats.chisquare参数。默认值为ddof=0,这不是您所拥有的:在您的情况下为len(rate) - 2,在我的上述情况下为len(rate) - 4
为了适合您(固定了t0x0),结果为p = 0.902。在所有参数均未设置的情况下,结果为0.999887(即再次为1)。



奖励:将周期t0固定为365.25时输出:

A = -4.05218922e-02
C = 2.74772524e-01
x0 = 8.69008279e+01

p = 0.997


和绘制的拟合:

python - scipy.stats.chisquare没有提供输入数据所期望的结果-LMLPHP

关于python - scipy.stats.chisquare没有提供输入数据所期望的结果,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/55618534/

10-16 07:31