import umontreal.iro.lecuyer.probdist.*; import umontreal.iro.lecuyer.util.*; /** * This class implements the algorithm NI2b wich uses the derivative * and simple integration. */ public class NI2b extends NortaInitDisc { private int m; /* Number of subintervals for the integration = max. number of iterations (also named m in the paper, paragraph "Method NI2" of section 3).*/ private double delta; /* Small positive parameter to make sure that rho_m is not too close to 1 or -1; (also named delta in the paper, paragraph "Method NI2" of section 3)*/ /** * Constructor of the class NI2b with the target rank correlation rX, * the two discrete marginals dist1 and dist2, * the parameter for truncation tr, the specific parameters m and * delta which correspond to m and delta in the paper * (paragraph "Method NI2" of section 3). */ public NI2b(double rX, DiscreteDistributionInt dist1, DiscreteDistributionInt dist2, double tr, int m, double delta) { super(rX, dist1, dist2, tr); this.m = m; this.delta = delta; computeParams(); } /** * Computes and returns the correlation rho_Z using the algorithm NI2b. */ public double computeCorr () { // x, y and c oefficients used for quadratic interpolation double[] x = new double[3]; double[] y = new double[3]; double[] c = new double[3]; double xtemp = 0.0, temp = 0.0; /* Values of rho and the recursive quantity I_k, given in paragraph "Method NI2" of section 3 in the paper,2 iterations before the last one. */ double xold, iold; /* Values of rho and I_k at one iteration before the last one.*/ double xnew, inew; // rho and I_k at the last iteration. double dold, dmid, dnew; /* Values of the derivative function g' at points xold, xmid (xold+h) and xnew. They correspond to g'(rho_0+2kh-2h), g'(rho_0+2kh-h) and g'(rho_0+2kh) in the formula of I_k in the paper (paragraph "Method NI2" of section 3). */ double h = (1 - delta) / (2 * m); /* Step size for the integration-grid spacing (2*h). It corresponds to h given in the third paragraph of section 4 in the paper */ double b = 0.0; // The returned solution. double rho1 = 0.0; // The initial guess. double intg1 = mu1 * mu2; // Computes g_r(rho1). double gr = rX * sd1 * sd2 + mu1 * mu2; /* Target value; integ(rho) = gr is equivalent to rho = the solution. */ // Pre-compute constants double hd3 = h / 3; double h2 = 2 * h; if (intg1 == gr) return rho1; if (0 < rX && rX < 1) { // Do search between 0 and rho_m=1-delta xold = rho1; dold = deriv(xold); iold = intg1; for (int i = 1; i <= m ;i++) { // Begin the search dmid = deriv(xold + h); xnew = xold + h2; dnew = deriv(xnew); inew = iold + hd3 * (dold + 4 * dmid + dnew); if (inew >= gr) { // The root is in current bracketing interval // Compute the parameters of quadratic interpolation x[0] = xtemp; x[1] = xold; x[2] = xnew; y[0] = temp; y[1] = iold; y[2] = inew; Misc.interpol(2, x, y, c); b = (c[2] * (xtemp + xold) - c[1] + Math.sqrt((c[1] - c[2] * (xtemp + xold)) * (c[1] - c[2] * (xtemp + xold)) - 4 * c[2] * (c[0] - c[1] * xtemp + c[2] * xtemp * xold - gr))) / (2 * c[2]); return b; } xtemp = xold; temp = iold; xold = xnew; dold = dnew ; iold = inew; } // Integration up to 1-delta did not bracket root // return 1-delta/2 ( = midpoint of current bracketing interval) b = 1 - delta / 2; } if ( -1 < rX && rX < 0) { // Do search between rho_m=-1+delta and 0 xold = rho1; dold = deriv(xold); iold = intg1; for (int i = 1; i <= m ;i++) { // Begin the search dmid = deriv(xold - h); xnew = xold - h2; dnew = deriv(xnew); inew = iold - hd3 * (dold + 4 * dmid + dnew); if (inew <= gr) { // The root is in current bracketing interval // Compute the parameters of quadratic interpolation x[0] = xnew; x[1] = xold; x[2] = xtemp; y[0] = inew; y[1] = iold; y[2] = temp; Misc.interpol(2, x, y, c); b = (c[2] * (xnew + xold) - c[1] + Math.sqrt((c[1] - c[2] * (xnew + xold)) * (c[1] - c[2] * (xnew + xold)) - 4 * c[2] * (c[0] - c[1] * xnew + c[2] * xnew * xold - gr))) / (2 * c[2]); return b; } xtemp = xold; temp = iold; xold = xnew; dold = dnew ; iold = inew; } // Integration up to -1+delta did not bracket root b = -1 + delta / 2; // return 1-delta/2 // ( = midpoint of current bracketing interval ) } return b; } /** * To display the inputs */ public String toString() { String desc = super.toString(); desc += "m : " + m + "\n"; desc += "delta : " + delta + "\n"; return desc; } }