Compute probit (inverse cdf) for ExpGamma distributions.

@param p Probability
@param k k, alpha aka. "shape" parameter
@param theta Theta = 1.0/Beta aka. "scaling" parameter
@param shift Shift parameter
@return Probit for ExpGamma distribution
186552
public static double quantile(double p, double k, double theta, double shift) {
    return FastMath.log(GammaDistribution.quantile(p, k, 1)) / theta + shift;
  }
868044
def cal_p_bm3(v, k, p_ref=0.0):
    
    vvr = v / k[0]
    p = (p_ref - 0.5 * (3. * k[1] - 5. * p_ref) * (1. - vvr**(-2. / 3.)) +
         9. / 8. * k[1] * (k[2] - 4. + 35. / 9. * p_ref / k[1]) *
         (1. - vvr**(-2. / 3.))**2.) * vvr**(-5. / 3.)
    return p
203797
def _kp(a, b):
    
    if a.shape != b.shape or a.shape[-1] != 1:
        raise(ValueError)
    N = a.shape[0]
    # take the outer product over the last two axes, then reshape:
    return np.einsum('ijk,ilk->ijkl', a, b).reshape(N, -1, 1)
1345849
protected void createBulge( int x1 , double p , boolean byAngle ) {
        double a11 = diag[x1];
        double a22 = diag[x1+1];
        double a12 = off[x1];
        double a23 = off[x1+1];

        if( byAngle ) {
            c = Math.cos(p);
            s = Math.sin(p);

            c2 = c*c;
            s2 = s*s;
            cs = c*s;
        } else {
            computeRotation(a11-p, a12);
        }

        // multiply the rotator on the top left.
        diag[x1]   = c2*a11 + 2.0*cs*a12 + s2*a22;
        diag[x1+1] = c2*a22 - 2.0*cs*a12 + s2*a11;
        off[x1]    = a12*(c2-s2) + cs*(a22 - a11);
        off[x1+1]  = c*a23;
        bulge = s*a23;

        if( Q != null )
            updateQ(x1,x1+1,c,s);
    }
591041
public static double klDivergence(double[] p, double[] q) {
        if (p.length != q.length) {
            throw new IllegalStateException("The length of p and q must be the same.");
        }
        double delta = 1e-8;
        if (!Multinomials.isMultinomial(p, delta)) {
            throw new IllegalStateException("p is not a multinomial");
        }
        if (!Multinomials.isMultinomial(q, delta)) {
            throw new IllegalStateException("q is not a multinomial");
        }
        
        double kl = 0.0;
        for (int i=0; i<p.length; i++) {
            if (p[i] == 0.0 || q[i] == 0.0) {
                continue;
            }
            kl += p[i] * FastMath.log(p[i] / q[i]);
        }
        return kl;
    }
572475
def pre_factor_kkt(Q, G, A):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    try:
        Q_LU = btrifact_hack(Q)
    except:
        raise RuntimeError("""
qpth Error: Cannot perform LU factorization on Q.
Please make sure that your Q matrix is PSD and has
a non-zero diagonal.
""")

    # S = [ A Q^{-1} A^T        A Q^{-1} G^T          ]
    #     [ G Q^{-1} A^T        G Q^{-1} G^T + D^{-1} ]
    #
    # We compute a partial LU decomposition of the S matrix
    # that can be completed once D^{-1} is known.
    # See the 'Block LU factorization' part of our website
    # for more details.

    G_invQ_GT = torch.bmm(G, G.transpose(1, 2).btrisolve(*Q_LU))
    R = G_invQ_GT.clone()
    S_LU_pivots = torch.IntTensor(range(1, 1 + neq + nineq)).unsqueeze(0) \
        .repeat(nBatch, 1).type_as(Q).int()
    if neq > 0:
        invQ_AT = A.transpose(1, 2).btrisolve(*Q_LU)
        A_invQ_AT = torch.bmm(A, invQ_AT)
        G_invQ_AT = torch.bmm(G, invQ_AT)

        LU_A_invQ_AT = btrifact_hack(A_invQ_AT)
        P_A_invQ_AT, L_A_invQ_AT, U_A_invQ_AT = torch.btriunpack(*LU_A_invQ_AT)
        P_A_invQ_AT = P_A_invQ_AT.type_as(A_invQ_AT)

        S_LU_11 = LU_A_invQ_AT[0]
        U_A_invQ_AT_inv = (P_A_invQ_AT.bmm(L_A_invQ_AT)
                           ).btrisolve(*LU_A_invQ_AT)
        S_LU_21 = G_invQ_AT.bmm(U_A_invQ_AT_inv)
        T = G_invQ_AT.transpose(1, 2).btrisolve(*LU_A_invQ_AT)
        S_LU_12 = U_A_invQ_AT.bmm(T)
        S_LU_22 = torch.zeros(nBatch, nineq, nineq).type_as(Q)
        S_LU_data = torch.cat((torch.cat((S_LU_11, S_LU_12), 2),
                               torch.cat((S_LU_21, S_LU_22), 2)),
                              1)
        S_LU_pivots[:, :neq] = LU_A_invQ_AT[1]

        R -= G_invQ_AT.bmm(T)
    else:
        S_LU_data = torch.zeros(nBatch, nineq, nineq).type_as(Q)

    S_LU = [S_LU_data, S_LU_pivots]
    return Q_LU, S_LU, R
1150709
public static double pdf(double x, double k, double theta, double shift) {
    if(x == Double.POSITIVE_INFINITY || x == Double.NEGATIVE_INFINITY) {
      return 0.;
    }
    if(x != x) {
      return Double.NaN;
    }
    x = (x - shift) * theta;
    final double ex = FastMath.exp(x);
    return ex < Double.POSITIVE_INFINITY ? FastMath.exp(k * x - ex - GammaDistribution.logGamma(k)) * theta : 0.;
  }
554801
public double randInverseCDF() {
        final double a0 = 2.50662823884;
        final double a1 = -18.61500062529;
        final double a2 = 41.39119773534;
        final double a3 = -25.44106049637;
        final double b0 = -8.47351093090;
        final double b1 = 23.08336743743;
        final double b2 = -21.06224101826;
        final double b3 = 3.13082909833;
        final double c0 = 0.3374754822726147;
        final double c1 = 0.9761690190917186;
        final double c2 = 0.1607979714918209;
        final double c3 = 0.0276438810333863;
        final double c4 = 0.0038405729373609;
        final double c5 = 0.0003951896511919;
        final double c6 = 0.0000321767881768;
        final double c7 = 0.0000002888167364;
        final double c8 = 0.0000003960315187;

        double y, r, x;

        double u = Math.random();
        while (u == 0.0) {
            u = Math.random();
        }
        
        y = u - 0.5;

        if (Math.abs(y) < 0.42) {
            r = y * y;
            x = y * (((a3 * r + a2) * r + a1) * r + a0)
                    / ((((b3 * r + b2) * r + b1) * r + b0) * r + 1);

        } else {
            r = u;
            if (y > 0) {
                r = 1 - u;
            }
            r = Math.log(-Math.log(r));
            x = c0 + r * (c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r * (c6 + r * (c7 + r * c8)))))));
            if (y < 0) {
                x = -(x);
            }
        }
        
        return mu + sigma * x;
    }
397403
private static long specialNPrimeB(final long numSamplesI, final double p, final double delta) {
    final double q, numSamplesF, oneMinusDelta;
    double tot, curTerm;
    long m;
    assertTrue(numSamplesI >= 1);
    assertTrue((0.0 < p) && (p < 1.0));
    assertTrue((0.0 < delta) && (delta < 1.0));
    q = 1.0 - p;
    oneMinusDelta = 1.0 - delta;
    numSamplesF = numSamplesI;
    curTerm = Math.pow(p, numSamplesF);  // curTerm = posteriorProbability (k, k, p)
    assertTrue(curTerm > 1e-100); // sanity check for non-use of logarithms
    tot = curTerm;
    m = numSamplesI;
    while (tot < oneMinusDelta) {
      curTerm = (curTerm * q * (m)) / ((m + 1) - numSamplesI);
      tot += curTerm;
      m += 1;
    }
    return (m); // don't need to back up
  }
1038308
function(p,phi,theta) {
        var r = 0;
        r += Math.pow(mathUtils.sin(this.m[0] * theta), this.m[1]);
        r += Math.pow(mathUtils.cos(this.m[2] * theta), this.m[3]);
        r += Math.pow(mathUtils.sin(this.m[4] * phi), this.m[5]);
        r += Math.pow(mathUtils.cos(this.m[6] * phi), this.m[7]);

        var sinTheta = mathUtils.sin(theta);
        p.x = r * sinTheta * mathUtils.cos(phi);
        p.y = r * mathUtils.cos(theta);
        p.z = r * sinTheta * mathUtils.sin(phi);
        return p;
    }
79347
def K(self):
        """"""
        K1 = np.float_power(pi, 0.5 * self.n_min * self.n_min)
        K1 /= (
            np.float_power(2, 0.5 * self.n_min * self._n_max)
            * self._mgamma(0.5 * self._n_max, self.n_min)
            * self._mgamma(0.5 * self.n_min, self.n_min)
        )

        K2 = np.float_power(
            2, self.alpha * self.size + 0.5 * self.size * (self.size + 1)
        )
        for i in xrange(self.size):
            K2 *= gamma(self.alpha + i + 1)
        return K1 * K2
592334
private double pVap( double T, double P ) {
        double A = 6.1121 * (1.0007 + 3.46E-6 * P);
        double b = 17.502;
        double c = 240.97;
        double e = A * exp(b * T / (c + T));
        return e;
    }
517035
public static double bernoulliCdf(int k, double p) {
        if(p<0) {
            throw new IllegalArgumentException("The probability p can't be negative.");
        }
        
        double probabilitySum=0.0;
        if(k<0) {
            
        }
        else if(k<1) { //aka k==0
            probabilitySum=(1-p);
        }
        else { //k>=1 aka k==1
            probabilitySum=1.0;
        }
        
        return probabilitySum;
    }
866588
public static double uniformCdf(int k, int n) {
        if(k<0 || n<1) {
            throw new IllegalArgumentException("All the parameters must be positive and n larger than 1.");
        }
        k = Math.min(k, n);
        
        double probabilitySum = k*uniform(n);
        
        return probabilitySum;
    }
12168
public static double binomial(int k, double p, int n) {
        if(k<0 ||  p<0 || n<1) {
            throw new IllegalArgumentException("All the parameters must be positive and n larger than 1.");
        }
        
        k = Math.min(k, n); 
        
        /*
        //Slow and can't handle large numbers
        $probability=StatsUtilities::combination($n,$k)*pow($p,$k)*pow(1-$p,$n-$k);
        */

        //fast and can handle large numbers
        //Cdf(k)-Cdf(k-1)
        double probability = approxBinomialCdf(k,p,n); 
        if(k>0) {
            probability -= approxBinomialCdf(k-1,p,n);
        }
        
        return probability;
    }
156191
@Reference(title = "Fast and accurate computation of binomial probabilities", //
      authors = "C. Loader", booktitle = "", //
      url = "http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf", //
      bibkey = "web/Loader00")
  public static double pmf(double x, int n, double p) {
    // Invalid values
    if(x < 0 || x > n) {
      return 0.;
    }
    // Extreme probabilities
    if(p <= 0.) {
      return x == 0 ? 1. : 0.;
    }
    if(p >= 1.) {
      return x == n ? 1. : 0.;
    }
    final double q = 1 - p;
    // FIXME: check for x to be integer, return 0 otherwise?

    // Extreme values of x
    if(x == 0) {
      if(p < .1) {
        return FastMath.exp(-devianceTerm(n, n * q) - n * p);
      }
      else {
        return FastMath.exp(n * FastMath.log(q));
      }
    }
    if(x == n) {
      if(p > .9) {
        return FastMath.exp(-devianceTerm(n, n * p) - n * q);
      }
      else {
        return FastMath.exp(n * FastMath.log(p));
      }
    }
    final double lc = stirlingError(n) - stirlingError(x) - stirlingError(n - x) - devianceTerm(x, n * p) - devianceTerm(n - x, n * q);
    final double f = (MathUtil.TWOPI * x * (n - x)) / n;
    return FastMath.exp(lc) / FastMath.sqrt(f);
  }
225817
@Override
    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
        /*
         * IMPLEMENTATION NOTES
         * --------------------
         * Where applicable, use is made of the one-sided Chebyshev inequality
         * to bracket the root. This inequality states that
         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
         * mu: mean, sig: standard deviation. Equivalently
         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
         * F(mu + k * sig) >= k^2 / (1 + k^2).
         *
         * For k = sqrt(p / (1 - p)), we find
         * F(mu + k * sig) >= p,
         * and (mu + k * sig) is an upper-bound for the root.
         *
         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
         * F(mu - k * sig) <= 1 / (1 + k^2).
         *
         * For k = sqrt((1 - p) / p), we find
         * F(mu - k * sig) <= p,
         * and (mu - k * sig) is a lower-bound for the root.
         *
         * In cases where the Chebyshev inequality does not apply, geometric
         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
         * the root.
         */
        if (p < 0.0 || p > 1.0) {
            throw new OutOfRangeException(p, 0, 1);
        }

        double lowerBound = getSupportLowerBound();
        if (p == 0.0) {
            return lowerBound;
        }

        double upperBound = getSupportUpperBound();
        if (p == 1.0) {
            return upperBound;
        }

        final double mu = getNumericalMean();
        final double sig = FastMath.sqrt(getNumericalVariance());
        final boolean chebyshevApplies;
        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || Double.isInfinite(sig) || Double.isNaN(sig));

        if (lowerBound == Double.NEGATIVE_INFINITY) {
            if (chebyshevApplies) {
                lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
            } else {
                lowerBound = -1.0;
                while (cumulativeProbability(lowerBound) >= p) {
                    lowerBound *= 2.0;
                }
            }
        }

        if (upperBound == Double.POSITIVE_INFINITY) {
            if (chebyshevApplies) {
                upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
            } else {
                upperBound = 1.0;
                while (cumulativeProbability(upperBound) < p) {
                    upperBound *= 2.0;
                }
            }
        }

        final UnivariateFunction toSolve = new UnivariateFunction() {

            public double value(final double x) {
                return cumulativeProbability(x) - p;
            }
        };

        double x = UnivariateSolverUtils.solve(toSolve, lowerBound, upperBound, getSolverAbsoluteAccuracy());

        if (!isSupportConnected()) {
            /* Test for plateau. */
            final double dx = getSolverAbsoluteAccuracy();
            if (x - dx >= getSupportLowerBound()) {
                double px = cumulativeProbability(x);
                if (cumulativeProbability(x - dx) == px) {
                    upperBound = x;
                    while (upperBound - lowerBound > dx) {
                        final double midPoint = 0.5 * (lowerBound + upperBound);
                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint;
                        } else {
                            upperBound = midPoint;
                        }
                    }
                    return upperBound;
                }
            }
        }
        return x;
    }