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