Monthly Archive for May, 2010

Vensim doesn’t have a function for the cumulative normal distribution, but it’s easy to implement via a macro. I used to use a polynomial cited in Numerical Recipes (error function, Ch. 6.2):

:MACRO: NCDF(x)

NCDF = 1-Complementary Normal CDF

~	dmnl

~		|

Complementary Normal CDF=
ERFCy/2
~	dmnl
~		|

ERFCy = IF THEN ELSE(y>=0,ans,2-ans)
~	dmnl
~	http://www.library.cornell.edu/nr/bookcpdf/c6-2.pdf
|

y = x/sqrt(2)
~	dmnl
~		|

ans=t*exp(-z*z-1.26551+t*(1.00002+t*(0.374092+t*(0.0967842+
t*(-0.186288+t*(0.278868+t*(-1.1352+t*(1.48852+
t*(-0.822152+t*0.170873)))))))))
~	dmnl
~		|

t=1/(1+0.5*z)
~	dmnl
~		|

z = ABS(y)
~	dmnl
~		|

:END OF MACRO:

I recently discovered a better approximation here, from algorithm 26.2.17 in Abromowitz and Stegun, Handbook of Mathematical Functions:

:MACRO: NCDF2(x)

NCDF2 =  IF THEN ELSE(x >= 0,

(1 - c * exp( -x * x / 2 ) * t *

( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 )),  ( c * exp( -x * x / 2 ) * t *

( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ))

)

~     dmnl

~     From http://www.sitmo.com/doc/Calculating_the_Cumulative_Normal_Distribution

Implements algorithm 26.2.17 from Abromowitz and Stegun, Handbook of Mathematical \

Functions. It has a maximum absolute error of 7.5e^-8.

http://www.math.sfu.ca/

|

c  =  0.398942

~     dmnl

~           |

t = IF THEN ELSE( x >= 0, 1/(1+p*x), 1/(1-p*x))

~     dmnl

~           |

b5 =  1.33027

~     dmnl

~           |

b4 = -1.82126

~     dmnl

~           |

b3 =  1.78148

~     dmnl

~           |

b2 = -0.356564

~     dmnl

~           |

b1 =  0.319382

~     dmnl

~           |

p  =  0.231642

~     dmnl

~           |

:END OF MACRO:

In advanced Vensim versions, paste the macro into the header of your model (View>As Text). Otherwise, you can implement the equations inside the macro directly in your model.