package org.geogebra.common.kernel.statistics; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.DecompositionSolver; import org.apache.commons.math3.linear.QRDecomposition; import org.apache.commons.math3.linear.RealMatrix; import org.geogebra.common.kernel.geos.GeoElement; import org.geogebra.common.kernel.geos.GeoList; import org.geogebra.common.kernel.geos.GeoPoint; import org.geogebra.common.util.debug.Log; import edu.umd.cs.findbugs.annotations.SuppressFBWarnings; /* GeoGebra - Dynamic Mathematics for Everyone http://www.geogebra.org This file is part of GeoGebra. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation. */ /** * <pre> * <h3>RegressionMath:</h3> * * RegressionMath is a library of sums, determinants and parameter calculations * used by the FitXxx[GeoList]:GeoFunction commands. * * Might be problems if callers are running in separate threads. Is this a * problem? * * * @author Hans-Petter Ulven * @version 20.02.10 Start 24.04.08 Update 15.11.08: public * det22(...),...,det44(...) for use in FitSin() and FitLogistic() * Update 27.01.09: doPolyN():boolean getPar():double[] to serve the * extending of FitPoly[...] to degree>=5 Matrix operations based on * Jama. Update 20.02.2010: Changed from JaMa to Apache matrix library * See doPolyN() Got rid of r=corrcoff(), not used. Exists as separate * command now. * * <ul> * <b>--- Interface: ---</b> * <li>RegressionMath(GeoList) * <li>det33(...), det44(...) determinants. (Faster than Gauss for n<5) * <li>doLinReg(),doQuadReg(),doCubicReg(),doQuartRet() * <li>doExpReg(),doLogReg(),doPowReg() * <li>getSigmaX(),getSigmaX2(),[getSigmaX3()..getSigmaX5() * <li>getSigmaY(),[getSigmaY2(),getSigmaXY(),getSigmaX2Y(),getSigmaX3Y * (),getSigmaX4Y()] * <li>getP1(),getP2(),getP3(),getP4(),getP5() //Parameters for * regression function * <li>getR() //regression coeffisient. (When users get regression, * they will certainly ask for this...) * </ul> */ public final class RegressionMath { public final static int LINEAR = 1, QUAD = 2, CUBIC = 3, QUART = 4, EXP = 5, LOG = 6, POW = 7; // / --- Properties --- /// private boolean error = false; // private int regtype = LINEAR; //Default private double // r, //Reg-coeff p1, p2, p3, p4, p5, // Parameters sigmax, sigmax2, sigmax3, sigmax4, // Sums of x,x^2,... sigmax5, sigmax6, sigmax7, sigmax8, sigmay, sigmay2, sigmaxy, sigmax2y, sigmax3y, sigmax4y; private GeoList geolist; private double[] xlist; private double[] ylist; private int size; // 27.01.09: private double[][] marray, yarray; // For (M_T*M)*Par=(M_T*Y) private double[] pararray; // Parameter array // / --- Interface --- /// public final double getP1() { return p1; } public final double getP2() { return p2; } public final double getP3() { return p3; } public final double getP4() { return p4; } public final double getP5() { return p5; } // public final double getR() {return r; } public final double getSigmaX() { return sigmax; } public final double getSigmaX2() { return sigmax2; } public final double getSigmaY() { return sigmay; } public final double getSigmaY2() { return sigmay2; } public final double getSigmaXy() { return sigmaxy; } // 27.01.09: /** Returns array with calculated parameters */ public final double[] getPar() { return pararray; } /** Does the Polynom regression for degree > 4 */ public final boolean doPolyN(GeoList gl, int degree) { error = false; geolist = gl; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } try { /* * Old Jama version: long time=System.currentTimeMillis(); * makeMatrixArrays(degree); //make marray and yarray Matrix M=new * Matrix(marray); Matrix Y=new Matrix(yarray); Matrix * Par=M.solve(Y); //Par.print(3,3); * pararray=Par.getRowPackedCopy(); * System.out.println(System.currentTimeMillis()-time); */ makeMatrixArrays(degree); // make marray and yarray RealMatrix M = new Array2DRowRealMatrix(marray, false); DecompositionSolver solver = new QRDecomposition(M).getSolver(); // time=System.currentTimeMillis(); RealMatrix Y = new Array2DRowRealMatrix(yarray, false); RealMatrix P = solver.solve(Y); pararray = P.getColumn(0); // System.out.println(System.currentTimeMillis()-time); // diff(pararray,par); } catch (Throwable t) { Log.debug(t.toString()); error = true; } // try-catch. ToDo: A bit more fine-grained error-handling... return !error; }// doPolyN() public final boolean doLinear(GeoList gl) { error = false; geolist = gl; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } doSums(LINEAR); // calculate neccessary sigmas if (error) { return false; } double n = det22(1.0d * size, sigmax, sigmax, sigmax2); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det22(sigmay, sigmax, sigmaxy, sigmax2) / n; p2 = det22(size, sigmay, sigmax, sigmaxy) / n; // r=corrCoeff(); return true; }// doLinearReg(GeoList) public final boolean doQuad(GeoList gl) { error = false; geolist = gl; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } doSums(QUAD); // calculate neccessary sigmas if (error) { return false; } double n = det33(1.0d * size, sigmax, sigmax2, sigmax, sigmax2, sigmax3, sigmax2, sigmax3, sigmax4); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det33(sigmay, sigmax, sigmax2, sigmaxy, sigmax2, sigmax3, sigmax2y, sigmax3, sigmax4) / n; p2 = det33(1.0d * size, sigmay, sigmax2, sigmax, sigmaxy, sigmax3, sigmax2, sigmax2y, sigmax4) / n; p3 = det33(1.0d * size, sigmax, sigmay, sigmax, sigmax2, sigmaxy, sigmax2, sigmax3, sigmax2y) / n; // r=0.0d; // Not useful return true; }// doQuad(Geolist) public final boolean doCubic(GeoList gl) { error = false; geolist = gl; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } doSums(CUBIC); // calculate neccessary sigmas if (error) { return false; } double n = det44(1.0d * size, sigmax, sigmax2, sigmax3, sigmax, sigmax2, sigmax3, sigmax4, sigmax2, sigmax3, sigmax4, sigmax5, sigmax3, sigmax4, sigmax5, sigmax6); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det44(sigmay, sigmax, sigmax2, sigmax3, sigmaxy, sigmax2, sigmax3, sigmax4, sigmax2y, sigmax3, sigmax4, sigmax5, sigmax3y, sigmax4, sigmax5, sigmax6) / n; p2 = det44(size, sigmay, sigmax2, sigmax3, sigmax, sigmaxy, sigmax3, sigmax4, sigmax2, sigmax2y, sigmax4, sigmax5, sigmax3, sigmax3y, sigmax5, sigmax6) / n; p3 = det44(size, sigmax, sigmay, sigmax3, sigmax, sigmax2, sigmaxy, sigmax4, sigmax2, sigmax3, sigmax2y, sigmax5, sigmax3, sigmax4, sigmax3y, sigmax6) / n; p4 = det44(size, sigmax, sigmax2, sigmay, sigmax, sigmax2, sigmax3, sigmaxy, sigmax2, sigmax3, sigmax4, sigmax2y, sigmax3, sigmax4, sigmax5, sigmax3y) / n; // r=0.0d; // Not useful return true; }// doCubic(Geolist) public final boolean doQuart(GeoList gl) { error = false; geolist = gl; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } doSums(QUART); // calculate neccessary sigmas if (error) { return false; } double n = det55(1.0d * size, sigmax, sigmax2, sigmax3, sigmax4, sigmax, sigmax2, sigmax3, sigmax4, sigmax5, sigmax2, sigmax3, sigmax4, sigmax5, sigmax6, sigmax3, sigmax4, sigmax5, sigmax6, sigmax7, sigmax4, sigmax5, sigmax6, sigmax7, sigmax8); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det55(sigmay, sigmax, sigmax2, sigmax3, sigmax4, sigmaxy, sigmax2, sigmax3, sigmax4, sigmax5, sigmax2y, sigmax3, sigmax4, sigmax5, sigmax6, sigmax3y, sigmax4, sigmax5, sigmax6, sigmax7, sigmax4y, sigmax5, sigmax6, sigmax7, sigmax8) / n; p2 = det55(1.0d * size, sigmay, sigmax2, sigmax3, sigmax4, sigmax, sigmaxy, sigmax3, sigmax4, sigmax5, sigmax2, sigmax2y, sigmax4, sigmax5, sigmax6, sigmax3, sigmax3y, sigmax5, sigmax6, sigmax7, sigmax4, sigmax4y, sigmax6, sigmax7, sigmax8) / n; p3 = det55(1.0d * size, sigmax, sigmay, sigmax3, sigmax4, sigmax, sigmax2, sigmaxy, sigmax4, sigmax5, sigmax2, sigmax3, sigmax2y, sigmax5, sigmax6, sigmax3, sigmax4, sigmax3y, sigmax6, sigmax7, sigmax4, sigmax5, sigmax4y, sigmax7, sigmax8) / n; p4 = det55(1.0d * size, sigmax, sigmax2, sigmay, sigmax4, sigmax, sigmax2, sigmax3, sigmaxy, sigmax5, sigmax2, sigmax3, sigmax4, sigmax2y, sigmax6, sigmax3, sigmax4, sigmax5, sigmax3y, sigmax7, sigmax4, sigmax5, sigmax6, sigmax4y, sigmax8) / n; p5 = det55(1.0d * size, sigmax, sigmax2, sigmax3, sigmay, sigmax, sigmax2, sigmax3, sigmax4, sigmaxy, sigmax2, sigmax3, sigmax4, sigmax5, sigmax2y, sigmax3, sigmax4, sigmax5, sigmax6, sigmax3y, sigmax4, sigmax5, sigmax6, sigmax7, sigmax4y) / n; // r=0.0d; // Not useful return true; }// doQuart(Geolist) public final boolean doExp(GeoList gl) { error = false; geolist = gl; double y = 0.0d; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } if (size == 0) { return false; } double ySign = ylist[0] < 0 ? -1 : 1; // Transform y->ln(y): for (int i = 0; i < size; i++) { y = ylist[i] * ySign; if (y < 0.0d) { // log(minus)! return false; } ylist[i] = Math.log(y); } // for all y doSums(LINEAR); // calculate neccessary sigmas if (error) { return false; } double n = det22(1.0d * size, sigmax, sigmax, sigmax2); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det22(sigmay, sigmax, sigmaxy, sigmax2) / n; p2 = det22(size, sigmay, sigmax, sigmaxy) / n; // transform back: p1 = Math.exp(p1) * ySign; // r=corrCoeff(); return true; }// doExp(GeoList) public final boolean doLog(GeoList gl) { error = false; geolist = gl; double x = 0.0d; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } // Transform x->ln(x): for (int i = 0; i < size; i++) { x = xlist[i]; if (x < 0.0d) { // log(minus)! return false; } xlist[i] = Math.log(xlist[i]); } // for all x doSums(LINEAR); // calculate neccessary sigmas if (error) { return false; } double n = det22(1.0d * size, sigmax, sigmax, sigmax2); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det22(sigmay, sigmax, sigmaxy, sigmax2) / n; p2 = det22(size, sigmay, sigmax, sigmaxy) / n; // No transformation of p1 or p2 neccessary // r=corrCoeff(); return true; }// doLog(GeoList) public final boolean doPow(GeoList gl) { error = false; geolist = gl; double x, y; size = geolist.size(); getPoints(); // getPoints from geolist if (error) { return false; } if (size == 0) { return false; } double ySign = ylist[0] < 0 ? -1 : 1; // Transform y->ln(y) and x->ln(x): for (int i = 0; i < size; i++) { x = xlist[i]; y = ylist[i] * ySign; if ((x < 0.0d) || (y < 0.0d)) { // log(minus)! return false; } xlist[i] = Math.log(x); ylist[i] = Math.log(y); } // for all x doSums(LINEAR); // calculate neccessary sigmas if (error) { return false; } double n = det22(1.0d * size, sigmax, sigmax, sigmax2); if (Math.abs(n - 0.0d) < 1.0E-15d) { return false; } p1 = det22(sigmay, sigmax, sigmaxy, sigmax2) / n; p2 = det22(size, sigmay, sigmax, sigmaxy) / n; // Transform back: p1 = Math.exp(p1) * ySign; // r=corrCoeff(); return true; }// doPow(GeoList) public final static double det22( // 15.11.08: public for FitSin and // FitLogisticc double a11, double a12, double a21, double a22) { return a11 * a22 - a21 * a12; }// det22() public final static double det33( // 15.11.08: public for FitSin and FitLogisticc double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33) { return a11 * (a22 * a33 - a32 * a23) - a12 * (a21 * a33 - a31 * a23) + a13 * (a21 * a32 - a31 * a22); }// det33() public final static double det44( // 15.11.08: public for FitSin and FitLogisticc double a11, double a12, double a13, double a14, double a21, double a22, double a23, double a24, double a31, double a32, double a33, double a34, double a41, double a42, double a43, double a44) { return a11 * a22 * a33 * a44 - a11 * a22 * a34 * a43 - a11 * a23 * a32 * a44 + a11 * a23 * a42 * a34 + a11 * a32 * a24 * a43 - a11 * a24 * a33 * a42 - a12 * a21 * a33 * a44 + a12 * a21 * a34 * a43 + a12 * a31 * a23 * a44 - a12 * a31 * a24 * a43 - a12 * a23 * a41 * a34 + a12 * a41 * a24 * a33 + a21 * a13 * a32 * a44 - a21 * a13 * a42 * a34 - a21 * a14 * a32 * a43 + a21 * a14 * a33 * a42 - a13 * a22 * a31 * a44 + a13 * a22 * a41 * a34 + a13 * a31 * a24 * a42 - a13 * a32 * a41 * a24 + a22 * a31 * a14 * a43 - a22 * a14 * a41 * a33 - a31 * a14 * a23 * a42 + a14 * a23 * a32 * a41; }// det44() public final static double det55( // 15.11.08: public for FitSin and FitLogisticc double a11, double a12, double a13, double a14, double a15, double a21, double a22, double a23, double a24, double a25, double a31, double a32, double a33, double a34, double a35, double a41, double a42, double a43, double a44, double a45, double a51, double a52, double a53, double a54, double a55) { return a11 * a22 * a33 * a44 * a55 - a11 * a22 * a33 * a45 * a54 - a11 * a22 * a34 * a43 * a55 + a11 * a22 * a34 * a53 * a45 + a11 * a22 * a43 * a35 * a54 - a11 * a22 * a35 * a44 * a53 - a11 * a23 * a32 * a44 * a55 + a11 * a23 * a32 * a45 * a54 + a11 * a23 * a42 * a34 * a55 - a11 * a23 * a42 * a35 * a54 - a11 * a23 * a34 * a52 * a45 + a11 * a23 * a52 * a35 * a44 + a11 * a32 * a24 * a43 * a55 - a11 * a32 * a24 * a53 * a45 - a11 * a32 * a25 * a43 * a54 + a11 * a32 * a25 * a44 * a53 - a11 * a24 * a33 * a42 * a55 + a11 * a24 * a33 * a52 * a45 + a11 * a24 * a42 * a35 * a53 - a11 * a24 * a43 * a52 * a35 + a11 * a33 * a42 * a25 * a54 - a11 * a33 * a25 * a52 * a44 - a11 * a42 * a25 * a34 * a53 + a11 * a25 * a34 * a43 * a52 - a12 * a21 * a33 * a44 * a55 + a12 * a21 * a33 * a45 * a54 + a12 * a21 * a34 * a43 * a55 - a12 * a21 * a34 * a53 * a45 - a12 * a21 * a43 * a35 * a54 + a12 * a21 * a35 * a44 * a53 + a12 * a31 * a23 * a44 * a55 - a12 * a31 * a23 * a45 * a54 - a12 * a31 * a24 * a43 * a55 + a12 * a31 * a24 * a53 * a45 + a12 * a31 * a25 * a43 * a54 - a12 * a31 * a25 * a44 * a53 - a12 * a23 * a41 * a34 * a55 + a12 * a23 * a41 * a35 * a54 + a12 * a23 * a51 * a34 * a45 - a12 * a23 * a51 * a35 * a44 + a12 * a41 * a24 * a33 * a55 - a12 * a41 * a24 * a35 * a53 - a12 * a41 * a33 * a25 * a54 + a12 * a41 * a25 * a34 * a53 - a12 * a24 * a33 * a51 * a45 + a12 * a24 * a51 * a43 * a35 + a12 * a33 * a51 * a25 * a44 - a12 * a51 * a25 * a34 * a43 + a21 * a13 * a32 * a44 * a55 - a21 * a13 * a32 * a45 * a54 - a21 * a13 * a42 * a34 * a55 + a21 * a13 * a42 * a35 * a54 + a21 * a13 * a34 * a52 * a45 - a21 * a13 * a52 * a35 * a44 - a21 * a14 * a32 * a43 * a55 + a21 * a14 * a32 * a53 * a45 + a21 * a14 * a33 * a42 * a55 - a21 * a14 * a33 * a52 * a45 - a21 * a14 * a42 * a35 * a53 + a21 * a14 * a43 * a52 * a35 + a21 * a32 * a15 * a43 * a54 - a21 * a32 * a15 * a44 * a53 - a21 * a15 * a33 * a42 * a54 + a21 * a15 * a33 * a52 * a44 + a21 * a15 * a42 * a34 * a53 - a21 * a15 * a34 * a43 * a52 - a13 * a22 * a31 * a44 * a55 + a13 * a22 * a31 * a45 * a54 + a13 * a22 * a41 * a34 * a55 - a13 * a22 * a41 * a35 * a54 - a13 * a22 * a51 * a34 * a45 + a13 * a22 * a51 * a35 * a44 + a13 * a31 * a24 * a42 * a55 - a13 * a31 * a24 * a52 * a45 - a13 * a31 * a42 * a25 * a54 + a13 * a31 * a25 * a52 * a44 - a13 * a32 * a41 * a24 * a55 + a13 * a32 * a41 * a25 * a54 + a13 * a32 * a24 * a51 * a45 - a13 * a32 * a51 * a25 * a44 + a13 * a41 * a24 * a52 * a35 - a13 * a41 * a25 * a34 * a52 - a13 * a24 * a42 * a51 * a35 + a13 * a42 * a51 * a25 * a34 + a22 * a31 * a14 * a43 * a55 - a22 * a31 * a14 * a53 * a45 - a22 * a31 * a15 * a43 * a54 + a22 * a31 * a15 * a44 * a53 - a22 * a14 * a41 * a33 * a55 + a22 * a14 * a41 * a35 * a53 + a22 * a14 * a33 * a51 * a45 - a22 * a14 * a51 * a43 * a35 + a22 * a41 * a15 * a33 * a54 - a22 * a41 * a15 * a34 * a53 - a22 * a15 * a33 * a51 * a44 + a22 * a15 * a51 * a34 * a43 - a31 * a14 * a23 * a42 * a55 + a31 * a14 * a23 * a52 * a45 + a31 * a14 * a42 * a25 * a53 - a31 * a14 * a25 * a43 * a52 + a31 * a23 * a15 * a42 * a54 - a31 * a23 * a15 * a52 * a44 - a31 * a15 * a24 * a42 * a53 + a31 * a15 * a24 * a43 * a52 + a14 * a23 * a32 * a41 * a55 - a14 * a23 * a32 * a51 * a45 - a14 * a23 * a41 * a52 * a35 + a14 * a23 * a42 * a51 * a35 - a14 * a32 * a41 * a25 * a53 + a14 * a32 * a51 * a25 * a43 + a14 * a41 * a33 * a25 * a52 - a14 * a33 * a42 * a51 * a25 - a23 * a32 * a41 * a15 * a54 + a23 * a32 * a15 * a51 * a44 + a23 * a41 * a15 * a34 * a52 - a23 * a15 * a42 * a51 * a34 + a32 * a41 * a15 * a24 * a53 - a32 * a15 * a24 * a51 * a43 - a41 * a15 * a24 * a33 * a52 + a15 * a24 * a33 * a42 * a51; }// det55 // / --- Private --- /// /* Do whatever sums neccessary */ @SuppressFBWarnings({ "SF_SWITCH_FALLTHROUGH", "missing break is deliberate" }) private final void doSums(int degree) { // do whatever sums neccessary double x, y, xx, xy; sigmax = sigmax2 = sigmax3 = sigmax4 = sigmax5 = sigmax6 = sigmax7 = sigmax8 = sigmaxy = sigmax2y = sigmax3y = sigmax4y = sigmay = sigmay2 = 0.0d; for (int i = 0; i < size; i++) { x = xlist[i]; y = ylist[i]; xx = x * x; xy = x * y; // save some calculations switch (degree) { // fall through switch default: // do nothing break; case 4: sigmax4y += x * xx * xy; sigmax7 += xx * xx * xx * x; sigmax8 += xx * xx * xx * xx; //$FALL-THROUGH$ case 3: sigmax3y += xx * xy; sigmax5 += xx * xx * x; sigmax6 += xx * xx * xx; //$FALL-THROUGH$ case 2: sigmax2y += xx * y; sigmax3 += xx * x; sigmax4 += xx * xx; //$FALL-THROUGH$ case 1: sigmay2 += y * y; // r only for linear } // always do these sigmax += x; sigmax2 += xx; sigmaxy += xy; sigmay += y; } } /* Get points to local array */ private final void getPoints() { // double x,y; double xy[] = new double[2]; GeoElement geoelement; // GeoPoint geopoint; xlist = new double[size]; ylist = new double[size]; for (int i = 0; i < size; i++) { geoelement = geolist.get(i); if (geoelement instanceof GeoPoint) { ((GeoPoint) geoelement).getInhomCoords(xy); xlist[i] = xy[0]; ylist[i] = xy[1]; } else { error = true; xlist[i] = 0.0d; ylist[i] = 0.0d; } // if } // for all points }// getPoints() // Make M with 1,x,x^2,... , and Y with y1,y2,.. for all datapoints private final void makeMatrixArrays(int degree) { marray = new double[size][degree + 1]; yarray = new double[size][1]; // column matrix for (int i = 0; i < size; i++) { // Y: yarray[i][0] = ylist[i]; // M: for (int j = 0; j < (degree + 1); j++) { marray[i][j] = Math.pow(xlist[i], j); } // for j (all degrees =columns in marray) } // for i (all datapoints = rows in marray, cols in yarray) }// makeMatrixArrays() // / --- DEBUG --- /// !!! Remember to comment out calls before distribution // !!! }// class RegressionMath