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