:

## Generating Samples from Custom Probability Distributions (II)

Maple 14

This is the second post in a four-part series that started with this post: Generating Samples from Custom Probability Distributions (I). We ended that post with the question of how we could improve the algorithm demonstrated there. (The third and fourth posts can be found at Generating Samples from Custom Probability Distributions (III) and Generating Samples from Custom Probability Distributions (IV).)

Typically the most expensive operation in this algorithm is checking whether the -coordinate is greater than , because we have to evaluate at every point. This is not such a big deal for the example function we are using here, but if would, for example, involve a number of special function evaluations or would require running a numerical integration, then this can be a great time saver. One way around evaluating all the time is if we have a good lower bound for that is much faster to evaluate; then we can compare the -coordinate to the lower bound first, and accept it already if it is smaller than the lower bound. Otherwise, the sample value is 'on probation' and needs to be compared to itself. For the upper bound, we can do something similar but better - if we have a tight upper bound for that satisfies the requirements for , then we can use it instead of . This reduces the amount of red surface area in this graph in the previous post; we need to generate fewer points that are to be rejected anyway.

For the upper bound and our new function , we will take this: Then is =   , so looks as follows: This distribution is available in Maple as the distribution. We generate the initial sample from this distribution.   (2)

For the lower bound, we can use  We see that and are indeed lower and upper bounds for . Now we generate the -coordinates of the points in .  (5)

The code that performs the actual accept/reject strategy looks as follows.      (6)

This is much better: for about 70% of the points, we did not even have to evaluate , and only about 1/8 of all points need to be rejected in the end. But we can still do much better. There will be more improvements in the next instalment of this post! 