diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -5,3 +5,6 @@ fixed configure.hs to work for local install updated to reflect exportation of check from Data.Packed.Development in hmatrix-0.9.3.0 + +0.1.03: + added Numeric.GSL.Regression.Linear diff --git a/hstatistics.cabal b/hstatistics.cabal --- a/hstatistics.cabal +++ b/hstatistics.cabal @@ -1,5 +1,5 @@ Name: hstatistics -Version: 0.1.0.2 +Version: 0.1.0.3 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010 @@ -10,7 +10,7 @@ Synopsis: Statistics Description: Purely functional interface for statistics based on hmatrix and GSL Category: Math -tested-with: GHC ==6.10.4 +tested-with: GHC ==6.12.1 cabal-version: >=1.2 @@ -35,12 +35,14 @@ Numeric.GSL.Distribution.Continuous Numeric.GSL.Distribution.Discrete Numeric.GSL.Distribution.Common + Numeric.GSL.Fitting.Linear other-modules: C-sources: lib/Numeric/GSL/statistics-aux.c lib/Numeric/GSL/sort-aux.c lib/Numeric/GSL/histogram-aux.c lib/Numeric/GSL/permutation-aux.c lib/Numeric/GSL/distribution-aux.c + lib/Numeric/GSL/fitting-aux.c ghc-prof-options: -auto diff --git a/lib/Numeric/GSL/Fitting/Linear.hs b/lib/Numeric/GSL/Fitting/Linear.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Fitting/Linear.hs @@ -0,0 +1,88 @@ +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.GSL.Fitting.Linear +-- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : haskell.vivian.mcphail gmail com +-- Stability : provisional +-- Portability : uses ffi +-- +-- GSL linear regression functions +-- +----------------------------------------------------------------------------- + +module Numeric.GSL.Fitting.Linear ( + linear, linear_w + ) where + +----------------------------------------------------------------------------- + +import Data.Packed.Vector +import Data.Packed.Matrix +import Data.Packed.Development + +--import Numeric.LinearAlgebra.Linear + +--import Control.Monad(when) + +import Foreign +import Foreign.ForeignPtr +--import Foreign.Marshal.Alloc(alloca) +import Foreign.C.Types(CInt) +--import Foreign.C.String(newCString,peekCString) + +--import GHC.ForeignPtr (mallocPlainForeignPtrBytes) + +--import GHC.Base +--import GHC.IOBase + +--import Prelude hiding(reverse) + +----------------------------------------------------------------------------- + +-- | fits the model Y = C X +linear :: Vector Double -- ^ x data + -> Vector Double -- ^ y data + -> (Vector Double,Vector Double,Matrix Double,Matrix Double,Matrix Double,Matrix Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq) +linear x y = unsafePerformIO $ do + let s = dim x + c0 <- createVector 1 + c1 <- createVector 1 + cov00 <- createMatrix RowMajor s s + cov01 <- createMatrix RowMajor s s + cov11 <- createMatrix RowMajor s s + chi_sq <- createMatrix RowMajor s s + app8 fitting_linear vec x vec y vec c0 vec c1 mat cov00 mat cov01 mat cov11 mat chi_sq "linear" + return (c0,c1,cov00,cov01,cov11,chi_sq) + + +----------------------------------------------------------------------------- + +foreign import ccall "fitting-aux.h linear" fitting_linear :: CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt + +----------------------------------------------------------------------------- + +-- | fits the model Y = C X, with x data weighted +linear_w :: Vector Double -- ^ x data + -> Vector Double -- ^ x weights + -> Vector Double -- ^ y data + -> (Vector Double,Vector Double,Matrix Double,Matrix Double,Matrix Double,Matrix Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq) +linear_w x w y = unsafePerformIO $ do + let s = dim x + c0 <- createVector 1 + c1 <- createVector 1 + cov00 <- createMatrix RowMajor s s + cov01 <- createMatrix RowMajor s s + cov11 <- createMatrix RowMajor s s + chi_sq <- createMatrix RowMajor s s + app9 fitting_linear_w vec x vec w vec y vec c0 vec c1 mat cov00 mat cov01 mat cov11 mat chi_sq "linear_w" + return (c0,c1,cov00,cov01,cov11,chi_sq) + + +----------------------------------------------------------------------------- + +foreign import ccall "fitting-aux.h linear_weighted" fitting_linear_w :: CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt + +----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/fitting-aux.c b/lib/Numeric/GSL/fitting-aux.c new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/fitting-aux.c @@ -0,0 +1,28 @@ +#include + +int linear(int xs, const double* x, int ys, const double* y, + int c0s, double* c0, int c1s, double* c1, + int cov00x, int cov00y, double* cov00, + int cov01x, int cov01y, double* cov01, + int cov11x, int cov11y, double* cov11, + int chi_sqx, int chi_sqy, double* chi_sq) +{ + if (xs != ys || xs != c0s || xs != c1s) return 2000; //BAD_SIZE + + return gsl_fit_linear(x,1,y,1,xs,c0,c1,cov00,cov01,cov11,chi_sq); +} + +int linear_weighted(int xs, const double* x, int ws, const double * w, + int ys, const double* y, + int c0s, double* c0, int c1s, double* c1, + int cov00x, int cov00y, double* cov00, + int cov01x, int cov01y, double* cov01, + int cov11x, int cov11y, double* cov11, + int chi_sqx, int chi_sqy, double* chi_sq) +{ + if (xs != ys || xs != ws || xs != c0s || xs != c1s) return 2000; //BAD_SIZE + + return gsl_fit_wlinear(x,1,w,1,y,1,xs,c0,c1,cov00,cov01,cov11,chi_sq); +} + +