[racket] adaptive simpson integration

On 01/21/2012 10:59 AM, Matthias Felleisen wrote:
>
> I have coded up the adaptive Simpson integration algorithm below,
> including a contract. The algorithm has been mentioned recently on
> 'users'. Any comments from the numerically inclined? -- Matthias

I have a way to visualize integration and more test cases. Put this at
the end of your code for now, and run it:


(require plot plot/utils (only-in plot/plot2d/kde kde))

;; integrates by summing the areas of a piecewise-linear approximation
;; based directly on Riemann's integral definition
(define ((naive-integrate f [steps 30]) x-min x-max)
   (define step-size (/ (- x-max x-min) (- steps 1)))
   (define xs (build-list steps (λ (n) (+ x-min (* n step-size)))))
   (define ys (map f xs))
   (for/fold ([sum 0] [funs empty] [mins empty] [maxs empty]
                      ) ([x1 xs] [x2 (rest xs)] [y1 ys] [y2 (rest ys)])
     (define dx (- x2 x1))
     (define dy (- y2 y1))
     (define y-avg (* 1/2 (+ y1 y2)))
     (define new-sum (+ sum (* dx y-avg)))
     (define new-fun (λ (x) (+ (* (/ (- x x1) dx) dy) y1)))
     (values new-sum
             (cons new-fun funs)
             (cons x1 mins)
             (cons x2 maxs))))

(define (do-integration-test f x-min x-max)
   (define simpson-f-area (time ((adaptive-simpson f) x-min x-max)))
   (define-values (naive-f-area f-funs f-mins f-maxs)
     (time ((naive-integrate f) x-min x-max)))
   (values
    simpson-f-area
    naive-f-area
    (plot (list (for/list ([f f-funs] [x-min f-mins] [x-max f-maxs]
                                      [color (in-cycle '(0 3))])
                  (function-interval (λ (x) 0) f x-min x-max
                                     #:color color))
                (function f #:style 'long-dash #:width 2))
          #:x-min x-min #:x-max x-max)))

(printf "step~n")
(printf "~v (actual value)~n" (+ (- 1 (cos 1)) 1))
(do-integration-test step 0 2)

(printf "sin~n")
(printf "~v (actual value)~n" (- 1 (cos 2)))
(do-integration-test sin 0 2)


(define-values (est-density x-min x-max)
   (let ([xs  (build-list 50 (λ _ (random)))])
     (kde xs 0.1)))

(printf "estimated density~n")
(printf "~v (actual value)~n" 1.0)
(do-integration-test est-density x-min x-max)


The code compares your Simpson implementation with a naive
implementation that sums the areas of uniformly spaced linear
approximations. It also plots the function and the linear
approximations. The `naive-integrate' function demonstrates how you
might instrument your Simpson code to do the same.

On "step" and "sin", your Simpson implementation is kicking naive's
butt, beating it by 4-5 decimal digits in only 2x-4x the time. Nice.

The "estimated density" case is intriguing. I'm a little surprised that
the naive integration algorithm does so well (generally about 10 digits 
correct), but that probably has some central-limit-y explanation. I'm 
most surprised that the Simpson implementation fails so spectacularly on 
what looks like a smooth function, getting only one digit right and 
taking 50x the time the naive algorithm takes.

My best guess: The estimated density may actually not be smooth enough.
(If you click and drag repeatedly to zoom way into part of an estimated
density plot, you'll eventually find places it isn't differentiable or 
continuous.) An estimated density is a sum of Gaussians, each placed 
over a sample point. It doesn't sum each Gaussian at each point, just 
the ones that contribute "enough" height to the function at that point. 
It looks smooth at reasonable scales, but it isn't.

So I tried setting the Simpson epsilon parameter to something low, but 
it didn't seem to help with accuracy, only speed. Curiouser and curiouser.

Finally, here are some test cases that fail because of division by zero:

(check-equal? ((adaptive-simpson (λ (x) 0.0)) 0.0 1.0)
               0.0)

(check-equal? ((adaptive-simpson (λ (x) 0)) 0 1)
               0)

Neil T

>
>
>
>
> #lang racket
>
> (require rackunit)
>
> ;;
> -----------------------------------------------------------------------------
>
>
;; interface, with contract
>
> ;; Simpson's rule for approximating an integral (define (simpson f L
> R) (* (/ (- R L) 6) (+ (f L) (* 4 (f (mid L R))) (f R))))
>
> (define (mid L R) (/ (+ L R) 2.))
>
> ;; adaptive Simpson limits approximate error to 15 ε (define
> (simpson-error-estimate f L R ε) [define M (mid L R)] (<  (abs (+
> (simpson f L M) (+ simpson f M R) (- (simpson f L R)))) (* 15 ε)))
>
> (provide/contract [adaptive-simpson (->i ((f (->  real? real?)))
> (#:epsilon (ε real?)) (r (f ε) (->i ((L real?) (R  (L) (and/c real?
> (>/c L)))) () (r real?) ;; adds six function calls to f, which seems
> tolerable #:post (L R) (simpson-error-estimate f L R ε))))])
>
> ;;
> -----------------------------------------------------------------------------
>
>
;; implementation
>
> (define (adaptive-simpson f #:epsilon [ε .000000001]) (lambda (L R)
> (define f at L (f L)) (define f at R (f R)) (define-values (M f at M whole)
> (simpson-1call-to-f f L f at L R f at R)) (asr f L f at L R f at R ε whole M
> f at M)))
>
> ;; computationally efficient: 2 function calls per step (define (asr
> f L f at L R f at R ε whole M f at M) (define-values (leftM  f at leftM  left*)
> (simpson-1call-to-f f L f at L M f at M)) (define-values (rightM f at rightM
> right*) (simpson-1call-to-f f M f at M R f at R)) (define delta* (- (+
> left* right*) whole)) (cond [(<= (abs delta*) (* 15 ε)) (+ left*
> right* (/ delta* whole 15))] [else (define epsilon1 (/ ε 2)) (+ (asr
> f L f at L M f at M epsilon1 left*  leftM  f at leftM) (asr f M f at M R f at R
> epsilon1 right* rightM f at rightM))]))
>
> (define (simpson-1call-to-f f L f at L R f at R) (define M (mid L R))
> (define f at M (f M)) (values M f at M (* (/ (abs (- R L)) 6) (+ f at L (* 4
> f at M) f at R))))
>
> ;; simplistic prototype: many calls to f per step (define (asr.v0 f L
> R ε whole) (define M  (mid L R)) (define left*  (simpson f L M))
> (define right* (simpson f M R)) (define delta* (- (+ left* right*)
> whole)) (cond [(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta*
> whole 15))] [else (define epsilon1 (/ ε 2)) (+ (asr.v0 f L M epsilon1
> left*) (asr.v0 f M R epsilon1 right*))]))
>
> ;;
> -----------------------------------------------------------------------------
>
>
;; examples
> (define int-const=5  (adaptive-simpson (lambda (x) 5))) (check-equal?
> (int-const=5 0 1) 5.0) (check-equal? (int-const=5 0 10) 50.0)
>
> (define int-identity (adaptive-simpson values)) (check-equal?
> (int-identity 0 1) .5)
>
> (define step (lambda (x) (if (<  x 1) (sin x) 1))) (define int-step
> (adaptive-simpson step))
>
> (int-step 0 2)
>
>
> ____________________ Racket Users list:
> http://lists.racket-lang.org/users