我有多个时间序列,包括观察值作为行(索引时间戳,n> 5000)和变量作为列(n = 419)。我在时间序列中选择了N百分比,然后打电话给groupby逐个分组。我想要的是平均值,标准差,然后是每年的95%置信区间。我可以通过下面的代码获得平均值和std,但我需要调用单独的引导程序函数以获得每年和每个组的95%CI:引导熊猫上的groupby对象
这是对分组数据看起来像什么的一瞥:(2013年和28列共有86行,数据始于20世纪70年代)。我需要为每个分组中的每一列使用'bootsrap'。
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
` 这里是我的GROUPBY和引导功能:
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]
一行。我怀疑我需要比bootstrap中的数据字段使用“分组”更具体,但我不确定如何执行此操作。我是否需要将此应用为lambda函数?或者也许用for循环?任何建议将不胜感激。
看完此页:Pandas, compute many means with bootstrap confidence intervals for plotting并尝试使用scikit bootstrap函数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])
如果这有效,我会为33年的每一年的28列中的每一列分配lo和hi,而我有一个有序的数组,似乎没有任何实际价值......这是yrly
其中包含日志转换组,每年通过手段启动与上面的数组不同,说唱的CIs接近这些数字。
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