How do I transform an Inverse CDF coded for a standard deviation of one to a different STDDEV?

Question

I am doing rejection sampling as part of an algorithm for numerical integration using an inverse cumulative distribution function (CDF).

I found two possible implementations:

1. In C and Java (and others, but not C#, so I would need to translate):

http://home.online.no/~pjacklam/notes/invnorm/

1. In C#, there is StatisticFormula.InverseNormalDistribution (in System.Windows.Forms.DataVisualization.Charting). Given Microsoft's track record with a bad NORMINV implementation in Excel many years ago, I am leery.

Both of these functions assume a mean of zero (which is what I will be using) and a standard deviation of one. How do I transform the input and/or output of this function so that the standard deviation is different from one?

I know with the Gaussian distribution:

f(x,mean,sd) = (1/(sd*sqrt(2*pi)))exp(-0.5((x-mean)/sd)^2)

So that if I have g(x) = f(x,0,1), that f(x,mean,sd) = (1/sd)*g(x/sd).

The transformation for the Gaussian to a different Std Dev is easy. Can I do something similar for the inverse CDF?

Answer 1

When all else fails, try MS Excel. I computed the inverse CDF using NORMINV (which lets you supply a standard deviation) and compared it to the value returned using NORMSINV (which assumes a std dev of one). The happy result was that you can just multiply the result obtained using a standard deviation of one by the desired standard deviation.

NORMSINV(x) * SD = NORMINV(x,0,SD)

BTW: I tried to use the Microsoft library, but without success. You can only call the functions in the context of a Chart, which is an unacceptable overhead. I ended up porting Jeremy Lea's C implementation to C#. It is cited the link I give in the question, http://home.online.no/~pjacklam/notes/invnorm/

I either made a mistake in my porting, or the C version has an error, because I get good accuracy for negative inputs, but when x > 6 accuracy degrades rapidly and is gone by x = 7.

Answer 2

Because of the broken link with the implementation in the answer by Paul Chernoch I had to look for my own answer. I found an implementation of R's version of the inverse CDF with a parametrized mean and standard deviation here: Standard Normal Distribution z-value function in C#

To prevent the risk of another broken link I am listing my copy of the method below.

        /// <summary>
    /// Quantile function (Inverse CDF) for the normal distribution.
    /// </summary>
    /// <param name="p">Probability.</param>
    /// <param name="mu">Mean of normal distribution.</param>
    /// <param name="sigma">Standard deviation of normal distribution.</param>
    /// <param name="isLowerTail">If true, probability is P[X <= x], otherwise P[X > x].</param>
    /// <param name="isLogValues">If true, probabilities are given as log(p).</param>
    /// <returns>P[X <= x] where x ~ N(mu,sigma^2)</returns>
    /// <remarks>See https://svn.r-project.org/R/trunk/src/nmath/qnorm.c </remarks>
    public static double QNorm(double p, double mu, double sigma, bool isLowerTail, bool isLogValues)
    {
        if (double.IsNaN(p) || double.IsNaN(mu) || double.IsNaN(sigma)) return (p + mu + sigma);
        double ans;
        bool isBoundaryCase = R_Q_P01_boundaries(p, double.NegativeInfinity, double.PositiveInfinity, isLowerTail, isLogValues, out ans);
        if (isBoundaryCase) return (ans);
        if (sigma < 0) return (double.NaN);
        if (sigma == 0) return (mu);

        double p_ = R_DT_qIv(p, isLowerTail, isLogValues);
        double q = p_ - 0.5;
        double r, val;

        if (Math.Abs(q) <= 0.425)  // 0.075 <= p <= 0.925
        {
            r = .180625 - q * q;
            val = q * (((((((r * 2509.0809287301226727 +
                       33430.575583588128105) * r + 67265.770927008700853) * r +
                     45921.953931549871457) * r + 13731.693765509461125) * r +
                   1971.5909503065514427) * r + 133.14166789178437745) * r +
                 3.387132872796366608)
            / (((((((r * 5226.495278852854561 +
                     28729.085735721942674) * r + 39307.89580009271061) * r +
                   21213.794301586595867) * r + 5394.1960214247511077) * r +
                 687.1870074920579083) * r + 42.313330701600911252) * r + 1.0);
        }
        else
        {
            r = q > 0 ? R_DT_CIv(p, isLowerTail, isLogValues) : p_;
            r = Math.Sqrt(-((isLogValues && ((isLowerTail && q <= 0) || (!isLowerTail && q > 0))) ? p : Math.Log(r)));

            if (r <= 5)              // <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11
            {
                r -= 1.6;
                val = (((((((r * 7.7454501427834140764e-4 +
                        .0227238449892691845833) * r + .24178072517745061177) *
                      r + 1.27045825245236838258) * r +
                     3.64784832476320460504) * r + 5.7694972214606914055) *
                   r + 4.6303378461565452959) * r +
                  1.42343711074968357734)
                 / (((((((r *
                          1.05075007164441684324e-9 + 5.475938084995344946e-4) *
                         r + .0151986665636164571966) * r +
                        .14810397642748007459) * r + .68976733498510000455) *
                      r + 1.6763848301838038494) * r +
                     2.05319162663775882187) * r + 1.0);
            }
            else                     // very close to  0 or 1 
            {
                r -= 5.0;
                val = (((((((r * 2.01033439929228813265e-7 +
                        2.71155556874348757815e-5) * r +
                       .0012426609473880784386) * r + .026532189526576123093) *
                     r + .29656057182850489123) * r +
                    1.7848265399172913358) * r + 5.4637849111641143699) *
                  r + 6.6579046435011037772)
                 / (((((((r *
                          2.04426310338993978564e-15 + 1.4215117583164458887e-7) *
                         r + 1.8463183175100546818e-5) * r +
                        7.868691311456132591e-4) * r + .0148753612908506148525)
                      * r + .13692988092273580531) * r +
                     .59983220655588793769) * r + 1.0);
            }
            if (q < 0.0) val = -val;
        }

        return (mu + sigma * val);
    }

    private static bool R_Q_P01_boundaries(double p, double left, double right, bool isLowerTail, bool isLogValues, out double ans)
    {
        if (isLogValues)
        {
            if (p > 0.0)
            {
                ans = double.NaN;
                return (true);
            }
            if (p == 0.0)
            {
                ans = isLowerTail ? right : left;
                return (true);
            }
            if (p == double.NegativeInfinity)
            {
                ans = isLowerTail ? left : right;
                return (true);
            }
        }
        else
        {
            if (p < 0.0 || p > 1.0)
            {
                ans = double.NaN;
                return (true);
            }
            if (p == 0.0)
            {
                ans = isLowerTail ? left : right;
                return (true);
            }
            if (p == 1.0)
            {
                ans = isLowerTail ? right : left;
                return (true);
            }
        }
        ans = double.NaN;
        return (false);
    }

    private static double R_DT_qIv(double p, bool isLowerTail, bool isLogValues)
    {
        return (isLogValues ? (isLowerTail ? Math.Exp(p) : -ExpM1(p)) : R_D_Lval(p, isLowerTail));
    }

    private static double R_DT_CIv(double p, bool isLowerTail, bool isLogValues)
    {
        return (isLogValues ? (isLowerTail ? -ExpM1(p) : Math.Exp(p)) : R_D_Cval(p, isLowerTail));
    }

    private static double R_D_Lval(double p, bool isLowerTail)
    {
        return isLowerTail ? p : 0.5 - p + 0.5;
    }

    private static double R_D_Cval(double p, bool isLowerTail)
    {
        return isLowerTail ? 0.5 - p + 0.5 : p;
    }
    private static double ExpM1(double x)
    {
        if (Math.Abs(x) < 1e-5)
            return x + 0.5 * x * x;
        else
            return Math.Exp(x) - 1.0;
    }

I ran some first very rudimentary tests of the results of this method against the implementation of NORM.INV in Excel 2016. The accuracy appears to be identical up to a precision of 10^-7.