4
$\begingroup$

From this post (95069) I can see that this question has been been given a workaround for symmetric PDFs and that the bug was eventually addressed. I checked and the fix is still in place for my current version 11.1.0.

My issue is that I am working with the CDF of the Generalised Pareto Distribution with is not symmetric and I am getting a similar issue.

With the CDF

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

I created a distribution function

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

It passes the basic checks including

gpdDist[μ, ξ, σ] /. ProbabilityDistribution -> Integrate
1

Creating an instance of this distribution and sampling with RandomVariate leads to values outside the support of the distribution.

dist = gpdDist[.5 10^6, .5, 1 10^6];
dist /. ProbabilityDistribution -> NIntegrate
Quantile[dist, 0.00000000001]
1.
500000.

The minimum value of the distribution is 500,000 == μ. However, RandomVariate routinely returns values far, far below μ. You may have to evaluate the lines more than once to see it occur.

Min@RandomVariate[dist, 10000]
-2.38204*10^10

and

BoxWhiskerChart[RandomVariate[dist, 10000], "Outliers"]

Mathematica graphics

Have I missed something? If not are there any workarounds?

$\endgroup$
5
$\begingroup$

I think this is a precision issue. Here is a function that computes sets of 10000 random variates with different random seeds:

minVariate[dist_, prec_] := Table[
    SeedRandom[n];
    Min @ RandomVariate[dist, 10000, WorkingPrecision->prec],
    {n, 100}
]

Let's see how many minimum variates are greater than 0 (there should be 100):

dist = gpdDist[.5 10^6, .5, 10^6];
Total @ UnitStep @ minVariate[dist, MachinePrecision]

29

Only 29 when using MachinePrecision. Now, let's increase the working precision (I also use Quiet since the distribution dist uses machine numbers):

Quiet @ Total @ UnitStep @ minVariate[dist, 20]
Quiet @ Total @ UnitStep @ minVariate[dist, 25]
Quiet @ Total @ UnitStep @ minVariate[dist, 30]

41

96

100

At a working precision of 30, the occurrence of spurious out of range variates has significantly decreased.

$\endgroup$
  • $\begingroup$ (1+) May you briefly elaborate on why my machine precision of 15.9546 is not enough. Is it due Log or InverseFunction being used for the default method in the pseudo-random number generator? $\endgroup$ – Edmund Jul 7 '17 at 2:26
  • 1
    $\begingroup$ @Edmund Actually, I don't know how RandomVariate works, sorry. The value of the out of range variate was about 10^16 less than the min value, which is why I suspected that precision played a role. $\endgroup$ – Carl Woll Jul 7 '17 at 2:45
6
$\begingroup$

Here's a workaround:

rGenPareto[μ_, ξ_, σ_, n_] := Module[{p}, p = RandomReal[{0, 1}, n];
  If[ξ == 0, μ + σ Log[1/(1 - p)], μ+((-1+(1-p)^-ξ) σ)/ξ]]

Min[rGenPareto[.5 10^6, .5, 1 10^6, 10000000]]
(* 500000.2263487703` *)

How to construct this? One can choose a uniform random number between 0 and 1 and then equate that to the associated value of the CDF to obtain a random sample from that distribution.

The generalized Pareto described above has two different equations depending on if $\xi=0$ or $\xi\neq 0$. Using Solve gets the associated inverse functions:

(* ξ != 0 *)
x /. FullSimplify[Solve[1 - (1 + ((x - μ)*ξ)/σ)^(-ξ^(-1)) == p, x]][[1]]
(* μ+((-1+(1-p)^-ξ) σ)/ξ *)

(* ξ == 0 *)
x /. (Solve[1 - E^(-((x - μ)/σ)) == p, x] /. C[1] -> 0)[[1]]

(* μ+σ Log[1/(1-p)] *)
$\endgroup$
  • $\begingroup$ (+1) By first principles? Do you have a source link as I was trying to find this by searching for "derive generalised pareto PDF" in Google but came up empty. $\endgroup$ – Edmund Jul 7 '17 at 2:27
  • $\begingroup$ Yes. I added what I hope is an explanation. $\endgroup$ – JimB Jul 7 '17 at 2:42
1
$\begingroup$

Is this satisfactory (using test values in OP):

pd = ProbabilityDistribution[
  D[genParetoCDF[0.5 10^6, 0.5, 1, x], x], {x, -Infinity, Infinity}]
Show[Histogram[rv = RandomVariate[pd, 1000], Automatic, "CDF"], 
 Plot[genParetoCDF[0.5 10^6, 0.5, 1, x], {x, 500000, 500008}]]
BoxWhiskerChart[rv]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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