如何在y误差不对称时计算线性拟合斜率的误差

Question

我有一组数据，可以把它们画成x轴和y轴的图。y轴上的数据有不对称的误差，也就是说，

我想用一个线性函数来拟合这些数据。在Python中，我可以用很多方法来进行这个拟合，但它们都有一个共同的问题，就是如何得到拟合参数的误差。

这和这个问题很相关，但因为我的数据在y轴上有不对称的误差，所以使用scipy.optimize.curve_fit来计算斜率的误差就变得不可能（或者说不太简单）。

那么，当y轴的误差条是不对称的时候，我该如何计算线性拟合的斜率误差呢？

有没有什么Python函数可以做到这一点？

Answer 1

有一篇由Barlow等人于2004年发表的论文，讨论了如何处理带有不对称误差条的变量的平均值计算。你可以考虑使用这些方法。

我通常采用的简单方法是，从一个叫做“分裂正态分布”的分布中多次抽取这个变量的值，并记录下最适合的多项式参数。然后，我会计算每个多项式参数的中位数和上下误差条（分别来自中位数的第84百分位和第16百分位）。

下面的代码是在Python 2.7.9中实现这个过程的。代码里有一个函数用来计算分裂正态值，还有用来从百分位计算误差和拟合多项式的功能。

希望这对你有帮助。

#! /bin/python

from random import choice, gauss
from numpy import polyfit

def split_normal(mus, sigmas_u68, sigmas_l68):
    """
    RET: A split-normal value.
    """

    split_normal = []

    for mu, sigma_u68, sigma_l68 in zip(mus, sigmas_u68, sigmas_l68):

        sigma = choice([sigma_u68, -sigma_l68])
        g = abs(gauss(0.0, 1.0)) * sigma + mu
        split_normal.append(g)

    return split_normal

def errors_84_16(x):
    """
    RET: 1-sigma upper/lower error bars from the 84/16th percentile
         from the median.
    """

    n = len(x)

    index_med = n / 2 # median.
    index_84 = int(round(n * 0.84135)) # 84th percentile from median.
    index_16 = int(round(n * 0.15865))

    x_sorted = sorted(x)
    x_med = x_sorted[index_med]
    x_u68 = x_sorted[index_84] - x_med # 1-sigma upper error.
    x_l68 = x_med - x_sorted[index_16] # 1-sigma lower error.

    return x_med, x_u68, x_l68

def assymetric_polyfit(x, y, y_u68, y_l68, n_mc=500):
    """
    DES: Solves y = a + b * x for assymentric y error bars.
    RET: [a, a_u68, a_l68, b, b_u68, b_l68].
    """

    a_mc = []
    b_mc = []

    for i in xrange(0, n_mc):
        y_mc = split_normal(y, y_u68, y_l68)
        pars = polyfit(x, y_mc, 2)
        a_mc.append(pars[2])
        b_mc.append(pars[1])


    a, a_u68, a_l68 = errors_84_16(a_mc)
    b, b_u68, b_l68 = errors_84_16(b_mc)

    return a, a_u68, a_l68, b, b_u68, b_l68

def example():
    """
    """

    x = [1.0, 2.0, 3.0, 4.0, 5.0]
    y = [5.0, 8.0, 11.0, 14.0, 17.0] # 2 + 3x
    y_u68 = [0.5, 0.5, 0.5, 0.5, 0.5]
    y_l68 = [1.5, 1.5, 1.5, 1.5, 1.5]

    pars = assymetric_polyfit(x, y, y_u68, y_l68)
    print(pars)

if __name__ == '__main__':

    example()