Monte Carlo Sampling in Lisp

creativelimits.net:programming/lisp-numerics/monte-carlo/

Monte Carlo Sampling in Lisp

The advantage of CL is that it provides an interactive environment for computation. The Monte Carlo method, which relies on repeated realizations of a random process, is a useful tool for exploring stochastic processes.

The ability to create macros in lisp will allow us to evaluate an arbitrary piece of code that contains a stochastic element. Rather than editing the source and recompiling every time we want to evaluate a new Monte Carlo model we can just write a new function and evaluate it with the macro.

A simple Monte Carlo macro

    (defmacro monte-carlo (code-to-run n)
    "A macro for taking code which has a random component and running the
    same piece of code n times and collecting the results."
    (let ((g (gensym)))
     `(let ((,g nil))
       (do ((i 0 (+ i 1)))
          ((>= i ,n) ,g)
	   (setf ,g
	          (cons
		  ,code-to-run
		  ,g))))))

This macro could be defined with the loop command and the collect keyword.

Example: Computing the expectation of a random product and sum

The goal is to compute the expectation of E(X) * E(Y), where X and Y are deviates from a uniform distribution on the interval [0,1]. The function (random 1.0) will generate a deviate from a uniform distribution on the interval [0,1].

A simple function for computing the mean of a list of values needs to be defined.

    
    (defun mean (nlist) (/ (reduce #'+ nlist) (length nlist)))

    CL-USER> (mean (monte-carlo (* (random 1.0) (random 1.0)) 10))
    0.28285187
    CL-USER> (mean (monte-carlo (* (random 1.0) (random 1.0)) 100))
    0.24429739
    CL-USER> (mean (monte-carlo (* (random 1.0) (random 1.0)) 1000))
    0.24217983
    CL-USER> (mean (monte-carlo (* (random 1.0) (random 1.0)) 10000))

As easy as it was to compute the expectation of the product of two random variables the sum can be computed.

    CL-USER> (mean (monte-carlo (+ (random 1.0) (random 1.0)) 10))
    1.0999103
    CL-USER> (mean (monte-carlo (+ (random 1.0) (random 1.0)) 100))
    1.043648
    CL-USER> (mean (monte-carlo (+ (random 1.0) (random 1.0)) 1000))
    1.0064476
    CL-USER> (mean (monte-carlo (+ (random 1.0) (random 1.0)) 10000))
    0.9966