Skew logistic distribution

Question

I am trying to program a skew logistic distribution function to return the PDF, the CDF, and (later) the inverse CDF. Although there are several forms of skew logistic, this is for perhaps the most widely known, a variation of the skew normal distribution introduced by Azzalini and covered by Saralees Nadarajah. Aside from a formal paper on Springer, it's covered in "Focus on African Diaspora Mathematics By Toka Diagana". A Google search will turn up the book and the skew logistic portion starts at the bottom of page 41.

I am comfortable with the PDF function (included below), but I am struggling with the CDF, primarily because I am not quite familiar with the math.

Here's the skew logistic CDF taken from the formal paper on Springer (Nadarajah, S. The skew logistic distribution, AStA Adv Stat Anal (2009) 93: 187–203):

$$F(x)=\left\{\begin{array}{lr}1-2\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}\binom{-1}{j}\binom{-2}{k}\frac{1}{1+\lambda j +k}\exp\left[-\frac{(1+\lambda j + k)x}{\beta}\right]&\text{if } x>0,\\ 2\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}\binom{-1}{j}\binom{-2}{k}\frac{1}{1+\lambda+\lambda j +k}\exp\left[\frac{(1+\lambda+\lambda j+k)x}{\beta}\right]&\text{if }x<0.\end{array}\right.$$

I understand everything except this portion:

$$\sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \binom{-1}{j}\binom{-2}{k}$$

My current approach is to simply start with the value, step up and down by 0.0001 to sum the PDF in each direction, and then divide the area below the value by the total area.

However, this strikes me as extremely inefficient and inaccurate past five or six decimal places. I have to believe a better function could be written based off the CDF formula Nadarajah gives.

Once I get the SkewLogisticCDF done properly, I'd like to also have a SkewLogisticCDFInverse as well.

// Skew Logistic PDF
double SkewLogisticPDF(double value, double location, double scale, double shape)
{

double x = exp((-(value - location)) / scale);
double y = exp((-(value - location)) * shape / scale);

x = 2 * x / (scale * pow(1 + x, 2) * (1 + y));

return x;
}


// Skew Logistic CDF
double SkewLogisticCDF(double value, double location, double scale, double shape)
{
const double step = 0.00001;
const double precision = 0.0000001;

double areaBelowValue = 0;
double areaAboveValue = 0;

double currentValue = 0;
double currentPDF = 0;

currentValue = value;
do
{
    currentPDF = SkewLogisticPDF(currentValue, location, scale, shape);
    areaBelowValue += currentPDF;

    currentValue -= step;

} while (currentPDF >= precision);

currentValue = value;
do
{
    currentPDF = SkewLogisticPDF(currentValue, location, scale, shape);
    areaAboveValue += currentPDF;

    currentValue += step;

} while (currentPDF >= precision);

double x = areaBelowValue / (areaBelowValue + areaAboveValue);

return x;
}

Answer 1

OK, after some digging I realized I was actually familiar with the notation below after all:

$$\sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \binom{-1}{j}\binom{-2}{k}$$

However I've only seen it in use for combinations where negative numbers are not allowed. I initially thought it must represent some other mathematical function, but I found a generalized version that easily allows for negative numbers. For those interested, it's covered under the Wikipedia article on "Binomial Coefficient" in the section on "Generalization and connection to the binomial series."

Below is the C++ function I came up with. It essentially starts the CDF off with zero and then adds to the CDF by increasing j and k from zero toward infinity until the change falls below the threshold given in the precision constant. Here the precision is set to 0.00000001, giving probably at least six decimal places of accuracy, possibly seven, which is good enough for my purposes. Decreasing the precision further would return even more accurate results at the cost of speed, but that might be necessary for a different application.

Some quick notes:

-- The net result of this function: $$\sum_{j=0}^{\infty}\binom{-1}{j}$$ is to simply return 1 if j is even and -1 if j is odd. I've substituted it with a much simpler jMult = ((j % 2) != 1) ? 1 : -1;

-- The net result of this function: $$\sum_{k=0}^{\infty}\binom{-2}{k}$$ is to simply return k + 1 if k is even and -(k + 1) if k is odd. I've substituted it with a much simpler kMult = -((k % 2) * 2 - 1)*(k + 1);

-- Although Nadarajah shows two different formulas in his paper, one for value greater than zero and one for value less than zero, I've written this function so that it can accept a positive or negative value and return the correct number with only one formula.

-- Also, Nadarajah doesn't mention this but I realized that when the shape is negative, there are j and k combinations that end up dividing by zero. However, the correct result is given by finding x using the negative of the shape and value and returning 1 - x as the result.

-- I accomplished both of the above by introducing the distance, distanceMultiplier, shapeTerm, and extraShapeTerm variables and adjusting Nadarajah's formula to use those.

-- I also use kLastSignificant and kPrevious variables to track how far out on k the previous iteration of j went before it fell below the precision threshold. I found that for larger values of j, the initial k would be below the precision threshold but later values of k would be above the threshold. Here I wanted to ensure the function captured those later values by going out at least as far on k as the previous iteration of j did.

[EDIT] -- I've updated the function to use Nadarajah's correct formula for the case of when the value is exactly zero. The psi function comes from http://people.sc.fsu.edu/~jburkardt/cpp_src/cdflib/cdflib.cpp

There are still two areas of improvement I can identify:

-- The function is extremely fast except when the value is close or equal to zero. For example, if the value is 0.0001, the j and k tend to stretch out for extremely long periods before they fall below the precision threshold.

-- For very large values, the j and k initially start below the precision threshold and the function results in a one or a zero, which are not possible in the logistic distribution. To fix this, I've added a portion at the end of the function to simply return either the precision or (1 - precision) in those cases.

In the end, the skew logistic CDF really needs an approximation formula that would easily address each of those issues and provide an even faster function at the cost of a little accuracy. It's how most other functions like this (i.e., the Normal CDF) are calculated and it would make the inverse much easier to calculate (and more practical to program).

// Skew Logistic CDF
double SkewLogisticCDF(double value, double location, double scale, double shape)
{
const double precision = 0.00000001;

double x = 0;

if (shape == 0)
{
    x = 1 / (1 + exp(-(value - location) / scale));
}
else
{

    // Evaluate the distance and the shape
    // If the shape is negative, then distance *= -1 and shape *= -1
    double shapeTerm = (shape > 0) ? shape : -shape;
    double distance = (shape > 0) ? (value - location) : (location - value);

    double distanceTermMultuplier = (distance < 0) ? 1 : -1;
    double extraShapeTerm = (distance < 0) ? shapeTerm : 0;

    int j = 0;
    int k = 0;

    double result;
    int kPrevious = 0;
    int kLastSignificant = 0;

    // distance != 0
    if (distance != 0)
    {
        // j
        do
        {

            // Returns 1 if j is even, -1 if j is odd
            // This is the ( -1 )
            //             (  j ) of the formula
            // Can also achieve using (j % 2)*-2 + 1;
            // Went with conditional operator to make quicker
            int jMult = ((j % 2) != 1) ? 1 : -1;

            k = 0;
            kLastSignificant = 0;
            // k
            do
            {

                // Returns (k + 1) if k is even, -(k + 1) if k is odd
                // This is the ( -2 )
                //             (  k ) of the formula
                // Can also achieve using -((k % 2) * 2 - 1)*(k + 1);
                // Went with conditional operator to make quicker
                int kMult = -((k % 2) * 2 - 1)*(k + 1);

                result = jMult * kMult * 1 / (1 + extraShapeTerm + shapeTerm * j + k) * exp((1 + extraShapeTerm + shapeTerm * j + k) * distanceTermMultuplier * distance / scale);

                x += result;

                kLastSignificant = (abs(result) >= precision) ? k : kLastSignificant;
                ++k;

            } while (abs(result) >= precision || k <= kPrevious);

            kPrevious = kLastSignificant;

            ++j;

        } while (kLastSignificant > 0);

        x *= 2;
        x = (distance < 0) ? x : 1 - x;

    }

    // distance == 0
    else
    {
        // j
        do
        {

            // Returns 1 if j is even, -1 if j is odd
            // This is the ( -1 )
            //             (  j ) of the formula
            // Can also achieve using (j % 2)*-2 + 1;
            // Went with conditional operator to make quicker
            int jMult = ((j % 2) != 1) ? 1 : -1;

            double psiAlpha = 1 + shapeTerm + shapeTerm * j;
            double psiOne = (1 + psiAlpha) / 2;
            double psiTwo = psiAlpha / 2;

            psiOne = psi(psiOne);
            psiTwo = psi(psiTwo);

            result = jMult * (0.5 - shapeTerm * (1 + j) * 0.5 * (psiOne - psiTwo));

            x += result;

            ++j;

        } while (abs(result) >= precision);

        x *= 2;

    }
}

// Reverse result if shape is negative
x = (shape > 0) ? x : 1 - x;

x = (x == 0) ? precision : x;
x = (x == 1) ? (1 - precision) : x;

return x;

}