Thanks Devon for the general code. It is always useful to have the general
case, even if it is often much slower.
The code John referenced converts to qnorm in the following J script. It
can be shortened and improved, but runs quite fast. The script provides
functions with the conventions for the names the same as those used in R.
NB. ======================================================================
Note 'Normal Distribution'
The following are functions based on the standard normal
distribution.
dnorm the density
pnorm the cumulative distribution function
qnorm the quantile function - inverse of pnorm
rnorm random deviates
The functions dnorm, pnorm and qnorm all take an array
argument. The function rnorm has a single argument, the number of
random deviates required.
These functions use definitions from the J library, Ewart Shaw,
Brian Schott, Ric Sherlock and the programming forum. qnorm uses the
algorithm from P.J. Acklam reported by John Randall.
The present implementation does not provide for any iterative improvement
as noted in http://home.online.no/~pjacklam/notes/invnorm/
)
dnorm =: (% %: 2p1) * ^@:(_0.5 * *:)
erf =: (*&(%:4p_1)%^@:*:)*[:1 H. 1.5*: NB. A&S 7.1.21 (right)
pnorm =: [:-:1:+[:erf%&(%:2) NB. A&S 26.2.29 (solved for P)
NB. Inverse of Normal Distribution Function
qnorm =: 3 : 0
polyn =: (|[EMAIL PROTECTED]) p. ] NB. function to reverse coefficients and evaluate
polynomials
a =. _3.969683028665376e01 2.209460984245205e02 _2.759285104469687e02
a =. a, 1.383577518672690e02 _3.066479806614716e01 2.506628277459239e00
b =. _5.447609879822406e01 1.615858368580409e02 _1.556989798598866e02
b =. b, 6.680131188771972e01 _1.328068155288572e01 1
c =. _7.784894002430293e_03 _3.223964580411365e_01 _2.400758277161838e00
c =. c, _2.549732539343734e00 4.374664141464968e00 2.938163982698783e00
d =. 7.784695709041462e_03 3.224671290700398e_01 2.445134137142996e00
d =. d, 3.754408661907416e00 1
NB. Define break-points.
p_low =. 0.02425
p_high =. 1 - p_low
NB. Rational approximation for lower region.
if. (0 < y) *. y < p_low do.
q =. %: -2*^. y
v =. (c polyn q) % d polyn q
end.
NB. Rational approximation for central region.
if. (p_low <: y) *. y <: p_high do.
q =. y - 0.5
r =. q*q
v =. q * (a polyn r) % b polyn r
end.
NB. Rational approximation for upper region.
if. (p_high < y) *. y < 1 do.
q =. %: -2* ^. 1-y
v =. -(c polyn q) % d polyn q
end.
v
)
rand01 =: [EMAIL PROTECTED] 0:
NB. from normalrand3 by Brian Schott
rnorm =: 3 : 0
n =. >. -: y
a=. %: - +: ^. rand01 n
b=. +: o. rand01 n
r1 =. a * 2 o. b
r2 =. a * 1 o. b
r1,r2
)
NB. ====================================================================
----- Original Message -----
From: "Devon McCormick" <[EMAIL PROTECTED]>
To: "Programming forum" <[email protected]>
Sent: Tuesday, September 09, 2008 3:13 PM
Subject: Re: [Jprogramming] Inverse of Normal CDF
I might point out that the version I submitted provides the framework for a
general inverse of a cdf, unlike the undoubtedly superior, but very
arbitrary-coefficienty version, to which John refers.
On 9/8/08, John Randall <[EMAIL PROTECTED]> wrote:
Sherlock, Ric wrote:
> Can anyone provide a definition for its inverse?
There is an algorithm and some implementations at
http://home.online.no/~pjacklam/notes/invnorm/
Best wishes,
John
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
--
Devon McCormick, CFA
^me^ at acm.
org is my
preferred e-mail
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
