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