Forum Navigation:

magazine

FORUMS > Programming and Software Forum < refresh >
Topic Title: c++ source code for Microsoft Excel's NORMINV() function?
Created On Wed May 24, 06 08:08 AM
Topic View:

View thread in raw text format


sjoo
Member

Posts: 109
Joined: Mar 2003

Wed May 24, 06 08:08 AM
User is offline View users profile

Hello everyone,

I 'm wondering if the c++ source code that is equivalent to Excel's NORMINV() function is available anywhere?

As you know, NORMINV(RAND(),Mu,Sigma) is used to simulate a value from the normal distribution in Microsoft Excel.

I'm hoping to exactly replicate the NORMINV() function.

Thanks very much.

Best regards.

sjoo


Edited: Wed May 24, 06 at 08:20 AM by sjoo
 
Reply
   
Quote
   
Top
   
Bottom
     



MattF
Senior Member

Posts: 892
Joined: Mar 2003

Wed May 24, 06 10:07 AM
User is offline View users profile

Don't bother because Excel's implementation is known to be garbage.

Read this: Numerics in Excel

Also here is a relevant paper by Graeme West - who I believe is also on wilmott but I've forgotten his forum name
AccurateCumNorm
 
Reply
   
Quote
   
Top
   
Bottom
     



DominicConnor
Senior Member

Posts: 11470
Joined: Jul 2002

Thu May 25, 06 02:48 PM
User is offline View users profile

I'm interested why you would want to exactly replicate it ?


-------------------------
Random Walkers: Meetup for 87 Quants in a Pub
 
Reply
   
Quote
   
Top
   
Bottom
     



Rufus
Member

Posts: 188
Joined: Jan 2002

Thu May 25, 06 04:19 PM
User is offline


Are you

(i) trying to generate normally distributed numbers;
(ii) reverse engineer excel norminv().

If (i) then see numerical recipes in c chapter seven verse 2

If (ii) then you must have a good reason!
 
Reply
   
Quote
   
Top
   
Bottom
     



sjoo
Member

Posts: 109
Joined: Mar 2003

Thu May 25, 06 11:47 PM
User is offline View users profile

I 'm intending to generate a return path using NORMINV(RAND(), mean, volatility).

With the simulated returns, I'll make functions to calculate performance & risk measures.

And I got a source code yesterday. I post the code.


Best ragards,
sjoo

/*
 *     Compute the quantile function for the normal distribution.
 *
 *     For small to moderate probabilities, algorithm referenced
 *     below is used to obtain an initial approximation which is
 *     polished with a final Newton step.
 *
 *     For very large arguments, an algorithm of Wichura is used.
 *
 *  REFERENCE
 *
 *     Beasley, J. D. and S. G. Springer (1977).
 *     Algorithm AS 111: The percentage points of the normal distribution,
 *     Applied Statistics, 26, 118-121.
 *
 *      Wichura, M.J. (1988).
 *      Algorithm AS 241: The Percentage Points of the Normal Distribution.
 *      Applied Statistics, 37, 477-484.
 */
#include <stdio.h>
#include <math.h>

#define R_D_Cval(p) (lower_tail ? (1 - (p)) : (p))  /*  1 - p */

#define R_DT_CIv(p) (log_p ? (lower_tail ? -expm1(p) : exp(p)) \
                             : R_D_Cval(p))
#define R_D_Lval(p)   (lower_tail ? (p) : (1 - (p)))
#define ML_NEGINF   ((-1.0) / 0.0)
#define R_D__0  (log_p ? ML_NEGINF : 0.)        /* 0 */
#define R_D__1  (log_p ? 0. : 1.)           /* 1 */
#define ML_POSINF       (1.0 / 0.0)

#define R_DT_0  (lower_tail ? R_D__0 : R_D__1)      /* 0 */
#define R_DT_1  (lower_tail ? R_D__1 : R_D__0)      /* 1 */
#define R_Q_P01_check(p)            \
         if ((log_p  && p > 0) ||            \
                        (!log_p && (p < 0 || p > 1)) )      \
    return 0;
#define R_DT_qIv(p) (log_p ? (lower_tail ? exp(p) : - expm1(p)) \
                             : R_D_Lval(p))

#define DBL_EPSILON 0.0000001

double expm1(double x)
{
    double y, a = fabs(x);

    if (a < DBL_EPSILON) return x;
    if (a > 0.697) return exp(x) - 1;  /* negligible cancellation */

    if (a > 1e-8)
    y = exp(x) - 1;
    else /* Taylor expansion, more accurate in this range */
    y = (x / 2 + 1) * x;

    /* Newton step for solving   log(1 + y) = x   for y : */
    /* WARNING: does not work for y ~ -1: bug in 1.5.0 */
    y -= (1 + y) * (log1p (y) - x);
    return y;
}

                                                          
double qnorm(double p, double mu, double sigma)
{
    double p_, q, r, val;
     int lower_tail = 1;
     int log_p = 0;
    if (p == R_DT_0)     return ML_NEGINF;
    if (p == R_DT_1)     return ML_POSINF;
    R_Q_P01_check(p);

    if(sigma  < 0)     return 0;
    if(sigma == 0)     return mu;

    p_ = R_DT_qIv(p);/* real lower_tail prob. p */
    q = p_ - 0.5;

/*-- use AS 241 --- */
/* double ppnd16_(double *p, long *ifault)*/
/*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3

        Produces the normal deviate Z corresponding to a given lower
        tail area of P; Z is accurate to about 1 part in 10**16.

        (original fortran code used PARAMETER(..) for the coefficients
         and provided hash codes for checking them...)
*/
    if (fabs(q) <= .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.);
    }
    else { /* closer than 0.075 from {0,1} boundary */

     /* r = min(p, 1-p) < 0.075 */
     if (q > 0)
         r = R_DT_CIv(p);/* 1-p */
     else
         r = p_;/* = R_DT_Iv(p) ^=  p */

     r = sqrt(- ((log_p &&
               ((lower_tail && q <= 0) || (!lower_tail && q > 0))) ?
              p : /* else */ log(r)));
        /* r = sqrt(-log(r))  <==>  min(p, 1-p) = exp( - r^2 ) */

        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.);
        }
        else { /* very close to  0 or 1 */
            r += -5.;
            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.);
        }

     if(q < 0.0)
         val = -val;
        /* return (q >= 0.)? r : -r ;*/
    }
    return mu + sigma * val;
}


int main(void) {

     int i =0;
     double mu = 100;
     double sigma = 10;
     
     for(i=0;i<100;i++) {
          printf("qnorm:%10.10f  \n",qnorm((double)i/100,mu,sigma));
     }
     return 0;

}


Edited: Fri May 26, 06 at 12:13 AM by sjoo

 
Reply
   
Quote
   
Top
   
Bottom
     



fuez
Junior Member

Posts: 16
Joined: Sep 2004

Wed May 31, 06 03:46 PM
User is offline

is this a source using the method "make it look big" ?
i am using the following function to get normal-distributed random numbers (sorry for not telling the source...i don't know it anymore):

Another possibility is the Box-Mueller-Method (but have no c++ sorce code for that).

#include <cmath>
double norminv(double q) {
     if(q == .5)
          return 0;

     q = 1.0 - q;

     double p = (q > 0.0 && q < 0.5) ? q : (1.0 - q);
     double t = sqrt(log(1.0 / pow(p, 2.0)));

     double c0 = 2.515517;
     double c1 = 0.802853;
     double c2 = 0.010328;

     double d1 = 1.432788;
     double d2 = 0.189269;
     double d3 = 0.001308;

     double x = t - (c0 + c1 * t + c2 * pow(t, 2.0)) /
                    (1.0 + d1 * t + d2 * pow(t, 2.0) + d3 * pow(t, 3.0));
    
     if(q > .5)
          x *= -1.0;

     return x;
}

 
Reply
   
Quote
   
Top
   
Bottom
     



richmondjamesp
Junior Member

Posts: 2
Joined: Jul 2006

Wed Jul 19, 06 02:38 PM
User is offline

fuez,

Thanks a bunch for this contribution. I converted it to VB and it must be 10 or 20 times faster than the Excel 2003 NormInv() I was using (which I've read lots about but that's another story - you'd think the folks at MS could get this right).

One minor issue I had was that when a probability of 1 is passed into the function, a divide by zero error occurs.

I'm not positive if it was a conversion error on my part or if a probability of 1 is inherently undefined (based on the normal distribution). I am using the Excel Rnd without parameters to generate the probability.

Anyhow, the converted code with my change is below. I don't need huge precision for what I'm using this for (portfolio withdrawal monte carlo simulation) and the simulation seems to behave the same with the new function as it did with NormInv from Excel 2003 so I think it's ok.

Any thoughts on this assumption would be appreciated.


Thanks again for the post...

Jim

Function NormInv(Probability As Double, Mu As Double, Sigma As Double)


Dim x As Double
Dim p As Double
Dim c0 As Double, c1 As Double, c2 As Double
Dim d1 As Double, d2 As Double, d3 As Double
Dim t As Double
Dim q As Double

q = Probability
If (q = 0.5) Then
NormInv = Mu
Else
q = 1# - q

If ((q > 0) And (q < 0.5)) Then
p = q
Else
If (q = 1) Then
p = 1 - 0.9999999 ' JPR - attempt to fix divide by zero below, what is NormInv(1,x,y)?
Else
p = 1# - q
End If
End If

t = Sqr(Log(1# / (p * p)))

c0 = 2.515517
c1 = 0.802853
c2 = 0.010328

d1 = 1.432788
d2 = 0.189269
d3 = 0.001308

x = t - (c0 + c1 * t + c2 * (t * t)) / (1# + d1 * t + d2 * (t * t) + d3 * (t ^ 3))

If (q > 0.5) Then
x = -1# * x
End If
End If

NormInv = (x * Sigma) + Mu

End Function


Edited: Wed Jul 19, 06 at 02:59 PM by richmondjamesp
 
Reply
   
Quote
   
Top
   
Bottom
     



fuez
Junior Member

Posts: 16
Joined: Sep 2004

Thu Jul 20, 06 09:29 AM
User is offline

Thank _you_ for converting it to visual basic ;-)
To your divisonbyzero-error: the function works with a range of ]0,1[. Ever tried Excels Norminv() function using a probability of 0 or 1? Well, you don't get a Division byZero, but an other error message.

Ok and as we are now in Visual-Basic, i post the BoxMuller-Algorithm too (again smaller) ;-)
Don't blame me for using "GOTO"! But it is (much) faster than a loop.

Function DrawNormal()
'returns normal distributed random number
'Method: Box-Müller
Dim Rand1 As Double, Rand2 As Double, S1 As Double, S2 As Double, x1 As Double, x2 As Double
start:
Rand1 = 2 * Rnd - 1
Rand2 = 2 * Rnd - 1
S1 = Rand1 ^ 2 + Rand2 ^ 2
If S1 > 1 Then GoTo start
S2 = Sqr(-2 * Log(S1) / S1)
x1 = Rand1 * S2
'x2 = Rand2 * S2 ' Not necessary to calculate.
'return one or the other
DrawNormal = x1
End Function

 
Reply
   
Quote
   
Top
   
Bottom
     



richmondjamesp
Junior Member

Posts: 2
Joined: Jul 2006

Fri Jul 21, 06 02:44 AM
User is offline

Fuez,

Cool, now he comes up with a name-brand algorithm. Thanks.

I compiled it in and it seemed a tad slower (~15%) than the other one. I didn't count the operations and I'm not that worried about it at this point. The simulation results seemed decent for both.

I'm definitely jazzed by the speed I get in the simulation now that I don't have Excel's NormInv() weighing things down (I needed two in an inner most loop). I've had more completed sim runs in the past day than I had in a week or two before that.

Also, on the error I had before, it was a probability of 0 that was causing the crash, not 1 as I said before (I missed the q=1-p statement)

In any case, just to warn other green horns - Excel 2003 Rnd returns 0,1[ rather than ]0,1[

Jim
 
Reply
   
Quote
   
Top
   
Bottom
     



goablin
Junior Member

Posts: 9
Joined: May 2007

Fri Jun 22, 07 07:18 AM
User is offline

http://home.online.no/~pjacklam/notes/invnorm/index.html#The_algorithm

-------------------------
Understanding comes with experience...
 
Reply
   
Quote
   
Top
   
Bottom
     



rwinston
Member

Posts: 52
Joined: Feb 2007

Sat Jun 23, 07 10:38 PM
User is offline

Another source (I'm sure there are plenty) might be the R repository:

http://www.google.com/codesearch?hl=en&q=show:fVJtn3LL6pw:Romt0EWtO38:6lgBJTgcS4U&sa=N&ct=rd&cs_p=http://cran.r-project.org/src/base/R-2/R-2.2.1.tar.gz&cs_f=R-2.2.1/src/nmath/dnorm.c
 
Reply
   
Quote
   
Top
   
Bottom
     



tkh
Member

Posts: 39
Joined: Dec 2005

Thu Apr 23, 09 09:26 AM
User is offline


Very, very sorry to bring up this old,old thread, but c++ norminv managed to match with excel's norminv when line

y -= (1 + y) * (log1p (y) - x);

was corrected to

y -= (1 + y) * (log(1+ y) - x);

in sjoo's post.

 
Reply
   
Quote
   
Top
   
Bottom
     



FaridMoussaoui
Member

Posts: 171
Joined: Jun 2008

Fri Apr 24, 09 10:52 AM
User is offline View users profile

please, find some references about the accuracy in Excel.

As you can remark by reading these papers, the conclusion is "stay away" from Excel.

best regards,

F.

March 2008 articles of "Computational Statistics and Data Analysis"

http://dx.doi.org/10.1016/j.csda.2008.03.004
http://dx.doi.org/10.1016/j.csda.2008.03.005
http://dx.doi.org/10.1016/j.csda.2008.03.006
http://dx.doi.org/10.1016/j.csda.2008.03.007
http://dx.doi.org/10.1016/j.csda.2008.03.008


Edited: Fri Apr 24, 09 at 10:53 AM by FaridMoussaoui
 
Reply
   
Quote
   
Top
   
Bottom
     



Actuary321
Member

Posts: 62
Joined: Nov 2009

Fri Mar 19, 10 05:50 PM
User is offline

This one is accurate to within 10-9

http://home.online.no/~pjacklam/notes/invnorm/
 
Reply
   
Quote
   
Top
   
Bottom
     



Actuary321
Member

Posts: 62
Joined: Nov 2009

Fri Mar 19, 10 05:50 PM
User is offline

This one is accurate to within 10-9

http://home.online.no/~pjacklam/notes/invnorm/
 
Reply
   
Quote
   
Top
   
Bottom
     



Cuchulainn
Senior Member

Posts: 36897
Joined: Jul 2004

Sat Mar 20, 10 10:38 AM
User is offline View users profile

My own personal preference is to use a peer-reviewed de facto standard such as here

then ws functions

-------------------------
www.datasimfinancial.com
 
Reply
   
Quote
   
Top
   
Bottom
     



AVt
Senior Member

Posts: 1051
Joined: Dec 2001

Sat Mar 20, 10 08:36 PM
User is offline View users profile

As already said by others it depends on what you are looking for.

If you want code to produce normal distributed numbers, then you
can use Box-Muller or better, search at Google or Wikipedia for
algorithms doing just that.

If you want to reverse engineer Excel (versions before 2010): if I
remember correctly experiments for the cdfN then it seemed that
Excel uses one of the (oldish) formulae in Abramowitz & Stegun.
So it is natural to guess they used that source for the inverse
as well, formula 26.2.23 may be the candidate to be checked by
shooting enough test values.

Note that the latter intention will become no longer necessary
if you have Excel 2010 in mind.
 
Reply
   
Quote
   
Top
   
Bottom
     



Actuary321
Member

Posts: 62
Joined: Nov 2009

Mon Mar 22, 10 02:02 PM
User is offline

Does Boost have NormInv()?

I wasnt able to find this on the website.
 
Reply
   
Quote
   
Top
   
Bottom
     



AVt
Senior Member

Posts: 1051
Joined: Dec 2001

Mon Mar 22, 10 02:10 PM
User is offline View users profile

---> quantile
 
Reply
   
Quote
   
Top
   
Bottom
     

View thread in raw text format
FORUMS > Programming and Software Forum < refresh >

Forum Navigation:

© All material, including contents and design, copyright Wilmott Electronic Media Limited - FuseTalk 4.01 © 1999-2014 FuseTalk Inc. Terms & Conditions