Commit e803849a authored by liam's avatar liam

Add linear least squares fitting. Define generic functions so that

warnings are suppressed.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3274 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 0a78898d
This diff is collapsed.
;********************************************************
; file: permutation.lisp
; description: Permutations
; date: Sun Mar 26 2006 - 11:51
; author: Liam M. Healy
; modified: Sun Jun 4 2006 - 21:02
;********************************************************
;;; $Id: $
;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-01-21 11:12:01EST permutation.lisp>
;; $Id: $
(in-package :gsl)
......@@ -33,7 +29,7 @@
(((pointer permutation) :pointer) ((first indices) :size))
:type :method
:c-return :size
:documentation "The ith element of the permutation.")
:documentation "The ith element of the permutation.") ; FDL
;;;;****************************************************************************
;;;; Setting values
......@@ -44,9 +40,9 @@
(((pointer permutation) :pointer))
:type :method
:c-return :void
:documentation
"Initialize the permutation @var{p} to the identity, i.e.
@math{(0,1,2,@dots{},n-1)}.")
:documentation ; FDL
"Initialize the permutation p to the identity, i.e.
(0,1,2,@dots{},n-1).")
(defun-gsl copy ((destination gsl-permutation) (source gsl-permutation))
"gsl_permutation_memcpy"
......@@ -54,18 +50,17 @@
((pointer source) gsl-permutation-c))
:type :method
:invalidate (destination)
:documentation
"Copy the elements of the permutation @var{src} into the
permutation @var{dest}. The two permutations must have the same size.")
:documentation ; FDL
"Copy the elements of the permutation source into the
permutation destination. The two permutations must have the same size.")
(defun-gsl swap-elements ((p gsl-permutation) i j)
"gsl_permutation_swap"
(((pointer p) gsl-permutation-c) (i :size) (j :size))
:type :method
:invalidate (p)
:documentation
"Exchanges the @var{i}-th and @var{j}-th elements of the
permutation @var{p}.")
:documentation ; FDL
"Exchanges the ith and jth elements of the permutation p.")
;;;;****************************************************************************
;;;; Permutation properties
......@@ -74,26 +69,26 @@
(defun-gsl permutation-size (p)
"gsl_permutation_size"
(((pointer p) gsl-permutation-c))
:c-return :size
:c-return :size ; FDL
:documentation
"The size of the permutation @var{p}.")
"The size of the permutation p.")
(defun-gsl permutation-data (p)
"gsl_permutation_data"
(((pointer p) gsl-permutation-c))
:c-return :pointer
:documentation
:documentation ; FDL
"A pointer to the array of elements in the
permutation @var{p}.")
permutation p.")
(defun-gsl data-valid ((permutation gsl-permutation))
"gsl_permutation_valid"
(((pointer permutation) :pointer))
:type :method
:c-return :boolean
:documentation
"Check that the permutation @var{p} is valid. The @var{n}
elements should contain each of the numbers 0 to @math{@var{n}-1} once and only
:documentation ; FDL
"Check that the permutation p is valid. The n
elements should contain each of the numbers 0 to n-1 once and only
once.")
;;;;****************************************************************************
......@@ -105,26 +100,26 @@
(((pointer p) gsl-permutation-c))
:invalidate (p)
:c-return :void
:documentation
"Reverse the order of the elements of the permutation @var{p}.")
:documentation ; FDL
"Reverse the order of the elements of the permutation p.")
(defun-gsl permutation-inverse (inv p)
"gsl_permutation_inverse"
(((pointer inv) gsl-permutation-c) ((pointer p) gsl-permutation-c))
:invalidate (inv)
:documentation
"Find the inverse of the permutation @var{p}.")
:documentation ; FDL
"Find the inverse of the permutation p.")
(defun-gsl permutation-next (p)
"gsl_permutation_next"
(((pointer p) gsl-permutation-c))
:c-return :success-failure
:invalidate (p)
:documentation
"Advance the permutation @var{p} to the next permutation
:documentation ; FDL
"Advance the permutation p to the next permutation
in lexicographic order and return p and T. If no further
permutations are available, return p and NIL with
@var{p} unmodified. Starting with the identity permutation and
p unmodified. Starting with the identity permutation and
repeatedly applying this function will iterate through all possible
permutations of a given order.")
......@@ -133,11 +128,11 @@
(((pointer p) gsl-permutation-c))
:c-return :success-failure
:invalidate (p)
:documentation
"Step backwards from the permutation @var{p} to the
:documentation ; FDL
"Step backwards from the permutation p to the
previous permutation in lexicographic order, returning p and T.
If no previous permutation is available, return
p and NIL with @var{p} unmodified.")
p and NIL with p unmodified.")
;;;;****************************************************************************
;;;; Applying Permutations
......@@ -146,41 +141,47 @@
(defun-gsl permute (p data stride n)
"gsl_permute"
(((pointer p) gsl-permutation-c) (data :pointer) (stride :size) (n :size))
:documentation
"Apply the permutation @var{p} to the array @var{data} of
size @var{n} with stride @var{stride}.")
:documentation ; FDL
"Apply the permutation p to the array data of
size n with stride stride.")
(defun-gsl permute-inverse (p data stride n)
"gsl_permute_inverse"
(((pointer p) gsl-permutation-c) (data :pointer) (stride :size) (n :size))
:documentation
"Apply the inverse of the permutation @var{p} to the array @var{data} of
size @var{n} with stride @var{stride}.")
:documentation ; FDL
"Apply the inverse of the permutation p to the array data of
size n with stride.")
(export 'permute-vector)
(defgeneric permute-vector (permutation vector)
(:documentation ; FDL
"Apply the permutation p to the elements of the
vector v considered as a row-vector acted on by a permutation
matrix from the right, v' = v P. The jth column of the
permutation matrix P is given by the p_j-th column of the
identity matrix. The permutation p and the vector v must
have the same length."))
(defun-gsl-vdsfc permute-vector (p (v gsl-vector))
"gsl_permute_vector"
(((pointer p) gsl-permutation-c) ((pointer v) gsl-vector-c))
:invalidate (v)
:documentation
"Apply the permutation @var{p} to the elements of the
vector @var{v}, considered as a row-vector acted on by a permutation
matrix from the right, @math{v' = v P}. The @math{j}-th column of the
permutation matrix @math{P} is given by the @math{p_j}-th column of the
identity matrix. The permutation @var{p} and the vector @var{v} must
have the same length.")
:invalidate (v))
(export 'permute-vector-inverse)
(defgeneric permute-vector-inverse (permutation vector)
(:documentation ; FDL
"Apply the inverse of the permutation p to the
elements of the vector v considered as a row-vector acted on by
an inverse permutation matrix from the right, v' = v P^T. Note
that for permutation matrices the inverse is the same as the transpose.
The jth column of the permutation matrix P is given by
the p_j-th column of the identity matrix. The permutation p
and the vector v must have the same length."))
(defun-gsl-vdsfc permute-vector-inverse (p (v gsl-vector))
"gsl_permute_vector_inverse"
(((pointer p) gsl-permutation-c) ((pointer v) gsl-vector-c))
:invalidate (v)
:documentation
"Apply the inverse of the permutation @var{p} to the
elements of the vector @var{v}, considered as a row-vector acted on by
an inverse permutation matrix from the right, @math{v' = v P^T}. Note
that for permutation matrices the inverse is the same as the transpose.
The @math{j}-th column of the permutation matrix @math{P} is given by
the @math{p_j}-th column of the identity matrix. The permutation @var{p}
and the vector @var{v} must have the same length.")
:invalidate (v))
(defun-gsl permutation* (p pa pb)
"gsl_permutation_mul"
......@@ -188,10 +189,10 @@
((pointer pa) gsl-permutation-c)
((pointer pb) gsl-permutation-c))
:invalidate (p)
:documentation
"Combine the two permutations @var{pa} and @var{pb} into a
single permutation @var{p}, where @math{p = pa . pb}. The permutation
@var{p} is equivalent to applying @math{pb} first and then @var{pa}.")
:documentation ; FDL
"Combine the two permutations pa and pb into a
single permutation p where p = pa . pb. The permutation
p is equivalent to applying pb first and then pa.")
;;;;****************************************************************************
;;;; Permutations in cyclic form
......@@ -201,24 +202,24 @@
"gsl_permutation_linear_to_canonical"
(((pointer q) gsl-permutation-c) ((pointer p) gsl-permutation-c))
:invalidate (q)
:documentation
"Compute the canonical form of the permutation @var{p} and
stores it in the output argument @var{q}.")
:documentation ; FDL
"Compute the canonical form of the permutation p and
stores it in the output argument q.")
(defun-gsl canonical-to-linear (p q)
"gsl_permutation_canonical_to_linear"
(((pointer p) gsl-permutation-c) ((pointer q) gsl-permutation-c))
:invalidate (p)
:documentation
"Convert a permutation @var{q} in canonical form back into
linear form storing it in the output argument @var{p}.")
:documentation ; FDL
"Convert a permutation q in canonical form back into
linear form storing it in the output argument p.")
(defun-gsl inversions (p)
"gsl_permutation_inversions" (((pointer p) gsl-permutation-c))
:c-return :size
:documentation
:documentation ; FDL
"Count the number of inversions in the permutation
@var{p}. An inversion is any pair of elements that are not in order.
p. An inversion is any pair of elements that are not in order.
For example, the permutation 2031 has three inversions, corresponding to
the pairs (2,0) (2,1) and (3,1). The identity permutation has no
inversions.")
......@@ -226,16 +227,15 @@
(defun-gsl linear-cycles (p)
"gsl_permutation_linear_cycles" (((pointer p) gsl-permutation-c))
:c-return :size
:documentation
"Count the number of cycles in the permutation @var{p}, given in linear form.")
:documentation ; FDL
"Count the number of cycles in the permutation p, given in linear form.")
(defun-gsl canonical-cycles (p)
"gsl_permutation_canonical_cycles"
(((pointer p) gsl-permutation-c))
:c-return :size
:documentation
"Count the number of cycles in the permutation @var{q},
given in canonical form.")
:documentation ; FDL
"Count the number of cycles in the permutation q, given in canonical form.")
;;;;****************************************************************************
;;;; Examples and unit test
......
;********************************************************
; file: eigensystems.lisp
; description: Eigenvectors and eigenvalues
; date: Sun May 21 2006 - 19:52
; author: Liam M. Healy
; modified: Tue Jul 4 2006 - 23:09
;********************************************************
;;; $Id: $
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2008-01-21 11:30:57EST eigensystems.lisp>
;; $Id: $
(in-package :gsl)
......@@ -45,33 +41,38 @@
(defun-gsl eigen-symm-alloc (n)
"gsl_eigen_symm_alloc" ((n :size))
:documentation "Allocate a workspace for computing eigenvalues of
@var{n}-by-@var{n} real symmetric matrices. The size of the workspace
is @math{O(2n)}."
:documentation ; FDL
"Allocate a workspace for computing eigenvalues of
n-by-n real symmetric matrices. The size of the workspace
is O(2n)."
:c-return :pointer)
(defun-gsl eigen-symm-free (w)
"gsl_eigen_symm_free" ((w :pointer))
:c-return :void
:documentation "Free the memory associated with the workspace @var{w}.")
:documentation ; FDL
"Free the memory associated with the workspace w.")
(defun-gsl eigen-symmv-alloc (n)
"gsl_eigen_symmv_alloc" ((n :size))
:documentation "Allocate a workspace for computing eigenvalues and
eigenvectors of @var{n}-by-@var{n} real symmetric matrices. The size of
the workspace is @math{O(4n)}."
:documentation ; FDL
"Allocate a workspace for computing eigenvalues and
eigenvectors of n-by-n real symmetric matrices. The size of
the workspace is O(4n)."
:c-return :pointer)
(defun-gsl eigen-symmv-free (w)
"gsl_eigen_symmv_free" ((w :pointer))
:c-return :void
:documentation "Free the memory associated with the workspace @var{w}.")
:documentation ; FDL
"Free the memory associated with the workspace w.")
(defun-gsl eigenvalues-symmetric (A eval ws)
"gsl_eigen_symm" ((A gsl-matrix-c) (eval gsl-vector-c) (ws :pointer))
:documentation "Eigenvalues of the real symmetric matrix
@var{A}. Additional workspace of the appropriate size must be provided
in @var{w}. The diagonal and lower triangular part of @var{A} are
:documentation ; FDL
"Eigenvalues of the real symmetric matrix
A. Additional workspace of the appropriate size must be provided
in w. The diagonal and lower triangular part of A are
destroyed during the computation, but the strict upper triangular part
is not referenced. The eigenvalues are stored in the vector @var{eval}
and are unordered."
......@@ -81,13 +82,14 @@
(defun-gsl eigenvalues-eigenvectors-symmetric (A eval evec ws)
"gsl_eigen_symmv"
((A gsl-matrix-c) (eval gsl-vector-c) (evec gsl-matrix-c) (ws :pointer))
:documentation "The eigenvalues and eigenvectors of the real
symmetric matrix @var{A}. Additional workspace of the appropriate size
must be provided in @var{w}. The diagonal and lower triangular part of
@var{A} are destroyed during the computation, but the strict upper
:documentation ; FDL
"The eigenvalues and eigenvectors of the real
symmetric matrix A. Additional workspace of the appropriate size
must be provided in w. The diagonal and lower triangular part of
A are destroyed during the computation, but the strict upper
triangular part is not referenced. The eigenvalues are stored in the
vector @var{eval} and are unordered. The corresponding eigenvectors are
stored in the columns of the matrix @var{evec}. For example, the
vector eval and are unordered. The corresponding eigenvectors are
stored in the columns of the matrix evec. For example, the
eigenvector in the first column corresponds to the first eigenvalue.
The eigenvectors are guaranteed to be mutually orthogonal and normalised
to unit magnitude."
......@@ -100,23 +102,25 @@
(defun-gsl eigen-herm-alloc (n)
"gsl_eigen_herm_alloc" ((n :size))
:documentation "Allocate a workspace for computing eigenvalues of
@var{n}-by-@var{n} complex hermitian matrices. The size of the workspace
is @math{O(3n)}."
:documentation ; FDL
"Allocate a workspace for computing eigenvalues of
n-by-n complex hermitian matrices. The size of the workspace
is O(3n)."
:c-return :pointer)
(defun-gsl eigen-herm-free (w)
"gsl_eigen_herm_free" ((w :pointer))
:c-return :void
:documentation "Free the memory associated with the workspace @var{w}.")
:documentation ; FDL
"Free the memory associated with the workspace w.")
(defun-gsl eigenvalues-hermitian (A eval w)
"gsl_eigen_herm"
(((pointer A) gsl-matrix-c) ((pointer eval) gsl-vector-c) (w :pointer))
:documentation
:documentation ; FDL
"Compute the eigenvalues of the complex hermitian matrix
@var{A}. Additional workspace of the appropriate size must be provided
in @var{w}. The diagonal and lower triangular part of @var{A} are
A. Additional workspace of the appropriate size must be provided
in w. The diagonal and lower triangular part of A are
destroyed during the computation, but the strict upper triangular part
is not referenced. The imaginary parts of the diagonal are assumed to be
zero and are not referenced. The eigenvalues are stored in the vector
......@@ -126,29 +130,31 @@
(defun-gsl eigen-hermv-alloc (n)
"gsl_eigen_hermv_alloc" ((n :size))
:documentation "Allocate a workspace for computing eigenvalues and
eigenvectors of @var{n}-by-@var{n} complex hermitian matrices. The size of
the workspace is @math{O(5n)}."
:documentation ; FDL
"Allocate a workspace for computing eigenvalues and
eigenvectors of n-by-n complex hermitian matrices. The size of
the workspace is O(5n)."
:c-return :pointer)
(defun-gsl eigen-hermv-free (w)
"gsl_eigen_hermv_free" ((w :pointer))
:c-return :void
:documentation "Free the memory associated with the workspace @var{w}.")
:documentation ; FDL
"Free the memory associated with the workspace w.")
(defun-gsl eigenvalues-eigenvectors-hermitian (A eval evec w)
"gsl_eigen_hermv"
(((pointer A) gsl-matrix-c) ((pointer eval) gsl-vector-c)
((pointer evec) gsl-matrix-c) (w :pointer))
:documentation
:documentation ; FDL
"Compute the eigenvalues and eigenvectors of the complex
hermitian matrix @var{A}. Additional workspace of the appropriate size
must be provided in @var{w}. The diagonal and lower triangular part of
@var{A} are destroyed during the computation, but the strict upper
hermitian matrix A. Additional workspace of the appropriate size
must be provided in w. The diagonal and lower triangular part of
A are destroyed during the computation, but the strict upper
triangular part is not referenced. The imaginary parts of the diagonal
are assumed to be zero and are not referenced. The eigenvalues are
stored in the vector @var{eval} and are unordered. The corresponding
complex eigenvectors are stored in the columns of the matrix @var{evec}.
stored in the vector eval and are unordered. The corresponding
complex eigenvectors are stored in the columns of the matrix evec.
For example, the eigenvector in the first column corresponds to the
first eigenvalue. The eigenvectors are guaranteed to be mutually
orthogonal and normalised to unit magnitude."
......@@ -163,15 +169,19 @@
"gsl_eigen_sort_t from /usr/include/gsl/gsl_eigen.h."
:value-ascending :value-descending :absolute-ascending :absolute-descending)
(export 'sort-eigenvalues-eigenvectors)
(defgeneric sort-eigenvalues-eigenvectors (eval evec sort-type)
(:documentation ; FDL
"Simultaneously sort the eigenvalues stored in the vector
eval and the corresponding real eigenvectors stored in the columns
of the matrix evec into ascending or descending order according to
the value of the parameter sort-type: :value-ascending,
:value-descending, :absolute-ascending, :absolute-descending."))
(defun-gsl sort-eigenvalues-eigenvectors
((eval gsl-vector-double) evec sort-type)
"gsl_eigen_symmv_sort"
((eval gsl-vector-c) (evec gsl-matrix-c) (sort-type eigen-sort-type))
:documentation "Simultaneously sort the eigenvalues stored in the vector
@var{eval} and the corresponding real eigenvectors stored in the columns
of the matrix @var{evec} into ascending or descending order according to
the value of the parameter @var{sort_type}, :value-ascending,
:value-descending, :absolute-ascending, :absolute-descending."
:type :method
:invalidate (eval evec))
......@@ -179,10 +189,6 @@
((eval gsl-vector-complex) evec sort-type)
"gsl_eigen_hermv_sort"
((eval gsl-vector-c) (evec gsl-matrix-c) (sort-type eigen-sort-type))
:documentation "Simultaneously sort the eigenvalues stored in the vector
@var{eval} and the corresponding complex eigenvectors stored in the
columns of the matrix @var{evec} into ascending or descending order
according to the value of the parameter @var{sort_type}."
:type :method
:invalidate (eval evec))
......
;********************************************************
; file: mathematical.lisp
; description: Mathematical functions
; date: Wed Mar 8 2006 - 22:09
; author: Liam M. Healy
; modified: Mon May 29 2006 - 23:06
;********************************************************
;; Mathematical functions
;; Liam Healy, Wed Mar 8 2006 - 22:09
;; Time-stamp: <2008-01-21 18:01:09EST mathematical.lisp>
;; $Id: $
(in-package :gsl)
......@@ -83,11 +80,13 @@
:documentation
"exp(x)-1, computed in a way that is accurate for small x.")
(defun-gsl hypotenuse (x y)
(defun-gsl hypotenuse* (x y)
;; This is redundant; there is "gsl_sf_hypot_e" defined as
;; #'hyptoenuse.
"gsl_hypot" ((x :double) (y :double))
:c-return :double
:documentation
"sqrt{x^2 + y^2} computed in a way that avoids overflow.")
"The hypotenuse sqrt{x^2 + y^2} computed in a way that avoids overflow.")
;; Not clear why this function exists
(defun-gsl gsl-asinh (x)
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-01-20 23:18:16EST gsll.asd>
;; Time-stamp: <2008-01-21 17:58:44EST gsll.asd>
;; $Id: $
(asdf:defsystem "gsll"
......@@ -92,7 +92,7 @@
:components
((:file "rng-types")
(:file "generators" :depends-on (rng-types))
(:file "quasi" :depends-on (rng-types))
(:file "quasi" :depends-on (rng-types generators))
(:file "gaussian" :depends-on (rng-types))
(:file "gaussian-tail" :depends-on (rng-types))
(:file "gaussian-bivariate" :depends-on (rng-types))
......@@ -172,4 +172,5 @@
(:file "roots-one" :depends-on (init general))
(:file "minimization-one" :depends-on (init general))
(:file "roots-multi" :depends-on (init general data))
(:file "minimization-multi" :depends-on (init general data))))
(:file "minimization-multi" :depends-on (init general data))
(:file "linear-least-squares" :depends-on (init general data))))
;********************************************************
; file: probability-distribution.lisp
; description: Histogram probability distribution.
; date: Mon Jan 1 2007 - 17:51
; author: Liam M. Healy
; modified: Mon Jan 1 2007 - 22:22
;********************************************************
;;; $Id: $
;; Histogram probability distribution.
;; Liam Healy, Mon Jan 1 2007 - 17:51
;; Time-stamp: <2008-01-21 18:02:00EST probability-distribution.lisp>
;; $Id: $
(in-package :gsl)
......@@ -116,6 +112,3 @@
@math{sum[i] <= r < sum[i+1]} and @math{delta} is
@math{(r - sum[i])/(sum[i+1] - sum[i])}.")
(export 'sample)
(defun sample (pdf value)
(histo-1d2d pdf sample (value)))
......@@ -324,7 +324,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p>
<tr>
<td><a href="Least_002dSquares-Fitting.html">
Least Squares Fitting</a></td>
<td>Not started</td>
<td>Done</td>
</tr>
<tr>
<td><a href="Nonlinear-Least_002dSquares-Fitting.html">
......@@ -374,7 +374,7 @@ and expect a few failures on amd64 (four), none? on i386. Often, the
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
<small>
Time-stamp: <2008-01-20 21:29:20EST index.html>
Time-stamp: <2008-01-21 17:58:20EST index.html>
</small>
<!-- hhmts end -->
</div>
......
;; Load GSL
;; Liam Healy Sat Mar 4 2006 - 18:53
;; Time-stamp: <2008-01-16 18:50:15 liam init.lisp>
;; Time-stamp: <2008-01-21 12:33:11EST init.lisp>
;; $Id: $
(defpackage gsll
......@@ -21,7 +21,9 @@
;;; The following define :size, from cffi-unix which became cffi-net,
;;; which is apparently turning into something even bigger and more
;;; irrelevant.
;;; irrelevant. All we need is the definition of :size. This should
;;; be changed to something outside the keyword package to suppress
;;; warnings from CFFI.
(cffi:defctype :uint32 :unsigned-int)
(cffi:defctype :uint64 :unsigned-long)
......
;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
;; Time-stamp: <2008-01-21 17:55:19EST linear-least-squares.lisp>
;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Linear regression
;;;;****************************************************************************
;;; Error in GSL documentation or in code: vectors must be equal
;;; lengths but strides re not? Shouldn't (floor length stride) be n?
;;; Error in GSL documentation for gsl_fit_linear_est, "c00" instead
;;; of "cov00" etc. Last arg to gsl_fit_wmul is labelled sumsq but
;;; referred to as chisq.
(defun-gsl linear-fit (x y &optional (x-stride 1) (y-stride 1))
"gsl_fit_linear"
(((gsl-array x) :pointer) (x-stride :size)
((gsl-array y) :pointer) (y-stride :size)
((dim0 x) :size) (c0 :double) (c1 :double)
(cov00 :double) (cov01 :double) (cov11 :double)
(sumsq :double))
:documentation ; FDL
"Determine the best-fit linear regression coefficients
and returns as the first two values c0,c1
of the model Y = c_0 + c_1 X for the dataset
(x,y), two vectors of equal length with strides
x-stride and y-stride. The errors on y are assumed unknown so
the variance-covariance matrix for the
parameters (c0, c1) is estimated from the scatter of the
points around the best-fit line and returned as the third, fourth
and fifth values cov00, cov01, cov11.
The sum of squares of the residuals from the best-fit line is returned
as the final value.")
(defun-gsl weighted-linear-fit
(x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
"gsl_fit_wlinear"
(((gsl-array x) :pointer) (x-stride :size)
((gsl-array weight) :pointer) (weight-stride :size)
((gsl-array y) :pointer) (y-stride :size)
((dim0 x) :size) (c0 :double) (c1 :double)
(cov00 :double) (cov01 :double) (cov11 :double)
(chisq :double))
:documentation ; FDL
"Compute the best-fit linear regression coefficients
c0, c1 of the model Y = c_0 + c_1 X for the weighted
dataset (x, y), two vectors of equal length with strides
x-stride and y-stride, and return as the first two values.
The vector weight, of the same length
and stride w-stride, specifies the weight of each datapoint. The
weight is the reciprocal of the variance for each datapoint in y.
The covariance matrix for the parameters (c0, c1) is
computed using the weights and returned via the parameters
(cov00, cov01, c0v01) as the next three values. The weighted sum
of squares of the residuals from the best-fit line, \chi^2, is
returned in as the last value.")
(defun-gsl linear-estimate (x c0 c1 cov00 cov01 cov11)
"gsl_fit_linear_est"
((x :double) (c0 :double) (c1 :double)
(cov00 :double) (cov01 :double) (cov11 :double)
(y :double) (y-error :double))
:documentation ; FDL
"Use the best-fit linear regression coefficients
c0, c1 and their covariance