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        }
173
174        /**
175         * 平均、共分散を利用して対象ベクトルとの距離を測ります。
176         *
177         * @param mtx1 共分散行列
178         * @param vec1 距離を測りたいベクトル
179         * @param vec2 平均ベクトル
180         * @return マハラノビス距離
181         */
182        private double distance(final RealMatrix mtx1, final double vec1[], final double vec2[]) {
183                // マハラノビス距離の公式
184                // マハラノビス距離 = (v1-v2)*inv(m1)*t(v1-v2)
185                // inv():逆行列
186                // t():転置行列
187
188                // ※getDeterminantは行列式(正方行列に対して定義される量)を取得
189                // javaの処理上、v1.lengthが2以上の場合、1/(v1.length)が0になる。
190                // その結果、行列式を0乗になるので、detに1が設定される。
191                // この式はマハラノビス距離を求める公式にない為、不要な処理?
192                final double det = Math.pow((new LUDecomposition(mtx1).getDeterminant()), 1/(vec1.length));
193
194                double[] tempSub = new double[vec1.length];
195
196                // (x - y)を計算
197                for(int i=0; i < vec1.length; i++) {
198                        tempSub[i] = (vec1[i]-vec2[i]);
199                }
200
201                double[] temp = new double[vec1.length];
202
203                //      (x - y) * det 不要な処理?
204                for(int i=0; i < temp.length; i++) {
205                        temp[i] = tempSub[i]*det;
206                }
207
208                // m2: (x - y)を行列に変換
209                final RealMatrix m2 = new Array2DRowRealMatrix(new double[][] { temp });
210
211                // m3: m2 * 共分散行列の逆行列
212                final RealMatrix m3 = m2.multiply(new LUDecomposition(mtx1).getSolver().getInverse());
213
214                // m4: m3 * (x-y)の転置行列
215                final RealMatrix m4 = m3.multiply((new Array2DRowRealMatrix(new double[][] { temp })).transpose());
216
217                // m4の平方根を返す
218                return Math.sqrt(m4.getEntry(0, 0));
219        }
220
221        // *** ここまでが本体 ***
222
223        /**
224         * ここからテスト用mainメソッド。
225         *
226         * @param args ****************************************
227         */
228        public static void main( final String [] args ) {
229                // 幾何的には、これらの重心を中心とした楕円の中に入っているかどうかを判定
230                final double[][] data = {
231                  {2, 10},
232                  {4, 21},
233                  {6, 27},
234                  {8, 41},
235                  {10, 50}
236                };
237
238                final double[] test = {12, 50};
239                final double[] test2 = {12, 59};
240
241                final HybsMahalanobis rtn = new HybsMahalanobis(data);
242
243                System.out.println( java.util.Arrays.toString(rtn.getDataDistance()) );
244
245                System.out.println(rtn.check( test ));
246                System.out.println(rtn.check( test2 ));
247        }
248}
249