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