個人網站 Geek交流圈
視頻演示 B站視頻

Grubbs準則建模與分析

目錄

一、Grubbs準則定義及原理 1

1、 Grubbs定義 2
2、Grubbs原理 3
（1）具體流程： 3
（2）案例分析: 4

二、matlab程序設計及分析 7

1、matlab簡介與使用方法 7
（1）開發環境 8
2、設計思路及相關函數 8
（1）定義測量值與顯著度alpha 8
（2）畫圖函數 8
（3）計算 9
（4）確定g值 9
（5）比較異常值 10
（6）系統測試 11
3、完整源程序 13

三、C程序設計及分析 15

1、C開發環境及簡介 15
（1）開發環境 15
2、設計思路及相關函數 15
（1）求平均值函數 15
（2）求標準差函數 16
（3）排序函數 16
（4）輸出函數 16
（5）計算函數 17
（6）系統測試 18
3、完整源程序 19
四、實驗總結 23

一、Grubbs準則定義及原理

1、Grubbs定義
格拉布斯準則，屬于正態分布的分支，數學詞匯。 格拉布斯準則是以正態分布為前提的，理論上較嚴謹，使用也方便。 某個測量值的殘余誤差的絕對值 |Vi |>Gg，則判斷此值中有較大誤差，應以剔除，此即格拉布斯準則。
格拉布斯準則是以正態分布為前提的，理論上較嚴謹，使用也方便。
某個測量值的殘余誤差的絕對值 |Vi |>Gg，則判斷此值中有較大誤差，應以剔除，此即格拉布斯準則。
利用格拉布斯(Grubbs)準則進行處理：
根據誤差理論，要有效地剔除偶然誤差，一般要測量10次以上，兼顧到精度和響應速度，取15次為一個單位。
在取得的15個數據中，有些可能含有較大的誤差，需要對它們分檢，剔除可疑值，提高自適應速度。對可疑值的剔除有多種準則，如萊以達準則、肖維勒（Chauvenet）準則、格拉布斯（Grubbs）準則等。
以Grubbs準則為例，它認為若某測量值 xi對應的殘差Vi滿足下式 |Vi|=| xi-x|>=g(n，a)× σ(X) 時應將該數據舍去。式中，x為n次采集到的AD 值的平均值，=(∑xi)/n ；σ(X)為測量數據組的標準差。
由貝塞爾函數可得： σ(X)=[(∑Vi2 )/(n-1)]1/2；g(n， a)是取決于測量次數n和顯著性水平a (相當于犯“棄真” 錯誤的概率系數)，a通常取0.01或0.05。
通過查表可得：當 n=15時，a=0.05, g(n，a)=2.41。
把15次采集到的AD值存入一個數組中然后求平均值，計算殘差，求標準差σ(X)。
將殘差絕對值與2.41倍的標準差σ(X)比較。剔除可疑值以后，再求平均值，求出新的平均值以后，應再重復以上過程，驗證是否還有可疑值存在。據我們對測量裝置大量的實際測試結果看，這樣做沒有什么必要，因為一般只有第一遍即可達到要求。
然而這種方法也有它的不足, 利用Grubbs準則需要處理大量的數據，而在一般的工業現場測試設備中，儀表結構大多采用嵌入式結構，如AVR單片機。這些MCU程序空間和數據空間有限，若處理大量數據，難以滿足資源要求。而且，由于Grubbs準則要求MCU進行大量數據處理，使得系統降低了信號采集速率，影響實時性。
什么數據可能含有粗大誤差？ 在一段序列中，過高或過低的數，如下圖。

圖1-1 粗大誤差分析圖
2、Grubbs原理
（1）具體流程：
1、定義輸入數據 data[n]
確定數據個數 n
2、選取顯著度 a 一般是0.01或者0.05

3、通過查閱格拉布斯系數表，可以確定格拉布斯系數 ，由n和a確定。
4、將數據進行排序，由小到大，由大到小都可以。
5、求數學期望，也就是平均值
6、求標準差

若認為x(n)可疑，則有

如果 ，則該測量值含有粗大誤差，應該剔除。
具體怎么判斷呢？
因為含有粗大誤差的數是過高或過低的數，那么我們就需要判斷最大值與最小值。假如最大值含有粗大誤差，那么需要剔除最大值，然后重新計算平均值與標準差，再求出 ，繼續判斷，直到最大值沒有含粗大誤差為止。需注意最大值剔除后，原來的第二大值就成為的最大值，所以我們一直是判斷最大值。
（2）案例分析
1、定義n個數據 1,2 ,3 ,4 ,9, 8, 7, 6, 17 ,5,10
這里n = 11；
2、設置顯著度a = 0.01；（自己設定，一般為0.01或0.05）
需注意，在不同的a下，同一數據可能得到的結果不同，案例分析。
通過查表可知：
當n = 11，a = 0.01時， = 2.48；
當n = 11，a = 0.05時， = 2.23；
3、從小到大排序為： 1 2 3 4 5 6 7 8 9 10 17
平均值 為：6.546
標準差 為：4.293

4、計算g(max) 與g(min) ，計算公式如下

判斷最小值： =1.292 < = 2.48 不含粗大誤差
判斷最大值： =2.435 < = 2.48 不含粗大誤差
說明在顯著度a = 0.01時，上面的數據沒有含有粗大誤差。

當我們取g = 0.05時， = 2.23
平均值 為：6.546
標準差 為：4.293
判斷最小值： =1.292 < = 2.23 不含粗大誤差
判斷最大值： =2.435 > = 2.23 含有粗大誤差
所以將最大值剔除，也就是將17剔除，那么最大值變為10
重新求平均值與標準差
平均值 為：5.50
標準差 為：2.87
再求此時的 與
判斷最小值： =1.567 < = 2.23 不含粗大誤差
判斷最大值： =1.567 < = 2.23 不含粗大誤差
判斷結束，消除粗大誤差的數據為
1 2 3 4 9 8 7 6 5 10 （未排序）
1 2 3 4 5 6 7 8 9 10 （排序后）

由上例分析，同樣的數據，當顯著度a=0.01時，數據未含有粗大誤差；當顯著度a=0.05時，數據含有粗大誤差。這說明數據是否含有粗大誤差是相對的，與自己的初始條件有關。
不同的分析方法得到的結果可能不同。
下面是Grubbs準則的流程圖，需要了解邏輯，在程序設計的過程中，就是根據該邏輯進行設計。

圖1-2 Grubbs流程圖

圖1-3 格拉布斯系數表

二、matlab程序設計及分析

1、matlab簡介與使用方法
（1）開發環境
MATLAB是美國MathWorks公司出品的商業數學軟件，用于數據分析、無線通信、深度學習、圖像處理與計算機視覺、信號處理、量化金融與風險管理、機器人，控制系統等領域。
MATLAB是matrix&laboratory兩個詞的組合，意為矩陣工廠（矩陣實驗室），軟件主要面對科學計算、可視化以及交互式程序設計的高科技計算環境。它將數值分析、矩陣計算、科學數據可視化以及非線性動態系統的建模和仿真等諸多強大功能集成在一個易于使用的視窗環境中，為科學研究、工程設計以及必須進行有效數值計算的眾多科學領域提供了一種全面的解決方案，并在很大程度上擺脫了傳統非交互式程序設計語言（如C、Fortran）的編輯模式。
本例程使用MATLAB R2016a進行實驗。
2、設計思路及相關函數
（1）定義測量值與顯著度alpha

%定義序列 2. a=[ 1,2 ,3 ,4 ,9, 8, 7, 6, 17 ,5,10]; 3. %輸入顯著度alpha 4. alpha=input('(請輸入alpha的值0.01或0.05)'); 程序說明： a 就是一維數組，數據個數為11，alpha為顯著度，讓我們自己輸入，一般為0.01或0.05。

（2）畫圖函數

1. subplot(121);stem(a,'filled'); %subplot函數是在平鋪位置創建坐標區 2. xlabel('序列個數n'); %xlabel函數 為 x 軸添加標簽 3. ylabel('幅值'); %ylabel函數 為 y 軸添加標簽 4. title('原序列'); %添加標題 5. %axis([0 25 0 30]); 6. axis([0 15 -20 20]); %設置坐標軸范圍和縱橫比

程序說明：
subplot(mnp)參數說明
m表示是圖排成m行
n表示圖排成n列
p是指你現在要把曲線畫到figure中哪個圖上，最后一個如果是1表示是從左到右第一個位置
stem函數是繪制離散序列數據，以實心的方式畫出莖稈。
如圖2-1所示，對比程序中的文字說明，查看各函數的對應關系。

圖2-1 畫圖程序顯示
（3）計算

1. n=length(a); %得到序列a的長度 2. avg=mean(a); %求數組的均值 3. st=std(a,1); %求標準差 4. a_max = max(a); % a數組中的最大值 5. a_min = min(a); % a數組中的最小值 6. reg_max=abs((max(a)-avg)/st); %求最大值的Grubbs系數 7. reg_min=abs((min(a)-avg)/st); %求最小值的Grubbs系數

程序說明：
在matlab中，很多數學公式都已經封裝成為了函數，只需要了解函數的運用就行，具體函數內容看程序中的注釋，每一步都有說明。
（4）確定 g 值

1. %確定g值 2. T_01=[0 0 1.155 1.492 1.749 1.944 2.097 2.22 2.323 2.41 2.485 2.55 2.607]; 3. T_05=[0 0 1.153 1.463 1.672 1.822 1.938 2.032 2.11 2.176 2.234 2.285 2.331]; 4. switch (alpha) %由顯著度 alpha 確定T的取值， 5. case 0.05 6. g=T_05(n); 7. case 0.01 8. g=T_01(n); 9. otherwise %disp('輸入了錯誤的alpha值'); 程序說明：在我們輸入顯著度a過后，在這段程序運行就能得到 。

（5）比較異常值

%比較異常值 2. if reg_max>g %判斷最大值 3. flag_2=1; 4. abn=max(a); 5. disp('被剔除的數據為') 6. disp(abn); 7. a(a==max(a))=[]; 8. else if reg_min>g %判斷最小值 9. flag_2=1; 10. abn=min(a); 11. disp('被剔除的數據為') 12. disp(abn); 13. a(a==min(a))=[]; 14. else %沒有含有粗大誤差的數啦 15. flag_1=0; %標志清為0，退出while循環 16. %畫出處理好的序列圖像 17. subplot(122);stem(a,'filled'); 18. title('處理后的序列'); 19. xlabel('序列個數n'); 20. ylabel('幅值'); 21. axis([0 15 -20 20]); 22. end 23. end

程序說明：
先判斷最大值是否含有粗大誤差，再判斷最小值是否含有粗大誤差，當都沒有粗大誤差的時候，輸出數據。
（6）系統測試
圖2-2 matlab操作頁面

圖2-3 程序運行指令

圖2-4 a = 0.05數據顯示圖

圖2-5 a = 0.02數據顯示圖
由程序運行可知，我們需要輸入顯著度alpha，0.01或者0.05。當我們輸入0.05時，經程序分析，數據17與-15含有粗大誤差，剔除數據17與-15過后，再進行計算，當沒有粗大誤差時，最后顯示出處理后的序列。
3、完整源程序

1. %定義序列 2. a=[ -15,2 ,3 ,4 ,9, 8, 7, 6, 17 ,5,10]; 3. 4. %輸入顯著度alpha 5. alpha=input('(請輸入alpha的值0.01或0.05)'); 6. 7. %關于畫圖程序========================================================================== 8. subplot(121);stem(a,'filled'); %subplot函數是在平鋪位置創建坐標區 stem函數是繪制離散序列數據，以實心的方式畫出莖稈。 9. %subplot(mnp)參數說明 10. %m表示是圖排成m行，n表示圖排成n列 11. %p是指你現在要把曲線畫到figure中哪個圖上，最后一個如果是1表示是從左到右第一個位置 12. 13. xlabel('序列個數n'); %xlabel函數 為 x 軸添加標簽 14. ylabel('幅值'); %ylabel函數 為 y 軸添加標簽 15. title('原序列'); %添加標題 16. %axis([0 25 0 30]); 17. axis([0 15 -20 20]); %設置坐標軸范圍和縱橫比 18. 19. flag_1=1; %設置標志變量 20. flag_2=0; 21. 22. %計算 23. while flag_1==1 %while開始 24. n=length(a); %得到序列a的長度 25. avg=mean(a); %求數組的均值 26. st=std(a,1); %求標準差 27. a_max = max(a); % a數組中的最大值 28. a_min = min(a); % a數組中的最小值 29. reg_max=abs((max(a)-avg)/st); %求最大值的Grubbs系數 30. reg_min=abs((min(a)-avg)/st); %求最小值的Grubbs系數 31. 32. %確定g值 33. T_01=[0 0 1.155 1.492 1.749 1.944 2.097 2.22 2.323 2.41 2.485 2.55 2.607]; 34. T_05=[0 0 1.153 1.463 1.672 1.822 1.938 2.032 2.11 2.176 2.234 2.285 2.331]; 35. switch (alpha) %由顯著度 alpha 確定T的取值， 36. case 0.05 37. g=T_05(n); 38. case 0.01 39. g=T_01(n); 40. otherwise %disp('輸入了錯誤的alpha值'); 41. end %while循環結束 42. 43. %比較異常值 44. if reg_max>g %判斷最大值 45. flag_2=1; 46. abn=max(a); 47. disp('被剔除的數據為') 48. disp(abn); 49. a(a==max(a))=[]; 50. else if reg_min>g %判斷最小值 51. flag_2=1; 52. abn=min(a); 53. disp('被剔除的數據為') 54. disp(abn); 55. a(a==min(a))=[]; 56. else %沒有含有粗大誤差的數啦 57. flag_1=0; %標志清為0，退出while循環 58. %畫出處理好的序列 59. subplot(122);stem(a,'filled'); 60. title('處理后的序列'); 61. xlabel('序列個數n'); 62. ylabel('幅值'); 63. axis([0 15 -20 20]); 64. end 65. end 66. 67. if flag_2==0 68. disp('沒有異常數據'); 69. end; 70. 71. end;

三、C程序設計及分析

1、C開發環境及簡介
Microsoft Visual Studio（簡稱VS）是美國微軟公司的開發工具包系列產品。VS是一個基本完整的開發工具集，它包括了整個軟件生命周期中所需要的大部分工具，如UML工具、代碼管控工具、集成開發環境(IDE)等等。所寫的目標代碼適用于微軟支持的所有平臺，包括Microsoft Windows、Windows Mobile、Windows CE、.NET Framework、.NET Compact Framework和Microsoft Silverlight 及Windows Phone。
本例程使用Visual Studio 2019進行實驗。
2、設計思路及相關函數
因為標準C中并沒有那么豐富的數學函數，很多需要自己造，所以在本次實驗過程中平均值，標準差，排序這些都是自己寫的函數。雖然自己寫繁瑣了一點，但個人覺得能更好的掌握原理，對自己有很大幫助。

（1）求平均值函數

1. //求數學期望 2. float MathExp(float* data, int start, int end) 3. { 4. float ex = 0; 5. int count = 0; 6. for (int i = start; i <= end; i++) 7. { 8. ex += data[i]; 9. count++; 10. } 11. ex /= count; 12. return ex; 13. }

（2）求標準差函數

1. //求標準差 2. float StdDev(float* data, int start, int end) 3. { 4. float ex = MathExp(data, start, end); 5. float sd = 0; 6. int count = 0; 7. for (int i = start; i <= end; i++) 8. { 9. sd += pow(data[i] - ex, 2); 10. count++; 11. } 12. sd = sqrt(sd / count); 13. return sd; 14. }

（3）排序函數

1. //排序 2. void Sort(float* data, int n) 3. { 4. for (int i = 1; i < n; i++) 5. { 6. for (int j = 0; j < n - i; j++) 7. { 8. if (data[j] > data[j + 1]) 9. { 10. float temp = data[j]; 11. data[j] = data[j + 1]; 12. data[j + 1] = temp; 13. } 14. } 15. } 16. }

（4）輸出函數

1. //輸出數組 2. void Print(float* data, int n) 3. { 4. for (int i = 0; i < n; i++) 5. { 6. printf("%0.3lf ", data[i]); 7. } 8. printf("\n"); 9. }

（5）計算函數

1. float ex = MathExp(data, start, end); //求數學期望 2. float sd = StdDev(data, start, end); //求標準差 3. 4. printf("數學期望為：%0.3lf\n 標準差為：%0.3lf\nT: %0.3lf\n", ex, sd, T[n - 1]); 5. float g0_min = (fabs(data[start] - ex)) / sd; 6. if ((fabs(data[start] - ex)) / sd >= T[n - 1]) //fabs是求絕對值 7. { 8. printf("%lf 為粗大誤差，應剔除\n", data[start]); 9. n--; 10. start++; 11. flag = 1; 12. } 13. 14. float g0_max = (fabs(data[end] - ex)) / sd; 15. if ((fabs(data[end] - ex)) / sd >= T[n-1]) 16. { 17. printf("%0.3lf 為粗大誤差，應剔除\n", data[end]); 18. n--; 19. end--; 20. flag = 1; 21. } 22. if (!flag) 23. { 24. printf("余下數據無壞值\n余下數據為：\n"); 25. Print_Part(data, start, end); 26. printf("\n"); 27. }

（6）系統測試

圖3-1 VS2019操作頁面

圖3-2 a = 0.05運行顯示

圖3-1 a = 0.01運行顯示
3、完整源程序

1. #define _CRT_SECURE_NO_WARNINGS 2. #include <stdio.h> 3. #include <stdlib.h> 4. #include <string.h> 5. #include <math.h> 6. 7. //求數學期望 8. float MathExp(float* data, int start, int end) 9. { 10. float ex = 0; 11. int count = 0; 12. for (int i = start; i <= end; i++) 13. { 14. ex += data[i]; 15. count++; 16. } 17. ex /= count; 18. return ex; 19. } 20. 21. //求標準差 22. float StdDev(float* data, int start, int end) 23. { 24. float ex = MathExp(data, start, end); 25. float sd = 0; 26. int count = 0; 27. for (int i = start; i <= end; i++) 28. { 29. sd += pow(data[i] - ex, 2); 30. count++; 31. } 32. sd = sqrt(sd / count); 33. return sd; 34. } 35. 36. 37. //排序 38. void Sort(float* data, int n) 39. { 40. for (int i = 1; i < n; i++) 41. { 42. for (int j = 0; j < n - i; j++) 43. { 44. if (data[j] > data[j + 1]) 45. { 46. float temp = data[j]; 47. data[j] = data[j + 1]; 48. data[j + 1] = temp; 49. } 50. } 51. } 52. } 53. 54. //輸出數組 55. void Print(float* data, int n) 56. { 57. for (int i = 0; i < n; i++) 58. { 59. printf("%0.3lf ", data[i]); 60. } 61. printf("\n"); 62. } 63. 64. //輸出部分數組 65. void Print_Part(float* data, int start, int end) 66. { 67. for (int i = start; i <= end; i++) 68. { 69. printf("%0.3lf ", data[i]); 70. } 71. printf("\n"); 72. 73. } 74. 75. int main() 76. { 77. /* 78. 本次測試的數據個數為11 顯著度設為0.01 79. 如果要測試其他數據，直接修改即可 80. 26.773, 50, 26.774, 26.778, 26.771, 26.780, 26.772, 26.777, 26.775, 26.774, 26.774 81. */ 82. float data[11] = { 1,2 ,3 ,4 ,9, 8, 7, 6, 17 ,5,10 }; 83. float T_05[11] = { 0, 0, 1.15, 1.46, 1.67, 1.82, 1.94, 2.03, 2.11, 2.18, 2.23 }; 84. float T_01[11] = { 0, 0, 1.15, 1.49, 1.75, 1.94, 2.10, 2.22, 2.32, 2.41, 2.48 }; 85. float T[11] = { 0 }; 86. char a_flag[5] = { 0 }; 87. 88. printf("請輸入顯著度a（0.01或0.05）："); 89. scanf("%s", &a_flag); 90. if (!strcmp(a_flag,"0.05")) //相同返回0 91. { 92. for (int i = 0; i < 11; i++) 93. { 94. T[i] = T_05[i]; 95. } 96. } 97. else if (!strcmp(a_flag, "0.01")) //相同返回0 98. { 99. for (int i = 0; i < 11; i++) 100. { 101. T[i] = T_01[i]; 102. } 103. } 104. else 105. { 106. printf("輸入參數錯誤，默認為a = 0.05\n"); 107. for (int i = 0; i < 11; i++) 108. { 109. T[i] = T_05[i]; 110. } 111. } 112. 113. int n = 11; 114. printf("本次處理的測量數據為: \n"); 115. Print(data, n); 116. printf("\n"); 117. 118. printf("從小到大排序: \n"); 119. Sort(data, n); 120. Print(data, n); 121. printf("\n"); 122. 123. int start = 0; 124. int end = n - 1; 125. char flag = 1; 126. while (flag) 127. { 128. flag = 0; 129. float ex = MathExp(data, start, end); //求數學期望 130. float sd = StdDev(data, start, end); //求標準差 131. 132. printf("數學期望為：%0.3lf\n 標準差為：%0.3lf\nT: %0.3lf\n", ex, sd, T[n - 1]); 133. float g0_min = (fabs(data[start] - ex)) / sd; 134. if ((fabs(data[start] - ex)) / sd >= T[n - 1]) //fabs是求絕對值 135. { 136. printf("%lf 為粗大誤差，應剔除\n", data[start]); 137. n--; 138. start++; 139. flag = 1; 140. } 141. 142. float g0_max = (fabs(data[end] - ex)) / sd; 143. if ((fabs(data[end] - ex)) / sd >= T[n-1]) 144. { 145. printf("%0.3lf 為粗大誤差，應剔除\n", data[end]); 146. n--; 147. end--; 148. flag = 1; 149. } 150. if (!flag) 151. { 152. printf("余下數據無壞值\n余下數據為：\n"); 153. Print_Part(data, start, end); 154. printf("\n"); 155. } 156. } 157. return 0; 158. }

四、實驗總結

• 通過建模分析，不僅讓我對Grubbs準則有了更深刻的認識，而且對matlab、C編程更加熟練。其實Grubbs原理并不是很難，但是建模與實際應用加大了難度，在編寫這篇論文的時候，反復去查閱資料與實際測試，因為要保證程序的健全性以及數據的準確性，總不能得到錯誤的運行結果吧。在這里我不得不吐槽一下，網上的一些資料真是誤人子弟，很多都是錯誤的，模棱兩可的也很多。

• 在實際編程過程中，用C的時間比較多，因為要自己寫些函數。后來想用C++實現一下的，因為有STL比較方便，后來覺得沒啥意思，覺得C和matlab是兩個最好的代表了，一個偏底層，一個偏應用。

• 最開始是想隨便實現下的，但是覺得既然做了，那就要做到最好，所以寫論文，刨析原理，不同的方法去建模分析，不僅要讓自己懂，也要讓別人懂。喜歡將復雜的知識用最簡易的話描述出來，而不是一堆學術用詞，讓不懂的人在里面反復打滾。

• 平時都對課程沒啥興趣，一般都是做幾個題，寫在本子上，念念理論，實在無聊，又沒啥進步。相反我更喜歡自由發揮的，做項目，寫論文，從寫方案，到設計與實現，每一步都能讓自己收獲巨大，這些都是課本上學不到的。

總結

以上是生活随笔為你收集整理的Grubbs准则建模与分析 C与Matlab实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。