diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2015, Alp Mestanogullari
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Alp Mestanogullari nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,8 @@
+regress
+=======
+
+[![Build Status](https://secure.travis-ci.org/alpmestan/regress.png?branch=master)](http://travis-ci.org/alpmestan/regress)
+
+Perform a linear or logistic regression using automatic differentiation ([ad](http://hackage.haskell.org/package/ad)).
+
+See the haddocks for `Numeric.Regression.Linear` and `Numeric.Regression.Logistic` for documentation and examples.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/regress.cabal b/regress.cabal
new file mode 100644
--- /dev/null
+++ b/regress.cabal
@@ -0,0 +1,33 @@
+name:                regress
+version:             0.1
+synopsis:            Linear and logistic regression through automatic differentiation
+description:
+  Linear and logistic regression through automatic differentiation
+  .
+  See "Numeric.Regression.Linear" and "Numeric.Regression.Logistic" for
+  docs and examples.
+
+homepage:            https://github.com/alpmestan/regress
+license:             BSD3
+license-file:        LICENSE
+author:              Alp Mestanogullari
+maintainer:          alpmestan@gmail.com
+copyright:           Alp Mestanogullari
+category:            Math
+build-type:          Simple
+cabal-version:       >=1.10
+extra-source-files:  README.md
+
+library
+  exposed-modules:
+      Numeric.Regression.Linear
+    , Numeric.Regression.Logistic
+  other-modules:
+    Numeric.Regression.Internal
+
+  build-depends:       base >=4.7 && <5,
+                       ad >= 4.2 && <4.3,
+                       vector >= 0.10 && <0.11
+  hs-source-dirs:      src
+  default-language:    Haskell2010
+  ghc-options:         -O2 -Wall
diff --git a/src/Numeric/Regression/Internal.hs b/src/Numeric/Regression/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Regression/Internal.hs
@@ -0,0 +1,21 @@
+module Numeric.Regression.Internal where
+
+import Control.Applicative
+import Data.Foldable
+import Data.Monoid
+
+data Acc a = Acc {-# UNPACK #-} !Int !a
+
+instance Monoid a => Monoid (Acc a) where
+  mempty = Acc 0 mempty
+
+  Acc m a `mappend` Acc n b = Acc (m + n) (a <> b)
+
+acc :: a -> Acc (Sum a)
+acc = Acc 1 . Sum
+
+dot :: (Applicative v, Foldable v, Num a)
+    => v a
+    -> v a
+    -> a
+dot x y = getSum . foldMap Sum $ liftA2 (*) x y
diff --git a/src/Numeric/Regression/Linear.hs b/src/Numeric/Regression/Linear.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Regression/Linear.hs
@@ -0,0 +1,88 @@
+module Numeric.Regression.Linear
+  (Model, compute, regress) where
+
+import Control.Applicative
+import Data.Foldable
+import Data.Monoid
+import Data.Traversable
+import Numeric.AD
+import Numeric.Regression.Internal
+
+-- | A model using the given @f@ to store parameters of type @a@.
+--   Can be thought of as some kind of vector throughough this
+--   package.
+type Model f a = f a
+
+-- | Compute the predicted value for
+--   the given model on the given observation
+compute :: (Applicative v, Foldable v, Num a)
+        => Model v a -- ^ theta vector, the model's parameters
+        -> v a       -- ^ @x@ vector, with the observed numbers
+        -> a         -- ^ predicted @y@ for this observation
+compute theta x = theta `dot` x
+{-# INLINE compute #-}
+
+-- | Cost function for a linear regression on a single observation
+cost :: (Applicative v, Foldable v, Floating a)
+     => Model v a -- ^ theta vector, the model's parameters
+     -> v a       -- ^ @x@ vector
+     -> a         -- ^ expected @y@ for the observation
+     -> a         -- ^ cost
+cost theta x y = 0.5 * (y - compute theta x) ^ (2 :: Int)
+{-# INLINE compute #-}
+
+-- | Cost function for a linear regression on a set of observations
+totalCost :: (Applicative v, Foldable v, Applicative f, Foldable f, Floating a)
+          => Model v a      -- ^ theta vector, the model's parameters
+          -> f a            -- ^ expected @y@ value for each observation
+          -> f (v a)        -- ^ input data for each observation
+          -> a              -- ^ total cost over all observations
+totalCost theta ys xs =
+  let Acc n (Sum s) = foldMap acc $ liftA2 (cost theta) xs ys
+  in s / fromIntegral n
+{-# INLINE totalCost #-}
+
+-- | Given some observed \"predictions\" @ys@, the corresponding
+--   input values @xs@ and initial values for the model's parameters @theta0@,
+--
+-- > regress ys xs theta0
+--
+-- returns a stream of values for the parameters that'll fit the data better
+-- and better.
+--
+-- Example:
+--
+-- @
+-- -- the theta we're approximating
+-- realtheta :: Model V.Vector Double
+-- realtheta = V.fromList [1.0, 2.0, 3.0]
+--
+-- -- let's start there and make 'regress'
+-- -- get values that better fit the input data
+-- theta0 :: Model V.Vector Double
+-- theta0 = V.fromList [0.2, 3.0, 2.23]
+--
+-- -- input data. (output value, vector of values for each input)
+-- ys_ex :: V.Vector Double
+-- xs_ex :: V.Vector (V.Vector Double)
+-- (ys_ex, xs_ex) = V.unzip . V.fromList $
+--   [ (3, V.fromList [0, 0, 1])
+--   , (1, V.fromList [1, 0, 0])
+--   , (2, V.fromList [0, 1, 0])
+--   , (6, V.fromList [1, 1, 1])
+--   ]
+--
+-- -- stream of increasingly accurate parameters
+-- thetaApproxs :: [Model V.Vector Double]
+-- thetaApproxs = learnAll ys_ex xs_ex theta0
+-- @
+regress :: (Traversable v, Applicative v, Foldable v, Applicative f, Foldable f, Ord a, Floating a)
+        => f a         -- ^ expected @y@ value for each observation
+        -> f (v a)     -- ^ input data for each observation
+        -> Model v a   -- ^ initial parameters for the model, from which we'll improve
+        -> [Model v a] -- ^ a stream of increasingly accurate values
+                       --   for the model's parameter to better fit the observations.
+regress ys xs t0 =
+  gradientDescent (\theta -> totalCost theta (fmap auto ys) (fmap (fmap auto) xs))
+                  t0
+{-# INLINE regress #-}
diff --git a/src/Numeric/Regression/Logistic.hs b/src/Numeric/Regression/Logistic.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Regression/Logistic.hs
@@ -0,0 +1,83 @@
+module Numeric.Regression.Logistic
+  (Model, regress) where
+
+import Control.Applicative
+import Data.Foldable
+import Data.Monoid
+import Data.Traversable
+import Numeric.AD
+import Numeric.Regression.Internal
+
+-- | A model using the given @f@ to store parameters of type @a@.
+--   Can be thought of as some kind of vector throughough this
+--   package.
+type Model f a = f a
+
+logit :: Floating a => a -> a
+logit x = 1 / (1 + exp (negate x))
+{-# INLINE logit #-}
+
+logLikelihood :: (Applicative v, Foldable v, Floating a)
+              => Model v a -- theta vector
+              -> a         -- y
+              -> v a       -- x vector (observation)
+              -> a
+logLikelihood theta y x =
+  y * log (logit z) + (1 - y) * log (1 - logit z)
+
+  where z = theta `dot` x
+{-# INLINE logLikelihood #-}
+
+totalLogLikelihood :: (Applicative v, Foldable v, Applicative f, Foldable f, Floating a)
+                   => a -- delta
+                   -> Model v a
+                   -> f a
+                   -> f (v a)
+                   -> a
+totalLogLikelihood delta theta ys xs =
+  (a - delta * b) / fromIntegral n
+
+  where Acc n (Sum a) = foldMap acc $ liftA2 (logLikelihood theta) ys xs
+        b = (/2) . getSum $ foldMap (\x -> Sum (x^(2::Int))) theta
+{-# INLINE totalLogLikelihood #-}
+
+-- | Given some observed \"predictions\" @ys@, the corresponding
+--   input values @xs@ and initial values for the model's parameters @theta0@,
+--
+-- > regress ys xs theta0
+--
+-- returns a stream of values for the parameters that'll fit the data better
+-- and better.
+--
+-- Example:
+--
+-- @
+-- ys_ex :: [Double]
+-- xs_ex :: [[Double]]
+-- (ys_ex, xs_ex) = unzip $
+--   [ (1, [1, 1])
+--   , (0, [-1, -2])
+--   , (1, [2, 5])
+--   , (0, [-1, 1])
+--   , (1, [2, -1])
+--   , (1, [1, -10])
+--   , (0, [-0.1, 30])
+--   ]
+--
+-- t0 :: [Double]
+-- t0 = [1, 0.1]
+--
+-- approxs' :: [Model [] Double]
+-- approxs' = learn 0.1 ys_ex xs_ex t0
+-- @
+regress :: (Traversable v, Applicative v, Foldable f, Applicative f, Ord a, Floating a)
+        => a           -- ^ learning rate
+        -> f a         -- ^ expect prediction for each observation
+        -> f (v a)     -- ^ input data for each observation
+        -> Model v a   -- ^ initial values for the model's parameters
+        -> [Model v a] -- ^ stream of increasingly accurate values for
+                       --   the model's parameters
+regress delta ys xs =
+  gradientAscent $ \theta ->
+    totalLogLikelihood (auto delta) theta (fmap auto ys) (fmap (fmap auto) xs)
+{-# INLINE regress #-}