import umontreal.iro.lecuyer.probdist.*; import umontreal.iro.lecuyer.util.*; /** * This class implements the algorithm NI1 based on Brent method for * root-finding. The root is found with accuracy = tolerance. */ public class NI1 extends NortaInitDisc { private double tolerance; /* Desired accuracy for the root-finder algorithm (epsilon in paragraph "Method NI1" of section 3 in the paper).*/ /** * Constructor of the class NI1 with the target rank correlation rX, * the two discrete marginals dist1 and dist2, * the parameter for truncation tr and the accuracy tolerance, * which corresponds to epsilon in the paper * (paragraph "Method NI1" of section 3). */ public NI1 (double rX, DiscreteDistributionInt dist1, DiscreteDistributionInt dist2, double tr, double tolerance) { super(rX, dist1, dist2, tr); this.tolerance = tolerance; computeParams(); } /** * Computes and returns the correlation rho_Z using the algorithm NI1. */ public double computeCorr () { final double ITMAX = 100; // Maximum number of iterations final double EPS = 1.0e-15; // Machine accuracy double b; /* Latest iterate and closest approximation of the root (the returned solution at the end). */ double a; // Previous iterate (previous value of b). double c; /* Previous or older iterate so that f(b) and f(c) have opposite signs (may coincide with a). */ double fa, fb, fc; /* Evaluations of the function f at points a, b and c. */ double tolerance1; /* Final tolerance wich involves the machine accuracy and the initial desired tolerance (epsilon). */ /** a, b, c, fa, fb, fc, and tolerance correspond to a, b, c, f(a) f(b), f(c) and epsilon, respectively, in the paper (paragraph "Method NI1" of section 3). */ double x1, x2; // Left and right endpoints of initial search interval double pp, q, rrr, s; /* Parameters to compute inverse quadratic interpolation. */ double xm; // Bisection point. double min1, min2; /* Criterias to check whether inverse quadratic interpolation can be performed or not. */ double e = 0.0, d = 0.0; /* Parameters to specify bounding interval and whether bisection or inverse quadratic interpolation was performed at one iteration before the last one. */ // Precompute constants. double cc = rX * sd1 * sd2; double ccc = cc + mu1 * mu2; if (rX == 0) return 0.0; if (rX > 0) { // Orient the search and initialize a, b, c x1 = 0.0; x2 = 1.0; a = x1; b = x2; c = x2; fa = - cc; fb = integ (b) - ccc; } else { x1 = -1.0; x2 = 0.0; a = x1; b = x2; c = x2; fa = integ (a) - ccc; fb = - cc; } fc = fb; for (int i = 1; i <= ITMAX; i++) { // Begin the search if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { // Rename a, b, c and adjust bounding interval d c = a; fc = fa; e = d = b - a; } if (Math.abs (fc) < Math.abs (fb)) { a = b; b = c; c = a; fa = fb; fb = fc; fc = fa; } // Convergence check tolerance1 = 2.0 * EPS * Math.abs (b) + 0.5 * tolerance; xm = 0.5 * (c - b); if (Math.abs (xm) <= tolerance1 || fb == 0.0) return b; if (Math.abs (e) >= tolerance1 && Math.abs (fa) > Math.abs (fb)) { s = fb / fa; // Attempt inverse quadratic interpolation if (a == c) { pp = 2.0 * xm * s; q = 1.0 - s; } else { q = fa / fc; rrr = fb / fc; pp = s * (2.0 * xm * q * (q - rrr) - (b - a) * (rrr - 1.0)); q = (q - 1.0) * (rrr - 1.0) * (s - 1.0); } if (pp > 0.0) q = -q; // Check whether in bounds pp = Math.abs (pp); min1 = 3.0 * xm * q - Math.abs (tolerance1 * q); min2 = Math.abs (e * q); if (2.0 * pp < (min1 < min2 ? min1 : min2)) { e = d; // Accept interpolation d = pp / q; } else { // Interpolation failed, use bisection d = xm; e = d; } } else { // Bounds decreasing too slowly, use bisection d = xm; e = d; } a = b; // a becomes the best trial fa = fb; if (Math.abs (d) > tolerance1) b += d; // Evaluate the new trial root else { if (xm > 0) b += Math.abs (tolerance1); else b += -Math.abs (tolerance1); } fb = integ (b) - ccc; } return b; } /** * To display the inputs. */ public String toString() { String desc = super.toString(); desc += "tolerance : " + tolerance + "\n"; return desc; } }