ピアースの基準

ロバスト統計学において、パースの基準は、データ セットから外れ値を 除去するための規則であり、ベンジャミン パースによって考案されました。

ピアースの基準によって除去された外れ値

外れ値の問題

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

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

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

ピアースの基準の用途

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

1852年、彼は研究者に外れ値を拒否すべきかどうかを判断するための最初の有意性検定法を発表しました(Peirce 1852、1878)。尤度比に基づくこの検定法は、そのような行動の賢明さについて国際的な議論を巻き起こしたという点で際立っています(Anscombe、1960、Rider、1933、Stigler、1973a)。

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

一連の観測において、誤差の限界を決定することが提案されている。この限界を超えると、そのような観測が複数存在する場合、その観測はすべて棄却される。この問題を解決するために提案されている原則は、提案された観測を保持することによって得られる誤差のシステムの確率が、それらを棄却することによって得られる誤差のシステムの確率に、それ以上の異常な観測を行わない確率を乗じた値よりも小さい場合、提案された観測は棄却されるべきであるというものである。[ 2 ]メートル{\displaystyle m}n{\displaystyle n}

ホーキンス[ 3 ]はこの基準の公式を提供している。

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

1852年から1867年まで、彼は米国沿岸測量局の経度測定責任者を務め、1867年から1874年までは同測量局の監督官を務めた。この間、彼のテストは、当時最も活発で数学的な傾向を持つこの統計組織の事務員全員によって一貫して採用された。[ 1 ]

ピアースの基準はウィリアム・ショーヴネの著書で議論されている。[ 2 ]

アプリケーション

パースの基準の応用例としては、2つの観測値間の回帰分析(例えば線形回帰)を行うために、観測ペアから品質の低いデータ点を除去することが挙げられます。パースの基準は観測データに依存せず(観測データの特性のみに依存する)、他のプロセスから独立して計算できる、非常に再現性の高いプロセスです。この特徴により、パースの基準は呼び出し関数として記述できるため、外れ値を識別するためのコンピュータアプリケーションに最適です。

以前の試み

1855年、B・A・グールドは、パースの方程式の値を表す値の表を作成することで、パースの基準をより簡単に適用できるようにしようとしました。[ 5 ]グールドのアルゴリズムとパースの基準の実際の適用の間には依然として乖離が残っています。

2003年、S.M.ロス(ニューヘイブン大学)は、グールドのアルゴリズム(現在は「パース法」と呼ばれている)を、新たな例題データセットとアルゴリズムの詳細な解説とともに再発表した。この手法は依然として参照表の使用に依存しており、この研究で更新された(パース基準表)。[ 6 ]

2008年、デンマークの地質学者K.トムセンが疑似コードを書こうとした。[ 7 ] このコードはグールドのアルゴリズムの枠組みを提供したが、ユーザーはピアースやグールドが報告した値を計算できなかった。

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

2013 年には、Gould のアルゴリズムの再検討と高度な Python プログラミング モジュール (numpy や scipy) の利用により、外れ値を識別するための二乗誤差のしきい値を計算できるようになりました。

Python実装

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

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

パースの基準は観測値、フィッティングパラメータ、残差誤差を入力として取らないため、出力をデータに再関連付けする必要があります。すべての二乗誤差の平均(つまり平均二乗誤差)に閾値二乗誤差(つまりこの関数の出力)を乗じることで、外れ値を識別するために使用されるデータ固有の閾値が得られます。

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

Pythonコード

#!/usr/bin/env python3インポートnumpyインポートscipy.specialdef peirce_dev ( N : int , n : int , m : int ) -> float : """Peirce の基準 グールドの方法論に基づく 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 ): # Lamda を計算します# (グールド方程式 A' の 1/(Nn) 乗根): ldiv = r_new ** n if ldiv == 0 : ldiv = 1.0e-6 Lamda = (( Q ** N ) / ( ldiv )) ** ( 1.0 / ( N - n )) # 計算しますx の 2 乗 (グールドの方程式 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 の 2 乗を使用して 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

Javaコード

org.apache.commons.math3.special.Erfをインポートしますパブリッククラス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 )) { // Lamda を計算します// (Gould の方程式 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 の 2 乗を計算します (Gould の方程式 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 { // x2 乗を使用して R(Gould の式 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 ; } }

R実装

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

Python版と同様に、「peirce_dev」関数によって返される二乗誤差(つまり)にモデル適合の平均二乗誤差を乗じて、二乗デルタ値(つまりΔ2)を取得する必要があります。Δ2を使用して、モデル適合の二乗誤差値を比較します。二乗誤差がΔ2を超える観測ペアは外れ値とみなされ、モデルから除外できます。識別された外れ値の数(Δ2とモデル適合の二乗誤差を比較)が想定される外れ値の数(つまりPeirceのn)よりも少なくなるまで、nの値を増加させていく反復子を作成する必要があります。 ×2{\displaystyle x^{2}}

Rコード

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 ) # # の対数グールド方程式 B: LnQN <- k * log ( k ) + ( N - k ) * log ( N - k ) - N * log ( N ) # # グールド方程式 D: R1 <- exp (( x ^ 2 - 1 ) / 2 ) * erfc ( x / sqrt ( 2 )) # # グールド方程式 A' をラムダ置換で R について解く: 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 ) } }リターンx }
peirce_dev <- function ( N , n , m ) { # N :: 観測値の総数# n :: 除去する外れ値の数# m :: モデルの未知数の数 (例: 回帰パラメータ) # # 観測値の数をチェック: if ( N > 1 ) { # Q (Gould 方程式 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 )) { # Lamda (Gould 方程式 A' の 1/(Nn) 乗根) を計算します: ldiv = Rnew ^ n if ( ldiv == 0 ) { ldiv = 1.0e-6 } Lamda = (( Q ^ N ) / ( ldiv )) ^ ( 1.0 / ( N - n )) # # xの2乗を計算します(Gouldの方程式C):x2 = 1.0 + ( N - m - n ) / n * ( 1.0 - Lamda ^ 2.0 ) # # x2が負になる場合は、ゼロに設定します:if ( x2 < 0 ) { x2 = 0 Rold = Rnew } else { # # xの2乗を使用してRを更新します(Gouldの方程式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 }

注記

  1. ^ a b S.M. Stigler, 「初期の州における数理統計」, The Annals of Statistics, vol. 6, no. 2, p. 246, 1978. オンラインで入手可能: https://www.jstor.org/stable/2958876
  2. ^ a b 『 Collected Writings of Peirce』(1982年版) 516ページの編集者注より引用。この引用はショーヴネ著『天文学の手引き』(2:558)を引用している。
  3. ^ DM Hawkins (1980). 「外れ値除去の初期の歴史」『外れ値の識別』(応用確率統計学モノグラフ)Chapman & Hall, 10ページ。
  4. ^ピアース(1878)
  5. ^ a b Gould, BA、「疑わしい観測結果の棄却に関するピアースの基準について、その適用を容易にするための表付き」、天文学ジャーナル、第83号、第4巻、第11号、pp.81–87、1855年。DOI: 10.1086/100480。
  6. ^ Ross, SM、「疑わしい実験データを排除するためのPeirceの基準」、Journal of Engineering Technology、第2巻、第2号、1~12頁、2003年。
  7. ^ Thomsen, K.、「トピック: Peirce の基準で使用するための計算表 - 1855 年と 2008 年」、The Math Forum @ Drexel、2008 年 10 月 5 日に投稿。2013 年 7 月 15 日にアクセス。
  8. ^ C. Dardis、「パッケージ: Peirce」、R-forge、オンラインでアクセス: https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/Peirce/Peirce-manual.pdf?root=peirce
  9. ^ C. Dardis, 「非正規分布外れ値の棄却に関するピアースの基準:適用範囲の定義」『Journal of Statistical Software』(未発表)。オンラインで入手可能: https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/Peirce/PeirceSub.pdf? root=peirce

参考文献