final Point2D.Double p = obs[i];
problem.addPoint(p.x, p.y);
}
// Direct solution (using simple regression).
final RealVector regress = new ArrayRealVector(problem.solve(), false);
// Dummy optimizer (to compute the chi-square).
final LeastSquaresProblem lsp = builder(problem).build();
final double[] init = { slope, offset };
// Get chi-square of the best parameters set for the given set of
// observations.
final double bestChi2N = getChi2N(lsp, regress);
final RealVector sigma = lsp.evaluate(regress).getSigma(1e-14);
// Monte-Carlo (generates a grid of parameters).
final int mcRepeat = MONTE_CARLO_RUNS;
final int gridSize = (int) FastMath.sqrt(mcRepeat);
// Parameters found for each of Monte-Carlo run.
// Index 0 = slope
// Index 1 = offset
// Index 2 = normalized chi2
final List<double[]> paramsAndChi2 = new ArrayList<double[]>(gridSize * gridSize);
final double slopeRange = 10 * sigma.getEntry(0);
final double offsetRange = 10 * sigma.getEntry(1);
final double minSlope = slope - 0.5 * slopeRange;
final double minOffset = offset - 0.5 * offsetRange;
final double deltaSlope = slopeRange/ gridSize;
final double deltaOffset = offsetRange / gridSize;
for (int i = 0; i < gridSize; i++) {
final double s = minSlope + i * deltaSlope;
for (int j = 0; j < gridSize; j++) {
final double o = minOffset + j * deltaOffset;
final double chi2N = getChi2N(lsp,
new ArrayRealVector(new double[] {s, o}, false));
paramsAndChi2.add(new double[] {s, o, chi2N});
}
}
// Output (for use with "gnuplot").
// Some info.
// For plotting separately sets of parameters that have a large chi2.
final double chi2NPlusOne = bestChi2N + 1;
int numLarger = 0;
final String lineFmt = "%+.10e %+.10e %.8e\n";
// Point with smallest chi-square.
System.out.printf(lineFmt, regress.getEntry(0), regress.getEntry(1), bestChi2N);
System.out.println(); // Empty line.
// Points within the confidence interval.
for (double[] d : paramsAndChi2) {
if (d[2] <= chi2NPlusOne) {
System.out.printf(lineFmt, d[0], d[1], d[2]);
}
}
System.out.println(); // Empty line.
// Points outside the confidence interval.
for (double[] d : paramsAndChi2) {
if (d[2] > chi2NPlusOne) {
++numLarger;
System.out.printf(lineFmt, d[0], d[1], d[2]);
}
}
System.out.println(); // Empty line.
System.out.println("# sigma=" + sigma.toString());
System.out.println("# " + numLarger + " sets filtered out");
}