beta.rkt
```#lang racket/base
(require racket/contract)
(provide/contract
(make-binomial (-> real? procedure?))
(make-beta (-> natural-number/c natural-number/c procedure?))
(make-cumulative-beta (-> natural-number/c natural-number/c procedure?))
(chebyshev-integral (-> natural-number/c natural-number/c procedure?)))

; Make a thunk which randomly returns #t with an expectation of x
(define (make-binomial x)
(lambda ()(if (< (random) x) #t #f)))

; Provide a beta function, a probability distribution function that describes how likely it is that
; the parameter x describes a process that produced t #t results and f #f results.
(define (make-beta t f)
(let [(norm ((chebyshev-integral t f) 1))]
(lambda (x)
(/ (* (expt x t)
(expt (- 1 x) f))
norm))))

; Provide a cumulative beta function, integral of beta from 0 to x.
(define (make-cumulative-beta t f)
(let* [(f (chebyshev-integral t f))
(norm (f 1))]
(lambda (x) (/ (f x) norm))))

; Integrate term x^p * (1-x)^q. To ensure that computation ends (q eventually reaches 0)
; q must be >= 0. To ensure that p doesn't hit 0 from below, p >= 0.
(define (chebyshev-integral p q)
(let* [(f (cond
[(= p 0) ; f = int (1-x)^q dx
(let [(q+1 (add1 q))]
(lambda (x)
(- (/ (expt (- 1 x) q+1)
q+1))))]
[(= q 0) ; f = int x^p dx
(let [(p+1 (add1 p))]
(lambda (x)
(/ (expt x p+1) p+1)))]
[else
; we do integration by parts, int u dv = uv - int v du
; let u =   (1-x)^q       dv = x^p dx
;    du = -q(1-x)^(q-1)dx  v = x^(p+1)/(p+1)
; this results in another integral of the same form,
; allowing an interesting recursion
(let [(p+1 (add1 p))]
(lambda (x)
(/ (+ (* (expt (- 1 x) q)
(expt x p+1))
(* q
; The tricky part - we incorporate the lambda
; returned by the function... into the lambda
; returned by the function!
(apply (chebyshev-integral p+1 (sub1 q))
(list x))))
p+1)))]))
(f0 (f 0))]
(lambda (x) (- (f x) f0))))
```