Exercise 3.5

This exercise has essentially two parts: the integration procedure itself, and then the means of applying it (which will also test that it functions correctly). I’ll discuss them separately.

Since we are using the Monte Carlo approach for integration, we first have to set up the parameters of a trial. We are given a hit/miss test, which is a procedure that returns true or false depending on whether the point is contained under the function or not. We pick random points using random-in-range with the range being x- and y-values of the bounding box. Then we run that through monte-carlo. The result will give us the fraction of points that were ‘hits’ in the defined region. As we test more points, the ratio of the valid points to the total points tested by the Monte Carlo trials approaches the ratio of the function’s area to the area of the bounding box. To get the estimated area of our function, we multiply the fractional ratio by the bounding box’s area, which comes from the rectangle defined by the corner points. Those steps will yield the desired integration result, and here is the code:

(define (estimate-integral P x1 x2 y1 y2 trials)
  (let ((total-area (* (- x2 x1) (- y2 y1)))
        (hit-fraction (monte-carlo trials (lambda () 
                                            (P (random-in-range x1 x2)
                                               (random-in-range y1 y2)
                                               )
                                            )
                                   )
                      )
        )
    (* hit-fraction total-area)
    )
  )

The application/testing portion is asking us to estimate the value of π by using the unit circle. Since our integrator simply measures the area given certain constraints, and π is part of the formula for the area of a circle, all we need to determine π is to run our integrator over a circle, and divide out the non-π portion of the area formula (i.e. radius squared). The bounding box can be set by adding and subtracting one radius length to the center of the circle, as that gives us a square that circumscribes the circle, enclosing it tightly.

We then need a proper predicate to test whether a point is within a circle. Given a point (x, y), we can determine whether it’s in a circle by checking to see if the point’s distance from the circle’s center is less than the radius of the circle. I wrote a generic procedure that can accept circles of any radius centered on other locations (though the variables are set external to the test itself). This does mean there is more computation involved for the test, but the procedure runs fast enough for me that the extra time is barely noticeable.

(define circleP (lambda (x y)
                  (< (+ (sqr (- x xc)) (sqr (- y yc))) r_squared) 
                  )
  )

This is a simple application of the distance formula, but to save time in the calculation, the value is compared to the square of the radius instead of taking the square root of the other side of the equation.

In testing, I added a few variations to see what effect it would have on the estimation. One thing that can change is the size of the circle. This doesn’t have a noticeable effect on how many trials it takes to get an accurate estimate. That make sense, since all that happens is both the estimated value and the expected value are scaled by some amount (though you would be limited at some point by the floating-point number representation).

The other variation was to change the size of the bounding box relative to the circle. This results in a typically less accurate estimate for the same number of trials, when compared to the tighter bounding box. The effect here is not too surprising, either. The Monte Carlo test grows more accurate as the proportion of hits approaches the actual proportion of the areas. With a lower chance to hit (because of a smaller ‘target’ area), the numbers involved in computing the ratio are smaller, and thus the result will be less accurate.

Exercise 3.6

Structurally, the random generator needs to first determine whether it is generating or in reset mode. It can be set up similarly to the bank account procedure of Section 3.1.1, with dispatch to locally defined functions depending on the argument. There is a slight difference in the dispatch method, however — rand might return either a procedure or a number, depending on the command.

In the back account example, the internal state variable was defined within make-account, which returned a procedure that could be called with the commands; those actions were limited solely to that account (i.e. procedure). We can use a similar approach here. The ‘internal state variable’ referred to in the text needs to be saved in some way if calls to rand are to work as specified.

One way would be to copy the bank account procedure, and create a pseudo-random number generator generator (and I think it’s good to shorten the first part to PRNG, so we can call it a PRNG Generator). Then rand would just send calls to itself to the generated PRNG. That is what I did in my first solution, but to demonstrate an alternative way of doing things, I included another approach (I call these ‘Method 1’ and ‘Method 2’ in the solution file, although Method 2 uses a Racket-specific routine). Only Method 1 is shown below.

(define (make-prng a b m s)
  (let ((random-seed s)
        )
    
    (define (rand-update)
      (set! random-seed (modulo (+ (* a random-seed) b) m))
      )
    
    (define (generate) 
      (rand-update)
      random-seed
      )  
    
    (define (reset new-value)
      (set! random-seed new-value) 
      )
    
    (define (dispatch arg)
    (cond
      ((eq? arg 'generate)  (generate))
      ((eq? arg 'reset) reset)
      (else
       (error "Unknown argument for rand :" arg)
       )
      )
      )
    dispatch
    )
  )

The definition of rand (in Method 1) resembles the definition of acc in the previous section — given a procedure make-prng that returns a procedure, rand is a named variable using that procedure’s result. Internally, calls to generate and reset are handled slightly differently, since one returns a value, and the other returns a procedure that accepts the new value; this is taken care of within the dispatch portion.

As for what is actually stored in the internal variable, it is the last random value, which can be sent to a procedure like rand-update that will generate the next value. My version has an internal rand-update, whose parameters are determined at the time make-prng is called. Those parameters are used by a LCG (Linear Congruential Generator), the method mentioned in the text’s footnote. While it’s not a great means for simulating randomness, it is fairly easy to implement and work with. Here’s a brief explanation of the formula, which also includes the requirements to make a ‘good’ sequence by choosing the right parameters.

The second method creates the PRNG in a different way. It turns out that Racket/Pretty Big actually has a generator already, and it’s called make-pseudo-random-generator. Calls to random (the built-in function) can also specify which PRNG to use by adding it as an optional final argument. Other Scheme implementations may have a similar procedure with an alternate name (MIT Scheme, for example, has make-random-state).

To reset the state of the PRNG in this method, I pass an entire PRNG in a known state as the reset value. It uses another procedure that copies the state of the PRNG without changing it (converting it to ‘vector’ form and then back). We can’t simply assign the PRNG argument directly, since otherwise any external call to the PRNG we passed would affect the rand function, and destroy the ‘known state’ it had when the reset was done. This does complicate how one chooses a new state in some situations, although the text implies more concern about repeating a sequence than about picking reset states. As long as the PRNG isn’t used elsewhere, it will remain in the state it was (thus allowing it to be passed again as a reset seed).

In Method 2, rand is the only procedure used (and thus it has no dispatch). It uses an external (global) variable, which rand relies on but can also modify. While that makes the code a tiny bit cleaner and simpler, the state of this version is not ‘local’ and confined. For random values, that may not be an issue, although it could lose the repeatibility we might want from a PRNG that the reset function took pains to preserve. However, it would not be that hard to continue relying on the built-in procedure and make another version, combining the built-in PRNG Generator with an approach using local state and dispatch.

Solutions 3.1.2