protected double createSample()

in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/StableSampler.java [626:689]

• 16 lines of code
• 5 McCabe index (conditional complexity)

        protected double createSample(double phi, double w) {
            // Here we use the formula provided by Weron for the 1-parameterization.
            // Note: Adding back zeta creates the 0-parameterization defined in Nolan (1998):
            // X ~ S0_alpha(s,beta,u0) with s=1, u0=0 for a standard random variable.
            // As alpha -> 1 the translation zeta to create the stable deviate
            // in the 0-parameterization is increasingly large as tan(pi/2) -> infinity.
            // The max translation is approximately 1e16.
            // Without this translation the stable deviate is in the 1-parameterization
            // and the function is not continuous with respect to alpha.
            // Due to the large zeta when alpha -> 1 the number of bits of the output variable
            // are very low due to cancellation.

            // As alpha -> 0 or 2 then zeta -> 0 and cancellation is not relevant.
            // The formula can be modified for infinite terms to compute a result for extreme
            // deviates u and w when the CMS formula fails.

            // Note the following term is subject to floating point error:
            // xi = atan(-zeta) / alpha
            // alphaPhiXi = alpha * (phi + xi)
            // This is required: cos(phi - alphaPhiXi) > 0 => phi - alphaPhiXi in (-pi/2, pi/2).
            // Thus we compute atan(-zeta) and use it to compute two terms:
            // [1] alpha * (phi + xi) = alpha * (phi + atan(-zeta) / alpha) = alpha * phi + atan(-zeta)
            // [2] phi - alpha * (phi + xi) = phi - alpha * phi - atan(-zeta) = (1-alpha) * phi - atan(-zeta)

            // Compute terms
            // Either term can be infinite or 0. Certain parameters compute 0 * inf.
            // t1=inf occurs alpha -> 0.
            // t1=0 occurs when beta = tan(-alpha * phi) / tan(alpha * pi / 2).
            // t2=inf occurs when w -> 0 and alpha -> 0.
            // t2=0 occurs when alpha -> 0 and phi -> pi/2.
            // Detect zeros and return as zeta.

            // Note sin(alpha * phi + atanZeta) is zero when:
            // alpha * phi = -atan(-zeta)
            // tan(-alpha * phi) = -zeta
            //                   = beta * tan(alpha * pi / 2)
            // Since |phi| < pi/2 this requires beta to have an opposite sign to phi
            // and a magnitude < 1. This is possible and in this case avoid a possible
            // 0 / 0 by setting the result as if term t1=0 and the result is zeta.
            double t1 = Math.sin(meps1 * phi + atanZeta);
            if (t1 == 0) {
                return zeta;
            }
            // Since cos(phi) is in (0, 1] this term will not create a
            // large magnitude to create t1 = 0.
            t1 /= Math.pow(Math.cos(phi), inv1mEps);

            // Iff Math.cos(eps * phi - atanZeta) is zero then 0 / 0 can occur if w=0.
            // Iff Math.cos(eps * phi - atanZeta) is below zero then NaN will occur
            // in the power function. These cases are avoided by phi=(-pi/2, pi/2) and direct
            // use of arctan(-zeta).
            final double t2 = Math.pow(Math.cos(eps * phi - atanZeta) / w, epsDiv1mEps);
            if (t2 == 0) {
                return zeta;
            }

            final double x = t1 * t2 * scale + zeta;

            // Check the bounds. Applies when alpha < 1 and beta = +/-1.
            if (x <= lower) {
                return lower;
            }
            return x < upper ? x : upper;
        }