# Random Variate generating strange results when using ProbabilityDistribution

Bug introduced in 10.1 and fixed in 10.3

I have created a probability distribution that follows a general normal distribution, given by:

$$\frac{b *e^{-\left(\frac{\sqrt{(x-u )^2}}{a }\right)^{b }}}{2 a \Gamma \left(\frac{1}{b }\right)}$$

When I define the Probability Distribution and generate Random values for it:

modelDrPD =ProbabilityDistribution[modelDr/.{β->4.2,α-> 0.39, μ->0.06}, {x, -∞, ∞}]
Histogram[RandomVariate[modelDrPD, 1000000], "FreedmanDiaconis"]


And I get the following histogram: The PDF integrates to 1, as needed... Why there is a second peak in the Histogram? It happens with many PDF I try

Thank you

Update: Here is my complete code:

modelDr = ( b /(2 a Gamma[1/b]) Exp[-(Sqrt[(x - u)^2]/a)^b])
modelDrPD =  ProbabilityDistribution[  modelDr /. {b -> 4.2, a -> 0.39, u -> 0.06}, {x, -Infinity,    Infinity}]
Histogram[RandomVariate[modelDrPD, 1000000], "FreedmanDiaconis"]

• I can't replicate the problem in 10.3 for Mac. What version are you using? – Andy Ross Jan 12 '16 at 14:05
• After adding in the definition modelDr = \[Beta] Exp[-(Abs[ x - \[Mu]]/\[Alpha])^\[Beta]]/(2 \[Alpha] Gamma[1/\[Beta]]);, it works fine on 10.2 Windows 7. What is your definition of modelDr? – JimB Jan 12 '16 at 15:42
• With 10.1 and @JimBaldwin definition I do get the strange plot along with a bunch of underflow warnings. Fabio you should include such warnings in your question if you got them. – george2079 Jan 12 '16 at 16:24
• It is RandomVariate throwing the errors. Making the parameters rational doesn't help. RandomVariate won't accept a WorkingPrecision option for some reason. – george2079 Jan 12 '16 at 16:39
• I've added the full code... I got no warnings... for some parameters values, everything goes fine... – Fábio Jan 12 '16 at 17:31

There seems to be a bug in version 10.1 that has been fixed in 10.3. You can always try writing your own random number generator. Here is a simple acceptance rejection method based on generalized Gaussian distributions as discussed here.

Here I use a very naive envelope, a uniform distribution over {mu - s*sd, mu + s*sd} where mu is the mean of your distribution, sd is the standard deviation and s is the number of standard deviations you would like to allow the envelope to encompass.

ar[a_, b_, mu_, k_, s_: 5] :=
Block[{bag = InternalBag[{0.}, 1], low, high, m, c, u, f, chunk},
(*bounds of envelope*)
{low, high} = {mu - s #, mu + s #} &[Sqrt[a^2*Gamma[3/b]/Gamma[1/b]]];

(*scaling factor to ensure envelope is above PDF*)
m = (0.5*b)/(a*Gamma[b^(-1)])*(high - low);

(*number of random variates to generate each pass optimized for
a minimal number of passes *)
chunk = Ceiling[k/(high - low)*m];

(*add accepted variates to a bag until enough have been collected*)
While[InternalBagLength[bag] < k,
c = RandomReal[{low, high}, chunk];
u = RandomReal[{0., 1.}, chunk];
f = b/(2*a*E^((a^(-1))^b*Abs[c - mu]^b)*Gamma[b^(-1)]);
InternalStuffBag[bag, Pick[c, UnitStep[f/(m/(high - low) ) - u], 1], 1];
];

(*return only the first k random variates*)
InternalBagPart[bag, All][[1 ;; k]]
]


This could probably be made faster by compilation or by choosing a better envelope but it is a good start and doesn't issue underflow messages.

modelDr = (b/(2 a Gamma[1/b]) Exp[-(Sqrt[(x - u)^2]/a)^b]);
modelDrPD =
ProbabilityDistribution[
modelDr /. {b -> 4.2, a -> 0.39, u -> 0.06}, {x, -Infinity,
Infinity}];

Show[Histogram[ar[.39, 4.2, .06, 1000000], "FreedmanDiaconis", "PDF"],
Plot[PDF[modelDrPD, x], {x, -1, 1}]] 