# Normal distribution function

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

Any ideas how to do this within Xojo?

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.

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.

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.

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 = (((((N6
z + N5)*z + N4)*z + N3)*z + N2)*z + N1)z + N0
d = ((((((M7
z + M6)*z + M5)*z + M4)*z + M3)*z + M2)z + M1)z + M0
c = e
n/d
else
f = z + 1.0/(z + 2.0/(z + 3.0/(z + 4.0/(z + 13.0/20.0))))
c = e/(RT2PI
f)
end if
end if

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

return c[/code]

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:

’ 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]``