diff --git a/README b/README
index 5260202d1c1b2277a1e1a016e545a5d6b63b00f9..b318a82eb91f70e74606a6c2e13f2b78d1b16bd1 100644
--- a/README
+++ b/README
@@ -13,6 +13,18 @@ Features:
 - This library is pure common lisp.  No C is used; the only requirement
 is that the GSL library must be loaded and accessible.
 
+---- policies:
+
+- All arrays and vectors that need to be passed to and from are
+  converted to or from GSL vectors or matrices; for GSL functions that
+  take/return raw C arrays (or pointers), a layer is wrapped
+  around them that do the conversion.
+
+-  Each defun-gsl corresponds to one C function.  Where multiple
+  such functions may be combined into a single regular or generic
+  function. 
+
+  
 ======================
 Needs work:
 See Thu Apr 20 2006 
@@ -45,3 +57,7 @@ This means that the following functions have been put aside:
 - #'polynomial-solve (requires vector_complex)
 - The index functions in sorting that are not returning permutations,
 because they need a vector_int.
+
+
+It would be nice to use e.g. vector-sap for direct vectors in
+e.g. SBCL.
diff --git a/data/data.lisp b/data/data.lisp
index e2be5a5ce395eb9483b7bb0353b97e3a387a129d..6a4bbb499e44af4aa37f58514136f9afec0150f3 100644
--- a/data/data.lisp
+++ b/data/data.lisp
@@ -3,7 +3,7 @@
 ; description: Using GSL storage.                        
 ; date:        Sun Mar 26 2006 - 16:32                   
 ; author:      Liam M. Healy                             
-; modified:    Sun Dec  3 2006 - 10:24
+; modified:    Sun Dec 31 2006 - 17:34
 ;********************************************************
 ;;; $Id: $
 
@@ -104,6 +104,7 @@
   '((FIXNUM . "_int")
     (SINGLE . "_float")
     (DOUBLE . "")
+    ;;(LONG-DOUBLE . "_long_double")
     (COMPLEX . "_complex")))
 
 (defmacro defdata
diff --git a/gsll.asd b/gsll.asd
index 10d3818e52ad7c1181145a90f3215ec2e380de0d..39a27adadc80f14f86bc9ccb3fcb6318e14d07ed 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Sun Dec  3 2006 - 21:58
+; modified:    Sun Dec 31 2006 - 22:07
 ;********************************************************
 ;;; $Id: $
 
@@ -133,4 +133,10 @@
    (:module statistics
 	    :depends-on (init data)
 	    :components
-	    ((:file "mean-variance")))))
+	    ((:file "mean-variance")
+	     (:file "absolute-deviation")
+	     (:file "higher-moments")
+	     (:file "autocorrelation")
+	     (:file "covariance")
+	     ;; minimum and maximum values provided in vector.lisp
+	     (:file "median-percentile")))))
diff --git a/init/interface.lisp b/init/interface.lisp
index 2b49944da6ecfcb1b5aacdffcd70d5ee073f25be..7ffa4facdce53c6e9032dd15fc05e88285635f00 100644
--- a/init/interface.lisp
+++ b/init/interface.lisp
@@ -3,7 +3,7 @@
 ; description: Macros to interface GSL functions.
 ; date:        Mon Mar  6 2006 - 22:35                   
 ; author:      Liam M. Healy
-; modified:    Fri Jul  7 2006 - 23:24
+; modified:    Sun Dec 31 2006 - 22:49
 ;********************************************************
 
 (in-package :gsl)
@@ -46,6 +46,9 @@
  ;;  (#:RET3504 (:DOUBLE (1+ KMAX))) (#:RET3505 (:DOUBLE (1+ KMAX)))
  ;;  (#:RET3506 :DOUBLE) (#:RET3507 :DOUBLE))
 
+(defun make-st (symbol type)
+  (list symbol type))
+
 (defun st-symbol (decl)
   (first decl))
 
@@ -169,7 +172,8 @@
 	  ,@(when invalidate `((cl-invalidate ,@invalidate)))
 	  ,@(when (or null-pointer-info (eq c-return :pointer))
 		  `((check-null-pointer ,cret-name
-		     ,@(or null-pointer-info '(:ENOMEM "No memory allocated")))))
+		     ,@(or null-pointer-info
+			   '(:ENOMEM "No memory allocated")))))
 	  ,@after
 	  (values
 	   ,@(case c-return
@@ -189,4 +193,22 @@
       (map-name ',(or index name) ,gsl-name)
       ,@(when export `((export ',name))))))
 
+;;;;****************************************************************************
+;;;; Macro defun-gsl 
+;;;;****************************************************************************
 
+(defmacro defun-optionals
+    (name arglist no-optional optionals &optional documentation)
+  "Define a function with and without optional arguments."
+  (let* ((optpos (position '&optional arglist))
+	 (mandatory-arglist (subseq arglist 0 optpos))
+	 (optional-arglist (subseq arglist (1+ optpos))))
+    `(defun ,name ,arglist
+       ,documentation
+       (if ,(first optional-arglist)
+	   (,(intern
+	      (concatenate 'string (string name) (string optionals)))
+	     ,@mandatory-arglist ,@optional-arglist)
+	   (,(intern
+	      (concatenate 'string (string name) (string no-optional)))
+	     ,@mandatory-arglist)))))
diff --git a/random/quasi.lisp b/random/quasi.lisp
index ddcdd76f049b4b75de9d37be83a1f52af1e88385..ed825e7030fad5dc999da60bd26cc7d4119ae45f 100644
--- a/random/quasi.lisp
+++ b/random/quasi.lisp
@@ -3,7 +3,7 @@
 ; description: Quasi-Random Number Generators            
 ; date:        Sun Jul 16 2006 - 15:54                   
 ; author:      Liam M. Healy                             
-; modified:    Mon Jul 17 2006 - 10:10
+; modified:    Sun Dec 31 2006 - 22:59
 ;********************************************************
 ;;; $Id: $
 
@@ -13,8 +13,12 @@
   ((dimension :initarg :dimension :accessor qr-dimension))
   (:documentation "A generator of quasi-random numbers."))
 
+(defparameter *default-quasi-random-number-generator* nil)
+
 (defun make-quasi-random-number-generator 
-    (dimension &optional (type *niederreiter2*) (generator t))
+    (dimension
+     &optional (type *default-quasi-random-number-generator*)
+     (generator t))
   "Make a random number generator; by default it is allocated on creation."
   (let ((instance
 	 (make-instance
@@ -110,6 +114,9 @@
      valid up to 12 dimensions."
   "gsl_qrng_niederreiter_2")
 
+(eval-when (:load-toplevel :execute)
+ (setf *default-quasi-random-number-generator* *niederreiter2*))
+
 (def-rng-type *sobol*
     "This generator uses the Sobol sequence described in Antonov, Saleev,
     @cite{USSR Comput. Maths. Math. Phys.} 19, 252 (1980). It is valid up to
diff --git a/statistics/absolute-deviation.lisp b/statistics/absolute-deviation.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..a3241ec99d1ae055cf7128dee5c36cee086833de
--- /dev/null
+++ b/statistics/absolute-deviation.lisp
@@ -0,0 +1,93 @@
+;********************************************************
+; file:        absolute-deviation.lisp                   
+; description: Absolute deviation                        
+; date:        Sun Dec 31 2006 - 13:19                   
+; author:      Liam M. Healy                             
+; modified:    Sun Dec 31 2006 - 21:47
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; To do: stride other than 1 when that information is availble from
+;;; the vector.
+
+(defun-gsl absolute-deviation-nom (data)
+  "gsl_stats_absdev"
+  (((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index absolute-deviation
+  :export nil
+  :documentation
+  "The absolute deviation from the mean of @var{data}.  The
+  absolute deviation from the mean is defined as
+  absdev  = (1/N) \sum |x_i - \Hat\mu|
+  where @math{x_i} are the elements of the dataset @var{data}.  The
+  absolute deviation from the mean provides a more robust measure of the
+  width of a distribution than the variance.  This function computes the
+  mean of @var{data} via a call to #'mean.")
+
+(defun-gsl absolute-deviation-m (data mean)
+  "gsl_stats_absdev_m"
+  (((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double))
+  :c-return :double
+  :index absolute-deviation
+  :export nil
+  :documentation
+  "The absolute deviation of the dataset @var{data}
+   relative to the given value of @var{mean},
+   absdev  = (1/N) \sum |x_i - mean|.
+   This function is useful if you have already computed the mean of
+   @var{data} (and want to avoid recomputing it), or wish to calculate the
+   absolute deviation relative to another value (such as zero, or the
+   median).")
+
+(export 'absolute-deviation)
+(defun-optionals absolute-deviation (data &optional mean)
+  -nom -m
+  "The absolute deviation from the mean of @var{data}.  The
+  absolute deviation from the mean is defined as
+  absdev  = (1/N) \sum |x_i - \Hat\mu|
+  where @math{x_i} are the elements of the dataset @var{data}.  The
+  absolute deviation from the mean provides a more robust measure of the
+  width of a distribution than the variance.  This function computes the
+  mean of @var{data} via a call to #'mean.")
+
+(defun-gsl weighted-absolute-deviation-nom (data weights)
+  "gsl_stats_wabsdev"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index weighted-absolute-deviation
+  :export nil)
+
+(defun-gsl weighted-absolute-deviation-m (data weights mean)
+  "gsl_stats_wabsdev_m"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double))
+  :c-return :double
+  :index weighted-absolute-deviation
+  :export nil)
+
+(export 'weighted-absolute-deviation)
+(defun-optionals weighted-absolute-deviation (data weights &optional mean)
+  -nom -m
+  "The weighted absolute deviation from the weighted
+   mean, defined as
+   absdev = (\sum w_i |x_i - \Hat\mu|) / (\sum w_i).")
+
+(lisp-unit:define-test absolute-deviation
+  (lisp-unit:assert-equal
+   '("0.618000000000d+01" "0.664777777778d+01" "0.618000000000d+01")
+   (with-data (vec vector-double 3)
+     (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
+     (with-data (weights vector-double 3)
+       (setf (data weights) #(3.0d0 1.0d0 2.0d0))
+       (let ((mean (mean vec)))
+	 (lisp-unit:fp-sequence
+	  (list
+	   (absolute-deviation vec)
+	   (weighted-absolute-deviation vec weights)
+	   (absolute-deviation vec mean))))))))
diff --git a/statistics/autocorrelation.lisp b/statistics/autocorrelation.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..da97ffb47e2b15348ce1d0cf17c6e99f94e73d35
--- /dev/null
+++ b/statistics/autocorrelation.lisp
@@ -0,0 +1,47 @@
+;********************************************************
+; file:        autocorrelation.lisp                   
+; description: Autocorrelation
+; date:        Sun Dec 31 2006 - 13:19                   
+; author:      Liam M. Healy                             
+; modified:    Sun Dec 31 2006 - 16:51
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; To do: stride other than 1 when that information is availble from
+;;; the vector.
+
+(defun-gsl autocorrelation-nom (data)
+  "gsl_stats_lag1_autocorrelation"
+  (((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index autocorrelation
+  :export nil)
+
+(defun-gsl autocorrelation-m (data mean)
+  "gsl_stats_lag1_autocorrelation_m"
+  (((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double))
+  :c-return :double
+  :index autocorrelation
+  :export nil)
+
+(export 'autocorrelation)
+(defun-optionals autocorrelation (data &optional mean)
+  -nom -m
+  "The lag-1 autocorrelation of the dataset @var{data}.
+  a_1 = {\sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i-1} - \Hat\mu)
+  \over
+  \sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i} - \Hat\mu)}.")
+
+(lisp-unit:define-test autocorrelation
+  (lisp-unit:assert-equal
+   '("-0.464636683425d-01" "-0.464636683425d-01")
+   (with-data (vec vector-double 3)
+     (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
+     (let ((mean (mean vec)))
+       (lisp-unit:fp-sequence
+	(list
+	 (autocorrelation vec)
+	 (autocorrelation vec mean)))))))
diff --git a/statistics/covariance.lisp b/statistics/covariance.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..2419427d3f540c1222356697d6c10d18cf590ed8
--- /dev/null
+++ b/statistics/covariance.lisp
@@ -0,0 +1,52 @@
+;********************************************************
+; file:        covariance.lisp                   
+; description: Covariance
+; date:        Sun Dec 31 2006 - 13:19                   
+; author:      Liam M. Healy                             
+; modified:    Sun Dec 31 2006 - 16:59
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; To do: stride other than 1 when that information is availble from
+;;; the vector.
+
+(defun-gsl covariance-nom (data1 data2)
+  "gsl_stats_covariance"
+  (((gsl-array data1) :pointer) (1 :int)
+   ((gsl-array data2) :pointer) (1 :int) ((dim0 data2) :size))
+  :c-return :double
+  :index covariance
+  :export nil)
+
+(defun-gsl covariance-m (data1 data2 mean1 mean2)
+  "gsl_stats_covariance_m"
+  (((gsl-array data1) :pointer) (1 :int)
+   ((gsl-array data2) :pointer) (1 :int) ((dim0 data2) :size)
+   (mean1 :double) (mean2 :double))
+  :c-return :double
+  :index covariance
+  :export nil)
+
+(export 'covariance)
+(defun-optionals covariance (data1 data2 &optional mean1 mean2)
+  -nom -m
+  "The covariance of the datasets @var{data1} and
+   @var{data2} which must both be of the same length @var{n}.
+   covar = {1 \over (n - 1)} \sum_{i = 1}^{n}
+      (x_{i} - \Hat x) (y_{i} - \Hat y).")
+
+(lisp-unit:define-test covariance
+  (lisp-unit:assert-equal
+   '("-0.293000000000d+00" "-0.293000000000d+00")
+   (with-data (vec1 vector-double 3)
+     (setf (data vec1) #(-3.21d0 1.0d0 12.8d0))
+     (with-data (vec2 vector-double 3)
+       (setf (data vec2) #(1.15d0 -1.0d0 0.5d0))
+       (let ((mean1 (mean vec1))
+	     (mean2 (mean vec2)))
+	 (lisp-unit:fp-sequence
+	  (list
+	   (covariance vec1 vec2)
+	   (covariance vec1 vec2 mean1 mean2))))))))
diff --git a/statistics/higher-moments.lisp b/statistics/higher-moments.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..cb94121ce39602fb24388f3e1721e8d4fdd9ba6e
--- /dev/null
+++ b/statistics/higher-moments.lisp
@@ -0,0 +1,132 @@
+;********************************************************
+; file:        higher-moments.lisp                       
+; description: Skewness and kurtosis.                    
+; date:        Sun Dec 31 2006 - 14:20                   
+; author:      Liam M. Healy                             
+; modified:    Sun Dec 31 2006 - 21:32
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+;;; To do: stride other than 1 when that information is availble from
+;;; the vector.
+
+(defun-gsl skewness-nomsd (data)
+  "gsl_stats_skew"
+  (((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index skewness
+  :export nil
+  :documentation
+  "The skewness of @var{data}, defined as
+  skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3
+  where @math{x_i} are the elements of the dataset @var{data}.
+  The skewness measures the asymmetry of the tails of a distribution.")
+
+(defun-gsl skewness-msd (data mean standard-deviation)
+  "gsl_stats_skew_m_sd"
+  (((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double) (standard-deviation :double))
+  :c-return :double
+  :index skewness
+  :export nil
+  :documentation
+  "The skewness of the dataset @var{data} using the
+   given values of the mean @var{mean} and standard deviation @var{sd},
+   skew = (1/N) \sum ((x_i - mean)/sd)^3
+   These functions are useful if you have already computed the mean and
+   standard deviation of @var{data} and want to avoid recomputing them.")
+
+(export 'skewness)
+(defun-optionals skewness
+    (data &optional mean standard-deviation)
+  -nomsd -msd
+  "The skewness of @var{data}, defined as
+  skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3
+  where @math{x_i} are the elements of the dataset @var{data}.
+  The skewness measures the asymmetry of the tails of a distribution.")
+
+(defun-gsl kurtosis-nomsd (data)
+  "gsl_stats_kurtosis"
+  (((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index kurtosis
+  :export nil)
+
+(defun-gsl kurtosis-msd (data mean standard-deviation)
+  "gsl_stats_kurtosis_m_sd"
+  (((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double) (standard-deviation :double))
+  :c-return :double
+  :index kurtosis
+  :export nil)
+
+(export 'kurtosis)
+(defun-optionals kurtosis
+    (data &optional mean standard-deviation)
+  -nomsd -msd
+  "The kurtosis of @var{data}, defined as
+   kurtosis = ((1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^4)  - 3
+   The kurtosis measures how sharply peaked a distribution is,
+   relative to its width.  The kurtosis is normalized to zero
+   for a gaussian distribution.")
+
+(defun-gsl weighted-skewness-nomsd (data weights)
+  "gsl_stats_wskew"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index weighted-skewness
+  :export nil)
+
+(defun-gsl weighted-skewness-msd (data weights mean standard-deviation)
+  "gsl_stats_wskew_m_sd"
+  (((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double) (standard-deviation :double))
+  :c-return :double
+  :index weighted-skewness
+  :export nil)
+
+(export 'weighted-skewness)
+(defun-optionals weighted-skewness
+    (data weights &optional mean standard-deviation)
+  -nomsd -msd
+  "The weighted skewness of the dataset.
+   skew = (\sum w_i ((x_i - xbar)/\sigma)^3) / (\sum w_i).")
+
+(defun-gsl weighted-kurtosis-nomsd (data weights)
+  "gsl_stats_wkurtosis"
+  (((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
+  :c-return :double
+  :index weighted-kurtosis
+  :export nil)
+
+(defun-gsl weighted-kurtosis-msd (data weights mean standard-deviation)
+  "gsl_stats_wkurtosis_m_sd"
+  (((gsl-array data) :pointer) (1 :int)
+   ((dim0 data) :size) (mean :double) (standard-deviation :double))
+  :c-return :double
+  :index weighted-kurtosis
+  :export nil)
+
+(export 'weighted-kurtosis)
+(defun-optionals weighted-kurtosis
+    (data weights &optional mean standard-deviation)
+  -nomsd -msd
+  "The weighted kurtosis of the dataset.
+   kurtosis = ((\sum w_i ((x_i - xbar)/sigma)^4) / (\sum w_i)) - 3.")
+
+(lisp-unit:define-test higher-moments
+  (lisp-unit:assert-equal
+   '("0.276511898399d+00" "0.276511898399d+00"
+     "-0.233333333333d+01" "-0.233333333333d+01")
+   (with-data (vec vector-double 3)
+     (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
+     (let* ((mean (mean vec))
+	    (sd (standard-deviation vec mean)))
+       (lisp-unit:fp-sequence
+	(list
+	 (skewness vec)
+	 (skewness vec mean sd)
+	 (kurtosis vec)
+	 (kurtosis vec mean sd)))))))
diff --git a/statistics/mean-variance.lisp b/statistics/mean-variance.lisp
index 608de09c0f32b77be9d08864b276913a63a08377..15f9add3fb33a7fa2be94b05095bca5a89ca93b6 100644
--- a/statistics/mean-variance.lisp
+++ b/statistics/mean-variance.lisp
@@ -3,7 +3,7 @@
 ; description: Mean, standard deviation, and variance    
 ; date:        Sat Dec  2 2006 - 22:15                   
 ; author:      Liam M. Healy                             
-; modified:    Sun Dec  3 2006 - 18:18
+; modified:    Sun Dec 31 2006 - 21:45
 ;********************************************************
 ;;; $Id: $
 
@@ -15,12 +15,18 @@
 (defmacro defun-gsl-stats (&rest args)
   "A defun-gsl for stats of double, single, and fixnum."
   (defun-gsl-all
+      ;;'(double single fixnum long)
       '(double single fixnum)
+      ;;'(:double :float :int :long-double)
       '(:double :float :int)
     "stats"
     'gsl-vector
     args))
 
+;;;;****************************************************************************
+;;;; Mean and weighted mean
+;;;;****************************************************************************
+
 (defgeneric mean (gsl-vector)
   (:documentation "The arithmetic mean of the vector.
    The arithmetic mean, or sample mean, is denoted by
@@ -32,7 +38,23 @@
   (((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
   :c-return :double)
 
-(defgeneric variance (gsl-vector)
+(defgeneric weighted-mean (gsl-vector weights)
+  (:documentation
+   "The weighted mean of the dataset, using the set of weights
+    The weighted mean is defined as
+    \Hat\mu = (\sum w_i x_i) / (\sum w_i)."))
+
+(defun-gsl-stats weighted-mean ((vector gsl-vector) weights)
+  "gsl_stats_wmean"
+  (((gsl-array weights) :pointer) (1 :int)
+    ((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
+  :c-return :double)
+
+;;;;****************************************************************************
+;;;; Variance
+;;;;****************************************************************************
+
+(defgeneric variance-nom (gsl-vector)
   (:documentation "The estimated, or @dfn{sample}, variance of
    @var{data}, a dataset of length @var{n} with stride @var{stride}.  The
    estimated variance is denoted by @math{\Hat\sigma^2} and is defined by
@@ -45,11 +67,13 @@
 
    This function computes the mean via a call to #'mean.  If
    you have already computed the mean then you can pass it directly to
-   #'variance-m"))
+   #'variance-m."))
 
-(defun-gsl-stats variance ((vector gsl-vector))
+(defun-gsl-stats variance-nom ((vector gsl-vector))
   "gsl_stats_variance"
   (((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
+  :index variance
+  :export nil
   :c-return :double)
 
 (defgeneric variance-m (gsl-vector mean)
@@ -60,14 +84,74 @@
   "gsl_stats_variance_m"
   (((gsl-array vector) :pointer) (1 :int)
    ((dim0 vector) :size) (mean :double))
+  :index variance
+  :export nil
   :c-return :double)
 
-(defgeneric standard-deviation (gsl-vector)
+(export 'variance)
+(defun-optionals variance (data &optional mean)
+  -nom -m
+  "The estimated, or @dfn{sample}, variance of
+   @var{data}, a dataset of length @var{n} with stride @var{stride}.  The
+   estimated variance is denoted by @math{\Hat\sigma^2} and is defined by
+   \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2
+   where @math{x_i} are the elements of the dataset @var{data}.  Note that
+   the normalization factor of @math{1/(N-1)} results from the derivation
+   of @math{\Hat\sigma^2} as an unbiased estimator of the population
+   variance @math{\sigma^2}.  For samples drawn from a gaussian distribution
+   the variance of @math{\Hat\sigma^2} itself is @math{2 \sigma^4 / N}.")
+
+;;;;****************************************************************************
+;;;; Weighted variance
+;;;;****************************************************************************
+
+(defgeneric weighted-variance-nom (gsl-vector weights)
+   (:documentation nil))
+
+(defun-gsl-stats weighted-variance-nom ((vector gsl-vector) weights)
+  "gsl_stats_wvariance"
+  (((gsl-array weights) :pointer) (1 :int)
+    ((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
+  :index weighted-variance
+  :export nil
+  :c-return :double)
+
+(defgeneric weighted-variance-m (gsl-vector weights mean)
+  (:documentation "Compute the weighted variance with the mean known to
+   avoid its recomputation."))
+
+(defun-gsl-stats weighted-variance-m ((vector gsl-vector) weights mean)
+  "gsl_stats_wvariance_m"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array vector) :pointer) (1 :int)
+   ((dim0 vector) :size) (mean :double))
+  :index variance
+  :export nil
+  :c-return :double)
+
+(export 'weighted-variance)
+(defun-optionals weighted-variance (data weights &optional mean)
+  -nom -m
+   "The estimated variance using weights
+   The estimated variance of a weighted dataset is defined as
+   \Hat\sigma^2 = ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) 
+                \sum w_i (x_i - \Hat\mu)^2
+   Note that this expression reduces to an unweighted variance with the
+   familiar @math{1/(N-1)} factor when there are @math{N} equal non-zero
+   weights.")
+
+;;;;****************************************************************************
+;;;; Standard deviation
+;;;;****************************************************************************
+
+(defgeneric standard-deviation-nom (gsl-vector)
   (:documentation "The standard deviation, square root of the variance."))
 
-(defun-gsl-stats standard-deviation ((vector gsl-vector))
+(defun-gsl-stats standard-deviation-nom ((vector gsl-vector))
   "gsl_stats_sd"
   (((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
+  :index standard-deviation
+  :export nil
   :c-return :double)
 
 (defgeneric standard-deviation-m (gsl-vector mean)
@@ -80,6 +164,45 @@
    ((dim0 vector) :size) (mean :double))
   :c-return :double)
 
+(export 'standard-deviation)
+(defun-optionals standard-deviation (data &optional mean)
+  -nom -m
+  "The standard deviation, square root of the variance.")
+
+;;;;****************************************************************************
+;;;; Weighted standard deviation
+;;;;****************************************************************************
+
+(defgeneric weighted-standard-deviation-nom (gsl-vector weights))
+
+(defun-gsl-stats weighted-standard-deviation-nom
+    ((vector gsl-vector) weights)
+  "gsl_stats_wsd"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
+  :index weighted-standard-deviation
+  :export nil
+  :c-return :double)
+
+(defgeneric weighted-standard-deviation-m (gsl-vector weights mean))
+
+(defun-gsl-stats weighted-standard-deviation-m
+    ((vector gsl-vector) weights mean)
+  "gsl_stats_wsd_m"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array vector) :pointer) (1 :int)
+   ((dim0 vector) :size) (mean :double))
+  :c-return :double)
+
+(export 'weighted-standard-deviation)
+(defun-optionals weighted-standard-deviation (data weights &optional mean)
+  -nom -m
+  "The standard deviation, square root of the variance.")
+
+;;;;****************************************************************************
+;;;; Variance with fixed mean
+;;;;****************************************************************************
+
 (defgeneric variance-with-fixed-mean (gsl-vector mean)
   (:documentation
    "An unbiased estimate of the variance of
@@ -108,37 +231,81 @@
    (mean :double))
   :c-return :double)
 
-;;; Examples and unit test
+;;;;****************************************************************************
+;;;; Weighted variance with fixed mean
+;;;;****************************************************************************
+
+(defgeneric weighted-variance-with-fixed-mean (gsl-vector weights mean)
+  (:documentation
+   "An unbiased estimate of the variance of weighted
+    dataset when the population mean @var{mean} of the underlying
+    distribution is known @emph{a priori}.  In this case the estimator for
+    the variance replaces the sample mean @math{\Hat\mu} by the known
+    population mean @math{\mu},
+    \Hat\sigma^2 = (\sum w_i (x_i - \mu)^2) / (\sum w_i)."))
+
+(defun-gsl-stats weighted-variance-with-fixed-mean
+    ((vector gsl-vector) weights mean)
+  "gsl_stats_wvariance_with_fixed_mean"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size)
+   (mean :double))
+  :c-return :double)
+
+(defgeneric weighted-standard-deviation-with-fixed-mean
+    (gsl-vector weights mean)
+  (:documentation
+   "The square root of the corresponding variance
+   function #'weighted-variance-with-fixed-mean."))
+
+(defun-gsl-stats weighted-standard-deviation-with-fixed-mean
+    ((vector gsl-vector) weights mean)
+  "gsl_stats_wsd_with_fixed_mean"
+  (((gsl-array weights) :pointer) (1 :int)
+   ((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size)
+   (mean :double))
+  :c-return :double)
+
+;;;;****************************************************************************
+;;;; Examples and unit test
+;;;;****************************************************************************
+
 (lisp-unit:define-test mean-variance
   (lisp-unit:assert-equal
-   '("0.353000000000d+01" "0.688807000000d+02" "0.688807000000d+02"
+   '("0.353000000000d+01" "0.282833333333d+01" "0.688807000000d+02"
+     "0.688807000000d+02" "0.849805863636d+02" "0.849805863636d+02"
      "0.829943974013d+01" "0.829943974013d+01" "0.461413666667d+02"
      "0.679274367739d+01")
    (with-data (vec vector-double 3)
      (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
+     (with-data (weights vector-double 3)
+       (setf (data weights) #(3.0d0 1.0d0 2.0d0))
+       (let ((mean (mean vec))
+	     (wmean (weighted-mean vec weights)))
+	 (lisp-unit:fp-sequence
+	  (list
+	   mean wmean
+	   (variance vec)
+	   (variance vec mean)
+	   (weighted-variance vec weights)
+	   (weighted-variance vec weights wmean)
+	   (standard-deviation vec)
+	   (standard-deviation vec mean)
+	   (variance-with-fixed-mean vec 4.0d0)
+	   (standard-deviation-with-fixed-mean vec 4.0d0)))))))
+  (lisp-unit:assert-equal
+   '("0.333333333333d+01" "0.253333333333d+02" "0.253333333333d+02"
+     "0.503322295685d+01" "0.503322295685d+01" "0.173333333333d+02"
+     "0.416333199893d+01")
+   (with-data (vec vector-fixnum 3)
+     (setf (data vec) #(8 4 -2))
      (let ((mean (mean vec)))
        (lisp-unit:fp-sequence
 	(list
 	 mean
 	 (variance vec)
-	 (variance-m vec mean)
+	 (variance vec mean)
 	 (standard-deviation vec)
-	 (standard-deviation-m vec mean)
+	 (standard-deviation vec mean)
 	 (variance-with-fixed-mean vec 4.0d0)
-	 (standard-deviation-with-fixed-mean vec 4.0d0)))))
-   (lisp-unit:assert-equal
-    '("0.333333333333d+01" "0.253333333333d+02" "0.253333333333d+02"
-      "0.503322295685d+01" "0.503322295685d+01" "0.173333333333d+02"
-      "0.416333199893d+01")
-    (with-data (vec vector-fixnum 3)
-      (setf (data vec) #(8 4 -2))
-      (let ((mean (mean vec)))
-	(lisp-unit:fp-sequence
-	 (list
-	  mean
-	  (variance vec)
-	  (variance-m vec mean)
-	  (standard-deviation vec)
-	  (standard-deviation-m vec mean)
-	  (variance-with-fixed-mean vec 4.0d0)
-	  (standard-deviation-with-fixed-mean vec 4.0d0))))))))
+	 (standard-deviation-with-fixed-mean vec 4.0d0)))))))
diff --git a/statistics/median-percentile.lisp b/statistics/median-percentile.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..86ff5dc6d2fcb08f033ef65fdb03868084f63415
--- /dev/null
+++ b/statistics/median-percentile.lisp
@@ -0,0 +1,67 @@
+;********************************************************
+; file:        median-percentile.lisp                   
+; description: Median and percentile
+; date:        Sun Dec 31 2006 - 13:19                   
+; author:      Liam M. Healy                             
+; modified:    Sun Dec 31 2006 - 22:07
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+;;; To do: stride other than 1 when that information is availble from
+;;; the vector.
+
+(defun-gsl median (sorted-data)
+  "gsl_stats_median_from_sorted_data"
+  (((gsl-array sorted-data) :pointer) (1 :int) ((dim0 sorted-data) :size))
+  :c-return :double
+  :documentation
+  "The median value of @var{sorted_data}, a dataset
+   of length @var{n} with stride @var{stride}.  The elements of the array
+   must be in ascending numerical order.  There are no checks to see
+   whether the data are sorted, so the function #'sort should
+   always be used first.
+   When the dataset has an odd number of elements the median is the value
+   of element @math{(n-1)/2}.  When the dataset has an even number of
+   elements the median is the mean of the two nearest middle values,
+   elements @math{(n-1)/2} and @math{n/2}.  Since the algorithm for
+   computing the median involves interpolation this function always returns
+   a floating-point number, even for integer data types.")
+
+(defun-gsl quantile (sorted-data fraction)
+  "gsl_stats_quantile_from_sorted_data"
+  (((gsl-array sorted-data) :pointer) (1 :int) ((dim0 sorted-data) :size)
+   (fraction :double))
+  :c-return :double
+  :documentation
+  "A quantile value of sorted-data, gsl-vector-double.  The
+   elements of the array must be in ascending numerical order.  The
+   quantile is determined by a fraction between 0 and 1.  For
+   example, to compute the value of the 75th percentile
+   'fraction should have the value 0.75.
+   There are no checks to see whether the data are sorted, so the function
+   #'sort should always be used first.
+   \hbox{quantile} = (1 - \delta) x_i + \delta x_{i+1}
+   where @math{i} is @code{floor}(@math{(n - 1)f}) and @math{\delta} is
+   @math{(n-1)f - i}.
+   Thus the minimum value of the array (@code{data[0*stride]}) is given by
+   'fraction equal to zero, the maximum value (@code{data[(n-1)*stride]}) is
+   given by 'fraction equal to one and the median value is given by 'fraction
+   equal to 0.5.  Since the algorithm for computing quantiles involves
+   interpolation this function always returns a floating-point number, even
+   for integer data types.")
+
+(lisp-unit:define-test median-percentile
+  (lisp-unit:assert-first-fp-equal
+   "0.100000000000d+01"
+   (with-data (vec vector-double 3)
+     (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
+     (median vec)))
+  (lisp-unit:assert-first-fp-equal
+   "0.185000000000d+01"
+   (with-data (vec vector-double 7)
+     (setf (data vec)
+	   #(-18.0d0 -12.0d0 -3.21d0 0.5d0 1.0d0 2.7d0 12.8d0))
+     (quantile vec 0.75d0))))
+