Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
branch: master
Fetching contributors…

Cannot retrieve contributors at this time

185 lines (148 sloc) 8.869 kb
-----------------------------------------------------------------------------
-- |
-- Module : Numeric.GSL.Fitting.Linear
-- Copyright : (c) A. V. H. McPhail 2010
-- License : BSD3
--
-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com
-- Stability : provisional
-- Portability : uses ffi
--
-- GSL linear regression functions
--
-- <http://www.gnu.org/software/gsl/manual/>
--
-----------------------------------------------------------------------------
module Numeric.GSL.Fitting.Linear (
linear, linear_w, linear_est,
multifit, multifit_w, multifit_est
) 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)
import System.IO.Unsafe(unsafePerformIO)
-----------------------------------------------------------------------------
-- | fits the model Y = C X
linear :: Vector Double -- ^ x data
-> Vector Double -- ^ y data
-> (Double,Double,Double,Double,Double,Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq)
linear x y = unsafePerformIO $ do
alloca $ \c0 ->
alloca $ \c1 ->
alloca $ \chi_sq ->
alloca $ \cov00 ->
alloca $ \cov01 ->
alloca $ \cov11 -> do
app2 (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) vec x vec y "linear"
c0' <- peek c0
c1' <- peek c1
cov00' <- peek cov00
cov01' <- peek cov01
cov11' <- peek cov11
chi_sq' <- peek chi_sq
return (c0',c1',cov00',cov01',cov11',chi_sq')
-----------------------------------------------------------------------------
foreign import ccall "fitting-aux.h linear" fitting_linear :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> Ptr Double -> 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
-> (Double,Double,Double,Double,Double,Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq)
linear_w x w y = unsafePerformIO $ do
alloca $ \c0 ->
alloca $ \c1 ->
alloca $ \chi_sq ->
alloca $ \cov00 ->
alloca $ \cov01 ->
alloca $ \cov11 -> do
app3 (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) vec x vec w vec y "linear_w"
c0' <- peek c0
c1' <- peek c1
cov00' <- peek cov00
cov01' <- peek cov01
cov11' <- peek cov11
chi_sq' <- peek chi_sq
return (c0',c1',cov00',cov01',cov11',chi_sq')
-----------------------------------------------------------------------------
foreign import ccall "fitting-aux.h linear_weighted" fitting_linear_w :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | computes the fitted function and standard deviation at the input point
linear_est :: Double -- ^ x data point
-> Double -- ^ c0
-> Double -- ^ c1
-> Double -- ^ cov00
-> Double -- ^ cov01
-> Double -- ^ cov11
-> (Double,Double) -- ^ (y,error)
linear_est x c0 c1 cov00 cov01 cov11 = unsafePerformIO $ do
alloca $ \y ->
alloca $ \e -> do
check "linear_est" $ fitting_linear_est x c0 c1 cov00 cov01 cov11 y e
y' <- peek y
e' <- peek e
return (y',e')
-----------------------------------------------------------------------------
foreign import ccall "fitting-aux.h linear_estimate" fitting_linear_est :: Double -> Double -> Double -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | fit the model Y = C X, with design matrix X
-- | X is a design matrix X_{ij} = x_j(i) with i observations and p predictors
-- | a polynomial would be X_{ij} = x_i^j
-- | a fourier series would be X_{ij} = sin (\omega_j x_i)
multifit :: Matrix Double -- ^ design matrix (X)
-> Vector Double -- ^ observations
-> (Vector Double,Matrix Double,Double) -- ^ (coefficients,covariance,chi_sq)
multifit x y = unsafePerformIO $ do
let p = cols x
cov <- createMatrix RowMajor p p
c <- createVector p
alloca$ \chi_sq -> do
app4 (fitting_multifit chi_sq) mat (cmat x) vec y vec c mat cov "multifit"
chi_sq' <- peek chi_sq
return (c,cov,chi_sq')
-----------------------------------------------------------------------------
foreign import ccall "fitting_aux.h multifit" fitting_multifit :: Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | fit the model Y = C X, with design matrix X, and x weighted
multifit_w :: Matrix Double -- ^ design matrix (X)
-> Vector Double -- ^ weights
-> Vector Double -- ^ observations
-> (Vector Double,Matrix Double,Double) -- ^ (coefficients,covariance,chi_sq)
multifit_w x w y = unsafePerformIO $ do
let p = cols x
cov <- createMatrix RowMajor p p
c <- createVector p
alloca$ \chi_sq -> do
app5 (fitting_multifit_w chi_sq) mat (cmat x) vec w vec y vec c mat cov "multifit"
chi_sq' <- peek chi_sq
return (c,cov,chi_sq')
-----------------------------------------------------------------------------
foreign import ccall "fitting_aux.h multifit_weighted" fitting_multifit_w :: Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | computes the fitted function and standard deviation at the input point
multifit_est :: Vector Double -- ^ input point
-> Vector Double -- ^ the coefficients
-> Matrix Double -- ^ the covariance matrix
-> (Double,Double) -- ^ (y,y_error)
multifit_est x c cov = unsafePerformIO $ do
alloca $ \y ->
alloca $ \e -> do
app3 (fitting_multifit_est y e) vec x vec c mat cov "multifit_estimate"
y' <- peek y
e' <- peek e
return (y',e')
-----------------------------------------------------------------------------
foreign import ccall "fitting_aux.h multifit_estimate" fitting_multifit_est :: Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
Jump to Line
Something went wrong with that request. Please try again.