:

## 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!