2
$\begingroup$

I would like to generate a random variates in Mathematica 9 with this distribution:

f[tetta_, x_] := ( 4 Csc[tetta/2]^4 Sin[tetta] Sin[2 x] (x Cos[x Sin[tetta/2]] -
     Csc[tetta/2] Sin[x Sin[tetta/2]])^2) / (2 Sin[2 x] (-1 - 2 x^2 + 2 x^4 + 2 x Sin[2 x]) +
     Sin[4 x]);

mydist[x_] := ProbabilityDistribution[f[t, x], {t, 0, Pi}, Assumptions -> x > 0];

And it's works well. But in some cases it's doesn't. For example:

RandomVariate[mydist[78]]
     (* 0.0354395 *)
RandomVariate[mydist[79]]

And for the last one Mathematica gives me: Hmm...

The same situation is for other parameter values: for 80.2 RandomVariate is working fine, but not for 80.3.

Could you help me to understand why is so happens and how is to solve that problem?

Thank you in advance.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Feb 17 '15 at 0:19
1
$\begingroup$

The issue is that your specification of the pdf, f has so many cycles in it for certain values of x that Mathematica is unable to verify that you have specified a valid pdf, because it loses too much precision in the integration.

I checked this using a simple Manipulate.

Manipulate[
 Plot[f[tt, m], {tt, 0, Pi}, PlotRange -> All, 
  Epilog ->  Text[ToString[NIntegrate[f[ttt, m], {ttt, 0, Pi}]], 
    Scaled[{0.7, 0.8}]]], {m, 1, 100}]

And for certain values, including 60.2, I got pictures like this:

enter image description here

with error messages about "catastrophic loss of precision".

One could try increasing the precision of the integration like this:

Manipulate[
 Plot[f[tt, m], {tt, 0, Pi}, PlotRange -> All, 
  Epilog -> 
   Text[ToString[
     NIntegrate[SetPrecision[f[ttt, m], 45], {ttt, 0, Pi}, 
      WorkingPrecision -> 40]], Scaled[{0.7, 0.8}]]], {m, 1, 100}]

enter image description here

You can see that this doesn't integrate to 1 and therefore isn't a valid probability density function for this value of $x$.

$\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.