@Override @Test(expected=SingularMatrixException.class) public void testNonInvertible() { /* * Overrides the method from parent class, since the default singularity * threshold (1e-14) does not trigger the expected exception. */ LinearProblem problem = new LinearProblem(new double[][] { { 1, 2, -3 }, { 2, 1, 3 }, { -3, 0, -9 } }, new double[] { 1, 1, 1 }); AbstractLeastSquaresOptimizer optimizer = createOptimizer(); PointVectorValuePair optimum = optimizer.optimize(100, problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 }); Assert.assertTrue(FastMath.sqrt(problem.target.length) * optimizer.getRMS() > 0.6); optimizer.computeCovariances(optimum.getPoint(), 1.5e-14); }
@Override @Test(expected=SingularMatrixException.class) public void testNonInvertible() { /* * Overrides the method from parent class, since the default singularity * threshold (1e-14) does not trigger the expected exception. */ LinearProblem problem = new LinearProblem(new double[][] { { 1, 2, -3 }, { 2, 1, 3 }, { -3, 0, -9 } }, new double[] { 1, 1, 1 }); AbstractLeastSquaresOptimizer optimizer = createOptimizer(); PointVectorValuePair optimum = optimizer.optimize(new MaxEval(100), problem.getModelFunction(), problem.getModelFunctionJacobian(), problem.getTarget(), new Weight(new double[] { 1, 1, 1 }), new InitialGuess(new double[] { 0, 0, 0 })); Assert.assertTrue(FastMath.sqrt(optimizer.getTargetSize()) * optimizer.getRMS() > 0.6); optimizer.computeCovariances(optimum.getPoint(), 1.5e-14); }
/** * This method updates the Kalman-Filter with the last known state/measurement of the observed value. * * @param xMeasured time step of measurement. * It is measured in full, i.e. as {@link Long}, seconds since midnight, January 1, 1970 UTC. * @param yMeasured Current measurement. * @return True if the update was successful. */ public boolean update(long xMeasured, double yMeasured) { boolean success = false; try { // Call predict(0), if no prediction was made since the last update // Reason: The Kalman-Filter needs a predict-update(correct)-cycle. if (!predictedSinceUpdate && lastUpdate != Constants.NO_PREDICTION) { predict(0); } filter.correct(new double[] {xMeasured, 0, yMeasured, 0 }); // When an older value is updated/corrected the attributes 'lastUpdated' and 'lastUpdate' do not change. if (lastUpdated < xMeasured) { lastUpdated = xMeasured; lastUpdate = yMeasured; } success = true; predictedSinceUpdate = false; lastMemUpdated(); } catch (NullArgumentException | DimensionMismatchException | SingularMatrixException e) { LogManager.getLogger(Kalman.class).error(e.getMessage(), e); } return success; }
@Override @Test(expected=SingularMatrixException.class) public void testNonInvertible() { /* * Overrides the method from parent class, since the default singularity * threshold (1e-14) does not trigger the expected exception. */ LinearProblem problem = new LinearProblem(new double[][] { { 1, 2, -3 }, { 2, 1, 3 }, { -3, 0, -9 } }, new double[] { 1, 1, 1 }); AbstractLeastSquaresOptimizer optimizer = createOptimizer(); optimizer.optimize(100, problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 }); Assert.assertTrue(FastMath.sqrt(problem.target.length) * optimizer.getRMS() > 0.6); optimizer.getCovariances(1.5e-14); }
private void outputResult() { output = new StringBuilder(); printDoubleArray("expectedPosition: ", expectedPosition); printDoubleArray("linear calculatedPosition: ", linearCalculatedPosition.toArray()); printDoubleArray("non-linear calculatedPosition: ", nonLinearOptimum.getPoint().toArray()); output.append("numberOfIterations: ").append(nonLinearOptimum.getIterations()).append("\n"); output.append("numberOfEvaluations: ").append(nonLinearOptimum.getEvaluations()).append("\n"); try { RealVector standardDeviation = nonLinearOptimum.getSigma(0); printDoubleArray("standardDeviation: ", standardDeviation.toArray()); output.append("Norm of deviation: ").append(standardDeviation.getNorm()).append("\n"); RealMatrix covarianceMatrix = nonLinearOptimum.getCovariances(0); output.append("covarianceMatrix: ").append(covarianceMatrix).append("\n"); } catch (SingularMatrixException e) { System.err.println(e.getMessage()); } System.out.println(output.toString()); }
/** * @param A a square matrix. * @return the inverse of A or null if A is non-square or singular. */ public static double[][] inverse(final double[][] A) { RealMatrix M = MatrixUtils.createRealMatrix(A); if (!M.isSquare()) return null; else { double[][] Ai = null; try { RealMatrix Mi = MatrixUtils.inverse(M); //new LUDecomposition(M).getSolver().getInverse(); Ai = Mi.getData(); } catch (SingularMatrixException e) {} return Ai; } }
public static double[] solve(final double[][] A, double[] b) { RealMatrix AA = MatrixUtils.createRealMatrix(A); RealVector bb = MatrixUtils.createRealVector(b); DecompositionSolver solver = new LUDecomposition(AA).getSolver(); double[] x = null; try { x = solver.solve(bb).toArray(); } catch (SingularMatrixException e) {} return x; }
@Nonnull public static RealMatrix inverse(@Nonnull final RealMatrix m, final boolean exact) throws SingularMatrixException { LUDecomposition LU = new LUDecomposition(m); DecompositionSolver solver = LU.getSolver(); final RealMatrix inv; if (exact || solver.isNonSingular()) { inv = solver.getInverse(); } else { SingularValueDecomposition SVD = new SingularValueDecomposition(m); inv = SVD.getSolver().getInverse(); } return inv; }
/** * L = A x R * * @return a matrix A that minimizes A x R - L */ @Nonnull public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R, final boolean exact) throws SingularMatrixException { LUDecomposition LU = new LUDecomposition(R); DecompositionSolver solver = LU.getSolver(); final RealMatrix A; if (exact || solver.isNonSingular()) { A = LU.getSolver().solve(L); } else { SingularValueDecomposition SVD = new SingularValueDecomposition(R); A = SVD.getSolver().solve(L); } return A; }
@Override @Test public void testNonInvertible() { try{ /* * Overrides the method from parent class, since the default singularity * threshold (1e-14) does not trigger the expected exception. */ LinearProblem problem = new LinearProblem(new double[][] { { 1, 2, -3 }, { 2, 1, 3 }, { -3, 0, -9 } }, new double[] { 1, 1, 1 }); final Optimum optimum = optimizer.optimize( problem.getBuilder().maxIterations(20).build()); //TODO check that it is a bad fit? Why the extra conditions? Assert.assertTrue(FastMath.sqrt(problem.getTarget().length) * optimum.getRMS() > 0.6); optimum.getCovariances(1.5e-14); fail(optimizer); }catch (SingularMatrixException e){ //expected } }
@Test(expected=SingularMatrixException.class) public void testQRSingular() { final RealMatrix a = MatrixUtils.createRealMatrix(new double[][] { { 1, 6, 4 }, { 2, 4, -1 }, { -1, 2, 5 } }); final RealVector b = new ArrayRealVector(new double[]{ 5, 6, 1 }); new QRDecomposition(a, 1.0e-15).getSolver().solve(b); }
/** Guess a start line using the last four results. * @param n number of cached results available * @param target target ground point * @return guessed start line */ private double guessStartLine(final int n, final Vector3D target) { try { // assume a linear model of the form: l = ax + by + cz + d final RealMatrix m = new Array2DRowRealMatrix(n, 4); final RealVector v = new ArrayRealVector(n); for (int i = 0; i < n; ++i) { m.setEntry(i, 0, cachedResults[i].getTarget().getX()); m.setEntry(i, 1, cachedResults[i].getTarget().getY()); m.setEntry(i, 2, cachedResults[i].getTarget().getZ()); m.setEntry(i, 3, 1.0); v.setEntry(i, cachedResults[i].getLine()); } final DecompositionSolver solver = new QRDecomposition(m, Precision.SAFE_MIN).getSolver(); final RealVector coefficients = solver.solve(v); // apply the linear model return target.getX() * coefficients.getEntry(0) + target.getY() * coefficients.getEntry(1) + target.getZ() * coefficients.getEntry(2) + coefficients.getEntry(3); } catch (SingularMatrixException sme) { // the points are not independent return Double.NaN; } }
/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p> * <p>Filters targets to populations where more than the minProportion lie in a single segment.</p> * <p>Returns null if no pass filtering. Please note that in these cases, * in the rest of this class, we use this to assume that a GMM is not a good model.</p> * * @param segments -- segments with segment mean in log2 copy ratio space * @param targets -- targets with a log2 copy ratio estimate * @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used * in the evaluation * @param numComponents -- number of components to use in the GMM. Usually, this is 2. * @return never {@code null}. Fitting result with indications whether it converged or was even attempted. */ private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments, final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents){ // For each target in a segment that contains enough targets, normalize the difference against the segment mean // and collapse the filtered targets into the copy ratio estimates. final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion); if (filteredTargetsSegDiff.size() < numComponents) { return new MixtureMultivariateNormalFitResult(null, false, false); } // Assume that Apache Commons wants data points in the first dimension. // Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API. final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2]; // Convert the filtered targets into 2d array (even if second dimension is length 1). The second dimension is // uncorrelated Gaussian. This is only to get around funny API in Apache Commons, which will throw an // exception if the length of the second dimension is < 2 final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED)); final NormalDistribution nd = new NormalDistribution(rng, 0, .1); for (int i = 0; i < filteredTargetsSegDiff.size(); i++) { filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i); filteredTargetsSegDiff2d[i][1] = nd.sample(); } final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents); final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d); try { multivariateNormalMixtureExpectationMaximization.fit(estimateEM0); } catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) { // We are done, we cannot make a fitting. We should return a result that we attempted a fitting, but it // did not converge. Include the model as it was when the exception was thrown. return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true); } return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true); }
private void updateInverseCovariance() { try { inverseCov = new LUDecomposition(cov).getSolver().getInverse(); } catch (SingularMatrixException e) { singularCovariances.inc(); inverseCov = new SingularValueDecomposition(cov).getSolver().getInverse(); } }
/** * Correct the current state estimate with an actual measurement. * * @param z * the measurement vector * @throws NullArgumentException * if the measurement vector is {@code null} * @throws DimensionMismatchException * if the dimension of the measurement vector does not fit * @throws SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) throws NullArgumentException, DimensionMismatchException, SingularMatrixException { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) - * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance) .multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // invert S // as the error covariance matrix is a symmetric positive // semi-definite matrix, we can use the cholesky decomposition DecompositionSolver solver = new CholeskyDecomposition(s).getSolver(); RealMatrix invertedS = solver.getInverse(); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); }
@Test(expected=SingularMatrixException.class) public void testNonInvertible() { QRDecomposition qr = new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)); final RealMatrix inv = qr.getSolver().getInverse(); }
@Override @Test(expected=SingularMatrixException.class) public void testNonInvertible() { /* * Overrides the method from parent class, since the default singularity * threshold (1e-14) does not trigger the expected exception. */ LinearProblem problem = new LinearProblem(new double[][] { { 1, 2, -3 }, { 2, 1, 3 }, { -3, 0, -9 } }, new double[] { 1, 1, 1 }); final LevenbergMarquardtOptimizer optimizer = createOptimizer() .withMaxEvaluations(100) .withMaxIterations(20) .withModelAndJacobian(problem.getModelFunction(), problem.getModelFunctionJacobian()) .withTarget(problem.getTarget()) .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) .withStartPoint(new double[] { 0, 0, 0 }); final double[] optimum = optimizer.optimize().getPoint(); Assert.assertTrue(FastMath.sqrt(optimizer.getTarget().length) * optimizer.computeRMS(optimum) > 0.6); optimizer.computeCovariances(optimum, 1.5e-14); }
/** * Returns a sample matrix from this distribution. * @return sampled matrix. */ public RealMatrix sample() { for (int i=0; i<100; i++) { try { RealMatrix A = sampleWishart(); RealMatrix result = new LUDecomposition(A).getSolver().getInverse(); LOG.log(Level.FINE, "Cov = {0}", result); return result; } catch (SingularMatrixException ex) { LOG.finer("Discarding singular matrix generated by the wishart distribution."); } } throw new RuntimeException("Unable to generate inverse wishart samples!"); }
/** * Correct the current state estimate with an actual measurement. * * @param z * the measurement vector * @throws NullArgumentException * if the measurement vector is {@code null} * @throws DimensionMismatchException * if the dimension of the measurement vector does not fit * @throws SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) throws NullArgumentException, DimensionMismatchException, SingularMatrixException { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance) .multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix .operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 // instead of calculating the inverse of S we can rearrange the formula, // and then solve the linear equation A x X = B with A = S', X = K' and // B = (H * P)' // K(k) * S = P(k)- * H' // S' * K(k)' = H * P(k)-' RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver() .solve(measurementMatrix.multiply(errorCovariance.transpose())) .transpose(); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain .getRowDimension()); errorCovariance = identity.subtract( kalmanGain.multiply(measurementMatrix)).multiply( errorCovariance); }
@Nonnull public static RealMatrix inverse(@Nonnull final RealMatrix m) throws SingularMatrixException { return inverse(m, true); }
/** * L = A x R * * @return a matrix A that minimizes A x R - L */ @Nonnull public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R) throws SingularMatrixException { return solve(L, R, true); }
/** * Creates a multivariate normal distribution with the given mean vector and * covariance matrix. * <br/> * The number of dimensions is equal to the length of the mean vector * and to the number of rows and columns of the covariance matrix. * It is frequently written as "p" in formulae. * * @param rng Random Number Generator. * @param means Vector of means. * @param covariances Covariance matrix. * @throws DimensionMismatchException if the arrays length are * inconsistent. * @throws SingularMatrixException if the eigenvalue decomposition cannot * be performed on the provided covariance matrix. * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is * negative. */ public MultivariateNormalDistribution(RandomGenerator rng, final double[] means, final double[][] covariances) throws SingularMatrixException, DimensionMismatchException, NonPositiveDefiniteMatrixException { super(rng, means.length); final int dim = means.length; if (covariances.length != dim) { throw new DimensionMismatchException(covariances.length, dim); } for (int i = 0; i < dim; i++) { if (dim != covariances[i].length) { throw new DimensionMismatchException(covariances[i].length, dim); } } this.means = MathArrays.copyOf(means); covarianceMatrix = new Array2DRowRealMatrix(covariances); // Covariance matrix eigen decomposition. final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix); // Compute and store the inverse. covarianceMatrixInverse = covMatDec.getSolver().getInverse(); // Compute and store the determinant. covarianceMatrixDeterminant = covMatDec.getDeterminant(); // Eigenvalues of the covariance matrix. final double[] covMatEigenvalues = covMatDec.getRealEigenvalues(); for (int i = 0; i < covMatEigenvalues.length; i++) { if (covMatEigenvalues[i] < 0) { throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0); } } // Matrix where each column is an eigenvector of the covariance matrix. final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim); for (int v = 0; v < dim; v++) { final double[] evec = covMatDec.getEigenvector(v).toArray(); covMatEigenvectors.setColumn(v, evec); } final RealMatrix tmpMatrix = covMatEigenvectors.transpose(); // Scale each eigenvector by the square root of its eigenvalue. for (int row = 0; row < dim; row++) { final double factor = FastMath.sqrt(covMatEigenvalues[row]); for (int col = 0; col < dim; col++) { tmpMatrix.multiplyEntry(row, col, factor); } } samplingMatrix = covMatEigenvectors.multiply(tmpMatrix); }
/** * Correct the current state estimate with an actual measurement. * * @param z * the measurement vector * @throws NullArgumentException * if the measurement vector is {@code null} * @throws DimensionMismatchException * if the dimension of the measurement vector does not fit * @throws SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) throws NullArgumentException, DimensionMismatchException, SingularMatrixException { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance) .multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 // instead of calculating the inverse of S we can rearrange the formula, // and then solve the linear equation A x X = B with A = S', X = K' and B = (H * P)' // K(k) * S = P(k)- * H' // S' * K(k)' = H * P(k)-' RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver() .solve(measurementMatrix.multiply(errorCovariance.transpose())) .transpose(); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); }
/** * Suppose y is the target to be checked, x is the impact to be tested, u is the universe * of all variables (including y, x and all other time series) and L is the given lag. The * question is then:<br> * <b>Does x Granger-cause y in the given universe u with the given lag L?</b><br> * <br> * The H0-test is then the following equation, which'll be solved via OLS: * <pre> * y_t = (y_t-1 ... y_y-L u0_t-1 ... u0_t-L ... un_t-1 ... un_t-L) * (beta_1 ... beta_(n*L))^T + (epsilon_1 ... epsilon_(n*L))^T * ... * y_t-d = (y_t-d-1 ... y_y-d-L u0_t-d-1 ... u0_t-d-L ... un_t-d-1 ... un_t-d-L) * (beta_1 ... beta_(n*L))^T + (epsilon_1 ... epsilon_(n*L))^T * which translates to * Y = U(without x) * BETA + EPSILON * </pre> * The H1-test is then the exact same as H0 but U contains also x. * * This translates to the following variable names: * <ul> * <li>Y => strippedY</li> * <li>X(without x) => laggedH0</li> * <li>X => laggedH1</li> * </ul> * * @param laggedH0Rows The number of rows in U without x (the laggedH0 matrix). This is the * d index in the above defined equation (which can also be called "number of observations"). * @param strippedY The Y vector that is cut off with the lag size. * @param laggedH0 The universe matrix without the x variable whose influence is to be tested. * The Matrix is also lagged with the given lag L. * @param variablesH0 The number of variables in the universe without the variable x. * @param laggedH1 The universe matrix with the x variable. The Matrix is also lagged with the * given lag L. * @param variablesH1 The number of variables in the universe with the variable x (naturally, * this should be variablesH0+1). * @return The computed GrangerCausalIndicator object, in case of success. * Null - If anything went wrong (e.g., the OLS can't compute the parameters due to a * SingularMatrixException) */ protected GrangerCausalIndicator performGranger(int laggedH0Rows, double[] strippedY, double[][] laggedH0, int variablesH0, double[][] laggedH1, int variablesH1) { try { OLSMultipleLinearRegression h0 = new OLSMultipleLinearRegression(); OLSMultipleLinearRegression h1 = new OLSMultipleLinearRegression(); h0.newSampleData(strippedY, laggedH0); // print(laggedH0); h1.newSampleData(strippedY, laggedH1); // print(laggedH1); double rs0[] = h0.estimateResiduals(); double rs1[] = h1.estimateResiduals(); // System.out.print("b = {"); // for(int i=0; i<b0.length;i++){System.out.print(FORMATTER.format(b0[i])+", ");}; // System.out.println(); // System.out.print("residuals = {"); // for(int i=0; i<rs0.length;i++){System.out.print(FORMATTER.format(rs0[i])+", ");}; // System.out.println(); double RSS0 = sqrSum(rs0); double RSS1 = sqrSum(rs1); int nbrParametersModel0 = this.lagSize * variablesH0; int nbrParametersModel1 = this.lagSize * variablesH1; // (RSS1-RSS2) / (p2-p1) // f = ------------------------- // RSS2 / (n-p2) // double ftest = ((RSS0 - RSS1)/this.lagSize) / (RSS1 / (laggedH0Rows - 2*this.lagSize - 1)); double ftest = ((RSS0 - RSS1) / (nbrParametersModel1-nbrParametersModel0)) / (RSS1 / (laggedH0Rows-nbrParametersModel1)); // System.out.println(RSS0 + " " + RSS1); // System.out.println("F-test " + ftest); FDistribution fDist = new FDistribution(this.lagSize, laggedH0Rows-2*this.lagSize-1); double pValue = 1.0 - fDist.cumulativeProbability(ftest); // System.out.println("P-value " + pValue); return new GrangerCausalIndicator(pValue, CRITICAL_VALUE, this.lagSize); } catch(SingularMatrixException smex) { smex.printStackTrace(); return null; } }
