;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; Includes information derived from the GSL (Gnu Scientific Library)
;;;   header files which are copyrighted and distributed under the GPL
;;;   by their respective authors.
;;;
;;; 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 2 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.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(libload "gsl/gsl-config")
(libload "gsl/gsl-idx")
(libload "gsl/statistics")

;; dummy function that adds GSL C header files in C file
(eval `(de dummy_stats () ,gsl-cpheader-cmd ()))

#? *** Statistical Functions
;; This is a set of simple statistical functions.
;; Most of them are high-level calls to GSL functions.
;; They operate on vectors of double.

#? (stat-mean <v>)
;; return the mean of elements of vector <v>.
(de stat-mean (v)
  ((-idx1- (-double-)) v)
  (let ((pv (idx-ptr v))
	(n (idx-dim v 0))
	(s (idx-modulo v 0)))
    (to-double #{ gsl_stats_mean( $pv, $s, $n) #})))

#? (stat-variance <v>)
;; return the variance of values in vector <v>.
(de stat-variance (v)
  ((-idx1- (-double-)) v)
  (let ((pv (idx-ptr v))
	(n (idx-dim v 0))
	(s (idx-modulo v 0)))
    (to-double #{ gsl_stats_variance( $pv, $s, $n) #})))

#? (stat-sd <v>)
;; return the standard deviation of values in vector <v>.
(de stat-sd (v)
  ((-idx1- (-double-)) v)
  (let ((pv (idx-ptr v))
	(n (idx-dim v 0))
	(s (idx-modulo v 0)))
    (to-double #{ gsl_stats_sd( $pv, $s, $n) #})))

#? (stat-skew <v>)
;; return the skew (measure asymmetry) of the distributions
;; of values in <v>.
(de stat-skew (v)
  ((-idx1- (-double-)) v)
  (let ((pv (idx-ptr v))
	(n (idx-dim v 0))
	(s (idx-modulo v 0)))
    (to-double #{ gsl_stats_skew( $pv, $s, $n) #})))

#? (stat-kurtosis <v>)
;; return the kurtosis (measure non-gaussianness) of the distributions
;; of values in <v>.
(de stat-kurtosis (v)
  ((-idx1- (-double-)) v)
  (let ((pv (idx-ptr v))
	(n (idx-dim v 0))
	(s (idx-modulo v 0)))
    (to-double #{ gsl_stats_kurtosis( $pv, $s, $n) #})))

#? (stat-covariance <v1> <v2>)
;; compute convariance of v1 and v2.
;; covariance = 1/(n - 1) sum_i (v1_i - [v1]) (v2_i - [v2])
(de stat-covariance (v1 v2)
  ((-idx1- (-double-)) v1 v2)
  (let ((pv1 (idx-ptr v1))
	(pv2 (idx-ptr v2))
	(n (idx-dim v1 0))
	(s1 (idx-modulo v1 0))
	(s2 (idx-modulo v2 0)))
    (when (<> n (idx-dim v1 0)) (error "vectors have different sizes"))
    (to-double #{ gsl_stats_covariance( $pv1, $s1, $pv2, $s2, $n) #})))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(dhc-make
 ()
 dummy_stats
 stat-mean
 stat-variance
 stat-sd
 stat-skew
 stat-kurtosis
 stat-covariance
 )



