;; @module stat.lsp
;; @description Basic statistics and plotting library
;; @version 3.0 - Eliminated plot functions and f-prob which now is built-in as prob-f
;; @version 3.1 - Documentation changes
;; @version 3.2 - Documentation, mention built-in 'stats' function since 10.4.2
;;
;; @author Lutz Mueller, 2001-2014
;; <h2>Functions for statistics</h2>
;; To use this module it has to be loaded at the beginning of the
;; program file:
;; <pre>
;; (load (append (env "NEWLISPDIR") "/modules/stat.lsp"))
;; ; or shorter
;; (module "stat.lsp")
;; </pre>
;; All functions work on integers and floats or a mix of both. <lists> are normal
;; LISP lists. <matrices> are lists of lists, one list for each row in the
;; two dimensional data matrix. See the function 'stat:matrix' on how to make matrices
;; from lists.
;;
;; In version 3.0 of 'stat.lsp' the usage for Gnuplot the stat:plot and stat:plotXY 
;; functions has been eliminated. Instead use the module 'plot.lsp' shipped with all binary and
;; source distributions. The F-distribution function 'f-prob' has also been eliminated
;; instead use one of new the built-in 'prob-f' or 'crit-f' functions.
;;
;; The documention contains only the call patterns. See the source for more
;; documentation.
;;
;; <h2>Summary of functions</h2><br>
;; <h3>General uni- and bi- variate statistics</h3>
;; <pre>
;; stat:sum      - sum of a vector of numbers (see also built-in stats since 10.4.2)
;; stat:mean     - arithmetik mean of a vector of numbers (see also built-in stats since 10.4.2)
;; stat:var      - estimated variance of numbers in a vector sample
;; stat:sdev     - estimated standard deviation of numbers in a vector (see also built-in stats since 10.4.2)
;; stat:sum-sq   - sum of squares of a data vector
;; stat:sum-xy   - sum of products of a two data vectors
;; stat:corr     - correlation coefficient between two vectors (built-in since 10.4.2)
;; stat:cov      - covariance of two number vectors
;; stat:sum-d2   - sum of squared differences of a vector from its mean
;; stat:sum-d2xy - sum of squared diffferences of two vectors
;; stat:regression - calculates the intecept and slope of a regression estimate
;; stat:fit      - return the fitted line using regression coefficients
;; stat:moments  - calulates 1st to 3rd moments from a vector of numbers
;; </pre>
;; <h3>Multi variate statistics</h3>
;; <pre>
;; stat:multiple-reg  - calculates a multiple regression
;; stat:cov-matrix    - calculates a covariance matrix
;; stat:corr-matrix   - calculates a correlation matrix
;; </pre>
;; <h3>Time series</h3>
;; <pre>
;; stat:smooth   - smoothes a vector of numbers
;; stat:lag      - calcultes a difference list with specified lag
;; stat:cumulate - cumulate a data vector
;; stat:power    - calculate the power spectrum of a time series
;; </pre>
;; <h3>Matrix and list utilities</h3>
;; <pre>
;; stat:matrix       - make a matrix from column vectors
;; stat:diagonal     - make a diagonal matrix
;; stat:get-diagonal - return the diagonal of a matrix in a vector
;; stat:mat-map      - map a binary function on to matrices
;; </pre>


(if (< (sys-info -2) 10111)
	(constant (global 'term) name))

;; @syntax (stat:corr <X> <Y>)
;; @param <X> A list of numbers.
;; @param <Y> A list of numbers.
;; @return Correlation coefficient of lists <X> and <Y>.
;; A 'corr' native function is built into newLISP since 10.4.2.

;; @syntax (stat:cov <X> <Y>)
;; @param <X> A list of numbers.
;; @param <Y> A list of numbers.
;; @return Covariance of data in lists <X> and <Y>

;; @syntax (stat:cov-matrix <X>)
;; @param <X> A matrix of numbers.
;; @return Covariance matrix of <X> with <N> rows and <k> columns.

;; @syntax (stat:corr-matrix <X>)
;; @param <X> A matrix of numbers.
;; @return Correlation matrix of <X> with <N> rows and <k> columns.

;; @syntax (stat:cumulate <X>)
;; @param <X> A list of numbers.
;; @return The cumulated list of <X>.

;; @syntax (stat:diagonal <item> <N>)
;; @param <item> The diagonal element.
;; @return A diagonal matrix of length <N> with <item> in the diagonal.

;; @syntax (stat:fit <X> <Y>)
;; @param <X> A list of numbers.
;; @param <Y> A list of numbers.
;; @return fitted line based on '(stat:regression X Y)'.

;; @syntax (stat:f-prob <F> <df1> <df2>)
;; @param <F> The variance ratio.
;; @param <df1> Degrees of freedom.
;; @param <df2> Degrees of freedom.
;; @return Probablity of F variance ratio for <df1>, <df2> degress of freedom.

;; @syntax (stat:get-diagonal <X>)
;; @param <X> An matrix filled with numbers.
;; @return A list from the diagonal elements of <X>.

;; @syntax (stat:lag <X> <n>)
;; @param <X> A list of numbers.
;; @param <n> Lag n.
;; @return A differenced list of <X> with a lag of <n>.
;; If the length of list <X> is <l> then the length of the resulting
;; differenced list is <l - n>.

;; @syntax (stat:mat-map <op> <A> <B>)
;; @return Matrix map, e.g. '(stat:mat-map + A B)'.
;; Used for adding and subtracting matrices.

;; @syntax (stat:matrix <C1> .... <CN>)
;; @param <C1> The first column list of values.
;; @param <CN> The Nth column list of values.
;; @return A matrix off  <1> to <N> columns <C>.

;; @syntax (stat:mean <X>) 
;; @param <X> A list of numbers.
;; @return The mean of data in list <X>.

;; @syntax (stat:moments <X>)
;; @param <X> A list of numbers.
;; @return Calculates all moments of list <X>.

;; @syntax (stat:multiple-reg <X> <offY>)
;; @param <X> A matrix of numbers.
;; @param <offY> Zero based offset into <Y>.
;; @return Multiple regression of vars in <X> onto <Y> at <offsetY>.

;; @syntax (stat:power <TS>)
;; @param <TS> A time series of numbers.
;; @return The power spectrum of a time series

;; @syntax (stat:regression <X> <Y>)
;; @param <X> A list of numbers.
;; @param <Y> A list of numbers.
;; returns <(b0 b1)> coefficients of regression <Y = b0 + b1*X>.

;; @syntax (stat:sdev <X>)
;; @param <X> A list of numbers.
;; @return Standard deviation of data in list <X>.

;; @syntax (stat:smooth <X> <alpha>)
;; @param <X> A list of numbers.
;; @param <alpha> Smoothing coefficient <0 &lt; alpha &lt; 1>.
;; @return Exponentially smoothed sequence in <X>.

;; @syntax (stat:sum <X>)
;; @param <X> A list of numbers,
;; @return Sum of data in list <X>.

;; @syntax (stat:sum-d2 <X>)
;; @param <X> A list of numbers.
;; @return Sum of squared diffs <(x - mean(X))^2> in list <X>.

;; @syntax (stat:sum-d2xy <X> <Y>)
;; @return Sum of squared differences <(x - y)^2> of elements in lists <X> and <Y>.

;; @syntax (stat:sum-sq <X>)              
;; @param <X> A list of numbers.
;; @return Sum of <x*x> data elements in list <X>.

;; @syntax (stat:sum-xy <X> <Y>)
;; @param <X> A list of numbers.
;; @param <Y> A list of numbers.
;; @return Sum of products <x*y> data elements in lists <X> and <Y>.

;; @syntax (stat:var <X>)
;; @param <X> A list of numbers.
;; @return The variance of the data in list <X>.



(context 'stat)

;-------------------  General uni and bi-variate statistics --------------------

; sum of a data vector X
(define (sum X)
	(apply add X))

; mean of a data vector X
(define (mean X)
	(div (sum X) (length X)))

; variance of a data vector X
(define (var X)
	(div (sum-d2 X) (sub (length X) 1)))

; standard deviation of a data vector X
(define (sdev X)
	(sqrt (var X)))

; sum of squares of a data vector X
(define (sum-sq X)
	(apply add (map mul X X)))

; sum of the product of data vectors X*Y
(define (sum-xy X Y)
	(apply add (map mul X Y)))

; covariance of data vectors X Y
(define (cov X Y)
	(sub (sum-xy X Y) (div (mul (sum X) (sum Y)) (length X))))

; sum of sqared differenses of X to mean of X
(define (sum-d2 X)
	(sub (sum-sq X) (div (mul (sum X) (sum X)) (length X))))

; Pearson r, product moment correlation of data vectors X and Y
(define (stat:corr X Y)
	(div (cov X Y) (sqrt (mul (sum-d2 X) (sum-d2 Y)))))

; regression Yest = b0 + b1*X calculates intercept b0 and slope b1
(define (regression X Y)
	(set 'b1 (div (cov X Y) (sum-d2 X)))
	(set 'b0 (sub (mean Y) (mul b1 (mean X))))
	(list b0 b1))

; fitted line using regression Y on X
(define (fit X Y, coeffs b0 b1)
	(set 'coeffs (regression X Y))
	(set 'b0 (first coeffs))
	(set 'b1 (last coeffs))
	(map (lambda (x) (add b0 (mul x b1))) X))

; sum of squared differences of X and Y
(define (sum-d2xy X Y)
  (apply add (map (lambda (x y) (mul (sub x y) (sub x y))) X Y)))


; moments of a vector of numbers
;
(define (moments vector, n median mean avg-dev std-dev var skew kurtosis dev sum)
  (set 'n (length vector))

  (set 'sum (apply add vector))
  (set 'mean (div sum n))

  (set 'avg-dev 0 'std-dev 0 'var 0 'skew 0 'kurtosis 0)

  (set 'dev (map sub vector (dup mean n)))
  (set 'avg-dev (div (apply add (map abs dev)) n))
  (set 'var (div (apply add (map mul dev dev)) (- n 1)))
  (set 'skew (apply add (map mul dev dev dev)))
  (set 'kurtosis (apply add (map  mul dev dev dev dev)))

  (set 'std-dev (sqrt var))

  (if (> var 0.0)
    (begin
	(set 'skew (div skew (mul n var std-dev)))
	(set 'kurtosis (sub (div kurtosis (mul n var var)) 3.0))))

  (sort vector)
  (set 'mid (/ n 2))

  (if (= 0 (% n 2))
	(set 'median (div (add (nth mid vector) (nth (- mid 1) vector)) 2))
	(set 'median (nth mid vector)))

  (list n median mean avg-dev std-dev var skew kurtosis)

;  (println (format "n:                  %d" n))
;  (println (format "median:             %f" median))
;  (println (format "mean:               %f" mean))
;  (println (format "average_deviation:  %f" avg-dev))
;  (println (format "standard_deviation: %f" std-dev))
;  (println (format "variance:           %f" var))
;  (println (format "skew:               %f" skew))
;  (println (format "kurtosis:           %f" kurtosis))
)

;-------------------------------- Time Series ----------------------------------

; expontial smoothing with 0 < alpha <= 1
(define (smooth lst alpha , previous slist)
  (set 'previous (first lst))
  (set 'slist '())
  (dolist (elmnt lst) 
   (set 'previous (add (mul alpha elmnt) (mul (sub 1 alpha) previous))) 
   (push previous slist))
  (reverse slist)) ; could be written shorter starting v.9.9.5
                   ; because push returns the modified list

;
; seasonal difference list with variable lag
; the resulting list is lag shorterm than the original
;
(define (lag lst n , sLst)
	(set 'sLst lst)
	(dotimes (i n) (pop lst))
	(set 'sLst (slice sLst 0 (length lst)))
	(map sub lst sLst))


;
; cumulate of a list
;
(define (cumulate lst , sc cum)
	(set 'sc 0.0)
	(set 'cum '())
	(dolist (x lst)
		(push (inc sc x) cum))
	(reverse cum)) ; could be written shorter after 9.9.5
	               ; because push returns the list

;
; power spectrum
;
; takes a rows by 2 columns matrix with real part in the first and 
; imagenary part the in the second column. If all numbers are real
; then the second column is just 0's.
;
; returns a matrix with two rows. First row contains power numbers
; and second row contains the respective frequencies
;

(define (power ts , lenOrg fts n n2 ps mid frqs)
	; remember original length
	(set 'lenOrg (length ts))
	; do discrete fourier transform
	(set 'fts (transpose (fft ts)))
	; calc power spectrum 
	(set 'n (length (transpose fts)))
	(set 'n2 (mul n n))
	(set 'ps (map (lambda (x y) (add (mul x x) (mul y y))) (nth 0 fts) (nth 1 fts)))
	(set 'ps (map (lambda (x) (mul (div x n2) 2)) ps))
	; the first and last are not multiplied by 2, divide them back
	; use deprecated nth-set in versions older than 9.9.02
	(setf (ps 0) (div (first ps) 2))
	(set 'mid (sub (div n 2) 1))
	(replace mid ps (div (nth mid ps) 2))
	; calc a vector with frequencies, adjusted for the new power-2 length
	; which came back from the FFT
	(set 'frqs (sequence 0 (- n 1)))
	(set 'frqs (map (lambda (x) (mul (div x n) lenOrg)) frqs))
	(transpose (matrix ps frqs)))



;------------------------- multivariate statistics -----------------------------

;
; multiple regression of variables in X onto one of variables in X, Y
; indicated by column offset offY
; 
; X is N rows by k columns, the column at offset offY is Y
;
; returns a matrix with two rows:
; first row is regression coefficients and multiple R: b0, b1, b2 ....., R
; second row is sum of squares: regression-SQ, error-SQ, total-SQ
; (the unused part of the second row is <tt>nil</tt> padded)
;
; the SQs can be used to calculate mean sqares for regression and error:
;
; regression-MSQ = regression-SQ / (k - 1)
; error-MSQ = error-sq / (n - k - 1)
;
; F-ratio = regression-MSQ / error-MSQ  
; with k and (n - k - 1) df degreees of freedom
;
;
;
;
(define (multiple-reg X offY ,  Y Ycoffs b b0 R2 Yest sqErr sqTotal sqReg sq d)
	(set 'covX (cov-matrix X))
	(set 'Y (extract-col X offY))
	; covX is the covariance matrix
	(pop covX offY)
	(set 'cvX (transpose covX))
	; the covariance matrix is reduced to cvX and the
	; extracted values put in Ycoffs
	(set 'Ycoffs (matrix (pop cvX offY)))
	; b contains the regression coefficients except for b0
	(set 'b (multiply (invert cvX) Ycoffs))
	; calculate multiple R2 as b'*b / sqTotal
	(set 'sqTotal (sum-d2 Y))
	(set 'R2 (div (first (first (multiply (transpose b) Ycoffs))) sqTotal))
	; estimate Y without b0
	(set 'Yest (multiply (reduce-col X offY) b))
	; calculate b0, d is the difference between Y and the Y estimate
	; b0 is the mean of differences between Y and Yest
	(set 'd (mat-map sub (matrix Y) Yest))
	(set 'b0 (mean (first (transpose d))))
	; estimate Y including b0
	(set 'Yest (mat-map add Yest (matrix (dup b0 (length Yest)))))
	; error sum of squares
	(set 'sqErr (sum-d2xy Y (first (transpose Yest))))
	; regression sum of squares
	(set 'sqReg (sub sqTotal sqErr))
	; make list b out of b0, b1, b2 ... sqrt(R2)
	(set 'b (append (list b0) (first (transpose b)) (list (sqrt R2))))
	; make list sq out of sqReg, sqErr and sqTotal
	(set 'sq (list sqReg sqErr sqTotal))
	; return matrix with two rows:
	(transpose (matrix b sq)))



;
; covariance matrix cov 
; 
; matrix x with N rows and k columns
; 
;
(define (cov-matrix X , XtX N I sumX sumX2)
	(set 'XtX (multiply (transpose X) X)) 
	(set 'N (length X))
	(set 'I (matrix (dup 1 N)))
	(set 'sumX (multiply (transpose X) I)) 
	(set 'sumX2 (multiply sumX (transpose sumX)))  
	(set 'sumX2 (multiply sumX2 (diagonal (div 1 N) (length sumX2))))
	(mat-map sub XtX sumX2))


;
; correlation matrix
;
; matrix X with N rows and k columns
;
;
(define (corr-matrix X , covX N d dd)
	(set 'covX (cov-matrix X))
	(set 'd (matrix (get-diagonal covX)))
	(set 'dd (multiply d (transpose d)))
	(set 'dd (map (lambda (z) (map sqrt z)) dd))
	(mat-map div covX dd))

	
;
; probablity of F variance ratio with degrees of freedom df1 df2
;
; 
(define (f-prob F df1 df2)
  (let (prob (mul 2 (betai (div df2 (add df2 (mul df1 F))) 
                           (mul 0.5 df2) 
                           (mul 0.5 df1)))) 
   (div (if (> prob 1) (sub 2 prob) prob) 2)))	

;----------------------------- utility functions -------------------------------

;
; make a matrix from 1 up to 16 lists
;
(define (matrix)
	(transpose (args)))


;
; make a diagonal matrix n by n and elmnt in the diagonal
; 
;
(define (diagonal elmnt n, m lst)
	(set 'm '())
	(dotimes (i n)
		(set 'lst (dup 0 n))
		; use deprecated nth-set in versions older than 9.9.02
		(if (< (sys-info -2) 9902)
			(nth-set (lst i) elmnt)
			(setf (lst i) elmnt))
		(push lst m))
	(reverse m)) ; make shorter starting v.9.9.5
		
;
; get the diagonal from a square matrix
;
(define (get-diagonal X , d x)
	(set 'd '())
	(dotimes (idx (length X))
		(push (nth idx (nth idx X)) d))
	(reverse d))

;
; matrix map
;
; e.g.: (mat-map sub A B) ;; for matrix subtraction
;
(define (mat-map op A B)
	(map (lambda (x y) (map op x y)) A B))

;
; reduce matrix by a column at offset
;
; returns the reduced matrix
;
(define (reduce-col matr off, X)
	(set 'matr (transpose matr))
	(pop matr off)
	(transpose matr))

;
; extract a column from a matrix
;
; returns the extracted column
;
(define (extract-col matr off, X)
	(pop (transpose matr) off))

;
; convert list to ascii lines terminated by CR-LF
; for storage in files usable by Gnuplot, R, Excel etc.
;
; example: 
;
; (write-file "MyData.txt" (list2ascii mydata-list))
;
(define (list2ascii lst)
	(append (join (map string lst) "\r\n") "\r\n"))

; eof