Path: csiph.com!x330-a1.tempe.blueboxinc.net!usenet.pasdenom.info!weretis.net!feeder4.news.weretis.net!eternal-september.org!feeder.eternal-september.org!.POSTED!not-for-mail
From: Joe Riel <joer@san.rr.com>
Newsgroups: comp.soft-sys.math.maple
Subject: Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs!
Date: Mon, 05 Dec 2011 10:18:55 -0800
Organization: A noiseless patient Spider
Lines: 53
Message-ID: <87ipluj300.fsf@san.rr.com>
References: <0db37bbb-8efc-409c-aa16-338c77d6b556@m10g2000vbc.googlegroups.com>
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Injection-Info: mx04.eternal-september.org; posting-host="fkgjSOZ+03w83W7QjZBmIg"; logging-data="19320"; mail-complaints-to="abuse@eternal-september.org";	posting-account="U2FsdGVkX18is44PJh4l8QJ/qViEJ7YO"
User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/23.3 (gnu/linux)
Cancel-Lock: sha1:/zuLGaOcrgLboGjiL+/KeNf9gNk= sha1:un33F68+YIM00BP0fHecKXjwhvs=
Xref: x330-a1.tempe.blueboxinc.net comp.soft-sys.math.maple:260

jdm <james.d.mclaughlin@gmail.com> writes:

> with(Statistics):
> Phiinverse := (x) -> Quantile(Normal(0, 1), x) :
>
> evalf(Phiinverse(1-1/2^24));
>
> and the result is 5.691523310-2.147372008*I
>
> Now, I can get a real-valued solution to this by using the equivalent
> functionality from stats instead of Statistics.
>
> with(stats):
> PhiinvClassic:= (x) -> statevalf[icdf, normald[0, 1]](x):
>
> evalf(PhiinvClassic(1-1/2^24));
>
> However, this is not a perfect solution either:
>
> evalf(PhiinvClassic(1-1/2^24.5));
>
> Float(undefined)
>
> This may seem like an odd thing to use Maple for, but it is relevant
> to some cryptanalytic research I'm doing, and I'd be interested to
> know if this is a known bug and if anyone knows how to work around it.
>
> (Anyone who knows whether Mathematica also has this problem, please
> post!)

Try the following,

with(Statistics):
Phiinverse := (x) -> Quantile(Normal(0, 1), x):

f1 := Phiinverse(x);
                                                       1/2
                     f1 := RootOf(-erf(_Z) - 1 + 2 x) 2

f2 := subsindets(f1
                 , 'specfunc(anything,RootOf)'
                 , rootof -> unapply('fsolve'(op(rootof), _Z), x)
                );
                                                            1/2
               f2 := (x -> fsolve(-erf(_Z) - 1 + 2 x, _Z)) 2


evalf(f2(1-1/2^24));
                                  5.294704084


-- 
Joe Riel