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