import umontreal.iro.lecuyer.probdist.*; /** * This class implements the algorithm NI3 (protected Newton-Raphson method). * The root is found with accuracy tolerance. */ public class NI3 extends NortaInitDisc { private double tolerance; /* Desired accuracy for the root-finder algorithm (epsilon in paragraph "Method NI3" of section 3 in paper).*/ /** * Constructor of the class NI3 with the target rank correlation rX, * the two discrete marginals dist1 and dist2, * the parameter for truncation tr and the specific parameter * tolerance which corresponds to epsilon in the paper * (paragraph "Method NI3" of section 3). */ public NI3 (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 NI3. */ public double computeCorr () { final double ITMAX = 100; // Maximum number of iterations. double xl, xh; /* Left and right endpoints of the root bracket at all iterations. */ double b = 0.0; // The returned solution. double f, df; // Function and its derivative evaluations. double dx; // Correction term. /** f, df, dx correspond to f(rho_k), f'(rho_k) and f(rho_k)/f'(rho_k), respectively, in the paper (paragraph "Method NI3"of section 3). */ double dxold; // The root correction at one iteration before the last one. double temp;// The root at one iteration before the last one. double ccc = rX * sd1 * sd2 + mu1 * mu2; // Precompute constant. if (rX == 0.0) return 0.0; if (rX > 0.0) { // Orient the search xl = 0.0; xh = 1.0; } else { xl = -1.0; xh = 0.0; } b = 2 * Math.sin (Math.PI * rX / 6); // Initial guess dxold = xh - xl; dx = dxold; f = integ (b) - ccc; df = deriv (b); for (int i = 1; i <= ITMAX; i++) { // Begin the search if ((((b - xh) * df - f) * ((b - xl) * df - f) > 0.0) || (Math.abs (2.0 * f) > Math.abs (dxold * df))) { // Do bisection if solution is out of range // or not decreasing fast enough dxold = dx; dx = 0.5 * (xh - xl); b = xl + dx; if (xl == b) // Change in root is negligible. return b; // Accept this root } else { dxold = dx; dx = f / df; temp = b; b -= dx; if (temp == b) return b; } if (Math.abs (dx) < tolerance) return b; // Convergence check f = integ (b) - ccc; df = deriv (b); if (f < 0.0) xl = b; // Maintain the brackets on the root else xh = b; } return b; } /** * To display the inputs. */ public String toString() { String desc = super.toString(); desc += "tolerance : " + tolerance + "\n"; return desc; } }