# RandomVariate sometimes fails for a custom distribution

From the following custom probability density function:

fhm[X_, Ao_, x_] := (X/Ao^X) x^(X - 1);


I want to get a random variable of $$10^6$$ points of values. So i do the following:

DD = ProbabilityDistribution[fhm[X, Ao, x], {x, 0, Ao}];
hm = RandomVariate[DD, 10^6];


After that hm for the same input parameters, sometimes creates a random variable and sometimes returns RandomVariate[] with the given input distribution.

I edited my question adding a full script below. The parameter that causes the issue is sigmas. For example for sigmas=0.1, i have no issue no matter how many times i restart mathematica and the script. Also note that X and Ao are always positive real numbers.

Nsamples = 10^6; (* Number of samples *)

f = 350 10^9; (*GHz to Hz*)
d = 60; (*distance m*)
sigmas = 0.2;

c = 299792458; (*m/s lightspeed*)
Gtdb = 55;
Grdb = 55;
Gt = 10^(Gtdb/10);
Gr = 10^(Grdb/10);
lamda[f_] := c/f; (*wavelength in m*)
heff = 1;
Aeff[Gr_, f_] := (Gr lamda[f]^2)/(4 π heff);
r[Gr_, f_] := Sqrt[(Aeff[Gr, f]/Pi) ];
Diam[Gr_, f_] := 2 r[Gr, f];
theta3dB[Gr_, f_] := 70 lamda[f]/Diam[Gr, f];
wd[d_, Gr_, f_] := d Tan[(theta3dB[Gr, f] (Pi/180))/2];
v[Gr_, f_, d_] := (Sqrt[Pi] r[Gr, f])/(Sqrt[2] wd[d, Gr, f]);
A0[Gr_, f_, d_] := Erf[v[Gr, f, d]]^2 ;
weq2[d_, Gr_, f_] := (wd[d, Gr, f])^2  (Sqrt[π] Erf[v[Gr, f, d]])/(
2 v[Gr, f, d] Exp[-(v[Gr, f, d])^2]);
Xi[d_, Gr_, f_] := weq2[d, Gr, f]/(4 (sigmas^2) );

X = Xi[d, Gr, f];
Ao = A0[Gr, f, d];

fhm[X_, Ao_, x_] := X/Ao^X x^(X - 1);(*PDF*)

In[28]:= N[X]
N[Ao]
Out[28]= 0.304121
Out[29]= 0.243095

In[30]:= Integrate[fhm[X, Ao, x], {x, 0, Ao}]
Out[30]= 1.

h = ProbabilityDistribution[fhm[X, Ao, x], {x, 0, Ao}]

hm = RandomVariate[h, Nsamples]
In[35]:= hm[[30]]

Out[35]= 0.0943069

In[36]:= hm[[1434]]

Out[36]= 0.0912831

In[34]:= hm1 = RandomVariate[h, Nsamples]

Out[34]= RandomVariate[
ProbabilityDistribution[
0.467567/\[FormalX]^0.695879, {\[FormalX], 0,
Erf[(21413747 Cot[(7 π^2)/(3600 10^(3/4))])/(
3000000000 2^(3/4) 5^(1/4) Sqrt[π])]^2}], 1000000]

• The function is returned with the given input distribution if one or more of the variables is not defined. As a sanity check, did you try also printing out the value of X and Ao at the same time? Perhaps you are clearing their values somewhere in the code...
– a20
May 27, 2019 at 12:53
• Yes, i checked it. X and Ao have numerical values May 27, 2019 at 13:15
• Works for me: i.stack.imgur.com/g6alu.png, i.stack.imgur.com/noHtY.png -- We need an example where it does not work, I guess. May 27, 2019 at 15:53
• I agree, we need an example, I do not see anything wrong with the code itself. I would guess it's either a bug, some computer issue (RAM overuse?), or that you are redefining/clearing the values of the parameters somewhere. What do you mean by that it works sometimes but sometimes not? Do you e.g. call the function at different places in your code, or do you rerun the same block of code, or do you have some loop over the function call?
– a20
May 27, 2019 at 17:46
• When using ProbabilityDistribution the constraints on the distribution's parameters should be included as Assumptions, i.e., DD[X_, Ao_] = ProbabilityDistribution[ fhm[X, Ao, x], {x, 0, Ao}, Assumptions -> Ao > 0 && X > 0]; These assumptions are then available to other built-in functions via DistributionParameterAssumptions[DD[X, Ao]] May 28, 2019 at 4:44

At the end, you have the following distribution:

h

ProbabilityDistribution[0.467567/\[FormalX]^0.695879, {\[FormalX], 0,
Erf[(21413747 Cot[(7 π^2)/(3600 10^(3/4))])/(
3000000000 2^(3/4) 5^(1/4) Sqrt[π])]^2}]


Note that the pdf (the first argument) is a machine number function. It seems that with a machine number argument, RandomVariate of the probability distribution occasionally fails. I think it's worth reporting this behavior to support. A simpler version that exhibits the same issue:

dist = N @ h

ProbabilityDistribution[0.467567/\[FormalX]^0.695879, {\[FormalX], 0., 0.243095}]


Code exhibiting issue:

Tally @ Table[ListQ @ RandomVariate[dist, 10], 100]


{{False, 61}, {True, 39}}

The output should consist of only True results.

Note that it is possible to workaround this issue by using exact arithmetic:

exact = Rationalize[dist, 0];
Tally @ Table[ListQ @ RandomVariate[exact, 10], 100]


{{True, 100}}

Or by using extended precision reals:

extended = dist /. r_Real :> SetPrecision[r, 20];
Tally @ Table[ListQ @ RandomVariate[extended, 10], 100]


{{True, 100}}