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