前言

numpy對(duì)python的意義非凡，在數(shù)據(jù)分析與機(jī)器學(xué)習(xí)領(lǐng)域?yàn)閜ython立下了汗馬功勞。現(xiàn)在用python搞數(shù)據(jù)分析或機(jī)器學(xué)習(xí)經(jīng)常使用的pandas、matplotlib、sklearn等庫(kù)，都需要基于numpy構(gòu)建。毫不夸張地說，沒有numpy，python今天在數(shù)據(jù)分析與機(jī)器學(xué)習(xí)領(lǐng)域只能是捉襟見肘。

什么是一門好的數(shù)據(jù)分析語(yǔ)言

數(shù)據(jù)分析面向的數(shù)據(jù)大多數(shù)是二維表。一門好的數(shù)據(jù)分析語(yǔ)言，首先需要能夠直接有個(gè)數(shù)據(jù)結(jié)構(gòu)存下這個(gè)二維表，然后要配上一套成熟的類SQL的數(shù)據(jù)操作接口，***要有一套好用的可視化工具。R語(yǔ)言就是一個(gè)極好的典范：用內(nèi)置的data.frame結(jié)構(gòu)做數(shù)據(jù)的存儲(chǔ);data.frame本身提供足夠強(qiáng)大的數(shù)據(jù)操作能力，另有dplyr、tidyr、data.table、plyr、reshape2等庫(kù)提供更好用更高效的數(shù)據(jù)操作能力;在繪圖上，除了基本的plot功能外，還提供了ggplot2這樣一套優(yōu)雅的繪圖語(yǔ)言，還通過htmlwidget庫(kù)與javascript各種繪圖庫(kù)建立了緊密的聯(lián)系，讓可視化的動(dòng)態(tài)展示效果更進(jìn)一步。Excel也是一個(gè)極好的例子，有單元格這種靈活的結(jié)構(gòu)為數(shù)據(jù)存儲(chǔ)做支撐，有大量的函數(shù)實(shí)現(xiàn)靈活的操作，也有強(qiáng)大的繪圖系統(tǒng)。

python目前在數(shù)據(jù)分析領(lǐng)域也已經(jīng)具備了相當(dāng)可觀的能力，包括pandas庫(kù)實(shí)現(xiàn)的DataFrame結(jié)構(gòu)，pandas本身提供的數(shù)據(jù)操作能力，matplotlib提供的數(shù)據(jù)可視化能力，而這一切都離不開numpy庫(kù)。

什么是一門好的機(jī)器學(xué)習(xí)語(yǔ)言

一般來講，一門好的機(jī)器學(xué)習(xí)語(yǔ)言在數(shù)據(jù)分析上也一定很吃得開，因?yàn)閿?shù)據(jù)分析往往是機(jī)器學(xué)習(xí)的基礎(chǔ)。但是機(jī)器學(xué)習(xí)的要求更高，因?yàn)樵谀Ｐ陀?xùn)練階段往往需要較為復(fù)雜的參數(shù)估計(jì)運(yùn)算，因此語(yǔ)言需要具備較強(qiáng)的科學(xué)計(jì)算能力。科學(xué)計(jì)算能力，最核心的就是矩陣運(yùn)算能力。關(guān)于矩陣運(yùn)算能力，這篇文章對(duì)各種語(yǔ)言有很好的比較。

如果沒有numpy，python內(nèi)部只能用list或array來表示矩陣。假如用list來表示[1,2,3]，由于list的元素可以是任何對(duì)象，因此list中所保存的是對(duì)象的指針，所以需要有3個(gè)指針和三個(gè)整數(shù)對(duì)象，比較浪費(fèi)內(nèi)存和CPU計(jì)算時(shí)間。python的array和list不同，它直接保存數(shù)值，和C語(yǔ)言的一維數(shù)組比較類似，但是不支持多維，表達(dá)形式很簡(jiǎn)陋，寫科學(xué)計(jì)算的算法很難受。numpy彌補(bǔ)了這些不足，其提供的ndarray是存儲(chǔ)單一數(shù)據(jù)類型的多維數(shù)組，且采用預(yù)編譯好的C語(yǔ)言代碼，性能上的表現(xiàn)也十分不錯(cuò)。

python***的機(jī)器學(xué)習(xí)庫(kù)sklearn構(gòu)建在numpy之上，提供了各種標(biāo)準(zhǔn)機(jī)器學(xué)習(xí)模型的訓(xùn)練與預(yù)測(cè)接口，其中模型訓(xùn)練接口的內(nèi)部實(shí)現(xiàn)是基于numpy庫(kù)實(shí)現(xiàn)的。比如很常見的線性回歸模型，參數(shù)估計(jì)調(diào)用的是numpy.linalg.lstsq函數(shù)。

numpy的核心結(jié)構(gòu)：ndarray

以下內(nèi)容摘錄自用Python做科學(xué)計(jì)算

a = np.array([[0,1,2],[3,4,5],[6,7,8]], dtype=np.float32) 

ndarray是numpy的核心數(shù)據(jù)結(jié)構(gòu)。我們來看一下ndarray如何在內(nèi)存中儲(chǔ)存的：關(guān)于數(shù)組的描述信息保存在一個(gè)數(shù)據(jù)結(jié)構(gòu)中，這個(gè)結(jié)構(gòu)引用兩個(gè)對(duì)象，一塊用于保存數(shù)據(jù)的存儲(chǔ)區(qū)域和一個(gè)用于描述元素類型的dtype對(duì)象。

 數(shù)據(jù)存儲(chǔ)區(qū)域保存著數(shù)組中所有元素的二進(jìn)制數(shù)據(jù)，dtype對(duì)象則知道如何將元素的二進(jìn)制數(shù)據(jù)轉(zhuǎn)換為可用的值。數(shù)組的維數(shù)、大小等信息都保存在ndarray數(shù)組對(duì)象的數(shù)據(jù)結(jié)構(gòu)中。

strides中保存的是當(dāng)每個(gè)軸的下標(biāo)增加1時(shí)，數(shù)據(jù)存儲(chǔ)區(qū)中的指針?biāo)黾拥淖止?jié)數(shù)。例如圖中的strides為12,4，即第0軸的下標(biāo)增加1時(shí)，數(shù)據(jù)的地址增加12個(gè)字節(jié)：即a[1,0]的地址比a[0,0]的地址要高12個(gè)字節(jié)，正好是3個(gè)單精度浮點(diǎn)數(shù)的總字節(jié)數(shù);第1軸下標(biāo)增加1時(shí)，數(shù)據(jù)的地址增加4個(gè)字節(jié)，正好是單精度浮點(diǎn)數(shù)的字節(jié)數(shù)。

以下內(nèi)容總結(jié)自Numpy官方文檔Numpy basics

關(guān)于ndarray的索引方式，有以下幾個(gè)重點(diǎn)需要記住：

• 雖然x[0,2] = x0，但是前者效率比后者高，因?yàn)楹笳咴趹?yīng)用***個(gè)索引后需要先創(chuàng)建一個(gè)temporary array，然后再應(yīng)用第二個(gè)索引，***找到目標(biāo)值。
• 分片操作不會(huì)引發(fā)copy操作，而是創(chuàng)建原ndarray的view;他們所指向的內(nèi)存是同一片區(qū)域，無論是修改原ndarray還是修改view，都會(huì)同時(shí)改變二者的值。
• index array和boolean index返回的是copy，不是view。

關(guān)于上面列舉的分片操作不會(huì)引發(fā)copy操作，我們來進(jìn)一步探討一下。先看一下numpy的例子：

 再來看一下R的例子：

 可以看到numpy和R在矩陣的分片操作有不同的設(shè)計(jì)理念：在R里分片操作會(huì)引起數(shù)據(jù)的復(fù)制，在numpy里不會(huì)。事實(shí)上，R的設(shè)計(jì)理念很多時(shí)候可以用一句話來概括：copy on modify，一旦對(duì)數(shù)據(jù)有修改就會(huì)引起內(nèi)存上的復(fù)制操作，這個(gè)操作要花不少時(shí)間，因此經(jīng)常會(huì)聽到人們抱怨R費(fèi)內(nèi)存且速度慢。所以，我們可以看到numpy在處理這件事情上明顯要用心很多，根據(jù)場(chǎng)景設(shè)計(jì)了不同的策略，不是簡(jiǎn)單地采用R的一刀切方式。當(dāng)然，這也帶來了一些學(xué)習(xí)成本，需要對(duì)numpy足夠熟悉才能避免踩坑。R社區(qū)里對(duì)copy on modify的哲學(xué)也有詬病并在努力改變，比如同是data.frame操作庫(kù)的data.table和dplyr，data.table性能比dplyr高很多，部分原因也是data.table規(guī)避了copy on modify的方式。

Structured Array

根據(jù)numpy的官方文檔，定義結(jié)構(gòu)化數(shù)組有四種方式。本文采用字典方法，通過定義一個(gè)dtype對(duì)象實(shí)現(xiàn)，需要指定的鍵值有names和formats。

persontype = np.dtype({ 
        'names': ['name', 'age', 'weight'],  
        'formats': ['S32', 'i', 'f'] 
    }) 
a = np.array([("Zhang", 32, 75.5), ("Wang", 24, 65.2)], dtype=persontype)  

我們用IPython的計(jì)時(shí)函數(shù)看一下提取數(shù)據(jù)的效率： 

%timeit a[1] 
%timeit a['name'] 
%timeit a[1]['name'] 
%timeit a['name'][1]  

輸出結(jié)果如下：

The slowest run took 46.83 times longer than the fastest. This could mean that an intermediate result is being cached. 
1000000 loops, best of 3: 153 ns per loop 
The slowest run took 34.34 times longer than the fastest. This could mean that an intermediate result is being cached. 
10000000 loops, best of 3: 174 ns per loop 
The slowest run took 13.00 times longer than the fastest. This could mean that an intermediate result is being cached. 
1000000 loops, best of 3: 1.08 µs per loop 
The slowest run took 9.84 times longer than the fastest. This could mean that an intermediate result is being cached. 
1000000 loops, best of 3: 412 ns per loop  

從上面的結(jié)果，我們發(fā)現(xiàn)，獲取相同的數(shù)據(jù)有多種操作，不同的操作性能差別很大。我做了一個(gè)推測(cè)，純粹是瞎猜：numpy在建立結(jié)構(gòu)化數(shù)組時(shí)，將整個(gè)結(jié)構(gòu)體連續(xù)存儲(chǔ)在一起，即按行存儲(chǔ)，因此a[1]的速度最快;但是為了保證提取列的效率，對(duì)a['name']建立了索引，因此a['name']的效率也很高;但是這個(gè)索引只對(duì)整個(gè)a起作用，如果輸入只有a的一部分，仍然需要遍歷整個(gè)a，去提取出對(duì)應(yīng)的數(shù)據(jù)，因此a[1]['name']比a['name'][1]的效率差很多。

實(shí)例

基于numpy過濾抖動(dòng)與填補(bǔ)

時(shí)間序列數(shù)據(jù)經(jīng)常會(huì)發(fā)現(xiàn)兩種情況：一種是抖得特別厲害，說明數(shù)據(jù)不穩(wěn)定不可信，支撐這個(gè)結(jié)果的數(shù)據(jù)量不夠;另一種是一動(dòng)不動(dòng)的一條直線，這往往是算法填充出來的默認(rèn)值，不是實(shí)際值。這些數(shù)據(jù)對(duì)于挖掘來說是噪音，應(yīng)該過濾掉。我們使用numpy來完成這個(gè)任務(wù)。抖動(dòng)的特點(diǎn)是頻繁跳動(dòng)，即一階差分有很多值絕對(duì)值比0大很多，那么我們將這些跳動(dòng)的點(diǎn)抓出來，統(tǒng)計(jì)下這些點(diǎn)之間的區(qū)間長(zhǎng)度，如果區(qū)間長(zhǎng)度過小，認(rèn)為是抖動(dòng)過多。填補(bǔ)的特點(diǎn)是數(shù)值長(zhǎng)期不變，即一階差分有很多值為0，那么我們統(tǒng)計(jì)一下連續(xù)為0的區(qū)間長(zhǎng)度分布，如果區(qū)間長(zhǎng)度過長(zhǎng)，比如連續(xù)填補(bǔ)了1小時(shí)，或者出現(xiàn)多個(gè)填補(bǔ)了30分鐘的區(qū)間，我們認(rèn)為是填補(bǔ)過多。

我們需要對(duì)跳點(diǎn)進(jìn)行定義：一階差分的絕對(duì)值超過dev_thresh，一階差分/max(基準(zhǔn)1，基準(zhǔn)2)的絕對(duì)值超過ratio_thresh。 

def jump(speed_array, dev_thresh, ratio_thresh): 
  diff_array = np.diff(speed_array, axis=0) 
  diff_array = diff_array.astype(np.float64) 
  ratio_array = diff_array/np.maxium(speed_array[:-1], speed_array[1:]) 
  ret_array = np.zeros(diff_array.size, dtype=np.int8) 
  for i in range(diff_array.size): 
    if abs(diff_array[i]) > diff_thresh and abs(ratio_array[i]) > ratio_thresh: 
      ret_array[i] = 1 
  return ret_array 
 
def interval(jump_array): 
  jump_idx = np.array([0] + [i for i,x in enumerate(jump_array) if x != 0] + [jump_array.size]) 
  interval_size = np.diff(jump_idx) 
  return interval_size 
 
def is_jump_too_much(interval_size): 
  flag = 0 
  if np.mean(interval_size) <= 10 or np.max(interval_size) <= 30: 
      flag = 1 
  return flag 
 
def is_fill_too_much(interval_size): 
  flag = 0 
  bin_array = np.bincount(interval_size) 
  if ( len(bin_array) >= 30 or  
       ( len(bin_array) >= 11 and np.sum(bin_array[10:]) >= 4 ) or  
       ( len(bin_array) >= 7 and np.sum(bin_array[6:]) >= 20 ) 
      ): 
      flag = 1 
  return flag  

基于numpy的局部趨勢(shì)擬合

用線性回歸可以得到時(shí)間序列的趨勢(shì)。 

def get_ts_trend(ts_array): 
  x = np.arange(0, len(ts_array), 1) 
  y = ts_array 
  A = np.vstack([x, np.ones(len(x))]).T 
  m, c = np.linalg.lstsq(A, y)[0] 
  return m  

堵點(diǎn)判別

交通數(shù)據(jù)比較復(fù)雜，不純粹是時(shí)間序列問題，而是時(shí)空數(shù)據(jù)，需要同時(shí)考慮時(shí)間關(guān)系和空間關(guān)系。本節(jié)介紹一個(gè)經(jīng)典特征的提取：堵點(diǎn)判別。

假設(shè)我們空間上有5個(gè)link，上游2個(gè)，自身1個(gè)，下游2個(gè);觀察5個(gè)時(shí)間點(diǎn)的擁堵狀態(tài)。判斷當(dāng)前l(fā)ink是不是堵點(diǎn)——即自身是拓?fù)渲?**個(gè)發(fā)生擁堵的點(diǎn);發(fā)生擁堵后，擁堵是擴(kuò)散的。 

def detect_congest_point(congest_array): 
    first_congest_flag = False 
    disperse_congest_flag = True 
    idx = np.where(congest_array == 1) 
    if idx[1][0] == congest_array.shape[1]/2: 
        first_congest_flag = True 
     
    disperse_dict = {} 
    for k in range(len(idx[0])): 
        if disperse_dict.has_key(idx[0][k]): 
            disperse_dict[idx[0][k]].append(idx[1][k]) 
        else: 
            disperse_dict[idx[0][k]] = [idx[1][k]] 
    sorted_disperse_list = sorted(disperse_dict.iteritems(), key=lambda d:d[0]) 
    for i in range(1, len(sorted_disperse_list)): 
        if not set(sorted_disperse_list[i-1][1]) <= set(sorted_disperse_list[i][1]): 
            disperse_congest_flag = False 
     
    return first_congest_flag and disperse_congest_flag  

關(guān)于作者：丹追兵：數(shù)據(jù)分析師一枚，編程語(yǔ)言python和R，使用Spark、Hadoop、Storm、ODPS。本文出自丹追兵的pytrafficR專欄，轉(zhuǎn)載請(qǐng)注明作者與出處：https://segmentfault.com/blog...