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