#lang racket/base
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; oneg-noncentral-ts-dist-en.rkt
;;; (rev 0.01 2024.07.07 english)
;;; Development: Enrique Comer-Barragan
;;;              Lab. CEMATi (cemati (dot) org)
;;;              TecNM/ITTijuana
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; ▼For usage examples please see line 226
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; Purpose: **basic implementation**
;;;   of the noncentral 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)

(provide make-pdf-ncts
         make-cdf-ncts
         make-inv-cdf-ncts
         pdf-ncts
         cdf-ncts
         inv-cdf-ncts)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; auxiliary local functions

;; regularized lower incomplete
;;  beta function
(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
;; 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)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Making of the PDF for noncentral t-Student Dist.
;; Formula source:
;; (https:)//pure.manchester.ac.uk/ws/portalfiles/
;;           portal/78743532/studentt.pdf
(define (make-pdf-ncts ν μ)
  (λ(x)
    (define (c x)
      (/ (* (exp (* -1/2 (sqr μ)))
            (expt ν (/ ν 2)))
         (* (sqrt pi)
            (expt (+ ν (sqr x))
                  (* 1/2 (+ ν 1)))
            (gamma (/ ν 2)))))
    (define (kth-term k)
      (/ (* (gamma (/ (+ ν k 1)
                      2))
            (expt μ k)
            (expt 2 (/ k 2))
            (expt x k))
         (* (gamma (+ k 1))
            (expt (+ ν (sqr x))
                  (/ k 2)))))
    ; Aproximation of inf summation
    ; first implementation:
    ; recursive & untyped
    (define (Sm m Sk k epsi)
      (define Snext
        (+ Sk (kth-term (add1 k))))
      #;(printf "Sk[~a]=~a~n"
              k Sk)
      (cond
        ((< (abs (- Snext Sk))
            epsi)
         Snext)
        ((>= k m) Snext)
        (else
          (Sm m Snext (add1 k) epsi))))
    ; current summation precision 1e-8
    (* (c x) (Sm 40
                 (kth-term 0) 0
                 1e-8))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; cumulative distribution function
;; with ν degrees of freedom
;; and μ non central parameter
;; Formula source:
;; (https:)//en.wikipedia.org/wiki/Noncentral_t-distribution
;; (current first formulation,
;;          the 2nd, doesn't work as expected)
(define (make-cdf-ncts ν μ)
  (define M
    (exp (* -1/2 (sqr μ))))
  (λ(x)
    (define (jth-term j x μ)
      (define ej (expt (* 1/2 (sqr μ))
                       j))
      (define pj
        (* (/ (factorial j))
           ej))
      (define qj
        (* (/ μ
              (* (sqrt 2)
                 (gamma (+ j 3/2))))
           ej))
      (define y (/ (sqr x)
                   (+ (sqr x) ν)))
      (+ (* pj
           (beta-inc-reg (+ j 1/2)
                       (/ ν 2)
                       y))
         (* qj
            (beta-inc-reg (+ j 1)
                          (/ ν 2)
                          y))))
    ; Aproximation of inf summation
    ; first implementation:
    ; recursive & untyped
    (define (Sm m Sj j epsi)
      (define Snext
        (+ Sj (jth-term (add1 j) x μ)))
      #;(printf "Sj[~a]=~a~n"
              j Sj)
      (cond
        ((< (abs (- Snext Sj))
            epsi) Snext)
        ((>= j m) Snext)
        (else
          (Sm m Snext (add1 j) epsi))))
    ; current summation precision 1e-8
    (define (G x μ)
      (+ (cdf (normal-dist 0 1)
              (- μ))
         (* 1/2 M (Sm 40
                      (jth-term 0 x μ)
                      0
                      1e-8))))
    (if (not (negative? x))
        (G x μ)
        (- 1 (G x (- μ))))))

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; pdf with parameters ν, μ
;; for t as main argument
(define (pdf-ncts ν μ t)
  ((make-pdf-ncts ν μ ) t))

;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; cdf with parameters ν, μ
;; for t as main argument
(define (cdf-ncts ν μ t)
  ((make-cdf-ncts ν μ) t))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; inv-cdf-ncts with parameters ν, μ
;; for t as main argument
;; note: r should satisfy
;;       0 < r < 1
(define (inv-cdf-ncts ν μ r)
  ((make-inv-cdf-ncts ν μ) r))

;;;;;;;;;;;;;;;;;;;;;;;;
;;; ▼Replace, comment or modify the
;;;   example function calls in
;;;   line 255, or write from line 255 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 
;;;   lines:
;;;   (pdf-ncts 4 2 3.5) ; f_{4,2}(3.5)
;;;   (cdf-ncts 4 2 3.5) ; F_{4,2}(3.5)
;;;   (inv-cdf-ncts 4 2 0.75) ; Fi_{4,2}(0.75)
;;; Return the numbers:
;;;   0.13858583179931994
;;;   0.7930150994699615
;;;   3.218631407678812 
;;; corresponding to the
;;; pdf, cdf, inv-cdf of the noncentral t-Student
;;; distribution with  ν=4, μ=2
;;; degrees of freedom.
;;; Compare with Wolfram alpha call:
;;; PDF[NoncentralStudentTDistribution[4,2],3.5]
;;; and WxMaxima code:
;;; load("distrib");
;;; cdf_noncentral_student_t (3.5,4,2)
;;; quantile_noncentral_student_t(0.75,4,2)
;;;;;;;;;;;;;;;;;;;;;;
(pdf-ncts 4 2 3.5) 
(cdf-ncts 4 2 3.5)
(inv-cdf-ncts 4 2 0.75)