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        /**
208         * ここからテスト用mainメソッド 。
209         *
210         * @param args 引数
211         */
212        public static void main( final String[] args ) {
213                // 3つの分類で分ける
214                final double[][] train_X = {
215                                {-2.0, 2.0}
216                                ,{-2.1, 1.9}
217                                ,{-1.8, 2.1}
218                                ,{0.0, 0.0}
219                                ,{0.2, -0.2}
220                                ,{-0.1, 0.1}
221                                ,{2.0, -2.0}
222                                ,{2.2, -2.1}
223                                ,{1.9, -2.0}
224                };
225
226                final int[][] train_Y = {
227                                {1, 0, 0}
228                                ,{1, 0, 0}
229                                ,{1, 0, 0}
230                                ,{0, 1, 0}
231                                ,{0, 1, 0}
232                                ,{0, 1, 0}
233                                ,{0, 0, 1}
234                                ,{0, 0, 1}
235                                ,{0, 0, 1}
236                };
237
238                 // test data
239                final double[][] test_X = {
240                                {-2.5, 2.0}
241                                ,{0.1, -0.1}
242                                ,{1.5,-2.5}
243                };
244
245                final double[][] test_Y = new double[test_X.length][train_Y[0].length];
246
247                final HybsLogisticRegression hlr = new HybsLogisticRegression( train_X, train_Y, 0.1, 500, 1 );
248
249                // テスト
250                // このデータでは2番目の条件には入りにくい?
251                for(int i=0; i<test_X.length; i++) {
252                         test_Y[i] = hlr.predict(test_X[i]);
253                         System.out.print( Arrays.toString(test_Y[i]) );
254                }
255        }
256}
257