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] ;           // 係数(0次、1次、2次)
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)=Σ(x^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         * 係数(0次、1次、2次)の順にセットした配列を返します。
093         *
094         * @return 係数の配列
095         */
096        @Override       // HybsRegression
097        public double[] getCoefficient() {
098                return Arrays.copyOf( cnst,cnst.length );
099        }
100
101        /**
102         * 決定係数の取得。
103         * @return 決定係数
104         */
105        @Override       // HybsRegression
106        public double getRSquare() {
107                return rsquare;
108        }
109
110        /**
111         * c2*x^2 + c1*x + c0を計算。
112         *
113         * @param in_x 必要な大きさの変数配列
114         * @return 計算結果
115         */
116        @Override       // HybsRegression
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