    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                // ３つの分類で分ける
    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