Normal distribution function

General
1

I need to replicate the NORMDIST(X) function from Excel.

Any ideas how to do this within Xojo?

2

You do a bit Googling and find the Wikipedia page at http://en.wikipedia.org/wiki/Normal_distribution, which contains the formula for the normal distribution function. It’s a simple power function with “e”. The power function is “^” in Xojo.

3

Thanks Beatrix - I was reading that wikipedia page overnight but its been 25 years since I took a stats class so am a bit rusty.

The formula is actually NORMSDIST(x) - not NORMDIST as I originally posted.

Microsoft explain the function here - http://support.microsoft.com/kb/827369
“NORMSDIST(z) returns the probability that the observed value of a standard normal random variable will be less than or equal to z. A standard normal random variable has mean 0 and standard deviation 1 (and also variance 1 because variance = standard deviation squared).”

So the outcome for 0 is 50% probability.
I could always generate an array independently and do a lookup.

4

I made a very small project with the fucntion you need. Here it is: https://dl.dropboxusercontent.com/u/3800071/CDF_Demo.xojo_binary_project

No documentation (it’s so small it is not required) and NO ERROR checking at all.

Julen.

5

This distribution is tricky. There are several approaches. Some are slow, and some lose accuracy on the tails. The following code is good to +/- 7 standard deviations for the area under the standard normal probability distribution function. I converted it from fortran, and it matches the Excel results.

[code] dim RT2PI as double
dim SPLIT as double
dim N0 as double
dim N1 as double
dim N2 as double
dim N3 as double
dim N4 as double
dim N5 as double
dim N6 as double
dim M0 as double
dim M1 as double
dim M2 as double
dim M3 as double
dim M4 as double
dim M5 as double
dim M6 as double
dim M7 as double
dim c as double
dim e as double
dim n as double
dim d as double
dim f as double
dim z as double

RT2PI = (4.0*acos(0.0))^0.5
SPLIT = 7.07106781186547

 N0 = 220.206867912376 
 N1 = 221.213596169931
 N2 = 112.079291497871
 N3 = 33.912866078383 
 N4 = 6.37396220353165 
 N5 = 0.700383064443688 
 N6 = 0.0352624965998911
 M0 = 440.413735824752 
 M1 = 793.826512519948 
 M2 = 637.333633378831 
 M3 = 296.564248779674 
 M4 = 86.7807322029461 
 M5 = 16.064177579207 
 M6 = 1.75566716318264 
 M7 = 0.0883883476483184

z = abs(x)
c = 0.0
’
if z<=37.0 then
e = exp(-zz/2.0)
if(z<SPLIT) then
n = (((((N6z + N5)*z + N4)*z + N3)*z + N2)*z + N1)z + N0
d = ((((((M7z + M6)*z + M5)*z + M4)*z + M3)*z + M2)z + M1)z + M0
c = en/d
else
f = z + 1.0/(z + 2.0/(z + 3.0/(z + 4.0/(z + 13.0/20.0))))
c = e/(RT2PIf)
end if
end if

if z>0.0 then
c = 1-c
end if

return c[/code]

6

Below is the inverse of the above,

[code] 'Purpose:
’
’ R4_NORMAL_01_CDF_INVERSE inverts the standard normal CDF.
’
’ Discussion:
’
’ The result is accurate to about 1 part in 10**7.
’
’ Licensing:
’
’ This code is distributed under the GNU LGPL license.
’
’ Modified:
’
’ 27 December 2004
’
’ Author:
’
’ Original FORTRAN77 version by Michael Wichura.
’ C version by John Burkardt.
’
’ Reference:
’
’ Michael Wichura,
’ The Percentage Points of the Normal Distribution,
’ Algorithm AS 241,
’ Applied Statistics,
’ Volume 37, Number 3, pages 477-484, 1988.
’
’ Parameters:
’
’ Input, float P, the value of the cumulative probability densitity function.
’ 0 < P < 1. If P is outside this range, an “infinite” result is returned.
’
’ Output, float R4_NORMAL_01_CDF_INVERSE, the normal deviate value with the
’ property that the probability of a standard normal deviate being less than or
’ equal to this value is P.

dim a(4) as double
dim b(4) as double
dim c(4) as double
dim d(3) as double
dim e(4) as double
dim f(3) as double
dim const2 as double
dim q as double
dim r as double
dim split1 as double
dim split2 as Double
dim value as double
dim const1 as double

a(1) = 3.3871327179
a(2) = 50.434271938
a(3) = 159.29113202
a(4) = 59.109374720

b(1) = 1.0
b(2) = 17.895169469
b(3) = 78.757757664
b(4) = 67.187563600

c(1) = 1.4234372777
c(2) = 2.7568153900
c(3) = 1.3067284816
c(4) = 0.17023821103

d(1) = 1.0
d(2) = 0.73700164250
d(3) = 0.12021132975

e(1) = 6.6579051150
e(2) = 3.0812263860
e(3) = 0.42868294337
e(4) = 0.017337203997

f(1) = 1.0
f(2) = 0.24197894225
f(3) = 0.012258202635

const1 = 0.180625
const2 = 1.6
split1 = 0.425
split2 = 5.0

if ( p <= 0.0 ) then
    value = - 999999999
   return value

end if

if ( 1.0 <= p ) then
  value = - 999999999
   return value

end if

q = p - 0.5

if (abs( q ) <= split1 ) then
     r = const1 - q * q
    value = q * app.R4POLY_VALUE(4, a(), r) / app.R4POLY_VALUE(4, b(), r)
elseif ( q < 0 ) then
      r = p

else
r = 1.0 - p
end if

if ( r <= 0.0 ) then
value = - 999999999
return value
end if

  r = ( -log ( r ) )^0.5

  if ( r <= split2 ) then
         r = r - const2

      value = app.R4POLY_VALUE(4, c(), r) / app.R4POLY_VALUE( 3, d(), r )
  else
      r = r - split2
      value = app.R4POLY_VALUE( 4, e(), r ) /app.R4POLY_VALUE(3, f(), r )
  end if

  if ( q > 0.0 ) then
      value = - value
  end if

return value[/code]

You need the following function for the above to work.

R4POLY_Value(n as int8, bdA() as double, x as double)

[code] Dim value as double
Dim i as int8
value =bdA(1)
for i = 2 to n
value = value + bdA(i)*x^(i-1)
next

return value[/code]