diff --git a/init/defmfun.lisp b/init/defmfun.lisp
index 942578f52065d003b497ac9820ad301733c5af04..d2ff4c308f076d29fef659abddd9b3556013cab4 100644
--- a/init/defmfun.lisp
+++ b/init/defmfun.lisp
@@ -1,6 +1,6 @@
 ;; Macro for defining GSL functions.
 ;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
-;; Time-stamp: <2008-08-17 16:50:24EDT defmfun.lisp>
+;; Time-stamp: <2008-08-17 18:14:19EDT defmfun.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -102,12 +102,14 @@
 	 (with-defmfun-key-args key-args
 	   (if (eq index t)
 	       (setf index name))
-	   (flet ((mapnfn (gsl-name) `(map-name ',index ,gsl-name)))
+	   (flet ((mapnfn (gslnm) `(map-name ',index ,gslnm)))
 	     (append
 	      (when index
 		(if indexed-functions
 		    (mapcar #'mapnfn indexed-functions)
-		    (list (mapnfn gsl-name))))
+		    (if (listp gsl-name)
+			(mapcar #'mapnfn gsl-name)
+			(list (mapnfn gsl-name)))))
 	      (when export `((export ',name))))))))
     (if index-export
 	(if (symbolp (first definition))
@@ -264,7 +266,9 @@
 		     (mapcar (lambda (args)
 			       (actual-element-c-type eltype args))
 			     c-arguments)
-		     key-args)
+		     key-args
+		     'body-optional-arg
+		     t)
 		    (complete-definition
 		     defn
 		     (if (eq defn :method) (list nil name) name)
@@ -418,30 +422,31 @@
     (name arglist gsl-name c-arguments key-args)
   "Expand defmfun where there is an optional argument
    present, giving the choice between two different GSL functions."
-  (with-defmfun-key-args key-args
-    (let ((optpos (position '&optional arglist)))
-      (if optpos
-	  (let ((mandatory-arglist (subseq arglist 0 optpos))
-		(optional-arglist (subseq arglist (1+ optpos))))
-	    (remf key-args :documentation)
-	    (setf indexed-functions gsl-name)
-	    (let ((body
-		   `(if ,(first optional-arglist)
-			,(expand-defmfun-body 
-			   nil
-			   (append mandatory-arglist optional-arglist)
-			   (second gsl-name)
-			   (second c-arguments)
-			   key-args)
-			,(expand-defmfun-body
-			   nil
-			   mandatory-arglist
-			   (first gsl-name)
-			   (first c-arguments)
-			   key-args))))
-	      (if (defgeneric-method-p name)
-		  body
-		  `((defun ,name ,arglist ,documentation ,@body)))))))))
+  (complete-definition 'cl:defun name arglist gsl-name c-arguments key-args
+		       'body-optional-arg))
+
+(defun body-optional-arg
+    (name arglist gsl-name c-arguments key-args)
+  "Create the body of a defmfun with &optional in its arglist,
+   where the presence or absence of the optional argument(s)
+   selects one of two GSL functions."
+  (let ((optpos (position '&optional arglist)))
+    (if optpos
+	(let ((mandatory-arglist (subseq arglist 0 optpos))
+	      (optional-arglist (subseq arglist (1+ optpos))))
+	  `(if ,(first optional-arglist)
+	       ,(body-no-optional-arg 
+		 name
+		 (append mandatory-arglist optional-arglist)
+		 (second gsl-name)
+		 (second c-arguments)
+		 key-args)
+	       ,(body-no-optional-arg
+		 name
+		 mandatory-arglist
+		 (first gsl-name)
+		 (first c-arguments)
+		 key-args))))))
 
 ;;;;****************************************************************************
 ;;;; A single function
@@ -478,7 +483,9 @@
 		    (list val)))))
 	   c-arguments)))
 
-(defun complete-definition (definition name arglist gsl-name c-arguments key-args)
+(defun complete-definition
+    (definition name arglist gsl-name c-arguments key-args
+     &optional (body-maker 'body-no-optional-arg) mapdown)
   "A complete definition form, starting with defun, :method, or defmethod."
   (destructuring-bind
 	(&key documentation &allow-other-keys) key-args
@@ -489,11 +496,14 @@
        ,(declaration-form
 	 (cl-argument-types arglist c-arguments)
 	 (set-difference (arglist-plain-and-categories arglist nil)
-			 (c-arguments c-arguments)))
+			 (if mapdown
+			     (apply 'union
+				    (mapcar 'c-arguments c-arguments))
+			     (c-arguments c-arguments))))
        ,@(when documentation (list documentation))
-       ,(expand-defmfun-body name arglist gsl-name c-arguments key-args))))
+       ,(funcall body-maker name arglist gsl-name c-arguments key-args))))
 
-(defun expand-defmfun-body (name arglist gsl-name c-arguments key-args)
+(defun body-no-optional-arg (name arglist gsl-name c-arguments key-args)
   "Expand the body (computational part) of the defmfun."
   (with-defmfun-key-args key-args
     (let* ((cret-type (if (member c-return *special-c-return*)
diff --git a/statistics/mean-variance.lisp b/statistics/mean-variance.lisp
index 0ce6e89706e81c658e9e258754976b11f9a992bb..fd030fcbcb8fd4655e1c20a8700ebfb726d9fe39 100644
--- a/statistics/mean-variance.lisp
+++ b/statistics/mean-variance.lisp
@@ -1,6 +1,6 @@
 ;; Mean, standard deviation, and variance    
 ;; Liam Healy, Sat Dec  2 2006 - 22:15
-;; Time-stamp: <2008-03-09 19:21:47EDT mean-variance.lisp>
+;; Time-stamp: <2008-08-17 12:06:30EDT mean-variance.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -20,33 +20,55 @@
 ;;;; Mean and weighted mean
 ;;;;****************************************************************************
 
-(defgeneric mean (vector)
-  (:documentation			; FDL
+(defmfun mean ((vector vector))
+  ("gsl_stats" :type "_mean")
+  (((c-pointer vector) :pointer) (1 :int) ((dim0 vector) sizet))
+  :definition :generic
+  :element-types :no-complex
+  :c-return :double
+  :documentation			; FDL
    "The arithmetic mean of the vector.
    The arithmetic mean, or sample mean, is denoted by
-   \Hat\mu and defined as \Hat\mu = (1/N) \sum x_i.  Returns a double-float."))
-
-(defmfun-stats mean ((vector vector))
-  "gsl_stats_mean"
-  (((gsl-array vector) :pointer) (1 :int) ((dim0 vector) size))
-  :c-return :double)
-
-(defgeneric weighted-mean (vector weights)
-  (:documentation			; FDL
-   "The weighted mean of the dataset, using the set of weights
+   \Hat\mu and defined as \Hat\mu = (1/N) \sum x_i.  Returns a double-float.")
+
+(defmfun weighted-mean ((vector vector) (weights vector))
+  ("gsl_stats" :type "_wmean")
+  (((c-pointer weights) :pointer) (1 :int)
+   ((c-pointer vector) :pointer) (1 :int)
+   ((dim0 vector) sizet))
+  :definition :generic
+  :element-types :float
+  :c-return :element-c-type
+  :documentation			; FDL
+  "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)."))
-
-(defmfun-stats-ds weighted-mean ((vector vector) weights)
-  "gsl_stats_wmean"
-  (((gsl-array weights) :pointer) (1 :int)
-    ((gsl-array vector) :pointer) (1 :int) ((dim0 vector) size))
-  :c-return :double)
+    \Hat\mu = (\sum w_i x_i) / (\sum w_i).")
 
 ;;;;****************************************************************************
 ;;;; Variance
 ;;;;****************************************************************************
 
+(defmfun variance ((vector vector) &optional mean)
+  (("gsl_stats" :type "_variance")
+   ("gsl_stats" :type "_variance_m"))
+  ((((c-pointer vector) :pointer) (1 :int) ((dim0 vector) sizet))
+   (((c-pointer vector) :pointer) (1 :int)
+   ((dim0 vector) sizet) (mean :double)))
+  :definition :generic
+  :element-types :no-complex
+  :c-return :double
+  :documentation			; FDL
+  "The estimated, or sample, variance of data.  The
+   estimated variance is denoted by \Hat\sigma^2 and is defined by
+   \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2
+   where x_i are the elements of the dataset data.  Note that
+   the normalization factor of 1/(N-1) results from the derivation
+   of \Hat\sigma^2 as an unbiased estimator of the population
+   variance \sigma^2.  For samples drawn from a gaussian distribution
+   the variance of \Hat\sigma^2 itself is 2 \sigma^4 / N.
+   If the mean value is known, it may be supplied which will use more
+   efficient routines to compute the variance.")
+
 (defgeneric variance-nom (vector)
   (:documentation			; FDL
    "The estimated, or sample, variance of data.  The