Ambrus's solution is probably better than mine, but I did this all in J this
way:
NB.* inverseCND: inverse of cumulative norm. dist.; OK between 1e_10 &
1-1e_10.
inverseCND=: 3 : 0"0
if. (y>0.4999999998)*.y<0.500000002 do. 0 NB. Unstable near inflection
else.
f=: 13 : ((19j17":y),'-cnd y') NB. Function to solve
f newton^:18]estInitInvCND y
end.
NB.EG _2.32634787404084-:inverseCND 0.01 NB. 1% of distribution to the
left
NB.EG 1.64485362695147-:inverseCND 0.95 NB. 5% of distribution to the
right
)
ESTINVCND=: (i:6),:~1-cnd i:6 NB. Use to estimate inverse of cnd
estInitInvCND=: 3 : 0"0
whst=. 0 i.~/:y,0{ESTINVCND
ix=. (<:1{$ESTINVCND)<.0>.whst-1 0
-:+/ix{1{ESTINVCND
)
NB.* newton: Newton's iterative method to find zero of fnc given estimate.
newton=: 1 : '- u % u D. 1'
NB.* erf: error function
erf =: (*&(%:4p_1)%^@:*:)*[:1 H. 1.5*: NB. A&S 7.1.21 (right)
NB.* cnd: cumulative normal distribution
cnd =: [:-:1:+[:erf%&(%:2) NB. A&S 26.2.29 (solved for P)
explainFns=: 0 : 0
[from http://www.jsoftware.com/jwiki/Doc/Articles/Play193 by Eugene
McDonnell]
We haven't ended quite yet. Perhaps you remember the article by Ewart Shaw
in Vector 18.4, in which he defined the error function erf using J's
hypergeometric conjunction:
erf =: (*&(%:4p_1)%^@:*:)*[:1 H. 1.5*: NB. A&S 7.1.21 (right)
and then defined the cumulative distribution function of the normal
distribution by:
cnd =: [:-:1:+[:erf%&(%:2) NB. A&S 26.2.29 (solved for P)
)
On Mon, Sep 8, 2008 at 1:28 AM, Zsbán Ambrus <[EMAIL PROTECTED]> wrote:
> On Mon, Sep 8, 2008 at 03:27, Sherlock, Ric <[EMAIL PROTECTED]>
> wrote:
> > I'd like to have a J implementation of the probit function (the inverse
> CDF (aka quantile function) of the standard normal distribution.
>
> Sure, just use the implementation in GSL
> (http://www.gnu.org/software/gsl/). I think this works:
>
> probit =: 'libgsl.so gsl_cdf_ugaussian_Pinv > d d' & (15!:0"_ 0)
> probit (%~i.)10
> __ _1.28155 _0.841621 _0.524401 _0.253347 0 0.253347 0.524401 0.841621
> 1.28155
>
> Ambrus
> ----------------------------------------------------------------------
> 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
