技術標籤：數學建模Pythonpython資料分析資料探勘數學建模

資料處理——拉伊達法則去除異常值（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('寫入成功')