Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.soft-sys.math.maple > #259 > unrolled thread

Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs!

Back to article view | Back to comp.soft-sys.math.maple

Contents

  Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! jdm <james.d.mclaughlin@gmail.com> - 2011-12-05 09:09 -0800
    Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! Joe Riel <joer@san.rr.com> - 2011-12-05 10:18 -0800
      Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! jdm <james.d.mclaughlin@gmail.com> - 2011-12-05 12:08 -0800
        Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! Joe Riel <joer@san.rr.com> - 2011-12-05 12:51 -0800
          Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! jdm <james.d.mclaughlin@gmail.com> - 2011-12-05 13:28 -0800
            Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! jdm <james.d.mclaughlin@gmail.com> - 2011-12-05 13:46 -0800
              Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! jdm <james.d.mclaughlin@gmail.com> - 2011-12-05 14:43 -0800
                Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! Joe Riel <joer@san.rr.com> - 2011-12-05 15:32 -0800
                  Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! jdm <james.d.mclaughlin@gmail.com> - 2011-12-05 16:04 -0800
                    Re: Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs! Herman Rubin <hrubin@skew.stat.purdue.edu> - 2011-12-07 02:06 +0000

#259 — Inverse CDF of standard normal is (according to Maple 15) complex valued for some positive inputs!

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!)

[toc] | [next] | [standalone]

#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

[toc] | [prev] | [next] | [standalone]

#261

Thanks! Do you know how I can get this to work with plot? When I try

> plot(f2(2^(-a)), a = 1 .. 56);

I get "Error, (in fsolve) a is in the equation, and is not solved
for".

James.

> 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

[toc] | [prev] | [next] | [standalone]

#262

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

> Thanks! Do you know how I can get this to work with plot? When I try
>
>> plot(f2(2^(-a)), a = 1 .. 56);
>
> I get "Error, (in fsolve) a is in the equation, and is not solved
> for".

Try plot('f2(2^(-a))', a=1..56);

I see, though, that there is a better solution on MaplePrimes.


>
> James.
>
>> 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
>

-- 
Joe Riel

[toc] | [prev] | [next] | [standalone]

#263

I can assure you that the issue is so far not solved by the discussion
on MaplePrimes. The most recent post I've made there explains that
Phiinverse was intended only as part of a more complicated formula
that I was trying to plot, and which I can't easily differentiate (and
going by the chain rule, the derivative involves Phiinverse anyway). I
discovered the problem with Phiinverse while trying to isolate the
root cause of the problem, and I have to get it to work perfectly for
the range a=1..256 at some point.

> Try plot('f2(2^(-a))', a=1..56);

I've already tried this - and the graph still goes wrong for the
higher values of a.

[toc] | [prev] | [next] | [standalone]

#264

mapleprimes.com is a bit haywire - so I'll need to add here that one
poster's suggestion of Digits:=18 has helped a little - using this in
the range a=1..56, the graphs contain discontinuities and correct
values, but no obviously incorrect values as they did before. I had to
increase this for higher values of a; digits:=40 has so far been okay-
ish in the same way for a up to 80. Hopefully increasing Digits will
suffice up to a=256.

Thanks to that poster - I can't seem to thank you on MaplePrimes as
posts (and parts thereof) keep on vanishing and reappearing.

[toc] | [prev] | [next] | [standalone]

#265

> Try plot('f2(2^(-a))', a=1..56);

I didn't notice those two apostrophes earlier - this does work a *lot*
better! So far so good...

[toc] | [prev] | [next] | [standalone]

#266

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

>> Try plot('f2(2^(-a))', a=1..56);
>
> I didn't notice those two apostrophes earlier - this does work a *lot*
> better! So far so good...

You will have to increase Digits to get to 256.

-- 
Joe Riel

[toc] | [prev] | [next] | [standalone]

#267

> You will have to increase Digits to get to 256.

True - that said, increasing from 100 to 170 didn't seem to offer much
improvement. I think I'll have to decide whether the gaps in the
graphs right now are something I can live with.

Thanks.

[toc] | [prev] | [next] | [standalone]

#268

On 2011-12-06, jdm <james.d.mclaughlin@gmail.com> wrote:
>> You will have to increase Digits to get to 256.

> True - that said, increasing from 100 to 170 didn't seem to offer much
> improvement. I think I'll have to decide whether the gaps in the
> graphs right now are something I can live with.

> Thanks.

One rarely gets means of evaluating functions in the tails of
distributions from packages which have not anticipated the problem.
One can do quite well if one thinks, and computers do not.

The problem is to find the solution of

	int(exp(-t^2/2) , t=x..infinity) = sqrt(2*Pi)*a.

Let us call b the right side, and take logarithms.  Using
Newton's method to solve that equation works well, as the
derivative of the logarithm of the integral is -exp(-t^2/2)/int.

One can do even better, but this will work well.
--
This address is for information only.  I do not claim that these views
are those of the Statistics Department or of Purdue University.
Herman Rubin, Department of Statistics, Purdue University
hrubin@stat.purdue.edu         Phone: (765)494-6054   FAX: (765)494-0558

[toc] | [prev] | [standalone]

Back to top | Article view | comp.soft-sys.math.maple

csiph-web