This is the fourth and final part of a blog post, the first three parts of which can be found here: Generating Samples from Custom Probability Distributions (I), Generating Samples from Custom Probability Distributions (II), and Generating Samples from Custom Probability Distributions (III). We ended the last part with a promise that we would get the efficiency of the random number generation up by a substantial factor from the methods discussed earlier by using a clever technique to subdivide the intervals of definition of the piecewise linear functions defining upper and lower bounds for the probability density function.

Maple uses a beautiful trick for this: it starts with the configuration as we have constructed it in the previous part of this post, and whenever it has to reject a sample value (a red point), or whenever it finds a "blue point" where the lower bound is at least 10% less than the value of , then it subdivides the interval containing that point, at that point. In this way, it is exactly in the areas where many red points would occur that Maple is most likely to adapt the shape of . The reason to start with the interval subdivided into subintervals so that  is all convex or all concave on each subinterval is clear as well now: then there is a clear solution to how to derive an upper and lower bound if such a subinterval is to be split up.

The actual code that Maple uses is written in C++, but for analysis we can reproduce something similar here. This means that there will be substantially more programming than in the versions of the code that we have seen before. Since this code is mostly for analysis and illustration of the ideas, we will sometimes use steps that are not maximally efficient or applicable in all circumstances, but that work for our example case. For example, we assume that the function  provided is twice differentiable on the closed interval that is specified.

We will put the code into a module, so that we don't run into issues with global variable names and the like. The module is setup for a certain distribution by an initial call, then further calls are made to generate the values one by one. (This as opposed to Maple's actual  ?Statistics/Sample call, which generates many values at once.) The intervals are represented by an Array of Records, which have the following keys:

The intervals are initially in order, but whenever an interval is split, the new interval is appended at the end of the Array.

Selection of the interval is done by a procedure called , which is updated after the initial creation of all the interval records and whenever a new interval is added.

Finally,  generates one sample value for the whole distribution. It lets the  select an interval, then asks that interval to select an -value using its  procedure. Then, a corresponding -value is selected, using the interval's  procedure. This is compared to the ; if it is lower, then the -value is returned after the point  is recorded in an Array . If it is higher than the lower bound (so we have a point 'on probation'), then we first split the interval if necessary. Finally, the point is returned if indeed it is lower than the actual value of  there (and  is recorded in an Array ), otherwise the process is started anew (after recording  in an Array ). Without further ado, here is the code:

generator := module()
local
    dispatcher, f, fprime, uniform01,
    createInterval, splitInterval, generateChords, generateTangents, 
    resetDispatcher;
    
export
    intervals, accepted, probation, rejected, setup, generateValue;
    
uses
    ST = Statistics;
        
    generateChords := proc(xleft, xmiddle, xright, yleft, ymiddle, yright)
    local x;
        return unapply(piecewise(
            x < xmiddle, ((x - xleft) * ymiddle +
                          (xmiddle - x) * yleft) / (xmiddle - xleft),
            ((x - xmiddle) * yright +
             (xright - x) * ymiddle) / (xright - xmiddle)), x);
    end proc;
    
    generateTangents := proc(xleft, xmiddle, xright, yleft, yright, 
                             derivativeLeft, derivativeRight)
    local x;
        return unapply(piecewise(
            x < xmiddle, yleft + derivativeLeft * (x - xleft),
            yright + derivativeRight * (x - xright)), x);
    end proc;
    
    setup := proc(ff :: appliable, 
                  a :: numeric, 
                  b :: numeric, $)
    local
        areas,
        fdoubleprime,
        i,
        lastpoint,
        derivativeLeft,
        pointTable,
        points,
        derivativeRight,
        xleft,
        xright,
        yleft,
        yright;
        
        if not a < b then
            error "invalid input, %1 expects its second argument, a, to be " 
            "less than its third argument, but received %2 and %3", 
            procname, a, b;
        end if;
        
        f := ff;
        fprime := D(f);
        fdoubleprime := D(fprime);
        
        uniform01 := ST:-Sample(ST:-RandomVariable(':-Uniform'(0, 1)));
        accepted, probation, rejected := Array(1..0), Array(1..0), Array(1..0);

        # Find all inflection points. We set them as keys in a table,
        # then gather them at the end.
        pointTable := table([b = 1]);
        lastpoint := a;
        
        while lastpoint <> FAIL and lastpoint < b do
            pointTable[lastpoint] := 1;
            lastpoint := RootFinding:-NextZero(fdoubleprime, lastpoint);
        end do;
        
        points := sort([indices(pointTable, ':-nolist')]);
        
        # Now set up the intervals in between these points. An
        # interval will be a record with these keys:
        # * xleft; the left boundary of the interval.
        # * xright; the right boundary of the interval.
        # * yleft, yright; the corresponding values of f.
        # * derivativeLeft, derivativeRight; the values of f' at xleft
        # and xright.
        # * lowerbound; a procedure that returns the lower bound for
        # points in this interval.
        # * upperbound; a procedure that returns the upper bound for
        # points in this interval.
        # * sample; a procedure that generates one sample value from
        # this interval.
        # * area; the total area under the upper curve defined for
        # this interval.

        intervals := Array(1 .. nops(points) - 1);
        
        xleft := points[1];
        yleft := f(xleft);
        derivativeLeft := fprime(xleft);
        for i to nops(points) - 1 do
            xright := points[i+1];
            yright := f(xright);
            derivativeRight := fprime(xright);
            
            intervals[i] := createInterval(xleft, xright, yleft, yright, 
                                           derivativeLeft, derivativeRight);

            xleft := xright;
            yleft := yright;
            derivativeLeft := derivativeRight;
        end do;
        
        resetDispatcher();

        return NULL;
    end proc;
    
    createInterval := proc(xleft, xright, yleft, yright, 
                           derivativeLeft, derivativeRight)
    local
        area,
        areas,
        generators,
        leftRectangleHeight,
        leftTriangleHeight,
        lowerbound,
        rightRectangleHeight,
        rightTriangleHeight,
        sample,
        selector,
        upperbound,
        xmiddle,
        ymiddle;
        
        # The intersection of the two tangents:
        xmiddle := (yright - yleft + derivativeLeft * xleft - derivativeRight * 
                    xright) / (derivativeLeft - derivativeRight);
        ymiddle := f(xmiddle);
        
        if ymiddle < ((xmiddle - xleft) * yright + (xright - xmiddle) * yleft) / 
        (xright - xleft) then
            # Middle point is lower than chord: convex, so chords
            # define upper curve, tangents lower curve.
            lowerbound := generateTangents(xleft, xmiddle, xright, yleft, yright,
                                           derivativeLeft, derivativeRight);
            upperbound := generateChords(xleft, xmiddle, xright, 
                                         yleft, ymiddle, yright);
        else
            # Middle point is higher than chord: concave, so
            # tangents define upper curve, chords lower curve.
            lowerbound := generateChords(xleft, xmiddle, xright, 
                                         yleft, ymiddle, yright);
            upperbound := generateTangents(xleft, xmiddle, xright, yleft, yright,
                                           derivativeLeft, derivativeRight);
            
            # Instead of f at xmiddle, we use the intersection of
            # the tangents for ymiddle in the remainder of this
            # iteration:
            ymiddle := upperbound(xmiddle);
        end if;
        
        leftRectangleHeight := min(yleft, ymiddle);
        leftTriangleHeight := abs(yleft - ymiddle);
        rightRectangleHeight := min(ymiddle, yright);
        rightTriangleHeight := abs(yright - ymiddle);
        
        areas := [leftRectangleHeight * (xmiddle - xleft),
                  leftTriangleHeight * (xmiddle - xleft)/2,
                  rightRectangleHeight * (xright - xmiddle),
                  rightTriangleHeight * (xright - xmiddle)/2];
        
        generators := map(ST:-Sample @ ST:-RandomVariable, [
            ':-Uniform'(xleft, xmiddle),
            ':-Triangular'(
                xleft, xmiddle, `if`(yleft > ymiddle, xleft, xmiddle)),
            ':-Uniform'(xmiddle, xright),
            ':-Triangular'(
                xmiddle, xright, `if`(ymiddle > yright, xmiddle, xright))]);
        selector := ST:-Sample(ST:-RandomVariable(
            ':-ProbabilityTable'(areas / add(x, x in areas))));
        
        sample := proc()
        local i;
            i := trunc(selector(1)[1]);
            return generators[i](1)[1];
        end proc;
        
        area := add(x, x in areas);
        
        return Record(':-xleft' = xleft,
                      ':-xright' = xright,
                      ':-yleft' = yleft,
                      ':-yright' = yright,
                      ':-derivativeLeft' = derivativeLeft,
                      ':-derivativeRight' = derivativeRight,
                      ':-lowerbound' = lowerbound,
                      ':-upperbound' = upperbound,
                      ':-sample' = sample,
                      ':-area' = area);
    end proc;

    resetDispatcher := proc()
    local areas, i, sample;
        areas := [seq(i:-area, i in intervals)];
        sample := ST:-Sample(ST:-RandomVariable(':-ProbabilityTable'(
            areas / add(x, x in areas))));
        dispatcher := proc() 
            return trunc(sample(1)[1]);
        end proc;
    end proc;
    
    generateValue := proc()
    local 
        fx,
        i,
        interval,
        lowerbound,
        x,
        y;
        
        while true do
            i := dispatcher();
            interval := intervals[i];
            
            x := interval:-sample();
            y := uniform01(1)[1] * interval:-upperbound(x);
            lowerbound := interval:-lowerbound(x);
            
            if y < lowerbound then
                accepted(rtable_num_elems(accepted) + 1) := [x, y];
                return x;
            end if;
            
            fx := f(x);
            
            if y > fx or lowerbound < 0.9 * fx then
                splitInterval(i, x, fx);
            end if;
            
            if y < fx then
                probation(rtable_num_elems(probation) + 1) := [x, y];
                return x;
            end if;
            
            rejected(rtable_num_elems(rejected) + 1) := [x, y];
        end do;
    end proc;

    splitInterval := proc(i, x, fx)
    local
        xDerivative,
        oldInterval;

        xDerivative := fprime(x);
        oldInterval := intervals[i];

        intervals[i] := createInterval(
            oldInterval:-xleft, x, oldInterval:-yleft, fx,
            oldInterval:-derivativeLeft, xDerivative);
        intervals(rtable_num_elems(intervals) + 1) := createInterval(
            x, oldInterval:-xright, fx, oldInterval:-yright,
            xDerivative, oldInterval:-derivativeRight);

        resetDispatcher();
    end proc;
end module:

Let us try this out.

We see that less than 0.5% of the attempts to generate a point resulted in evaluating the function , and of the cases in which that was necessary, a bit more than half resulted in rejecting the point under consideration.

The algorithm has done a pretty interesting job in picking areas where it was necessary to subdivide the intervals. The intervals are somewhat shorter where we would expect it: the areas where we saw more than a few red and blue spots in our previous attempt.

The final piecewise linear approximation by the lower and upper bounds is quite good.