7
$\begingroup$

I've declared my own probability distribution as follows

dist[a_, b_] := ProbabilityDistribution[Cos[(b x)/2]^2 Sinc[a x]^2, {x, -400, 400}];

When I try to create a random variate using this probability distribution I only get data points in a range (-400, -245). 

data = RandomVariate[dist[200, 500], 10^4];
Histogram[data, 200, "ProbabilityDensity"]

Histogram of prob. dist.

Just to clarify, this isn't an issue with the histogram, when I look at data it does only contain numbers in the range seen above.

Improve this question
$\endgroup$
6
  • $\begingroup$ Your distribution is not define as it is called in the RandomVariate call. $\endgroup$
    – Andy Ross
    Commented Oct 2, 2013 at 12:57
  • $\begingroup$ Thank you Andy, I've edited the question so that it makes sense! $\endgroup$
    – mrkprc1
    Commented Oct 2, 2013 at 13:03
  • 2
    $\begingroup$ Also, don't use D as your distribution name, that is a built-in symbol. As a general rule it is always good to start your own function names in lower case to avoid such conflicts. $\endgroup$
    – Andy Ross
    Commented Oct 2, 2013 at 13:31
  • $\begingroup$ Ok, again this is a problem that comes from copy/pasting straight from mathematica notebook. My distribution is not called D, I'm using the curly D that I could not get on this page, I'll edit again to save confusion! $\endgroup$
    – mrkprc1
    Commented Oct 2, 2013 at 13:35
  • 1
    $\begingroup$ You can gain some insight into what happends by plotting the CDF of your not properly normalised pdf. The CDF ranges from 0 to the integral of the function, not 1. This is useful, you can check CDF[dist, xmax ] == 1 to verify a properly normalized distribution. Now the InverseCDF simply inverts the ill formed CDF, but truncates to the range 0-1. If you go through that exercise you can see the weird behavior makes some sense. $\endgroup$
    – george2079
    Commented Oct 2, 2013 at 16:09

1 Answer 1

Reset to default
10
$\begingroup$

A pdf is only well-defined if it integrates to unity over the domain of support. You have set up the domain of support as (-400, 400). You would have to check what values of $a$ and $b$ (if any) are appropriate for that domain of support. But, more to the point, your domain of support is ill-defined for the parameter values you have provided:

f = Cos[(b x)/2]^2 Sinc[a x]^2;

Integrate[ f /. {a -> 200, b -> 500}, {x, -400, 400}] // N

0.00785395

If you set up the domain of support on the real line {x,-Infinity, Infinity}:

Integrate[f /. {a -> 200, b -> 500}, {x, -Infinity, Infinity}]

Pi/400

... you just need to multiply your pdf by 400/Pi, and you are all set (given those parameter values for $a$ and $b$).

Improve this answer
$\endgroup$
4
  • $\begingroup$ Ok, so if I change my probability distribution to D[a_, b_, A_]:=ProbabilityDistribution[ (Cos[(b x)/2]^2 Sinc[a x]^2)/Integrate[Cos[(b x)/2]^2 Sinc[a x]^2, {x, -A, A}], {x, -A, A}]; should that then work since it guarantees normalisation for any range A? $\endgroup$
    – mrkprc1
    Commented Oct 2, 2013 at 13:30
  • $\begingroup$ Wow "The integral of the PDF over the distribution domain needs to be unity:" is burried in the docs under "properties and relations". Poor documentation considering how ill behaved the function is if you dont precisely meet that criteria. $\endgroup$
    – george2079
    Commented Oct 2, 2013 at 15:23
  • 3
    $\begingroup$ @george2079 Pardon me, but this shouldn't be in the docs at all, as this is elementary maths. According to the law of total probability, for a "set of pairwise disjoint events whose union is the entire sample space" the sum of all probabilities should be 1. (For continuous distributions, replace sum with integral.) $\endgroup$
    – István Zachar
    Commented Oct 2, 2013 at 21:04
  • 1
    $\begingroup$ It should be clearly documented that it doesn't check for validity. $\endgroup$
    – george2079
    Commented Oct 3, 2013 at 12:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.