數據處理——拉伊達法則去除異常值（Python實現）

背景：

題目出自2020年中國研究生數學建模競賽B題

代碼及附件
上傳時間：2020.12.24

1 數據采集

原始數據采集來自於中石化高橋石化實時數據庫（霍尼韋爾PHD）及LIMS實驗數據庫。其中操作變量數據來自於實時數據庫，采集時間為2017年4月至2020年5月，采集操作位點數共354個。2017年4月至2019年9月，數據采集頻次為3分鍾/次；2019年10月至2020年5月，數據采集頻次為6分鍾/次。原料、產品和催化劑數據來自於LIMS實驗數據庫，數據時間范圍為2017年4月至2020年5月。其中原料及產品的辛烷值是重要的建模變量，該數據采集頻次為每周2次。

2 數據整定

原始數據中，大部分變量數據正常，但每套裝置的數據均有部分位點存在問題：部分變量只含有部分時間段的數據，部分變量的數據全部為空值或部分數據為空值。因此對原始數據進行處理后才可以使用。數據處理方法如下：

（1）對於只含有部分時間點的位點，如果其殘缺數據較多，無法補充，將此類位點刪除；

（2）刪除325個樣本中數據全部為空值的位點；

（3）對於部分數據為空值的位點，空值處用其前后兩個小時數據的平均值代替；

（4）根據工藝要求與操作經驗，總結出原始數據變量的操作范圍，然后采用最大最小的限幅方法剔除一部分不在此范圍的樣本；

（5）根據拉依達准則（3σ准則）去除異常值。

3σ准則：設對被測量變量進行等精度測量，得到x1，x2，……，xn，算出其算術平均值x及剩余誤差vi=xi-x（i=1,2,...,n），並按貝塞爾公式算出標准誤差σ，若某個測量值xb的剩余誤差vb（1<=b<=n），滿足|vb|=|xb-x|>3σ，則認為xb是含有粗大誤差值的壞值，應予剔除。貝塞爾公式如下：

3 樣本確定

本題目標為降低S Zorb裝置產品辛烷值損失，故確定樣本的主要依據為樣品的辛烷值數據。由於辛烷值的測定數據相對於操作變量數據而言相對較少，而且辛烷值的測定往往滯后，因此確定某個樣本的方法為：以辛烷值數據測定的時間點為基准時間，取其前2個小時的操作變量數據的平均值作為對應辛烷值的操作變量數據。

4 Python代碼

• 讀取 excel表格數據
• 去除空值、殘缺數據過多（大於10個）的位點
• 殘缺數據少的位點 數據用同位點平均值替代
• 根據 拉伊達准則 去除異常值

import numpy as np
import pandas as pd
import xlrd
import xlwt
import openpyxl
import math
import matplotlib.pyplot as plt
from datetime import date, datetime

# 用來正常顯示中文標簽
plt.rcParams['font.sans-serif'] = ['SimHei']
# 用來正常顯示負號
plt.rcParams['axes.unicode_minus'] = False


# 定義正態分布函數
def pdf(x, ave, std):
    y = (1.0 / math.sqrt(2 * math.pi * std)) * np.exp(-(x - ave) ** 2 / (2 * std ** 2))
    return y


# 通過路徑，打開文件工作簿
#313
# ExcelFile = xlrd.open_workbook(r'F:\0.0\競賽\2020數學建模\1111\附件三：285號和313號樣本原始數據 1.xlsx')
#285
ExcelFile = xlrd.open_workbook(r'F:\0.0\競賽\2020數學建模\1111\313.xlsx')


# 通過sheet名查找打開工作表
sheet = ExcelFile.sheet_by_name('操作變量')
# 通過sheet索引查找打開工作表
# sheet=ExcelFile.sheet_by_index(4)

# 打印輸出當前工作表名，行數以及列數
print('當前工作表名，行數以及列數' + sheet.name, sheet.nrows, sheet.ncols)

# 獲取行或者列中部分元素的值sheet.row_values(a,b,c) sheet.col_values(a,b,c)
# a,b,c均從0開始計數，從b號元素開始到c-1號元素，不包含c號元素

# 確定 各列輸出列表
outX1 = []  #均值輸出
list_1 = []  #空 列表
Z1 = 0  # 計算多少列剔除異常值


for a in range(1,355):

    # rowsData=sheet.row_values(3,1,355)   #取第4列，第2號元素到第355號元素

    # 285 號
    colsData = sheet.col_values(a, 3, 43)  #取第2列，第4號元素到第43號元素


    #313 號
    #colsData = sheet.col_values(a, 43, 84)  # 取第2列，第4號元素到第43號元素

    # 打印輸出取出的數據元素，默認數據類型為list
    # print("取出的數據元素:")
    # print(colsData)
    # 將list類型的數據轉換成array類型
    data = np.array(colsData)

    # 計算數據的算術平均值、標准誤差以及剩余誤差
    arithmeticMean = data.mean()
    standardDeviation = data.std()
    residualError = abs(data - arithmeticMean)

    # 打印輸出均值和標准差
    # print('列算數平均值、標准誤差以及剩余誤差：')
    # print(arithmeticMean, standardDeviation, residualError)

    # 未測位點判斷
    if standardDeviation == 0:
        X1 = 0

        print('第' + str(a) + '列數據全為0，算術平均值：' + str(X1))
        #print('注：本列全為0，缺損值過多')
        continue

    else :

        #解決缺損值問題
        Y1 = 0         #計算缺損數據量個數
        for i in range(0,40):
            if colsData[i] == 0:
                #求缺損值的個數
                Y1 = Y1 + 1
        #缺損值不多時,用平均值替代
        if 0 < Y1 <= 10:
            arithmeticMean = arithmeticMean * 40 / (40 - Y1)
            print('第' + str(a) + '列有少量缺損值，已替換')
        elif  Y1 > 10:
            arithmeticMean = arithmeticMean * 40 / (40 - Y1)
            print('第' + str(a) + '列缺損值過多，已刪除此為點')



        if __name__ == '__main__':

            # plot histgram of its distribution
            x = np.sort(data)
            y = pdf(x, arithmeticMean, standardDeviation)


            # 去除異常值
            good_x = []
            outliers = []


            for i, j in zip(residualError, x):
                if i < 3 * standardDeviation:
                    good_x.append(j)
                else:
                    outliers.append(j)

            good_x = np.array(good_x)
            goodArithmeticMean = good_x.mean()

            goodStandardDeviation = good_x.std()
            good_y = pdf(good_x, goodArithmeticMean, goodStandardDeviation)





            #求解 剔除異常值之后的算術平均值
            if outliers != list_1:
                X1 = (arithmeticMean * 40 - sum(outliers))/(40 - len(outliers))
                print('第' + str(a) + '列剔除的異常值為:', outliers)
                print('第' + str(a) + '列數據的移除異常值之后的算數平均值：' + str(X1))

                Z1 = Z1 + 1
            else:
                X1 = arithmeticMean
                print('第' + str(a) + '列數據的算術平均值：' + str(X1))

# 將處理之后的數據打印輸出

    outX1.append(X1)

    # # 數據可視化
    # plt.plot(x, y, c='b', label=u'原始值')
    # plt.plot(good_x, good_y, c='r', label=u'去除異常值后數據')
    # plt.title('Normalization distribution curve  313號第' + str(a) + '列元素')
    #
    # if outliers != list_1:
    #     plt.ylabel(u"標准差太小")
    #     plt.xlabel(u'原始數據  剔除的異常值為:'+ str(outliers))
    # else :
    #     plt.xlabel(u"原始數據")
    #
    # plt.legend()
    # plt.show()

    # 數據可視化


    if outliers != list_1:
        plt.ylabel(u"標准差太小")
        plt.xlabel(u'原始數據  剔除的異常值為:'+ str(outliers))
        plt.plot(x, y, c='b', label=u'原始值')
        plt.plot(good_x, good_y, c='r', label=u'去除異常值后數據')
        plt.title('Normalization distribution curve  313號第' + str(a) + '列元素')
        plt.legend()
        plt.show()

print('處理之后的數據：' )
print(outX1)
print('共計' + str(Z1) + '列剔除了 異常值')



# #將數據寫入新文件
# #打開工作簿，進行操作
# wb = openpyxl.Workbook()
# sheet = wb.active
#
# # sheet['A1'] = 'hello'
# # print(sheet['A1'].value)
#
# #新創建頁
# ws1 = wb.create_sheet('285')
# #ws2 = wb.create_sheet('313')
# row = 0
#
# ws1.append(outX1)
# #ws2.append(outX1)
# #313號樣本給入 313  #285號樣本給入 285
#
# wb.save(filename='date.xlsx')
# print('寫入成功')