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. Ewart Shaw 1-Nov-2006 last modified 6-Nov-2006
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