#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))))