;; 
;;
;; probdist documentation for the Chicken Scheme module system.
;;
;; This program is free software: you can redistribute it and/or
;; modify it under the terms of the GNU General Public License as
;; published by the Free Software Foundation, either version 3 of the
;; License, or (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful, but
;; WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;; General Public License for more details.
;;
;; A full copy of the GPL license can be found at
;; <http://www.gnu.org/licenses/>.
;;
;;

(use eggdoc)

(define doc
  `((eggdoc:begin
     (name "probdist")
     (description "Probability density functions.")
     (author (url "http://chicken.wiki.br/users/ivan-raikov" "Ivan Raikov"))

     (history 
      (version "1.5" "Ported to Chicken 4.")
      (version "1.4" "Added some error checking in make-normal-pdf and a record printer for the normal-pdf datatype.")
      (version "1.3" "Build script updated for better cross-platform compatibility")
      (version "1.2" "More metafile fixes")
      (version "1.1" "Metafile fixes")
      (version "1.0" "Initial release"))

     (requires (url "datatype.html"     "datatype")
	       (url "syntax-case.html"  "syntax-case")
	       (url "srfi-4-utils.html" "srfi-4-utils")
	       (url "matrix-utils.html" "matrix-utils")
	       (url "blas.html"         "blas")
	       (url "atlas-lapack.html" "atlas-lapack"))

     (usage  (ul (li "(require-extension sampled-pdf)")
		 (li "(require-extension normal-pdf)")))

     (download "probdist.egg")

     (documentation
      (p "probdist is a library of routines to compute univariate or multivariate "
	 "normal probability function for a given mean and covariance; and to approximate "
	 "an arbitrary probability distribution as a weighted set of samples (sampled PDF). ")
      (p "This code is based on code in the " (tt "dysii") " project by Lawrence Murray. ")
      (p "The sampled PDF algorithm is based on the paper by Kitagawa: " 
	 (tt "_Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear State Space Models_,"
	     " Journal of Computational and Graphical Statistics, 1996, 5, 1-25") ". ")

       (subsection "Normal PDF"

		   ,(eggdoc:make-defsig 
		     'datatype "(define-datatype normal-pdf normal-pdf? (NormalPDF S mu sigma R sigmaInv sigmaDet ZI mu-zero? ))"
		     `(p "A representation of a normal PDF:"
			 (symbol-table
			  (describe "S"      "number of dimensions")
			  (describe "mu"     "expectation")
			  (describe "sigma"  "covariance")
			  (describe "R"      "upper triangular Cholesky decomposition of covariance matrix")
			  (describe "sigmaInv"  "inverse of covariance matrix")
			  (describe "sigmaDet"  "determinant of covariance matrix")
			  (describe "ZI"        ("the constant factor " (tt "(2*pi)^(-S/2) / prod(diag(R))")))
			  (describe "mu-zero?"  "true if expectation is zero"))))

		   
		   (procedure "make-normal-pdf:: S * MU * SIGMA -> NORMAL-PDF"
			      (p "Creates a new normal PDF object with the specified dimension, expectation, and "
				 "covariance. If the dimension " (tt "S") " is 1, then the expectation "
				 (tt "MU") " and the covariance " (tt "SIGMA") " must be scalars; "
				 "otherwise, they must be SRFI-4 " (tt "f64") " vectors. "
				 "In the latter case, " (tt "MU") " must be a vector of size " (tt "S") ", "
				 "and " (tt "SIGMA") " must be a matrix of size " (tt "S x S") ". "))
		   
		   (procedure "normal-pdf:density:: NORMAL-PDF * X -> DENSITY"
			      (p "Computes the density of the distribution at the given point " (tt "X") ". "
				 "If the dimension " (tt "S") " of the distribution is 1, then " (tt "X") 
				 " must be a scalar, otherwise it must be an SRFI-4 " (tt "f64") " vector "
				 "of length " (tt "S") ". "))
		   
		   (procedure "normal-pdf:sample:: NORMAL-PDF * X -> SAMPLE"
			      (p "Sample from the distribution. Let " (tt "X") " be a sample from a "
				 "standard normal distribution. Then the sample is " (tt "SAMPLE = MU + R*X")
				 ", where " (tt "R") " is the upper triangular Cholesky decomposition "
				 "of the covariance matrix. "))
				  
		   (procedure "normal-pdf:expectation:: NORMAL-PDF -> MU"
			      (p "Return the expectation value of the distribution (scalar or f64vector). "))

		   (procedure "normal-pdf:covariance:: NORMAL-PDF -> SIGMA"
			      (p "Return the covariance value of the distribution (scalar or f64vector). ")) )

       (subsection "Sampled PDF"

		   ,(eggdoc:make-defsig 
		     'datatype "(define-datatype sampled-pdf sampled-pdf? (SampledPDF N W WXS MU SIGMA))"
		     `(p "A representation of a weighted sampled PDF:"
			 (symbol-table
			  (describe "N"     "sample set size")
			  (describe "W"     "total weight of the sample set")
			  (describe "wxs"   "weighted samples (see below)")
			  (describe "mu"    "expectation")
			  (describe "sigma" "covariance"))))
			  
		   ,(eggdoc:make-defsig 
		     'datatype "(define-datatype weighted-samples weighted-samples? (WeightedSamples S senv make-senv))"
		     `(p "A representation of weighted samples:"
			 (symbol-table
			  (describe "S"      "number of dimensions")
			  (describe "senv"   ("an ordered dictionary data structure that follows the API of "
					      (url "rb-tree.html" "rb-tree")))
			  (describe "make-senv" "a procedure to create an empty samples dictionary"))))


		   (procedure "make-sampled-pdf:: S * MAKE-SENV * XS [* WS * XCAR * XCDR * XNULL?] -> SAMPLED-PDF"
			      (p "Creates a new sampled PDF object with the specified dimension, "
				 "samples environment creation procedure that follows the API of make-rb-tree, "
				 "and weighted samples. "
				 "The samples can be specified in one of two ways. If both " (tt "XS") " and "
				 (tt "WS") " are provided, then " (tt "XS") " is a sequence that contains the "
				 "samples, and " (tt "WS") " is a sequence that contains the corresponding weights. "
				 "If only " (tt "XS") " is provided, or " (tt "WS") " is false, then " (tt "XS")
				 " must consist of cons cells, where the car is the weight, and the cdr is the sample. "
				 "Optional arguments " (tt "XCAR XCDR XNULL?") " could be used for sequences other "
				 "than lists (e.g. streams). "))

		   (procedure "sampled-pdf:expectation:: SAMPLED-PDF -> MU"
			      (p "Return the expectation value of the distribution (scalar or f64vector). "))

		   (procedure "sampled-pdf:covariance:: SAMPLED-PDF -> SIGMA"
			      (p "Return the covariance value of the distribution (scalar or f64vector). "))

		   (procedure "sampled-pdf:find-sample:: U * SAMPLED-PDF -> X"
			      (p "Given a number " (tt "u in [0,1]") ", returns the sample " (tt "x(i)")
				 "such that " (tt "\\sum_{j=1}^{i-1} w_j < Wu <= \\sum_{j=1}^{i} w_j") ", "
				 "where " (tt "W") " is the total weight of the sample set. "))

		   (procedure "sampled-pdf:normalize:: SAMPLED-PDF -> SAMPLED-PDF"
			      (p "Normalizes the sample weights to sum to 1."))
		   
		   (procedure "sampled-pdf:resample:: MAKE-SENV * SAMPLED-PDF [ * M] -> SAMPLED-PDF"
			      (p "Resamples the distribution. This produces a new approximation "
				 "of the same distribution using a set of equally weighted sample "
				 "points. Sample points are selected using the deterministic "
				 "resampling method given in the appendix of Kitagawa. "
				 "Optional argument " (tt "M") " is the number of samples to take "
				 "for the new distribution. If not given, defaults to the number of "
				 "samples in the existing distribution.")) ))

     (examples (pre #<<EOF
;;
;; Multivariate test of normal-pdf.
;;
;; This test creates a 3-dimensional normal distribution with a
;; random mean and variance. It then samples from this distribution
;; and calculates the sample mean and variance for comparison.
;;

(require-extension srfi-1)
(require-extension srfi-4)
(require-extension srfi-4-utils)
(require-extension matrix-utils)
(require-extension blas)
(require-extension normal-pdf)
(require-extension random-mtzig)

(define-utility-matrices f64)

;; Number of samples 
(define N 100000)

;; Number of dimensions (sample size)
(define S 3)


(define (f64vector-scale a x)
  (let ((n (f64vector-length x)))
    (blas:dscal n a x)))

(define (f64vector-sum x y)
  (let ((n (f64vector-length x)))
    (blas:daxpy n 1.0 x y)))

(define (f64vector-sub x y)
  (let ((n (f64vector-length x)))
    (blas:daxpy n -1.0 y x)))

(define (f64vector-mul a b)
  (f64vector-map (lambda (v1 v2) (* v1 v2)) a b))


(define f64matrix-map (make-matrix-map f64vector-ref f64vector-set!))

(define (upper M N A)
  (f64matrix-map blas:RowMajor M N A (lambda (i j v) (if (>= j i) v 0.0))))

;;
;;  Computes the arithmetic mean of a list of f64 vectors using the
;;  recurrence relation mean_(n) = mean(n-1) + (v[n] -
;;  mean(n-1))/(n+1)
;;
(define (mean s vs)
    (let loop ((i 0) (vs vs) (mean (make-f64vector s 0.0)))
      (if (null? vs)  mean
	  (loop (fx+ 1 i) (cdr vs)
		(let ((d (f64vector-sub (car vs) mean)))
		  (f64vector-sum mean (f64vector-scale (/ 1 (+ 1 i)) d)))))))

(define (covariance s vs mean)
   (let* ((cov  (matrix-zeros S S)))
     (let ((ds (map (lambda (x) (f64vector-sub x mean)) vs)))
       (let ((n (fold (lambda (d n) 
			(blas:dsyr! blas:RowMajor blas:Upper s 1.0 d cov)
			(fx+ 1 n)) 0 ds)))
	 (f64vector-scale (/ 1 (- n 1)) cov)))))
			 
      
;; Given a vector of values in the range [0..1], scale all values in
;; the vector to the specified range.
(define (f64vector-urange x lo hi)
  (if (< lo hi) (let ((d (- hi lo)))
		  (f64vector-map (lambda (x) (+ lo (* d x))) x))))

(define rng  (random-mtzig:init))

;;  distribution parameters
(define mu     (f64vector-urange (random-mtzig:f64vector-randu! S rng) -5 5))
(define sigma  (let ((x (f64vector-urange (random-mtzig:f64vector-randu! (* S S) rng) 0 5)))
		 (blas:dgemm! blas:RowMajor blas:NoTrans blas:Trans S S S 
			      1.0 x x 0.0 (make-f64vector (* S S)))))


(define gpdf  (begin
		(print "mu = " mu)
		(print "sigma = " sigma)
		(make-normal-pdf S mu sigma)))

(define input (list-tabulate N (lambda (i) (random-mtzig:f64vector-randn! S rng))))

(define samples
  (let loop ((i 0) (input input) (samples (list)))
    (if (null? input)  samples
	(let ((samples (cons (normal-pdf:sample gpdf  (car input))  samples)))
	  (loop (fx+ 1 i) (cdr input) samples)))))

(if (<= N 100) (print "input = " input))
(if (<= N 100) (print "samples = " samples))

;; test for mean and covariance equality to one decimal place
(equal? (f64vector-map (lambda (x) (truncate (* 10 x))) mu)
	(f64vector-map (lambda (x) (truncate (* 10 x))) 
		       (mean S samples)))

(equal? (f64vector-map (lambda (x) (truncate (* 10 x))) (upper S S sigma))
	(f64vector-map (lambda (x) (truncate (* 10 x)))
		       (covariance S samples (mean S samples))))

EOF
))

     (license
      "Copyright 2007-2009 Ivan Raikov and the Okinawa Institute of Science and Technology.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

A full copy of the GPL license can be found at
<http://www.gnu.org/licenses/>."))))

(if (eggdoc->html doc) (void))