(Python)使用scikits bootstrap估计回归参数置信区间

Question

我刚开始尝试通过scikits提供一个很好的bootstrapping包:https: //github.com/cgevans/scikits-bootstrap

但是我在尝试从线性回归估计相关系数的置信区间时遇到了问题.返回的置信区间完全在原始统计量的范围之外.

这是代码:

import numpy as np
from scipy import stats
import bootstrap as boot

np.random.seed(0)
x     = np.arange(10)
y     = 10 + 1.5*x + 2*np.random.randn(10)
r0    = stats.linregress(x, y)[2]

def my_function(y):
    return stats.linregress(x, y)[2]

ci    = boot.ci(y, statfunction=my_function, alpha=0.05, n_samples=1000, method='pi')

这产生ci = [ - 0.605,0.644]的结果,但原始统计量为r0 = 0.894.

我在R中试过这个并且似乎在那里工作得很好:ci横跨r0如预期的那样.

请帮忙!

Answer 1

你能提供你的R代码吗？我有兴趣知道如何处理这个问题.

这里的问题是你只是传递y给boot.ci,但每次运行my_function时,它都会使用整个原始的x(注意缺少x输入到my_function).Bootstrapping将统计函数应用于重采样数据,因此如果您使用原始数据x和样本应用统计函数y,则会产生无意义的结果.这就是为什么BCA方法不会在所有的工作,实际上是:它不能应用于您的统计功能以刀切样,不具有相同数目的元素.

样本沿轴0(行)进行,因此如果要将多个1D数组传递给统计函数,可以使用多个列:可以使用xy = vstack((x,y)).T,然后使用从这些列获取数据的统计函数:

def my_function(xysample):
    return stats.linregress(xysample[:,0], xysample[:,1])[2]

或者,如果您想避免弄乱数据,可以定义一个对索引进行操作的函数,然后将索引传递给boot.ci:

def my_function2(i):
    return stats.linregress(x[i], y[i])[2]

boot.ci(np.arange(len(x)), statfunction=my_function2, alpha=0.05, n_samples=1000, method='pi')

请注意,在这两种情况下,BCA都可以工作,所以你也可以使用method ='bca',除非你确实想要使用百分比间隔; BCA几乎总是更好.

我确实意识到这两种方法都不太理想.老实说,我从来没有需要将这样的多个数组传递给我的统计功能,而且大多数人都可能使用mean它们作为他们的统计功能.我认为这里最好的想法可能是允许传递大小相等的[0]数组的列表,例如,boot.ci([x,y],...)然后同时对所有这些数组进行采样,并将它们全部作为单独的参数传递给statfunction.在那种情况下,你可以有一个my_function(x,y).我会看看我是否可以这样做,但是如果你能告诉我你的R代码,那就太好了,因为我想知道是否有更好的方法可以解决这个问题.

更新:

在最新版本的scikits.bootstrap(v0.3.1)中,可以提供一个数组元组,它们的样本将作为独立的参数传递给statfunction.此外,statfunction可以提供数组输出,并且将为输出中的每个点计算置信区间.因此,现在这很容易做到.以下将给出linregress的每个输出的置信区间:

cis = boot.ci( (x,y), statfunction=stats.linregress )

cis[:,2] 在这种情况下,将是所需的置信区间.