    001package org.opengion.penguin.math.statistics;
    002
    003import org.apache.commons.math3.stat.StatUtils;
    004import org.apache.commons.math3.linear.RealMatrix;
    005import org.apache.commons.math3.linear.Array2DRowRealMatrix;
    006import org.apache.commons.math3.linear.LUDecomposition;
    007import org.apache.commons.math3.stat.correlation.Covariance;
    008
    009/**
    010 * apache.commons.mathを利用した、マハラノビス距離関係の処理クラスです。
    011 *
    012 * 相関を考慮した距離が求まります。
    013 * 教師無し学習的に、異常値検知に利用可能です。
    014 * 閾値は95%区間の2.448がデフォルトです。(3なら99%)
    015 *
    016 * 「Juan Francisco Quesada-Brizuela」氏の距離計算PGを参照しています。
    017 * 学術的には様々な改良が提案されていますが、このクラスでは単純なマハラノビス距離を扱います。
    018 */
    019public class HybsMahalanobis {
    020
    021        private double[] dataDistance;                  // 元データの各マハラノビス距離
    022        private double[] average;                               // 平均
    023        private RealMatrix covariance;                  // 共分散
    024        private double limen=2.448;                             // 異常値検知をする際の閾値(初期値は95%信頼楕円)
    025
    026        /**
    027         * コンストラクタ。
    028         * 与えたデータマトリクスを元にマハラノビス距離を求めるための準備をします。
    029         * (平均と共分散を求めます）
    030         * 引数calcにtrueをセットすると各点のマハラノビス距離を計算します。
    031         *
    032         * データ = { { 90 ,60 }, { 70, 80 } }
    033         * のような形としてデータを与えます。
    034         *
    035         * @param matrix 値のデータ
    036         * @param calc 距離計算を行うかどうか
    037         */
    038        public HybsMahalanobis( final double[][] matrix, final boolean calc ) {
    039                // 一応元データをセットしておく
    040                final RealMatrix dataMatrix = new Array2DRowRealMatrix( matrix );
    041
    042                // 共分散行列を作成
    043                covariance = new Covariance(matrix).getCovarianceMatrix();
    044                //平均の配列を作成
    045                average = new double[matrix[0].length];
    046                for( int i=0; i<matrix[0].length; i++) {
    047                        average[i] = StatUtils.mean(dataMatrix.getColumn(i));
    048                }
    049
    050                if(calc) {
    051                        dataDistance = new double[matrix.length];
    052                        for( int i=0; i< matrix.length; i++ ) {
    053        //                      dataDistance[i] = distance( matrix[i] );
    054                                dataDistance[i] = distance( covariance,matrix[i],average );             // PMD:Overridable method 'distance' called during object construction
    055                        }
    056                        // 標準偏差、平均を取る場合
    057                        //double maxDst = StatUtils.max( dataDistance );                // 最大
    058                        //double vrDst = StatUtils.variance( dataDistance );    // 分散
    059                        //double shigma = Math.sqrt(vrDst);                                             // シグマ
    060                        //double meanDst = StatUtils.mean( dataDistance );              // 平均
    061                }
    062        }
    063
    064        /**
    065         * 距離計算がtrueの形の簡易版コンストラクタです。
    066         *
    067         * @param matrix 値データ
    068         */
    069        public HybsMahalanobis(final double[][] matrix) {
    070                this(matrix,true);
    071        }
    072
    073        /**
    074         * コンストラクタ。
    075         * 計算済みの共分散と平均、閾値を与えるパターン。
    076         *
    077         * @param covarianceData 共分散
    078         * @param averageData  平均配列
    079         */
    080        public HybsMahalanobis(final double[][] covarianceData, final double[] averageData) {
    081                this.covariance = new Array2DRowRealMatrix(covarianceData);
    082                this.average = averageData;
    083        }
    084
    085        /**
    086         * 平均配列を返します。
    087         *
    088         * @return 平均
    089         */
    090        public double[] getAverage() {
    091                return average;
    092        }
    093
    094        /**
    095         * 共分散配列を返します。
    096         *
    097         * @return 共分散
    098         */
    099        public double[][] getCovariance() {
    100                return covariance.getData();
    101        }
    102
    103        /**
    104         * 閾値を返します。
    105         *
    106         * @return 閾値
    107         */
    108        public double getLimen() {
    109                return limen;
    110        }
    111
    112        /**
    113         * 平均配列をセットします。
    114         *
    115         * @param ave 平均
    116         */
    117        public void setAverage( final double[] ave ) {
    118                this.average = ave;
    119        }
    120
    121        /**
    122         * 共分散配列をセットします。
    123         *
    124         * @param cvr 共分散
    125         */
    126        public void setCovariance( final double[][] cvr ) {
    127                this.covariance = new Array2DRowRealMatrix(cvr);
    128        }
    129
    130        /**
    131         * 閾値をセットします。
    132         * 距離の二乗がカイ2乗分布となるため、
    133         * 初期値は2.448で、95%区間を意味します。
    134         * 2が86%、3が99％です。
    135         *
    136         * @param lim 閾値
    137         */
    138        public void setLimen( final double lim ) {
    139                this.limen = lim;
    140        }
    141
    142        /**
    143         * コンストラクタで元データを与え、計算させた場合のマハラノビス距離の配列を返します。
    144         *
    145         * @return 各点のマハラノビス距離の配列
    146         */
    147        public double[] getDataDistance() {
    148                return dataDistance;
    149        }
    150
    151        /**
    152         * マハラノビス距離を計算します。
    153         *
    154         * @param vec 判定する点（ベクトル）
    155         * @return マハラノビス距離
    156         */
    157        public double distance( final double[] vec) {
    158                return distance( covariance, vec, average  );
    159        }
    160
    161        /**
    162         * 与えたベクトルが閾値を超えたマハラノビス距離かどうかを判定します。
    163         * 閾値以下ならtrue、超えている場合はfalseを返します。
    164         * （異常値判定）
    165         *
    166         * @param vec 判定する点（ベクトル）
    167         * @return 閾値以下かどうか
    168         */
    169        public boolean check( final double[] vec) {
    170//              final double dist =  distance( covariance, vec, average  );
    171//              return ( dist <= limen );
    172                return distance( covariance, vec, average ) <= limen ;                          // 6.9.7.0 (2018/05/14) PMD Useless parentheses.
    173        }
    174
    175        /**
    176         * 平均、共分散を利用して対象ベクトルとの距離を測ります。
    177         *
    178         * @og.rev 6.9.8.0 (2018/05/28) det を削除します。
    179         * @og.rev 6.9.9.0 (2018/08/20) ﾛｼﾞｯｸ修正ミス
    180         *
    181         * @param mtx1 共分散行列
    182         * @param vec1 距離を測りたいベクトル
    183         * @param vec2 平均ﾍﾞｸﾄﾙ
    184         * @return マハラノビス距離
    185         */
    186        private double distance(final RealMatrix mtx1, final double vec1[], final double vec2[]) {
    187                // ﾏﾊﾗﾉﾋﾞｽ距離の公式
    188                // ﾏﾊﾗﾉﾋﾞｽ距離 = (v1-v2)*inv(m1)*t(v1-v2)
    189                // inv():逆行列
    190                // t():転置行列
    191
    192                // ※getDeterminantは行列式(正方行列に対して定義される量)を取得
    193                // javaの処理上、v1.lengthが2以上の場合、1/(v1.length)が0になる。
    194                // その結果、行列式を0乗になるので、detに1が設定される。
    195                // この式はﾏﾊﾗﾉﾋﾞｽ距離を求める公式にない為、不要な処理？
    196//              final double det = Math.pow((new LUDecomposition(mtx1).getDeterminant()), 1/(vec1.length));
    197                // 6.9.8.0 (2018/05/28) det を削除します。
    198                // PMD で、1/(vec1.length) が指摘され、FindBugs で、整数同士の割り算を、double にキャストしている警告が出ます。
    199                // vec1 の配列が1の場合のみ有効にするなら、他の方法があるはずで、不要な処理？ というコメントとあわせて、
    200                // とりあえずコメントアウトしておきます。
    201        //      final double det = Math.pow( new LUDecomposition(mtx1).getDeterminant() , 1/ vec1.length );                     // 6.9.7.0 (2018/05/14) PMD Useless parentheses.
    202
    203                double[] temp = new double[vec1.length];
    204                // (x - y)を計算
    205                for(int i=0; i < vec1.length; i++) {
    206                        temp[i] = vec1[i]-vec2[i];                                              // 6.9.7.0 (2018/05/14) PMD Useless parentheses.
    207                }
    208
    209        //      double[] tempSub = new double[vec1.length];
    210
    211        //      // (x - y)を計算
    212        //      for(int i=0; i < vec1.length; i++) {
    213        //              tempSub[i] = vec1[i]-vec2[i];                                   // 6.9.7.0 (2018/05/14) PMD Useless parentheses.
    214        //      }
    215
    216        //      double[] temp = new double[vec1.length];
    217
    218        //      //      (x - y) * det 不要な処理？
    219        //      for(int i=0; i < temp.length; i++) {
    220        //              temp[i] = tempSub[i]*det;
    221        //      }
    222
    223                // m2: (x - y)を行列に変換
    224                final RealMatrix m2 = new Array2DRowRealMatrix( new double[][] { temp } );
    225
    226                // m3: m2 * 共分散行列の逆行列
    227                final RealMatrix m3 = m2.multiply( new LUDecomposition(mtx1).getSolver().getInverse() );
    228
    229                // m4: m3 * (x-y)の転置行列
    230//              final RealMatrix m4 = m3.multiply((new Array2DRowRealMatrix(new double[][] { temp })).transpose());
    231//              final RealMatrix m4 = m3.multiply( new Array2DRowRealMatrix( new double[][] { temp } ) ).transpose() ;                  // 6.9.7.0 (2018/05/14) PMD Useless parentheses.
    232                final RealMatrix m4 = m3.multiply( new Array2DRowRealMatrix( new double[][] { temp } ).transpose() ) ;                  // 6.9.9.0 (2018/08/20) ﾛｼﾞｯｸ修正ミス
    233
    234                // m4の平方根を返す
    235                return Math.sqrt(m4.getEntry(0, 0));
    236        }
    237
    238        // *** ここまでが本体 ***
    239
    240        /**
    241         * ここからテスト用mainメソッド。
    242         *
    243         * @param args ****************************************
    244         */
    245        public static void main( final String [] args ) {
    246                // 幾何的には、これらの重心を中心とした楕円の中に入っているかどうかを判定
    247                final double[][] data = {
    248                  {2, 10},
    249                  {4, 21},
    250                  {6, 27},
    251                  {8, 41},
    252                  {10, 50}
    253                };
    254
    255                final double[] test = {12, 50};
    256                final double[] test2 = {12, 59};
    257
    258                final HybsMahalanobis rtn = new HybsMahalanobis(data);
    259
    260                System.out.println( java.util.Arrays.toString(rtn.getDataDistance()) );
    261
    262                System.out.println(rtn.check( test ));
    263                System.out.println(rtn.check( test2 ));
    264        }
    265}
    266