本文在 Python 中用箱線圖、傅里葉變換、熵、自相關(guān)和 PCA 分析時(shí)間序列數(shù)據(jù)。數(shù)據(jù)可視化是任何數(shù)據(jù)相關(guān)項(xiàng)目中最重要的階段之一。根據(jù)數(shù)據(jù)可視化的對(duì)象有：

　　1.數(shù)據(jù)可視化報(bào)告結(jié)果。

　　2.數(shù)據(jù)可視化來(lái)分析數(shù)據(jù)，換句話說(shuō)，數(shù)據(jù)科學(xué)家內(nèi)部使用的可視化來(lái)提取有關(guān)數(shù)據(jù)的信息，然后實(shí)施模型。

　　本文主要關(guān)注后一種，因?yàn)樗忉屃艘恍┯兄诜治鰰r(shí)間序列數(shù)據(jù)的方法。

　　什么是時(shí)間序列？

　　基本數(shù)值時(shí)間序列是有序的、帶時(shí)間戳的觀測(cè)值(測(cè)量值)的集合，其中每個(gè)觀測(cè)值都是從同一測(cè)量過(guò)程中獲得的數(shù)值標(biāo)量。

　　什么是時(shí)間戳？

　　在我們將“時(shí)間”捕獲為數(shù)據(jù)點(diǎn)之前，我們不會(huì)深入探討需要精確定義的許多細(xì)節(jié)(準(zhǔn)確性、格式、日歷約定、時(shí)區(qū)等等)。我們將時(shí)間戳定義為具有所需精度的時(shí)間點(diǎn)的表示就足夠了。例如，這可能是根據(jù)某個(gè)日歷的日期約定(例如“08-06-2020”)，或者自 1970 年以來(lái)以整數(shù)表示的毫秒數(shù)(這實(shí)際上是 UNIX 紀(jì)元約定!)

　　Python類(lèi)庫(kù)

　　首先，這些是與 notebook 一起使用的庫(kù)。大多數(shù)代碼都圍繞 NumPy 和 Pandas庫(kù)，因?yàn)閿?shù)據(jù)主要以 Pandas Dataframe 表現(xiàn)的 NumPy 數(shù)組。

　　導(dǎo)入文件

　　下載數(shù)據(jù)后，運(yùn)行以下代碼將其導(dǎo)入。

　　正如所觀察到的，數(shù)據(jù)包含六個(gè)傳感器的傳感器數(shù)據(jù)、每個(gè)數(shù)據(jù)點(diǎn)的日期時(shí)間以及機(jī)器狀態(tài)。這是“BROKEN”、“NORMAL”或“RECOVERING”，但為了簡(jiǎn)化可視化，它被分組如下：

　　在任何編程語(yǔ)言中使用日期時(shí)間總是具有挑戰(zhàn)性的，Python 也不例外。盡管處理日期時(shí)間有多種方法，但這里使用函數(shù) pandas.to_datetime 將 datetime 列(讀取為字符串)轉(zhuǎn)換為時(shí)間戳。

　　數(shù)據(jù)預(yù)處理

　　在進(jìn)行可視化之前，分析了本次數(shù)據(jù)的重復(fù)值和缺失值。并且刪除重復(fù)項(xiàng)的函數(shù)：

def drop_duplicates(df: pd.DataFrame(), subset: list = ['DATE_TIME']) -> pd.DataFrame():

    df = df.drop_duplicates((subset))

    return df

　　填充缺失值的函數(shù)：

def fill_missing_date(df: pd.DataFrame(), column_datetime: str ='DATE_TIME'):

    print(f'輸入形狀: {df.shape}')

    data_s = df.drop([column_datetime], axis=1)

    datetime_s = df[column_datetime].astype(str)

    start_date = min(df[column_datetime])

    end_date = max(df[column_datetime])

    date_s = pd.date_range(start_date, end_date, freq="min").strftime('%Y-%m-%d %H:%M:%S')

    data_processed_s = []

    for date_val in date_s:

        pos = np.where(date_val == datetime_s)[0]        

        assert len(pos) in [0, 1]

        if len(pos) == 0:

            data = [date_val] + [0] * data_s.shape[1]

        elif len(pos) == 1:

            data = [date_val] + data_s.iloc[pos].values.tolist()[0]

        data_processed_s.append(data)

    df_processed = pd.DataFrame(data_processed_s, columns=[column_datetime] + data_s.columns.values.tolist())

    df_processed[column_datetime] = pd.to_datetime(df_processed[column_datetime])

    print(f'輸出形狀: {df_processed.shape}')

    return df_processed

　　這是預(yù)處理階段的整個(gè)管道。此外，數(shù)據(jù)分為輸入數(shù)據(jù)和輸出數(shù)據(jù)。

　　輸入形狀：(10081, 7)

　　輸出形狀：(10081, 2)

　　數(shù)據(jù)可視化

　　現(xiàn)在，準(zhǔn)備開(kāi)始數(shù)據(jù)可視化。這是傳感器數(shù)據(jù)和異常情況的圖。完整代碼可以在公眾號(hào):機(jī)器學(xué)習(xí)研習(xí)院 后臺(tái)回復(fù) 時(shí)間序列可視化獲取.

　　均值和標(biāo)準(zhǔn)

　　可以更好地總結(jié)數(shù)據(jù)隨時(shí)間變化的行為的最基本圖之一是均值標(biāo)準(zhǔn)圖，我們?cè)谄渲酗@示按時(shí)間范圍分組的均值和標(biāo)準(zhǔn)差。這主要有助于分析指定時(shí)間范圍內(nèi)的基線和噪聲。

df_data_hour = df_data.groupby(pd.Grouper(key='datetime', axis=0, freq='H')).mean()

df_labels_hour = df_labels.groupby(pd.Grouper(key='datetime', axis=0, freq='H')).sum()

df_rollmean = df_data_hour.resample(rule='D').mean()

df_rollstd = df_data_hour.resample(rule='D').std()

for name in df.columns:

    if name not in ['datetime', 'machine_status']:

        fig, axs = plt.subplots(1, 1, figsize=(15, 2))

        axs.plot(df_data_hour[name], color='blue', label='Original')

        axs.plot(df_rollmean[name], color='red', label='Rolling Mean')

        plt.plot(df_rollstd[name], color='black', label='Rolling Std' )

        axs.set_title(name)

        plt.legend()

        plt.show()

　　箱形圖

　　另一個(gè)有趣的圖表是通過(guò)箱線圖顯示的。箱線圖是一種通過(guò)四分位數(shù)以圖形方式顯示數(shù)值數(shù)據(jù)的局部性、擴(kuò)散性和偏度組的方法。有兩個(gè)主要框表示從第25個(gè)百分位數(shù)到第75個(gè)百分位數(shù)的數(shù)據(jù)，兩者之間用分布的中位數(shù)隔開(kāi)。除了盒子之外，還有從盒子延伸出來(lái)的晶須，表明上四分位和下四分位之外的變異性。與數(shù)據(jù)集其他部分顯著不同的異常值也被繪制為箱線圖上須之外的單獨(dú)點(diǎn)。

　　這一個(gè)類(lèi)似于平均和標(biāo)準(zhǔn)圖，因?yàn)樗砻鲾?shù)據(jù)的平穩(wěn)性。但是，它也可以顯示異常值，這有助于從視覺(jué)上檢測(cè)異常和數(shù)據(jù)之間的任何關(guān)系。

　　傅里葉變換

　　快速傅里葉變換(FFT)是一種計(jì)算序列離散傅里葉變換的算法。這種類(lèi)型的圖很有趣，因?yàn)樗翘幚頃r(shí)間序列時(shí)特征提取的主要方法之一。通常的做法不是用時(shí)間序列來(lái)訓(xùn)練模型，而是應(yīng)用傅里葉變換來(lái)提取頻率，然后訓(xùn)練模型。

　　為此，我們必須選擇一個(gè)滑動(dòng)窗口來(lái)計(jì)算FFT。滑動(dòng)窗口越寬，頻率數(shù)越高。缺點(diǎn)是您將得到更少的時(shí)間戳，從而丟失數(shù)據(jù)的時(shí)間分辨率。當(dāng)減小窗口的大小時(shí)，我們得到了相反的結(jié)果：更少的頻率但更高的時(shí)間分辨率。然后，窗口的大小應(yīng)該取決于任務(wù)。

　　FFT的滑動(dòng)窗口 對(duì)于如下圖所示，我選擇了一個(gè)包含64個(gè)數(shù)據(jù)的時(shí)間窗口。因此，頻率從1 - 32hz。

def fft(data, nwindow=64, freq = 32):

    ffts = []

    for i in range(0, len(data)-nwindow, nwindow//2):

        sliced = data[i:i+nwindow]

        fft = np.abs(np.fft.rfft(sliced*np.hamming(nwindow))[:freq])

        ffts.append(fft.tolist())

    ffts = np.array(ffts)

    return ffts

def data_plot(date_time, data, labels, ax):

    ax.plot(date_time, data)

    ax.set_xlim(date2num(np.min(date_time)), date2num(np.max(date_time)))

    axs_twinx = ax.twinx()

    axs_twinx.plot(date_time, labels, color='red')

    ax.set_ylabel('Label')

def fft_plot(ffts, ax):

    ax.imshow(np.flipud(np.rot90(ffts)), aspect='auto', cmap=matplotlib.cm.bwr, 

               norm=LogNorm(vmin=np.min(ffts), vmax=np.max(ffts)))

    ax.set_xlabel('Timestamp')

    ax.set_ylabel('Freq')

df_fourier = df_data.copy()

for name in df_boxplot.columns:

    if name not in ['datetime', 'date']:

        fig, axs = plt.subplots(2, 1, figsize=(15, 6))

        data = df_fourier[name].to_numpy()

        ffts = fft(data, nwindow=64, freq = 32)

        data_plot(df_fourier['datetime'], data, df_labels['machine_status'], axs[0])

        fft_plot(ffts, axs[1])

        axs[0].set_title(name)

        plt.show()

　　熵

　　可視化信息和熵是機(jī)器學(xué)習(xí)中的一個(gè)有用工具，因?yàn)樗鼈兪窃S多特征選擇、構(gòu)建決策樹(shù)和擬合分類(lèi)模型的基礎(chǔ)。

　　熵的計(jì)算如下：

　　歸一化頻率分布

　　最低熵是針對(duì)某一隨機(jī)變量計(jì)算的，該隨機(jī)變量的單個(gè)事件的概率為1.0，即確定性。一個(gè)隨機(jī)變量的最大熵是當(dāng)所有事件都是等可能的。

def entropy(data, nwindow=64, freq = 32):

    entropy_s = []

    for i in range(0, len(data)-nwindow, nwindow//2):

        sliced = data[i:i+nwindow]

        fft = np.abs(np.fft.rfft(sliced*np.hamming(nwindow))[:nwindow//2])

        p = fft / np.sum(fft)

        entropy = - np.sum(p * np.log(p))

        entropy_s.append(entropy)

    entropy_s = np.array(entropy_s)

    return entropy_s

def data_plot(date_time, data, labels, ax):

    ax.plot(date_time, data)

    axs_twinx = ax.twinx()

    axs_twinx.plot(date_time, labels, color='red')

    ax.set_xlabel('Value')

    ax.set_ylabel('Label')

def entropy_plot(data, ax):

    ax.plot(data, c='k')

    ax.set_xlabel('Timestamp')

    ax.set_ylabel('Entropy')

df_entropy = df_data.copy()

for name in df_boxplot.columns:

    if name not in ['datetime', 'date']:

        fig, axs = plt.subplots(2, 1, figsize=(15, 6))

        data = df_entropy[name].to_numpy()

        entropy_s = entropy(data, nwindow=64, freq = 32)

        data_plot(df_entropy['datetime'], data, df_labels['machine_status'], axs[0])

        entropy_plot(entropy_s, axs[1])

        axs[0].set_title(name)

        plt.show()

　　降維

　　當(dāng)有多個(gè)傳感器時(shí)，實(shí)現(xiàn)一種降維方法來(lái)獲得包含大部分信息的1、2或3個(gè)主要組件總是很有趣的。

　　對(duì)于這個(gè)例子，我實(shí)現(xiàn)了主成分分析(PCA)。這是計(jì)算主要組件并使用它們對(duì)數(shù)據(jù)進(jìn)行基礎(chǔ)更改的過(guò)程。

　　被解釋方差比率是每一個(gè)被選擇的組成部分的方差百分比。

　　對(duì)于第一個(gè)PCA組件，可以繪制數(shù)據(jù)，并直觀地檢查異常和時(shí)間序列之間是否存在關(guān)系。

pca = PCA(n_components=2)

principalComponents = pca.fit_transform(x)

principalDf = pd.DataFrame(data = principalComponents, columns = ['pc1', 'pc2'])

df_pca = df_data.copy()

df_pca['pca1'] = pd.Series(principalDf['pc1'].values, index=df.index)

df_pca['pca2'] = pd.Series(principalDf['pc2'].values, index=df.index)

print(df_pca.shape)

print(df_pca.head())

df_pca_hour = df_pca.groupby(pd.Grouper(key='datetime', axis=0, freq='H')).mean()

df_labels_hour = df_labels.groupby(pd.Grouper(key='datetime', axis=0, freq='H')).sum()

for name in df_pca.columns:

    if name in ['pca1', 'pca2']:

        fig, axs = plt.subplots(1, 1, figsize=(15, 2))

        axs.plot(df_pca_hour[name], color='blue')

        axs_twinx = axs.twinx()

        axs_twinx.plot(df_labels_hour['machine_status'], color='red')

        axs.set_title(name)

        plt.show()

　　自相關(guān)

　　最后，特別是對(duì)于預(yù)測(cè)任務(wù)，繪制數(shù)據(jù)的自相關(guān)性是很有趣的。這個(gè)表示給定的時(shí)間序列和它自己在連續(xù)時(shí)間間隔中的滯后版本之間的相似程度。

　　與自相關(guān)相關(guān)的是增強(qiáng)迪基-富勒統(tǒng)計(jì)檢驗(yàn)，用于檢驗(yàn)給定的時(shí)間序列是否平穩(wěn)。