4
$\begingroup$

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:

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"] 
$\endgroup$
  • $\begingroup$ I can't replicate the problem in 10.3 for Mac. What version are you using? $\endgroup$ – Andy Ross Jan 12 '16 at 14:05
  • $\begingroup$ 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? $\endgroup$ – JimB Jan 12 '16 at 15:42
  • $\begingroup$ 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. $\endgroup$ – george2079 Jan 12 '16 at 16:24
  • $\begingroup$ It is RandomVariate throwing the errors. Making the parameters rational doesn't help. RandomVariate won't accept a WorkingPrecision option for some reason. $\endgroup$ – george2079 Jan 12 '16 at 16:39
  • $\begingroup$ I've added the full code... I got no warnings... for some parameters values, everything goes fine... $\endgroup$ – Fábio Jan 12 '16 at 17:31
11
$\begingroup$

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

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.