001package org.opengion.penguin.math.statistics; 002 003import java.util.Arrays; 004import java.util.Collections; 005import java.util.List; 006 007/** 008 * 多項ロジスティック回帰の実装です。 009 * 確率的勾配降下法(SGD)を利用します。 010 * 011 * ロジスティック回帰はn次元の情報からどのグループに所属するかの予測値を得るための手法の一つです。 012 * 013 * 実装は 014 * http://nbviewer.jupyter.org/gist/mitmul/9283713 015 * https://yusugomori.com/projects/deep-learning/ 016 * を参考にしています。 017 */ 018public class HybsLogisticRegression { 019 private final int n_N; // データ個数 020 private final int n_in; // データ次元 021 private final int n_out; // ラベル種別数 022 023 // 写像変数ベクトル f(x) = Wx + b 024 private double[][] vW; 025 private double[] vb; 026 027 /** 028 * コンストラクタ。 029 * 030 * 学習もしてしまう。 031 * 032 * xはデータセット各行がn次元の説明変数となっている。 033 * trainはそれに対する{0,1,0},{1,0,0}のようなラベルを示すベクトルとなる。 034 * 学習率は通常、0.1程度を設定する。 035 * このロジックではループ毎に0.95をかけて徐々に学習率が下がるようにしている。 036 * 全データを利用すると時間がかかる場合があるので、確率的勾配降下法を利用しているが、 037 * 選択個数はデータに対する割合を与える。 038 * データ個数が少ない場合は1をセットすればよい。 039 * 040 * @param data データセット配列 041 * @param label データに対応したラベルを示す配列 042 * @param learning_rate 学習係数(0から1の間の数値) 043 * @param loop 学習のループ回数(ミニバッチを作る回数) 044 * @param minibatch_rate 全体に対するミニバッチの割合(0から1の間の数値) 045 * 046 */ 047 public HybsLogisticRegression(final double data[][], final int label[][], final double learning_rate ,final int loop, final double minibatch_rate ) { 048 // List<Integer> indexList; //シャッフル用 049 050 this.n_N = data.length; 051 this.n_in = data[0].length; 052 this.n_out = label[0].length; // ラベル種別 053 054 vW = new double[n_out][n_in]; 055 vb = new double[n_out]; 056 057 // 確率勾配に利用するための配列インデックス配列 058 final Integer[] random_index = new Integer[n_N]; //プリミティブ型だとasListできないため 059 for( int i=0; i<n_N; i++) { 060 random_index[i] = i; 061 } 062 final List<Integer> indexList = Arrays.asList( random_index ); 063 064 double localRate = learning_rate; 065 for(int epoch=0; epoch<loop; epoch++) { 066 Collections.shuffle( indexList ); 067 // random_index = indexList.toArray(new Integer[indexList.size()]); 068 069 //random_indexの先頭からn_N*minibatch_rate個のものを対象に学習をかける(ミニバッチ) 070 for(int i=0; i< n_N * minibatch_rate; i++) { 071 // final int idx = random_index[i]; 072 final int idx = indexList.get(i); 073 train(data[idx], label[idx], localRate); 074 } 075 localRate *= 0.95; //徐々に学習率を下げて振動を抑える。 076 } 077 } 078 079 /** 080 * データを与えて学習をさせます。 081 * パラメータの1行を与えています。 082 * 083 * 0/1のロジスティック回帰の場合は 084 * ラベルc(0or1)が各xに対して与えられている時 085 * s(x)=σ(Wx+b)=1/(1+ exp(-Wx-b))として、 086 * 確率の対数和L(W,b)の符号反転させたものの偏導関数 087 * ∂L/∂w=-∑x(c-s(x)) 088 * ∂L/∂b=-∑=(c-s(x)) 089 * が最小になるようなW,bの値をパラメータを変えながら求める。 090 * というのが実装になる。(=0を求められないため) 091 * 多次元の場合はシグモイド関数σ(x)の代わりにソフトマックス関数π(x)を利用して 092 * 拡張したものとなる。(以下はソフトマックス関数利用) 093 * 094 * @param in_x 1行分のデータ 095 * @param in_y xに対するラベル 096 * @param lr 学習率 097 * @return 差分配列 098 */ 099 private double[] train( final double[] in_x, final int[] in_y, final double lr ) { 100 final double[] p_y_given_x = new double[n_out]; 101 final double[] dy = new double[n_out]; 102 103 for(int i=0; i<n_out; i++) { 104 p_y_given_x[i] = 0; 105 for(int j=0; j<n_in; j++) { 106 p_y_given_x[i] += vW[i][j] * in_x[j]; 107 } 108 p_y_given_x[i] += vb[i]; 109 } 110 softmax( p_y_given_x ); 111 112 // 勾配の平均で更新? 113 for(int i=0; i<n_out; i++) { 114 dy[i] = in_y[i] - p_y_given_x[i]; 115 116 for(int j=0; j<n_in; j++) { 117 vW[i][j] += lr * dy[i] * in_x[j] / n_N; 118 } 119 120 vb[i] += lr * dy[i] / n_N; 121 } 122 123 return dy; 124 } 125 126 /** 127 * ソフトマックス関数。 128 * π(xi) = exp(xi)/Σexp(x) 129 * @param in_x 変数X 130 */ 131 private void softmax( final double[] in_x ) { 132 // double max = 0.0; 133 double sum = 0.0; 134 135 // for(int i=0; i<n_out; i++) { 136 // if(max < x[i]) { 137 // max = x[i]; 138 // } 139 // } 140 141 for(int i=0; i<n_out; i++) { 142 //x[i] = Math.exp(x[i] - max); // maxとの差分を取ると利点があるのか分からなかった 143 in_x[i] = Math.exp(in_x[i]); 144 sum += in_x[i]; 145 } 146 147 for(int i=0; i<n_out; i++) { 148 in_x[i] /= sum; 149 } 150 } 151 152 /** 153 * 写像式 Wx+b のW、係数ベクトル。 154 * @return 係数ベクトル 155 */ 156 public double[][] getW() { 157 return vW; 158 } 159 160 /** 161 * 写像式 Wx + bのb、バイアス。 162 * @return バイアスベクトル 163 */ 164 public double[] getB() { 165 return vb; 166 } 167 168 /** 169 * 出来た予測式に対して、データを入力してyを出力する。 170 * (yは各ラベルに対する確率分布となる) 171 * @param in_x 予測したいデータ 172 * @return 予測結果 173 */ 174 public double[] predict(final double[] in_x) { 175 final double[] out_y = new double[n_out]; 176 177 for(int i=0; i<n_out; i++) { 178 out_y[i] = 0.; 179 for(int j=0; j<n_in; j++) { 180 out_y[i] += vW[i][j] * in_x[j]; 181 } 182 out_y[i] += vb[i]; 183 } 184 185 softmax(out_y); 186 187 return out_y; 188 } 189 190 //************** ここまでが本体 ************** 191 /** 192 * ここからテスト用mainメソッド 。 193 * 194 * @param args 引数 195 */ 196 public static void main( final String[] args ) { 197 // 3つの分類で分ける 198 final double[][] train_X = { 199 {-2.0, 2.0} 200 ,{-2.1, 1.9} 201 ,{-1.8, 2.1} 202 ,{0.0, 0.0} 203 ,{0.2, -0.2} 204 ,{-0.1, 0.1} 205 ,{2.0, -2.0} 206 ,{2.2, -2.1} 207 ,{1.9, -2.0} 208 }; 209 210 final int[][] train_Y = { 211 {1, 0, 0} 212 ,{1, 0, 0} 213 ,{1, 0, 0} 214 ,{0, 1, 0} 215 ,{0, 1, 0} 216 ,{0, 1, 0} 217 ,{0, 0, 1} 218 ,{0, 0, 1} 219 ,{0, 0, 1} 220 }; 221 222 // test data 223 final double[][] test_X = { 224 {-2.5, 2.0} 225 ,{0.1, -0.1} 226 ,{1.5,-2.5} 227 }; 228 229 final double[][] test_Y = new double[test_X.length][train_Y[0].length]; 230 231 final HybsLogisticRegression hlr = new HybsLogisticRegression( train_X, train_Y, 0.1, 500, 1 ); 232 233 // テスト 234 // このデータでは2番目の条件には入りにくい? 235 for(int i=0; i<test_X.length; i++) { 236 test_Y[i] = hlr.predict(test_X[i]); 237 System.out.print( Arrays.toString(test_Y[i]) ); 238 } 239 } 240} 241