本文在 Python 中用箱線圖、傅里葉變換、熵、自相關和 PCA 分析時間序列數據。數據可視化是任何數據相關項目中最重要的階段之一。根據數據可視化的對象有：

　　1.數據可視化報告結果。

　　2.數據可視化來分析數據，換句話說，數據科學家內部使用的可視化來提取有關數據的信息，然后實施模型。

　　本文主要關注后一種，因為它解釋了一些有助于分析時間序列數據的方法。

　　什么是時間序列？

　　基本數值時間序列是有序的、帶時間戳的觀測值(測量值)的集合，其中每個觀測值都是從同一測量過程中獲得的數值標量。

　　什么是時間戳？

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

　　Python類庫

　　首先，這些是與 notebook 一起使用的庫。大多數代碼都圍繞 NumPy 和 Pandas庫，因為數據主要以 Pandas Dataframe 表現的 NumPy 數組。

　　導入文件

　　下載數據后，運行以下代碼將其導入。

　　正如所觀察到的，數據包含六個傳感器的傳感器數據、每個數據點的日期時間以及機器狀態。這是“BROKEN”、“NORMAL”或“RECOVERING”，但為了簡化可視化，它被分組如下：

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

　　數據預處理

　　在進行可視化之前，分析了本次數據的重復值和缺失值。并且刪除重復項的函數：

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

    df = df.drop_duplicates((subset))

    return df

　　填充缺失值的函數：

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

　　這是預處理階段的整個管道。此外，數據分為輸入數據和輸出數據。

　　輸入形狀：(10081, 7)

　　輸出形狀：(10081, 2)

　　數據可視化

　　現在，準備開始數據可視化。這是傳感器數據和異常情況的圖。完整代碼可以在公眾號:機器學習研習院 后臺回復 時間序列可視化獲取.

　　均值和標準

　　可以更好地總結數據隨時間變化的行為的最基本圖之一是均值標準圖，我們在其中顯示按時間范圍分組的均值和標準差。這主要有助于分析指定時間范圍內的基線和噪聲。

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()

　　箱形圖

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

　　這一個類似于平均和標準圖，因為它表明數據的平穩性。但是，它也可以顯示異常值，這有助于從視覺上檢測異常和數據之間的任何關系。

　　傅里葉變換

　　快速傅里葉變換(FFT)是一種計算序列離散傅里葉變換的算法。這種類型的圖很有趣，因為它是處理時間序列時特征提取的主要方法之一。通常的做法不是用時間序列來訓練模型，而是應用傅里葉變換來提取頻率，然后訓練模型。

　　為此，我們必須選擇一個滑動窗口來計算FFT。滑動窗口越寬，頻率數越高。缺點是您將得到更少的時間戳，從而丟失數據的時間分辨率。當減小窗口的大小時，我們得到了相反的結果：更少的頻率但更高的時間分辨率。然后，窗口的大小應該取決于任務。

　　FFT的滑動窗口 對于如下圖所示，我選擇了一個包含64個數據的時間窗口。因此，頻率從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()

　　熵

　　可視化信息和熵是機器學習中的一個有用工具，因為它們是許多特征選擇、構建決策樹和擬合分類模型的基礎。

　　熵的計算如下：

　　歸一化頻率分布

　　最低熵是針對某一隨機變量計算的，該隨機變量的單個事件的概率為1.0，即確定性。一個隨機變量的最大熵是當所有事件都是等可能的。

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()

　　降維

　　當有多個傳感器時，實現一種降維方法來獲得包含大部分信息的1、2或3個主要組件總是很有趣的。

　　對于這個例子，我實現了主成分分析(PCA)。這是計算主要組件并使用它們對數據進行基礎更改的過程。

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

　　對于第一個PCA組件，可以繪制數據，并直觀地檢查異常和時間序列之間是否存在關系。

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()

　　自相關

　　最后，特別是對于預測任務，繪制數據的自相關性是很有趣的。這個表示給定的時間序列和它自己在連續時間間隔中的滯后版本之間的相似程度。

　　與自相關相關的是增強迪基-富勒統計檢驗，用于檢驗給定的時間序列是否平穩。