8

我正在使用 scipy.odr 使用权重拟合数据，但我不知道如何获得拟合优度或 R 平方的度量。有没有人建议如何使用函数存储的输出来获得这个度量？

4

3 回答 3

8

的res_var属性Output是拟合优度的所谓缩减卡方值，这是拟合优度统计量的流行选择。不过，非线性拟合有些问题。您可以直接查看残差（对于out.delta残差X和残差）。正如链接论文中所建议的那样，实施交叉验证或引导方法来确定拟合优度，留给读者作为练习。out.epsY

于 2014-01-28T13:02:39.757 回答
2

ODR的输出给出了估计的参数beta以及这些参数的标准偏差sd_beta。继第页。ODRPACK 文档的第 76章，您可以将这些值转换为 t 统计量(beta - beta_0) / sd_beta，其中beta_0是您要测试显着性的数字（通常为零）。从那里，您可以使用 t 分布来获取 p 值。

这是一个工作示例：

import numpy as np
from scipy import stats, odr


def linear_func(B, x):
    """
    From https://docs.scipy.org/doc/scipy/reference/odr.html
    Linear function y = m*x + b
    """
    # B is a vector of the parameters.
    # x is an array of the current x values.
    # x is in the same format as the x passed to Data or RealData.
    #
    # Return an array in the same format as y passed to Data or RealData.
    return B[0] * x + B[1]


np.random.seed(0)
sigma_x = .1
sigma_y = .15
N = 100
x_star = np.linspace(0, 10, N)
x = np.random.normal(x_star, sigma_x, N)
# the true underlying function is y = 2*x_star + 1
y = np.random.normal(2*x_star + 1, sigma_y, N)

linear = odr.Model(linear_func)
dat = odr.Data(x, y, wd=1./sigma_x**2, we=1./sigma_y**2)
this_odr = odr.ODR(dat, linear, beta0=[1., 0.])
odr_out = this_odr.run()
# degrees of freedom are n_samples - n_parameters
df = N - 2  # equivalently, df = odr_out.iwork[10]
beta_0 = 0  # test if slope is significantly different from zero
t_stat = (odr_out.beta[0] - beta_0) / odr_out.sd_beta[0]  # t statistic for the slope parameter
p_val = stats.t.sf(np.abs(t_stat), df) * 2
print('Recovered equation: y={:3.2f}x + {:3.2f}, t={:3.2f}, p={:.2e}'.format(odr_out.beta[0], odr_out.beta[1], t_stat, p_val))

Recovered equation: y=2.00x + 1.01, t=239.63, p=1.76e-137

在使用这种方法处理非线性问题时需要注意的一点，来自同一个 ODRPACK 文档：

[Boggs and Rogers, 1990b] 中介绍了蒙特卡洛实验的结果，该实验检查了四种不同测量误差模型的线性化置信区间的准确性。这些结果表明，Δ 的置信区域和区间不如 β 的准确。

尽管存在潜在的不准确性，但协方差矩阵经常用于为非线性普通最小二乘和测量误差模型构建置信区域和区间，因为由此产生的区域和区间计算成本低、通常足够且从业者熟悉。然而，在使用这些区域和区间时必须小心，因为近似值的有效性将取决于模型的非线性、误差的方差和分布以及数据本身。当需要更可靠的区间和区域时，应使用其他更准确的方法。（参见，例如，[Bates and Watts, 1988]、[Donaldson and Schnabel, 1987] 和 [Efron, 1985]。）

于 2020-05-20T05:52:25.013 回答
1

正如 R. Ken 所提到的，残差的卡方或方差是更常用的拟合优度检验之一。ODR 将残差平方和存储在其中out.sum_square，您可以验证自己out.res_var = out.sum_square/degrees_freedom对应于通常所说的缩减卡方：即卡方检验结果除以其预期值。

至于线性回归中另一个非常流行的拟合优度估计量 R 平方及其调整版本，我们可以定义函数

import numpy as np

def R_squared(observed, predicted, uncertainty=1):
    """ Returns R square measure of goodness of fit for predicted model. """
    weight = 1./uncertainty
    return 1. - (np.var((observed - predicted)*weight) / np.var(observed*weight))

def adjusted_R(x, y, model, popt, unc=1):
    """
    Returns adjusted R squared test for optimal parameters popt calculated
    according to W-MN formula, other forms have different coefficients:
    Wherry/McNemar : (n - 1)/(n - p - 1)
    Wherry : (n - 1)/(n - p)
    Lord : (n + p - 1)/(n - p - 1)
    Stein : (n - 1)/(n - p - 1) * (n - 2)/(n - p - 2) * (n + 1)/n

    """
    # Assuming you have a model with ODR argument order f(beta, x)
    # otherwise if model is of the form f(x, a, b, c..) you could use
    # R = R_squared(y, model(x, *popt), uncertainty=unc)
    R = R_squared(y, model(popt, x), uncertainty=unc)
    n, p = len(y), len(popt)
    coefficient = (n - 1)/(n - p - 1)
    adj = 1 - (1 - R) * coefficient
    return adj, R

从 ODR 运行的输出中，您可以找到模型参数的最佳值，out.beta此时我们拥有计算 R 平方所需的一切。

from scipy import odr

def lin_model(beta, x):
    """
    Linear function y = m*x + q
    slope m, constant term/y-intercept q
    """
    return beta[0] * x + beta[1]

linear = odr.Model(lin_model)
data = odr.RealData(x, y, sx=sigma_x, sy=sigma_y)
init = odr.ODR(data, linear, beta0=[1, 1])
out = init.run()

adjusted_Rsq, Rsq = adjusted_R(x, y, lin_model, popt=out.beta)
于 2021-03-26T22:51:49.510 回答