python - numpy polyfit 没有根据卡方给出最佳线性拟合

Question

我必须找到适合某些数据的模型，然后计算其卡方值以查看它的拟合程度。所以我使用 numpy 的 polyfit 函数来找到拟合，然后计算卡方并绘制它。它看起来很合适，但是，我通过优化最小卡方来设法获得更好的拟合。我原以为 numpy 的 polyfit 会尝试做同样的事情，不是吗？

这是我的代码：

import numpy as np
from matplotlib import pyplot as plt
from scipy import stats

data = np.loadtxt('MW_Cepheids.dat', usecols=(1,2,3,4,5,6), dtype=float)

plx = data[:,0]                     # Parallax data
plx_err = data[:,1]                 # Error in plx.

P = data[:,2]                       # Period data
log_P = np.log10(P)                 # Logarithm of P.

m = data[:,3]                       # Apparent magnitude data

A = data[:,4]                       # Extinction data
A_err = data[:,5]                   # Error in ext.

dst = 1000/plx                      # Distance calculation
dst_err = 1000*(plx_err/plx**2)     # Error in dist. calc.


# Calculating Absolute Magnitude and its error
# This would correspond to y and the error in y
M = m - 5*np.log10(dst) + 5 - A
M_err = np.sqrt((5*dst_err/(np.log(10)*dst))**2 + A_err**2)


# Using numpy's polyfit function
a,b = np.polyfit(log_P, M, 1)

M_m = a*log_P + b

# Calculating chi squared from the polyfit result
chi2 = np.sum(((M-M_m)/(M_err))**2)
print(a, b, chi2)


def lin_fit (x, y, y_err):
    """ Cycles through probable values of gradient and intercept, calculating chi2 for each.
    Returns best values for a and b, ie when chi2 is smallest."""

    chi2_min = 999999
    a_range = np.arange(-2.5, -2.3, 0.001)
    b_range = np.arange(-1.7, -1.4, 0.001)

    for a in a_range:
        for b in b_range:
            y_m = a*x + b
            chi2 = np.sum(((y-y_m)/(y_err))**2)

            if chi2 < chi2_min:
                chi2_min = chi2
                a_fit = a
                b_fit = b
    return a_fit, b_fit, chi2_min


a, b, chi2 = lin_fit(log_P, M, M_err)
print(a, b, chi2)

数据在这个文件中：我们传输链接

我知道这里的简单选择是坚持使用另一种找到合适的方法，但这要麻烦得多，我想确保我正确计算 chi2 并理解为什么 polyfit 没有找到“最佳”拟合。

非常感谢！

标签： pythonnumpydata-analysis