在pandas的groupby对象上进行自助法

Question

我有多个时间序列数据，每一行是观察值（用时间戳来标记，超过5000行），每一列是变量（有419列）。我想从这些时间序列中选择一定比例的包含数据，然后按年份分组。我的目标是计算每一年的平均值、标准差，以及95%的置信区间。我可以用下面的代码轻松得到平均值和标准差，但我需要调用一个单独的自助法（bootstrap）函数来为每一年和每个组计算95%的置信区间：

这是分组后数据的一个快照：（2013年有86行数据，28列，数据从1970年代开始）。我需要对分组后的每一列在每一年中使用“自助法”。

for year, group in grouped:
print year
print group

2013
                  101        102        103        104       105        109     
2013-04-02    3162.84    4136.02   77124.56       0.00    973.18    9731.81   
2013-04-04    1033.81    5464.44   87283.30    3692.19   4282.94     295.37   
2013-04-04     640.75    4164.87  131033.14    2563.00   1121.31     961.12   
2013-04-10     246.87    4196.84   88380.57    4443.72    493.75    1234.37   
2013-04-13       0.00    8300.49  114291.42   10003.16    212.83    6385.00   

` 这里是我用于分组和自助法的函数：

def gbY_20pct(nm): # sort into 20% timeseries inclusion, groupby year, take mean for year
        nm1=nm.replace('0', np.nan) # remove 0 for logical count
        coun=nm1.count(axis=0,numeric_only=True)
        pct=(coun/len(nm1)) *100
        pCount=pct.loc[pct >= 20]
        nm1=nm.loc[:, pCount.index]
        grouped = nm1.groupby(nm1.index.map(lambda x: x.year))
        yrly=grouped.mean().astype(int)
        yrly_coun=grouped.count().astype(int)
        yrly_std=grouped.std().astype(int)
        yrly_max=grouped.max().astype(int)
        yrM1=yrly.join(yrly_std, lsuffix=' mean', rsuffix=' std', how='outer')
        yrM2=yrly_max.join(yrly_coun, lsuffix=' max', rsuffix=' count', how='outer')
        data=yrM1.join(yrM2, how='outer')
        return data

`

import numpy as np
import numpy.random as npr  
def bootstrap(data, num_samples, statistic, alpha):
    """Returns bootstrap estimate of 100.0*(1-alpha) CI for statistic."""
    n = len(data)
    idx = npr.randint(0, n, (num_samples, n))
    samples = data[idx]
    stat = np.sort(statistic(samples, 1))
    return (stat[int((alpha/2.0)*num_samples)],
            stat[int((1-alpha/2.0)*num_samples)])

为了测试代码，我手动调用它（分组已经定义，函数还没有结束）

from bootstrap import bootstrap
low, high = bootstrap(grouped, 100000, np.mean, 0.05)
Traceback (most recent call last):

  File "<ipython-input-49-cd362c7908d1>", line 1, in <module>
    low, high = bootstrap(grouped, 100000, np.mean, 0.05)

  File "bootstrap.py", line 14, in bootstrap

  File "C:\Users\ryan.morse\AppData\Local\Continuum\Anaconda\lib\site-packages\pandas\core\groupby.py", line 2991, in __getitem__
    bad_keys = list(set(key).difference(self.obj.columns))

TypeError: unhashable type: 'numpy.ndarray'

问题出在 samples = data[idx] 这一行。我怀疑在自助法中使用“分组”作为数据字段时需要更具体一些，但我不太确定该怎么做。我需要用lambda函数吗？或者用for循环？任何建议都非常感谢。

在查看了这个页面后：Pandas，计算多个均值并使用自助法置信区间进行绘图，并尝试了scikit自助法函数 https://scikits.appspot.com/bootstrap，我测试了上面定义的函数，发现它运行得快得多，结果也相当不错。

---

编辑：

我在想这样的做法可能有效，但我仍然无法正确写出语法：

groups=dict(list(grouped)) # this allows me to visualize the data and call values

for key, value in groups.iteritems():
low_i, high_i = bootstrap(groups.values(), 100000, np.mean, 0.05) 

Traceback (most recent call last):

  File "<ipython-input-36-7a8e261d656e>", line 2, in <module>
    low_i, high_i=bootstrap(groups.values(), 10000, np.mean, 0.05)

  File "<ipython-input-15-3ce4acd651dc>", line 7, in bootstrap
    samples = data[idx]

TypeError: only integer arrays with one element can be converted to an index

我不确定如何为自助法函数调用“数据”，以及如何遍历所有年份并保持所有年份的低值和高值（要么在同一个数据框中，要么在两个单独的数据框中）。

任何帮助都将不胜感激...

---

编辑2 我也可以很容易地使用lambda函数，但我似乎无法得到正确的输出：

for col, group in nm1.groupby(nm1.index.year):
    lo,hi=bootstrap(group,1000, np.mean, 0.05)

lo
Out[117]: 
array([ 0.05713616,  0.30724739,  0.39592714,  0.55113183,  0.68623155,
        0.69493923,  0.73513661,  0.84086099,  0.85882618,  0.86698939,
        0.99399694,  1.04415927,  1.06553914,  1.11306698,  1.15344871,
        1.27943327,  1.43275895,  1.81076036,  2.21647657,  2.37724615,
        2.39004626,  2.43154256,  2.89940325,  3.02234954,  3.30773642,
        3.96535146,  3.98973744,  4.38873853])

hi
Out[118]: 
array([ 0.20584822,  0.38832222,  0.42140066,  0.48615202,  0.59686031,
        0.67388397,  0.84269082,  0.84532503,  0.87078368,  0.9033272 ,
        0.90765817,  0.97523759,  0.99186096,  1.01668772,  1.06681722,
        1.18205259,  1.38524423,  1.79908484,  2.22314773,  2.33789105,
        2.5521743 ,  2.64242269,  2.88851233,  2.94387756,  3.44294791,
        3.63914938,  3.99185026,  4.36450246])

如果这能成功，我会得到每28列和33年对应的低值和高值，但实际上我得到的是一个有序的数字数组，似乎没有任何实际意义……这是yrly的一个片段，它包含每年的对数变换后的分组均值，自助法的置信区间应该接近这些数字，而不是上面的数组。

           101       102       103       104       105       109       135  
1978  3.416638  3.701268  3.828442  2.911944  2.687491  2.076515  1.232035   
1979  2.710939  3.172061  4.234109  1.666818  3.390646  1.355179  3.003813   
1980  2.652617  2.375495  3.316380  1.101594  2.220028  1.195449  1.998862   
1981  3.363424  3.485015  3.441784  2.242618  2.256745  1.719140  1.150454   
1982  2.791865  2.019883  4.093960  1.038226  2.106627  1.180935  2.456144   
1983  2.597307  2.213450  4.458691  1.274352  2.820910  1.705242  3.452762   
1984  3.042197  4.023952  3.816964  2.499883  2.445258  1.769485  1.690180   
1985  2.669850  2.162608  3.600731  1.400102  1.845218  1.234235  2.517108   
1986  3.597527  2.763436  2.790792  1.410343  2.116275  1.042812  1.528532

Answer 1

经过这些思考，我得出的结论是：

import scipy.stats
ci = grouped.aggregate(lambda x: scipy.stats.sem(x, ddof=1) * 1.96) #use mean +(-) ci to get 95% conf interval

其实我并不需要用引导法来处理数据，所以我可以直接根据平均值的标准误差，乘以正态分布的0.975分位数，来估算95%的置信区间。前提是这些数据是符合正态分布的（不过这又是另一个话题...）。

           101       102       103       104       105       109       135
1978  0.230630  0.191651  0.168648  0.282588  0.237939  0.288924  0.257476   
1979  0.192579  0.147305  0.120740  0.225826  0.145646  0.266530  0.199315   
1980  0.189258  0.195263  0.182756  0.166479  0.166401  0.172550  0.189483   
1981  0.200727  0.169663  0.184478  0.232392  0.198591  0.230457  0.194084   
1982  0.271740  0.267881  0.164450  0.248718  0.246636  0.260973  0.253430   
1983  0.253495  0.279114  0.116744  0.266888  0.236672  0.317195  0.155766