【问题标题】:SciPy LeastSq Goodness of Fit EstimatorSciPy LeastSq 拟合优度估计器
【发布时间】:2011-11-27 03:38:20
【问题描述】:

我有一个使用 SciPy 的 leastsq 函数拟合的数据表面。

我想对leastsq 返回后的合身质量进行一些估计。我希望这将作为函数的返回包含在内,但如果是这样,它似乎没有明确记录。

是否有这样的返回,或者,除此之外,某些函数我可以传递我的数据,返回的参数值和拟合函数会给我一个拟合质量的估计值(R^2 或类似的)?

谢谢!

【问题讨论】:

    标签: python scipy data-fitting


    【解决方案1】:

    如果你像这样打电话给leastsq

    import scipy.optimize
    p,cov,infodict,mesg,ier = optimize.leastsq(
            residuals,a_guess,args=(x,y),full_output=True)
    

    在哪里

    def residuals(a,x,y):
        return y-f(x,a)
    

    然后,使用给定hereR^2 的定义,

    ss_err=(infodict['fvec']**2).sum()
    ss_tot=((y-y.mean())**2).sum()
    rsquared=1-(ss_err/ss_tot)
    

    你问什么infodict['fvec']?这是残差数组:

    In [48]: optimize.leastsq?
    ...
          infodict -- a dictionary of optional outputs with the keys:
                      'fvec' : the function evaluated at the output
    

    例如:

    import scipy.optimize as optimize
    import numpy as np
    import collections
    import matplotlib.pyplot as plt
    
    x = np.array([821,576,473,377,326])
    y = np.array([255,235,208,166,157])
    
    def sigmoid(p,x):
        x0,y0,c,k=p
        y = c / (1 + np.exp(-k*(x-x0))) + y0
        return y
    
    def residuals(p,x,y):
        return y - sigmoid(p,x)
    
    Param=collections.namedtuple('Param','x0 y0 c k')
    p_guess=Param(x0=600,y0=200,c=100,k=0.01)
    p,cov,infodict,mesg,ier = optimize.leastsq(
        residuals,p_guess,args=(x,y),full_output=True)
    p=Param(*p)
    xp = np.linspace(100, 1600, 1500)
    print('''\
    x0 = {p.x0}
    y0 = {p.y0}
    c = {p.c}
    k = {p.k}
    '''.format(p=p))
    pxp=sigmoid(p,xp)
    
    # You could compute the residuals this way:
    resid=residuals(p,x,y)
    print(resid)
    # [ 0.76205302 -2.010142    2.60265297 -3.02849144  1.6739274 ]
    
    # But you don't have to compute `resid` -- `infodict['fvec']` already
    # contains the info.
    print(infodict['fvec'])
    # [ 0.76205302 -2.010142    2.60265297 -3.02849144  1.6739274 ]
    
    ss_err=(infodict['fvec']**2).sum()
    ss_tot=((y-y.mean())**2).sum()
    rsquared=1-(ss_err/ss_tot)
    print(rsquared)
    # 0.996768131959
    
    plt.plot(x, y, '.', xp, pxp, '-')
    plt.xlim(100,1000)
    plt.ylim(130,270)
    plt.xlabel('x')
    plt.ylabel('y',rotation='horizontal')
    plt.grid(True)
    plt.show()
    

    【讨论】:

    • 上述rsquared的定义是否正确?对于每个bin,residual^2的值应该除以variance^2,然后在最后对所有结果求和。
    • @user1016260 上面的例子没有任何错误,但是是的,应该这样考虑。我的问题是针对具有不同数量参数的不同模型执行此操作。基本上使用这种 R^2 方法,没有办法考虑参数的数量。使用类似减少卡方的东西,它确实如此。那不是最好用于模型吗?这可能不适合在这里讨论。
    • @astromax 您可能需要考虑类似 . Akaike Information Criterion (AIC)
    • @iamgin: leastsq 使用Levenberg–Marquardt algorithm。根据链接,“在只有一个最小值的情况下,一个不知情的标准猜测......会正常工作;在有多个最小值的情况下,算法只有在初始猜测已经有点接近最终解决方案时才会收敛。”
    • @iamgin:上面使用的 sigmoid 函数的形式为y = c / (1 + np.exp(-k*(x-x0))) + y0。当x 很大时,exp(...) 基本上为零,所以y 倾向于c + y0。当x 是巨大的负数时,exp(...) 是巨大的,c / (huge) 趋向于 0,所以y 趋向于y0。所以y0是下限——图中140左右,c + y0是上限——250左右。所以c是跨度或差异250 - 140 ~ 110。
    猜你喜欢
    • 2014-02-19
    • 1970-01-01
    • 2012-07-01
    • 2020-07-31
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多