scipy python中最小平方拟合的置信区间

发布于 2024-11-04 05:12:26 字数 55 浏览 5 评论 0 原文

如何在Python中计算最小二乘拟合（scipy.optimize.leastsq）的置信区间？

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

雪若未夕 2024-11-11 05:12:26

我会使用引导方法。
请参阅此处：http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html< /a>

噪声高斯的简单示例：

x = arange(-10, 10, 0.01)

# model function
def f(p):
    mu, s = p
    return exp(-(x-mu)**2/(2*s**2))

# create error function for dataset    
def fff(d):
    def ff(p):
        return d-f(p)
    return ff

# create noisy dataset from model
def noisy_data(p):
    return f(p)+normal(0,0.1,len(x))

# fit dataset to model with least squares    
def fit(d):
    ff = fff(d)
    p = leastsq(ff,[0,1])[0]
    return p

# bootstrap estimation        
def bootstrap(d):
    p0 = fit(d)
    residuals = f(p0)-d
    s_residuals = std(residuals)

    ps = []
    for i in range(1000):
        new_d = d+normal(0,s_residuals,len(d))
        ps.append(fit(new_d))

    ps = array(ps)
    mean_params = mean(ps,0)
    std_params = std(ps,0)

    return mean_params, std_params

data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)

print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996

I would use bootstrapping method.

See here: http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html

Simple example for noisy gaussian:

x = arange(-10, 10, 0.01)

# model function
def f(p):
    mu, s = p
    return exp(-(x-mu)**2/(2*s**2))

# create error function for dataset    
def fff(d):
    def ff(p):
        return d-f(p)
    return ff

# create noisy dataset from model
def noisy_data(p):
    return f(p)+normal(0,0.1,len(x))

# fit dataset to model with least squares    
def fit(d):
    ff = fff(d)
    p = leastsq(ff,[0,1])[0]
    return p

# bootstrap estimation        
def bootstrap(d):
    p0 = fit(d)
    residuals = f(p0)-d
    s_residuals = std(residuals)

    ps = []
    for i in range(1000):
        new_d = d+normal(0,s_residuals,len(d))
        ps.append(fit(new_d))

    ps = array(ps)
    mean_params = mean(ps,0)
    std_params = std(ps,0)

    return mean_params, std_params

data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)

print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996

回复收藏 0 原文

心作怪 2024-11-11 05:12:26

我不确定你所说的置信区间是什么意思。

一般来说，leastsq 对您尝试最小化的函数了解不多，因此它无法真正给出置信区间。然而，它确实返回了 Hessian 矩阵的估计，换句话说，就是二阶导数对多维问题的推广。

正如函数的文档字符串中所暗示的，您可以使用该信息以及残差（拟合解与实际数据之间的差异）来计算参数估计的协方差，这是置信区间的局部猜测。

请注意，这只是局部信息，我怀疑严格来说，只有当您的目标函数是严格凸的时，您才能得出结论。我没有关于该声明的任何证据或参考:)。

回复收藏 0 原文

烟凡古楼 2024-11-11 05:12:26

估计置信区间 (CI) 的最简单方法是将标准误差（标准偏差）乘以常数。要计算常数，您需要知道自由度 (DOF) 数以及要计算 CI 的置信水平。
以这种方式估计的 CI 有时称为渐近 CI。
您可以在 Motulsky & 的“使用线性和非线性回归将模型拟合到生物数据”中阅读更多相关信息。克里斯托普洛斯 (Google 图书）。同一本书（或非常相似）可以免费作为作者软件的手册。

您还可以阅读如何使用 C++ Boost.Math 库计算 CI。在此示例中，CI 是针对一个变量的分布计算的。在最小二乘拟合的情况下，DOF 不是 N-1，而是 NM，其中 M 是参数数量。在 Python 中做同样的事情应该很容易。

这是最简单的估计。我不知道zephyr提出的bootstrapping方法，但它可能比我写的方法更可靠。