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
(lambda (x)
(- (/ (expt (- 1 x) q+1)
q+1))))]
[(= q 0) ; f = int x^p dx
(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