#lang racket/base
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; oneg-t-student-en.rkt (rev 0.02 2024.07.03 english)
;;; Development: Enrique Comer-Barragan
;;;              Lab. CEMATi (cemati (dot) org)
;;;              TecNM/ITTijuana
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; ▼For usage examples please see line 139
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; Purpose: basic implementation
;;;   of the t-Student distribution
;;;   with functions for: pdf, cdf & inv-cdf.
;;; Expected usage:
;;;   mainly academic and educational
;;;   (for the current implementation).
;;; Observation: the code is not optimal
;;;              with respect of time or memory,
;;;              and currently is untyped.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; License: Copyright 2024 Enrique Comer-Barragan

;;;   Licensed under the Apache License, Version 2.0
;;;   (the "License"); you may not use this file except
;;;   in compliance with the License.

;;;   You may obtain a copy of the License at
;;;       (http:) //www.apache.org/licenses/LICENSE-2.0

;;;   Unless required by applicable law or agreed to
;;;   in writing, software distributed under the
;;;   License is distributed on an "AS IS" BASIS,
;;;   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
;;;   either express or implied.

;;;   See the License for the specific language
;;;   governing permissions and limitations under
;;;   the License.
;;;;;;;;;;;;;;;;;
  
(require math)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; auxiliary local functions
(define (beta-inc-reg a b x)
  (/ (beta-inc a b x)
     (beta a b)))

;; secant method (elemental)
;; Notice:
;; in the context of the expected
;; applications, we assume convergence.
;; That is: we applied to monotonous
;; increasing functions, of the cdf type,
;; in another context, the return
;; value may be #f.
(define (secant-method f x0 x1 epsi n
                   #:debug? (debug? #f))
  (define-values (f0 f1)
    (values (f x0) (f x1)))
  (define x2
    (if (not (= f0 f1))
        (- x1 (* f1 (/ (- x1 x0)
                       (- f1 f0))))
        #f))
  (when debug?
    (printf "[~a: ~a, ~a; ~a (~a : ~a)]~n"
            n x0 x1 x2 f0 f1))
  (cond
    ((not x2) #f) ; division by zero
    ((< (abs (- x2 x1)) epsi)
     x2)
    ((<= n 0) #f) ; non convergence
    (else (secant-method f x1 x2 epsi
                     (sub1 n)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; function for (centered) t-Student

;; probability density function
;; with ν degrees of freedom
(define (make-pdf-ts ν)
  (λ(t)
    (* (/ (* (sqrt ν)
             (beta 1/2 (/ ν 2))))
       (expt
        (+ 1 (/ (sqr t) ν))
        (* -1/2 (+ ν 1))))))

;; cumulative distribution function
;; with ν degrees of freedom
(define (make-cdf-ts ν)
  (define (cdf-ts t)
    (cond
      ((< t 0)
       (- 1 (cdf-ts (- t))))
      ((= t 0) 0.5)
      (else
       (- 1
          (* 1/2
             (beta-inc-reg
              (/ ν 2)
              1/2
              (/ ν
                 (+ (sqr t)
                    ν))))))))
  cdf-ts)

;; inverse cdf with ν degrees of freedom
;; notice that we use epsilon=1e-12
;;        and maximum number of
;;        iterations in the secant-method
;;        of 30.
(define (make-inv-cdf-ts ν)
  (λ(r)
    (define (f x)
      (- ((make-cdf-ts ν) x) r))
    (secant-method f 1/2 1/3 1e-12 30)))

;; pdf with parameter ν
;; for t as main argument
(define (pdf-ts v t)
  ((make-pdf-ts v) t))

;; cdf with parameter ν
;; for t as main argument
(define (cdf-ts v t)
  ((make-cdf-ts v) t))

;; inv-cdf with parameter ν
;; for r as main argument
;; note: r should satisfy
;;       0 < r < 1
(define (inv-cdf-ts v r)
  ((make-inv-cdf-ts v) r))

;;;;;;;;;;;;;;;;;;;;;;;;
;;; ▼Replace or modify the
;;;   example function call in
;;;   line 163, or write from line 162 on
;;;   the function(s) that you want
;;;   to evaluate in order.
;;;   (after that, you can press
;;;   the [RUN] button, in order to
;;;   evaluate your function(s).
;;;
;;; For example: the following three
;;;   lines:
;;;   (pdf-ts 21 2) ; f_{21}(2)
;;;   (cdf-ts 21 2) ; F_{21}(2)
;;;   (inv-cdf-ts 21 0.7) ; Fi_{21}(0.7)
;;; Return the numbers:
;;;   0.057918057638886446
;;;   0.9706999939756457 
;;;   0.5324551470144957
;;; corresponding to the
;;; pdf, cdf, inv-cdf of the t-Student
;;; distribution with 21
;;; degrees of freedom.
;;;;;;;;;;;;;;;;;;;;;;

(inv-cdf-ts 21 0.7)