    001package org.opengion.penguin.math.statistics;
    002
    003import java.util.Arrays;
    004
    005/**
    006 * 独自実装の二次回帰計算クラスです。
    007 * f(x) = c1x^2 + c2x + c3
    008 * の曲線を求めます。
    009 */
    010public class HybsSquadraticRegression implements HybsRegression {
    011        private final double[] cnst = new double[3] ;           // 係数(０次、１次、２次)
    012        private double rsquare;         // 決定係数 今のところ求めていない
    013
    014        /**
    015         * コンストラクタ。
    016         * 与えた二次元データを元に二次回帰を計算します。
    017         *
    018         * @param data xとyの組み合わせの配列
    019         */
    020        public HybsSquadraticRegression( final double[][] data ) {
    021                //二次回帰曲線を求めるが、これはapacheにはなさそうなので自前で計算する。
    022                train( data );
    023        }
    024
    025        /**
    026         * 係数計算
    027         *
    028         *      c3Σ＋c2Σx＋c1Σx^2＝Σy
    029         *      c3Σx＋c2Σ(x^2)＋c1Σx^3＝Σ(xy)
    030         *      c3Σ(x^2)＋c2Σ(x^3)＋c1Σ(x^4)＝Σ(ｘ^2*y)
    031         *      この三元連立方程式を解くことになる。
    032         *
    033         * @param data x,yの配列
    034         */
    035        private void train( final double[][] data ) {
    036                // xの二乗等の総和用
    037                final int data_n=data.length;;
    038                double sumx2    = 0;
    039                double sumx             = 0;
    040                double sumxy    = 0;
    041                double sumy             = 0;
    042                double sumx3    = 0;
    043                double sumx2y   = 0;
    044                double sumx4    = 0;
    045
    046                // まずは計算に使うための和を計算
    047                for( int i=0; i < data_n; i++ ) {
    048                        final double data_x     = data[i][0];
    049                        final double data_y     = data[i][1];
    050                        final double x2         = data_x*data_x;
    051
    052                        sumx    += data_x;
    053                        sumx2   += x2;
    054                        sumxy   += data_x * data_y;
    055                        sumy    += data_y;
    056                        sumx3   += x2 * data_x;
    057                        sumx2y  += x2 * data_y;
    058                        sumx4   += x2 * x2;
    059                }
    060
    061                // ガウス・ジョルダン法で係数計算
    062                final double diffx2 = sumx2 - sumx * sumx / data_n;
    063                final double diffxy = sumxy - sumx * sumy / data_n;
    064                final double diffx3 = sumx3 - sumx2 * sumx /data_n;
    065                final double diffx2y = sumx2y - sumx2 * sumy /data_n;
    066                final double diffx4 = sumx4 - sumx2 * sumx2 /data_n;
    067                final double diffd = diffx2 * diffx4 - diffx3 * diffx3;
    068
    069                cnst[2] = ( diffx2y * diffx2 - diffxy * diffx3 ) / diffd;
    070                cnst[1] = ( diffxy * diffx4 - diffx2y * diffx3 ) / diffd;
    071                cnst[0] = sumy/data_n - cnst[1]*sumx/ data_n - cnst[2]*sumx2/data_n;
    072
    073                rsquare = 0;            // 決定係数 今のところ求めていない
    074        }
    075
    076        /**
    077         * 決定係数の取得。
    078         * @return 決定係数
    079         */
    080        @Override
    081        public double getRSquare() {
    082                return rsquare;
    083        }
    084
    085        /**
    086         * 係数(０次、１次、２次)の順にセットした配列を返します。
    087         *
    088         * @return 係数の配列
    089         */
    090        @Override
    091        public double[] getCoefficient() {
    092                return Arrays.copyOf( cnst,cnst.length );
    093        }
    094
    095        /**
    096         * c2*x^2 + c1*x + c0を計算。
    097         *
    098         * @param in_x 必要な大きさの変数配列
    099         * @return 計算結果
    100         */
    101        @Override
    102        public double predict( final double... in_x ) {
    103                return cnst[2] * in_x[0] * in_x[0] + cnst[1] * in_x[0] + cnst[0];
    104        }
    105
    106        //************** ここまでが本体 **************
    107        /**
    108         * ここからテスト用mainメソッド 。
    109         *
    110         * @param args 引数
    111         */
    112        public static void main( final String[] args ) {
    113                final double[][] data = {{1, 2.3}, {2, 5.1}, {3, 9.1}, {4, 16.2}};
    114
    115                final HybsSquadraticRegression sr = new HybsSquadraticRegression(data);
    116
    117                final double[] cnst = sr.getCoefficient();
    118
    119                System.out.println(cnst[2]);
    120                System.out.println(cnst[1]);
    121                System.out.println(cnst[0]);
    122
    123                System.out.println(sr.predict( 5 ));
    124        }
    125}
    126