diff --git a/Statistics/LinearRegression.hs b/Statistics/LinearRegression.hs --- a/Statistics/LinearRegression.hs +++ b/Statistics/LinearRegression.hs @@ -1,6 +1,12 @@ {-# LANGUAGE BangPatterns #-} -module Statistics.LinearRegression (linearRegressionRSqr, linearRegression, correl, covar) where +module Statistics.LinearRegression ( + linearRegression, + linearRegressionRSqr, + linearRegressionTLS, + correl, + covar, + ) where import qualified Data.Vector.Unboxed as U import qualified Statistics.Sample as S @@ -9,13 +15,11 @@ -- | Covariance of two samples covar :: S.Sample -> S.Sample -> Double -covar xs ys = U.sum (U.zipWith (*) (U.map f1 xs) (U.map f2 ys)) / (n-1) +covar xs ys = U.sum (U.zipWith (*) (U.map (subtract m1) xs) (U.map (subtract m2) ys)) / (n-1) where !n = fromIntegral $ U.length xs !m1 = S.mean xs !m2 = S.mean ys - f1 = \x -> (x - m1) - f2 = \x -> (x - m2) {-# INLINE covar #-} @@ -35,7 +39,7 @@ linearRegressionRSqr :: S.Sample -> S.Sample -> (Double, Double, Double) linearRegressionRSqr xs ys = (alpha, beta, r*r) where - !c = U.sum (U.zipWith (*) (U.map (subtract m1) xs) (U.map (subtract m2) ys)) / (n-1) + !c = covar xs ys !r = c / (sx * sy) !m1 = S.mean xs !m2 = S.mean ys @@ -48,17 +52,24 @@ -- | Simple linear regression between 2 samples. -- Takes two vectors Y={yi} and X={xi} and returns --- (alpha, beta, r*r) such that Y = alpha + beta*X +-- (alpha, beta) such that Y = alpha + beta*X linearRegression :: S.Sample -> S.Sample -> (Double, Double) linearRegression xs ys = (alpha, beta) where - !c = U.sum (U.zipWith (*) (U.map (subtract m1) xs) (U.map (subtract m2) ys)) / (n-1) - !r = c / (sx * sy) + (alpha, beta, _) = linearRegressionRSqr xs ys +{-# INLINE linearRegression #-} + +-- | Total Least Squares (TLS) linear regression. +-- Assumes x-axis values (and not just y-axis values) are random variables and that both variables have similar distributions. +-- interface is the same as linearRegression. +linearRegressionTLS :: S.Sample -> S.Sample -> (Double,Double) +linearRegressionTLS xs ys = (alpha, beta) + where + !c = covar xs ys + !b = (S.varianceUnbiased xs - (S.varianceUnbiased ys)) / c !m1 = S.mean xs !m2 = S.mean ys - !sx = S.stdDev xs - !sy = S.stdDev ys - !n = fromIntegral $ U.length xs - !beta = r * sy / sx + !betas = [(-b - sqrt(b^2+4))/2,(-b + sqrt(b^2+4)) /2] + !beta = if c > 0 then maximum betas else minimum betas !alpha = m2 - beta * m1 -{-# INLINE linearRegression #-} +{-# INLINE linearRegressionTLS #-} diff --git a/statistics-linreg.cabal b/statistics-linreg.cabal --- a/statistics-linreg.cabal +++ b/statistics-linreg.cabal @@ -1,70 +1,44 @@ --- linreg.cabal auto-generated by cabal init. For additional options, --- see --- http://www.haskell.org/cabal/release/cabal-latest/doc/users-guide/authors.html#pkg-descr. --- The name of the package. Name: statistics-linreg - --- The package version. See the Haskell package versioning policy --- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for --- standards guiding when and how versions should be incremented. -Version: 0.2.1 - --- A short (one-line) description of the package. -Synopsis: Linear regression between two samples, based on the 'statistics' package - --- A longer description of the package. -Description: Provides a function to perform a linear regression between 2 samples, see the documentation of the linearRegression function. This library is based on the 'statistics' package. +Version: 0.2.2 +Synopsis: Linear regression between two samples, based on the 'statistics' package. +Description: Provides functions to perform a linear regression between 2 samples, see the documentation of the linearRegression functions. This library is based on the 'statistics' package. . - 0.2.*: added the r-squared version and improved the performances - --- URL for the project homepage or repository. + * 0.2.2: added the Total-Least-Squares version and made some refactoring to eliminate code duplication + . + * 0.2.1: added the r-squared version and improved the performances. + . + Code sample: + . + > import qualified Data.Vector.Unboxed as U + > + > test :: Int -> IO () + > test k = do + > let n = 10000000 + > let a = k*n + 1 + > let b = (k+1)*n + > let xs = U.fromList [a..b] + > let ys = U.map (\x -> x*100 + 2000) xs + > -- thus 100 and 2000 are the alpha and beta we want + > putStrLn "linearRegression:" + > print $ linearRegression xs ys + . + The r-squared and Total-Least-Squares versions work the same way. Homepage: http://github.com/alpmestan/statistics-linreg Bug-reports: https://github.com/alpmestan/statistics-linreg/issues - --- The license under which the package is released. License: MIT - --- The file containing the license text. License-file: LICENSE - --- The package author(s). -Author: Alp Mestanogullari - --- An email address to which users can send suggestions, bug reports, --- and patches. +Author: Alp Mestanogullari , Uri Barenholz Maintainer: Alp Mestanogullari - --- A copyright notice. -Copyright: 2010-2011 Alp Mestanogullari - --- Stability of the pakcage (experimental, provisional, stable...) +Copyright: 2010-2012 Alp Mestanogullari Stability: Experimental - Category: Math, Statistics - Build-type: Simple - --- Extra files to be distributed with the package, such as examples or --- a README. --- Extra-source-files: - --- Constraint on the version of Cabal needed to build this package. Cabal-version: >=1.6 - Library - -- Modules exported by the library. Exposed-modules: Statistics.LinearRegression - - -- Packages needed in order to build this package. Build-depends: statistics >= 0.5, vector >= 0.5, base >= 4 && < 5 - - Ghc-options: -funbox-strict-fields - -- Modules not exported by this package. - -- Other-modules: - - -- Extra tools (e.g. alex, hsc2hs, ...) needed to build the source. - -- Build-tools: + Ghc-options: -funbox-strict-fields -O2 source-repository head type: git