4
$\begingroup$

From this answer (19542) to my question I need to add the WorkingPrecision option to RandomVariate for that particular ProbabilityDistribution not to experience what I guess is an overflow off some sort and return incorrect values.

Another option given by WRI was to use Rationalize[ #, 0]& on each of the parameter values to effectively give them infinite precision such that the calculations would not experience any precision anomalies.

Is there a way to set the precision on the function such that it does all its calculations at the prescribe precision?

I tried by setting $MinPrecision and $MaxPrecision as in this answer (19542).

ClearAll[genParetoCDF]; 
genParetoCDF[μ_, ξ_, σ_, x_] := 
 Block[{$MinPrecision = 2 $MachinePrecision, $MaxPrecision = 2 $MachinePrecision},
  Piecewise[{
   {1 - (1 + ((x - μ)*ξ)/σ)^(-ξ^(-1)), 
     (x >= μ && ξ > 0) || (μ <= x <= μ - σ/ξ && ξ < 0)}, 
   {1 - E^(-((x - μ)/σ)), 
     x >= μ && ξ == 0}
  }]
 ]

However it still returns MachinePrecision values.

Precision@genParetoCDF[.2, .5 10^6, 100000, 1 10^6]
MachinePrecision

When used with RandomVariate without WorkingPrecision it still results in negatives (from what I think is the overflow).

With

ClearAll[gpdDist];
gpdDist[μ_, ξ_, σ_] := 
 ProbabilityDistribution[{"CDF", genParetoCDF[μ, ξ, σ, x]}, {x, μ, ∞}, 
  Assumptions -> {{μ, ξ, σ} ∈ Reals, σ > 0}]

Then (Also curious as to why this does not generate some sort of overflow error)

Min@RandomVariate[dist, 1000000]
-7.49163*10^6

Is it possible to set the precision on the function? I still need to be able to use the function symbolically as well. For example,

Assuming[{ξ > 0, σ > 0, {ξ, μ, σ} ∈ Reals},
 Limit[genParetoCDF[ξ, μ, σ, x], x -> ∞]
 ]

Any ideas?

$\endgroup$
1
  • $\begingroup$ It yields machine precision because you feed it machine precision arguments. Machine precision is infectious: almost any calculation that has machine precision as input will have machine precision as output unless you do something that explicitly demands otherwise. $\endgroup$
    – John Doty
    Commented Jul 11, 2017 at 17:02

2 Answers 2

1
$\begingroup$

Use SetPrecision[] to set the initial precision of your arguments:

ClearAll[genParetoCDF];
genParetoCDF[μ0_, ξ0_, σ0_, x0_] := 
  Module[{μ, ξ, σ, x},
   {μ, ξ, σ, x} = SetPrecision[{μ0, ξ0, σ0, x0}, 2 $MachinePrecision];
   Block[{$MinPrecision = 2 $MachinePrecision, $MaxPrecision = 2 $MachinePrecision},
    Piecewise[
       {{1 - (1 + ((x - μ)*ξ)/σ)^(-ξ^(-1)),
         (x >= μ && ξ > 0) || (μ <= x <= μ - σ/ξ && ξ < 0)},
        {1 - E^(-((x - μ)/σ)), x >= μ && ξ == 0}}]
    ]];

genParetoCDF[.2, .5 10^6, 100000, 1 10^6]
Precision[%]
(*
  0.000030849421087619473139643876639982
  31.9092
*)

You might consider controlling the precision of the inputs to gpdDist, too:

gpdDist[μ_, ξ_, σ_] := 
 ProbabilityDistribution[{"CDF", genParetoCDF[μ, ξ, σ, x]},
  {x, SetPrecision[μ, 2*$MachinePrecision], ∞}, 
  Assumptions -> {{μ, ξ, σ} ∈ Reals, σ > 0}];

Then I get results like this (no overflow, no variates outside domain):

dist = gpdDist[.5 10^6, .2, .1 10^6];

RandomVariate[dist, WorkingPrecision -> 2*$MachinePrecision]
(*  512877.53799123430547484091911214  *)

Min@RandomVariate[dist, 1000000, WorkingPrecision -> 2*$MachinePrecision]
(*  500000.01897993857486312711526885  *)

It also works symbolically (see caveat below about symbolic expressions and numerics):

Assuming[{ξ > 0, σ > 0, {ξ, μ, σ} ∈ Reals}, 
 Limit[genParetoCDF[ξ, μ, σ, x], x -> ∞]]
(*  Piecewise[{{1, μ >= 0}}, 0]  *)

You might want to consider issuing a warning if the input has a precision less than 2 $MachinePrecision. Using input less than the working precision should be considered poor practice, if not an error. Further, once you have the symbolic CDF function, the symbol-variables in it will no longer have the precision control, since such control applies only to numbers and by extension can be applied only in numeric routines. Code like RandomVariate[dist] will use machine precision inputs, and numeric error in this case will lead to bogus results. If you want the symbolic Piecewise CDF expression to use in calculations, then you have to give input values for the symbols of sufficient precision. (Actually, that would be my recommended practice in this case, as I currently understand it. I would probably get rid of the $MinPrecision/$MaxPrecision stuff. If you're trying to deal with unreliable users, then your code must catch any of their numeric input and check/fix its precision. This is hard to do, if they are allowed to use ReplaceAll and other built-in functions, bypassing your code, on the symbolic CDF.)

$\endgroup$
3
  • $\begingroup$ This does get genPareto returning the requested precision. However, it still suffers the overflow issue when used with RandomVariate, with dist = gpdDist[.5 10^6, .2, .1 10^6] for example. I can't use ?NumericQ as that will prevent me from using the function symbolically. $\endgroup$
    – Edmund
    Commented Jul 11, 2017 at 19:04
  • $\begingroup$ @Edmund Note that my answer works symbolically and allows arbitrary argument precision to override the default. Can't help you with RandomVariate, though. $\endgroup$
    – John Doty
    Commented Jul 11, 2017 at 19:56
  • $\begingroup$ @Edmund That's because you don't set the precision of gpdDIst and you get the same issue mentioned in chat -- it uses machine precision when you give m.p. input, and therefore so does dist and RandomVariate. The overflow/negative random variates has to do with rounding error, I suppose. I get no such problems with high precision inputs. $\endgroup$
    – Michael E2
    Commented Jul 12, 2017 at 1:06
2
$\begingroup$

Write a generic function, and then define special cases that recognize machine number arguments. Translate those to the necessary precision, and then use your generic function.

 genParetoCDF[μ_, ξ_, σ_, x_] := 
   Piecewise[{{1 - (1 + ((x - μ)*ξ)/σ)^(-ξ^(-1)), (x >= μ && ξ > 0) || (μ <= 
        x <= μ - σ/ξ && ξ < 0)}, {1 - 
     E^(-((x - μ)/σ)), x >= μ && ξ == 0}}]

genParetoCDF[μ_, ξ_, σ_, x_] := 
 genParetoCDF[
   SetPrecision[μ, 2 $MachinePrecision], ξ, σ, x] /; 
MachineNumberQ[μ]

genParetoCDF[μ_, ξ_, σ_, x_] := 
 genParetoCDF[μ, 
   SetPrecision[ξ, 2 $MachinePrecision], σ, x] /; 
  MachineNumberQ[ξ]

genParetoCDF[μ_, ξ_, σ_, x_] := 
 genParetoCDF[μ, ξ, 
   SetPrecision[σ, 2 $MachinePrecision], x] /; 
  MachineNumberQ[σ]

genParetoCDF[μ_, ξ_, σ_, x_] := 
 genParetoCDF[μ, ξ, σ, 
   SetPrecision[x, 2 $MachinePrecision]] /; MachineNumberQ[x]

Precision[genParetoCDF[.2, .5 10^6, 100000, 1 10^6]]
(* 31.8819 *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.