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