パースの基準

実数値の測定値を含むデータセットにおいて、疑わしい外れ値とは、他のほとんどのデータ値のクラスターから外れているように見える測定値のことです。位置の要約統計量として算術平均を用いる場合、これらの外れ値によって位置の推定値が大きく変化します。問題は、算術平均は外れ値が含まれると非常に敏感であるということです。統計用語で言えば、算術平均は頑健ではありません。

外れ値が存在する場合、統計学者には2つの選択肢があります。1つ目は、データセットから疑わしい外れ値を除外し、算術平均を用いて位置パラメータを推定する方法です。2つ目は、中央値などのロバストな統計量を用いる方法です。

ピアースの基準は、外れ値を除去するための統計的手法である。

統計学者であり統計史家でもあるスティーブン・M・スティグラーは、ベンジャミン・パースについて次のように書いています。^{[ 1 ]}

「1852年、彼は外れ値を棄却すべきかどうかを研究者に知らせるための最初の有意性検定を発表した（Peirce 1852, 1878）。尤度比に基づくこの検定は、そのような行動の妥当性について国際的な議論を巻き起こしたという点で特筆すべきものである（Anscombe 1960、Rider 1933、Stigler 1973a）。」

ピアースの基準は、ガウス分布の統計分析から導き出されたものです。外れ値を除去するための他の基準とは異なり、ピアースの方法は2つ以上の外れ値を特定するために適用できます。

「一連の決定において、${\displaystyle m}$観測誤差の限界、これを超えると、これほど大きな誤差を含むすべての観測は棄却される可能性があるが、観測数が十分であれば、${\displaystyle n}$このような観察。この問題を解決するために提案されている原則は、提案された観察を保持することによって得られるエラーシステムの確率が、それらを拒否することによって得られるエラーシステムの確率に、その数以上の異常な観察を行う確率を掛けたものよりも小さい場合に、提案された観察を拒否すべきであるというものである。^{[ 2 ]}

ホーキンス^{[ 3 ]}は基準の公式を示している。

ピアースの基準は、 1878年に米国沿岸測地測量局と改名された米国沿岸測量局で数十年にわたって使用されました^{[ 4 ]}。

「1852年から1867年まで、彼は米国沿岸測量局の経度測定責任者を務め、1867年から1874年まで同測量局の監督官を務めた。これらの期間中、彼のテストは、当時最も活動的で数学的傾向の強い統計組織であった同局のすべての事務員によって一貫して使用された。」^{[ 1 ]}

パースの基準はウィリアム・ショーヴネの著書で論じられている。^{[ 2 ]}

1855年、BA グールドは、パースの方程式の値を表す値の表を作成することで、パースの基準をより簡単に適用できるようにしようと試みた。^{[ 5 ]}グールドのアルゴリズムとパースの基準の実際的な適用との間には、依然として乖離が存在する。

2003年、SM Ross（ニューヘイブン大学）は、新しいサンプルデータセットとアルゴリズムの解説を用いて、グールドのアルゴリズム（現在は「パース法」と呼ばれている）を再提示した。この手法は依然としてルックアップテーブルの使用に依存しており、この研究ではルックアップテーブルが更新されている（パースの基準テーブル）。^{[ 6 ]}

2008年にデンマークの地質学者K.トムセンが擬似コードの作成を試みた。^{[ 7 ]} このコードはグールドのアルゴリズムの枠組みを提供したが、ユーザーはピアースまたはグールドが報告した値を計算することに成功しなかった。

2012年、C. Dardisは、外れ値除去の比較を含むさまざまな手法（Peirceの基準とChauvenet法）を備えたRパッケージ「Peirce」をリリースしました。Dardisと共同貢献者のSimon Mullerは、Thomsenの擬似コードを「findx」という関数に実装することに成功しました。そのコードは、以下のR実装のセクションで紹介します。Rパッケージの参考文献はオンラインで入手可能です^{[ 8 ]}。また、Rパッケージの結果に関する未発表のレビューも公開されています^{[ 9 ] 。}

2013年、グールドのアルゴリズムの再検討と高度なPythonプログラミングモジュール（numpyやscipyなど）の利用により、外れ値を特定するための二乗誤差の閾値を計算することが可能になった。

ピアースの基準を使用するには、まず入力値と戻り値を理解する必要があります。 回帰分析（またはデータへの曲線のフィッティング）では、残差誤差（またはフィッティングされた曲線と観測点の差）が発生します。したがって、各観測点には、フィッティングされた曲線に関連付けられた残差誤差があります。残差誤差を二乗（つまり、残差誤差を2乗）すると、残差誤差は正の値として表されます。二乗誤差が大きすぎる場合（つまり、観測が不十分な場合）、曲線フィッティングから得られる回帰パラメータ（たとえば、線形曲線の傾きと切片）に問題が発生する可能性があります。

パースの発案は、何が「大きすぎる」誤差を構成し、したがって「外れ値」として識別されるかを統計的に特定し、観測値と曲線との適合性を向上させるために観測値から除外するというものでした。K. トムセンは、計算を実行するには 3 つのパラメータが必要であることを特定しました。観測値のペアの数 (N)、除外する外れ値の数 (n)、および残差を取得するために曲線フィッティングで使用される回帰パラメータ (係数など) の数 (m) です。このプロセスの最終結果は、閾値 (二乗誤差) を計算することです。この閾値よりも二乗誤差が小さい観測値は保持され、この値よりも二乗誤差が大きい観測値は (外れ値として) 除外されます。

ピアースの基準は、観測値、適合パラメータ、残差誤差を入力として受け取らないため、出力をデータと再関連付ける必要があります。すべての二乗誤差の平均（すなわち、平均二乗誤差）を求め、それを閾値二乗誤差（すなわち、この関数の出力）で乗算すると、外れ値を識別するために使用されるデータ固有の閾値が得られます。

次の Python コードは、Gould 1855 の表 1 (m = 1) と表 2 (m = 2) の指定されたN (最初の列) とn (最上行) に対する x の二乗値を返します。 ^{[ 5 ]}ニュートン法による反復により、N 対 log Q (Gould、1855 の表 III) や x 対 log R (Peirce、1852 の表 III および Gould、1855 の表 IV) などのルックアップ テーブルは不要になりました。

#!/usr/bin/env python3 import numpy import scipy.specialdef peirce_dev ( N : int , n : int , m : int ) -> float : """Peirce の基準 Gould の方法論に基づく Peirce の基準を使用して 外れ値を識別するための二乗閾値誤差偏差を返します。 引数:  - int、観測の総数 (N)  - int、削除する外れ値の数 (n)  - int、モデルの未知数の数 (m) 戻り値:  float、二乗誤差閾値 (x2)  """ # float を入力変数に割り当てます: N = float ( N ) n = float ( n ) m = float ( m )# 観測数をチェック: if N > 1 : # Q (グールド方程式 B の N 乗根) を計算: Q = ( n ** ( n / N ) * ( N - n ) ** (( N - n ) / N )) / N # # R 値を初期化 (浮動小数点数として) r_new = 1.0 r_old = 0.0 # <- while ループを促すために必要# # R に収束するための反復を開始: while abs ( r_new - r_old ) > ( N * 2.0e-16 ): # ラムダを計算# (グールド方程式 A' の 1/(Nn) 乗根): ldiv = r_new ** n if ldiv == 0 : ldiv = 1.0e-6 Lamda = (( Q ** N ) / ( ldiv )) ** ( 1.0 / ( N - n )) # x の二乗を計算(グールドの式 C): x2 = 1.0 + ( N - m - n ) / n * ( 1.0 - Lamda ** 2.0 ) # x2 が負になった場合は 0 を返す: if x2 < 0 : x2 = 0.0 r_old = r_new else : # x の二乗を使用して R を更新する (グールドの式 D): r_old = r_new r_new = numpy . exp (( x2 - 1 ) / 2.0 ) * scipy . special . erfc ( numpy . sqrt ( x2 ) / numpy . sqrt ( 2.0 ) ) else : x2 = 0.0 return x2

import org.apache.commons.math3.special.Erf ;public class PierceCriterion {/**  * ピアースの基準 * <p> * グールドの方法論に基づくピアースの基準を使用して 外れ値を識別するための二乗閾値誤差偏差を返します。  * <p>  * 引数:  * - int、観測の総数 (N)  * - int、削除する外れ値の数 (n)  * - int、モデルの未知数の数 (m)  * 戻り値:  * float、二乗誤差閾値 (x2)  **/ public static final double peirce_dev ( double N , double n , double m ) { // 観測数をチェック: double x2 = 0.0 ; if ( N > 1 ) { // Q (グールドの方程式 B の N 乗根) を計算: double Q = ( Math . pow ( n , ( n / N )) * Math . pow (( N - n ), (( N - n ) / N ))) / N ;// Rの値を初期化します（浮動小数点数として）double r_new = 1.0 ; double r_old = 0.0 ; // <- whileループを実行するために必要// R に収束するための反復を開始します: while ( Math . abs ( r_new - r_old ) > ( N * 2.0e-16 )) { // ラムダを計算します// (グールド方程式 A' の 1 / (N - n) 乗根): double ldiv = Math . pow ( r_new , n ); if ( ldiv == 0 ) { ldiv = 1.0e-6 ; } double Lamda = Math . pow (( Math . pow ( Q , N ) / ( ldiv )), ( 1.0 / ( N - n ))); // x の二乗 (グールド方程式 C): x2 = 1.0 + ( N - m - n ) / n * ( 1.0 - Math . pow ( Lamda , 2.0 )); // x2 が負の値になった場合は 0 を返す: if ( x2 < 0 ) { x2 = 0.0 ; r_old = r_new ; } else { // x の二乗を使用して R(グールドの方程式 D) を更新する: r_old = r_new ; r_new = Math . exp (( x2 - 1 ) / 2.0 ) * Erf . erfc ( Math . sqrt ( x2 ) / Math . sqrt ( 2.0 )); } } } else { x2 = 0.0 ; } return x2 ; } }

トムセンのコードは、2012年にC. DardisとS. Mullerによって次の関数呼び出し「findx」に正常に書き込まれ、最大誤差偏差を返します。${\displaystyle x}$前のセクションで紹介した Python コードを補完するために、二乗最大誤差偏差を返す R 版の「peirce_dev」もここに示します。${\displaystyle x^{2}}$これら 2 つの関数は、"findx" 関数から返された値を二乗するか、"peirce_dev" 関数から返された値の平方根を取るかのいずれかによって、同等の値を返します。違いはエラー処理にあります。たとえば、"findx" 関数は無効なデータに対して NaN を返しますが、"peirce_dev" は 0 を返します (これにより、追加の NA 値処理なしで計算を続行できます)。また、"findx" 関数は、潜在的な外れ値の数が観測数に近づくにつれてエラー処理をサポートしていません (欠損値エラーと NaN 警告をスローします)。

Python バージョンと同様に、二乗誤差 (つまり、${\displaystyle x^{2}}$）「peirce_dev」関数によって返される値は、モデル適合の平均二乗誤差で乗算して二乗デルタ値（つまり、Δ2）を取得する必要があります。Δ2を使用して、モデル適合の二乗誤差値を比較します。二乗誤差がΔ2より大きい観測ペアは外れ値とみなされ、モデルから削除できます。外れ値の数（Δ2とモデル適合の二乗誤差を比較）が想定値（つまり、Peirceのn）より少なくなるまで、nの値を増やしてテストするイテレータを作成する必要があります。

findx <- function ( N , k , m ) { # K. Thomsen (2008) によるメソッド# C. Dardis と S. Muller (2012) によって作成# オンラインで入手可能: https://r-forge.r-project.org/R/?group_id=1473 # # 変数の定義: # N :: 観測値の数# k :: 除去する可能性のある外れ値の数# m :: 未知の量の数# # 相補誤差関数 erfc が必要です: erfc <- function ( x ) 2 * pnorm ( x * sqrt ( 2 ), lower.tail = FALSE ) # x <- 1 if (( N - m - k ) <= 0 ) { return ( NaN ) print ( NaN ) } else { x <- min ( x , sqrt (( N - m ) / k ) - 1e-10 ) # # Log ofグールドの方程式 B: LnQN <- k * log ( k ) + ( N - k ) * log ( N - k ) - N * log ( N ) # # グールドの方程式 D: R1 <- exp (( x ^ 2 - 1 ) / 2 ) * erfc ( x / sqrt ( 2 )) # # ラムダ置換を用いて R について解いたグールドの方程式 A': R2 <- exp ( ( LnQN - 0.5 * ( N - k ) * log (( N - m - k * x ^ 2 ) / ( N - m - k ))) ) / k ) # # 2 つの R 方程式を等式化します: R1d <- x * R1 - sqrt ( 2 / pi / exp ( 1 )) R2d <- x * ( N - k ) / ( N - m - k * x ^ 2 ) * R2 # # x を更新します: oldx <- x x <- oldx - ( R1 - R2 ) / ( R1d - R2d ) # # 収束するまでループします: while ( abs ( x - oldx ) >= N * 2e-16 ) { R1 <- exp (( x ^ 2 - 1 ) / 2 ) * erfc ( x / sqrt ( 2 )) R2 <- exp ( ( LnQN - 0.5 * ( N - k ) * log (( N - m - k * x ^ 2 ) / ( N - m - k )) ) / k ) R1d <- x * R1 - sqrt ( 2 / pi / exp ( 1 )) R2d <- x * ( N - k ) / ( N - m - k * x ^ 2 ) * R2 oldx <- x x <- oldx - ( R1 - R2 ) / ( R1d - R2d ) } }return ( x ) }

peirce_dev <- function ( N , n , m ) { # N :: 観測値の総数# n :: 除去する外れ値の数# m :: モデルの未知数 (回帰パラメータなど) の数# # 観測値の数を確認します: if ( N > 1 ) { # Q (グールド方程式 B の N 乗根) を計算します: Q = ( n ^ ( n / N ) * ( N - n ) ^ (( N - n ) / N )) / N # # R 値を初期化します: Rnew = 1.0 Rold = 0.0 # <- while ループを促すために必要# while ( abs ( Rnew - Rold ) > ( N * 2.0e-16 )) { # ラムダ (グールド方程式 A' の 1/(Nn) 乗根) を計算します: ldiv = Rnew ^ n if ( ldiv == 0 ) { ldiv = 1.0e-6 } Lamda = (( Q ^ N ) / ( ldiv )) ^ ( 1.0 / ( N - n )) # # x の二乗 (グールドの方程式 C) を計算します: x2 = 1.0 + ( N - m - n ) / n * ( 1.0 - Lamda ^ 2.0 ) # # x2 が負になった場合は、ゼロに設定します: if ( x2 < 0 ) { x2 = 0 Rold = Rnew } else { # # x の二乗を使用して R (グールドの方程式 D) を更新します: # 注: エラー関数 (erfc) は pnorm (Rbasic) に置き換えられます: # ソース: # http://stat.ethz.ch/R-manual/R-patched/library/stats/html/Normal.html Rold = Rnew Rnew =exp (( x2 -1 ) / 2.0 ) * ( 2 * pnorm ( sqrt ( x2 ) / sqrt ( 2 ) * sqrt ( 2 ), lower = FALSE )) } } } else { x2 = 0 } x2 }