Random Variate generating strange results when using ProbabilityDistribution

Question

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"] 

Answer 1

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 = Internal`Bag[{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[Internal`BagLength[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)]);
      Internal`StuffBag[bag, Pick[c, UnitStep[f/(m/(high - low) ) - u], 1], 1];
    ];

    (*return only the first k random variates*)
    Internal`BagPart[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}]]