Message Boards Message Boards

0
|
6289 Views
|
8 Replies
|
0 Total Likes
View groups...
Share
Share this post:

[?] Speed up RandomVariate using custom distribution?

Posted 7 years ago

I am trying to use the following code to obtain random samples according to custom distribution, but Mathematica is very often not able to give it in reasonable time. One random sample takes either between 0,015-0,05 seconds or it runs without the result for more than one minute at which point I stop the calculation.

nlm2[x_] := -602.989505885402` - 61.69262610012762` x + 1.2379005632470046` x^2 - 0.013278068455723385` x^3 + 0.00008027261407746941` x^4 - 2.592723712534534`* 10^-7 x^5 +3.4947822352866564` * 10^-10 x^6

Prob[DA_] :=Exp[-(nlm2[DA] + 1887.5)*1063.2151454471323`]/560.6805023451591`

Dist = ProbabilityDistribution[Prob[x], {x, 110, 145}]

RandomVariate[Dist]
POSTED BY: Michal Jex
8 Replies
Posted 7 years ago

Good. I wish I knew what the exact problem is. (I get the same symptoms for your density: first random sample seems to work but then stalls forever for any subsequent random samples.)

POSTED BY: Jim Baldwin
Posted 7 years ago

I suspect that the issue is that of "round off error" because of the large numbers involved (such as 145^6 = 9294114390625) when an attempt (internally) is made to match a uniform random number to an approximation of the cumulative distribution function. (Maybe scaling things using Horner's scheme might be helpful to reduce the round-off error.)

One approximation to get around this is to generate a list of values for the cumulative distribution function (with the associated values of DA) and then use interpolation. One can make the approximation more accurate but at the expense of more time creating the interpolation table. Here's an example:

cumulative = Table[{NIntegrate[Prob[x], {x, 110, t}], t}, {t, 110, 145, 1/100}];
rSample = Interpolation[cumulative];
Histogram[rSample[#] & /@ RandomReal[{0, 1}, 1000]]

Histogram of random sample

POSTED BY: Jim Baldwin
Posted 7 years ago

Thank you. Your suggestion, how to bypass the problem works.

POSTED BY: Michal Jex
Posted 7 years ago

It would be better to edit your original question to put in the missing characters and use the formatting tools. You might also want to add in what you mean by "too slow".

POSTED BY: Jim Baldwin
Posted 7 years ago

Thank you for the suggestion. I modified the original question.

POSTED BY: Michal Jex
Posted 7 years ago

(Comment removed after question was modified.)

POSTED BY: Jim Baldwin
Posted 7 years ago

I tried to plot the probability density and it seems to be working properly.enter image description here

POSTED BY: Michal Jex
Posted 7 years ago

The forum dropped stars between the numbers and hats so in the expression of the polynom there is 0.00008027261407746941 x^4 - 2.592723712534534 * 10^-7 x^5 +3.4947822352866564` * 10^-10 x^6

POSTED BY: Michal Jex
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract