2014-10-17 72 views
2

我有多個時間序列,包括觀察值作爲行(索引時間戳,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 

回答

0

畢竟,我想出的答案是:

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

原來我並不真的需要引導數據,這樣我就可以估算基於95%的置信區間對於均值的標準誤差乘以正態分佈的0.975分位數,假設數據是正態分佈的(但這是另一個問題......)。

  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