# Translating algorithms from FORTRAN etc.

I'm intending to rationalise my J code for various probability distributions. Several algorithms have been translated from other sources such as the mainly FORTRAN Applied Statistics algorithms at http://lib.stat.cmu.edu/apstat/. Translated algorithms must, at least initially, compromise between (1) mirroring the original presentation (making errors easier to spot) and (2) converting to 'proper' J.

In some cases an alternative would be to link to (e.g.) the GNU scientific library.

Here's an example (still under development) that some people may find useful. I'd be grateful for any comments.

## Inverse of the Standard Normal N(0,1) CDF

```NB. N(0,1) inverse cdf
NB. Translated from  http://lib.stat.cmu.edu/apstat/241
NB. =========================================================

NB. =========================================================
NB. Utilities (inserted for completeness)
NB.
NB.* vftxt    v  numeric vector from text string
NB.* ratpoly  c  rational polynomial approximation
NB.

WHITESPACE=. 9 10 13 32 { a.
vftxt=. [: ". e.&WHITESPACE {"0 1 ,.&' '

ratpoly=. 2 : 'm&p. % (1,n)&p.'

NB. =========================================================
NB. ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
NB.
NB. Produces the normal deviate Z corresponding to a given lower
NB. tail area of P; Z is accurate to about 1 part in 10^16.
NB.

SPLIT1=. 0.425 [ SPLIT2=. 5.0
CONST1=. 0.180625 [ CONST2=. 1.6

NB. Coefficients for P close to 0.5

A=. vftxt 0 : 0
3.3871328727963666080
1.3314166789178437745e2
1.9715909503065514427e3
1.3731693765509461125e4
4.5921953931549871457e4
6.7265770927008700853e4
3.3430575583588128105e4
2.5090809287301226727e3
)

B=. vftxt 0 : 0
4.2313330701600911252e1
6.8718700749205790830e2
5.3941960214247511077e3
2.1213794301586595867e4
3.9307895800092710610e4
2.8729085735721942674e4
5.2264952788528545610e3
)

ratAB=. A ratpoly B

NB. Coefficients for P not close to 0, 0.5 or 1.

C=. vftxt 0 : 0
1.42343711074968357734
4.63033784615654529590
5.76949722146069140550
3.64784832476320460504
1.27045825245236838258
2.41780725177450611770e_1
2.27238449892691845833e_2
7.74545014278341407640e_4
)

D=. vftxt 0 : 0
2.05319162663775882187
1.67638483018380384940
6.89767334985100004550e_1
1.48103976427480074590e_1
1.51986665636164571966e_2
5.47593808499534494600e_4
1.05075007164441684324e_9
)

ratCD=. C ratpoly D

NB. Coefficients for P near 0 or 1.

E=. vftxt 0 : 0
6.65790464350110377720
5.46378491116411436990
1.78482653991729133580
2.96560571828504891230e_1
2.65321895265761230930e_2
1.24266094738807843860e_3
2.71155556874348757815e_5
2.01033439929228813265e_7
)

F=. vftxt 0 : 0
5.99832206555887937690e_1
1.36929880922735805310e_1
1.48753612908506148525e_2
7.86869131145613259100e_4
1.84631831751005468180e_5
1.42151175831644588870e_7
2.04426310338993978564e_15
)

ratEF=. E ratpoly F

qfp=. -&0.5
test1=. SPLIT1 < |@qfp
r1fq=. CONST1 - *:
nd1=. (] * ratAB@r1fq)@qfp
r2fp=. [: %: [: - [: ^. ] <. -.
test2=. >&SPLIT2
nd2fr=. ratCD@-&CONST2
nd3fr=. ratEF@-&SPLIT2

nd=. nd1 ` (*@qfp * (nd2fr`nd3fr @. test2)@r2fp) @. test1

in01=. >&0 * 1 + <&1   NB. test for y in range
n01cdfinv=: (__"0 ` _: `nd @. in01)"0 f.```

## Example

```   20j14 ": n01cdfinv ,. 0 0.158655253931457 0.977249868051821 2.86651571879194e_7 1 10
__
_1.00000000000000
2.00000000000000
_5.00000000000000
_
_```

Contributed by EwartShaw

EwartShaw/N01CdfInv (last edited 2009-09-08 08:50:05 by OlegKobchenko)