diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,11 @@
+# Changelog for `welford`
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to the
+[Haskell Package Versioning Policy](https://pvp.haskell.org/).
+
+## Unreleased
+
+## 0.1.0.0 - YYYY-MM-DD
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Author name here (c) 2022
+
+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 Author name here 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,1 @@
+# Welford: Online mean and variance computation
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/src/Lib.hs b/src/Lib.hs
new file mode 100644
--- /dev/null
+++ b/src/Lib.hs
@@ -0,0 +1,6 @@
+module Lib
+    ( someFunc
+    ) where
+
+someFunc :: IO ()
+someFunc = putStrLn "someFunc"
diff --git a/src/Statistics/Sample/WelfordOnlineMeanVariance.hs b/src/Statistics/Sample/WelfordOnlineMeanVariance.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Sample/WelfordOnlineMeanVariance.hs
@@ -0,0 +1,167 @@
+{-# LANGUAGE DeriveAnyClass    #-}
+{-# LANGUAGE DeriveGeneric     #-}
+{-# LANGUAGE FlexibleInstances #-}
+module Statistics.Sample.WelfordOnlineMeanVariance
+  ( WelfordExistingAggregate(..)
+  , addValue
+  , finalize
+  , nextValue
+  , newWelfordAggregate
+  , Mean
+  , Variance
+  , SampleVariance
+  ) where
+
+import           Control.DeepSeq
+import           Data.Serialize
+import qualified Data.Vector          as VB
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Unboxed  as VU
+import           GHC.Generics
+
+type Mean a = a
+type Variance a = a
+type SampleVariance a = a
+
+
+-- | For the storage of required information.
+data WelfordExistingAggregate a
+  = WelfordExistingAggregateEmpty
+  | WelfordExistingAggregate
+      { welfordCount :: !Int
+      , welfordMean  :: !a
+      , welfordM2    :: !a
+      }
+  deriving (Eq, Show, Read, Generic, NFData, Serialize)
+
+
+-- | Create a new empty Aggreate for the calculation.
+newWelfordAggregate :: WelfordExistingAggregate a
+newWelfordAggregate = WelfordExistingAggregateEmpty
+
+
+-- | Class for all data strucutres that can be used to computer the Welford approximation. For instance, this can be used to compute the Welford algorithm on a `Vector`s of `Fractional`, while only
+-- requiring to handle one `WelfordExistingAggregate`.
+class WelfordOnline a where
+  plus :: a -> a -> a
+  minus :: a -> a -> a
+  multiply :: a -> a -> a
+  divideInt :: a -> Int -> a
+
+
+instance (WelfordOnline a) => WelfordOnline (VB.Vector a) where
+  plus          = VB.zipWith plus
+  {-# INLINE plus #-}
+  minus         = VB.zipWith minus
+  {-# INLINE minus #-}
+  multiply      = VB.zipWith multiply
+  {-# INLINE multiply #-}
+  divideInt x i = VB.map (`divideInt` i) x
+  {-# INLINE divideInt #-}
+
+instance (WelfordOnline a, VS.Storable a) => WelfordOnline (VS.Vector a) where
+  plus          = VS.zipWith plus
+  {-# INLINE plus #-}
+  minus         = VS.zipWith minus
+  {-# INLINE minus #-}
+  multiply      = VS.zipWith multiply
+  {-# INLINE multiply #-}
+  divideInt x i = VS.map (`divideInt` i) x
+  {-# INLINE divideInt #-}
+
+
+instance (WelfordOnline a, VU.Unbox a) => WelfordOnline (VU.Vector a) where
+  plus          = VU.zipWith plus
+  {-# INLINE plus #-}
+  minus         = VU.zipWith minus
+  {-# INLINE minus #-}
+  multiply      = VU.zipWith multiply
+  {-# INLINE multiply #-}
+  divideInt x i = VU.map (`divideInt` i) x
+  {-# INLINE divideInt #-}
+
+
+instance WelfordOnline Double where
+  plus = (+)
+  {-# INLINE plus #-}
+  minus = (-)
+  {-# INLINE minus #-}
+  multiply = (*)
+  {-# INLINE multiply #-}
+  divideInt x i = x / fromIntegral i
+  {-# INLINE divideInt #-}
+
+instance WelfordOnline Float where
+  plus = (+)
+  {-# INLINE plus #-}
+  minus = (-)
+  {-# INLINE minus #-}
+  multiply = (*)
+  {-# INLINE multiply #-}
+  divideInt x i = x / fromIntegral i
+  {-# INLINE divideInt #-}
+
+instance WelfordOnline Rational where
+  plus = (+)
+  {-# INLINE plus #-}
+  minus = (-)
+  {-# INLINE minus #-}
+  multiply = (*)
+  {-# INLINE multiply #-}
+  divideInt x i = x / fromIntegral i
+  {-# INLINE divideInt #-}
+
+
+-- | Add one value to the current aggregate.
+addValue :: (WelfordOnline a) => WelfordExistingAggregate a -> a -> WelfordExistingAggregate a
+addValue (WelfordExistingAggregate count mean m2) val =
+  let count' = count + 1
+      delta = val `minus` mean
+      mean' = mean `plus` (delta `divideInt` count')
+      delta2 = val `minus` mean'
+      m2' = m2 `plus` (delta `multiply` delta2)
+   in WelfordExistingAggregate count' mean' m2'
+addValue WelfordExistingAggregateEmpty val =
+  let count' = 1
+      delta' = val
+      mean' = delta' `divideInt` count'
+      delta2' = val `minus` mean'
+      m2' = delta' `multiply` delta2'
+   in WelfordExistingAggregate count' mean' m2'
+
+
+-- | Calculate mean, variance and sample variance from aggregate.
+finalize :: (WelfordOnline a) => WelfordExistingAggregate a -> (Mean a, Variance a, SampleVariance a)
+finalize (WelfordExistingAggregate count mean m2)
+  | count < 2 = (mean, m2, m2)
+  | otherwise = (mean, m2 `divideInt` count, m2 `divideInt` (count - 1))
+finalize WelfordExistingAggregateEmpty = error "finalize: Emtpy Welford Online Aggreate. Add data first!"
+
+-- | Add a new sample to the aggregate and compute mean and variances.
+nextValue :: (WelfordOnline a) => WelfordExistingAggregate a -> a -> (WelfordExistingAggregate a, (Mean a, Variance a, SampleVariance a))
+nextValue agg val =
+  let agg' = addValue agg val
+   in (agg', finalize agg')
+
+
+-- # For a new value newValue, compute the new count, new mean, the new M2.
+-- # mean accumulates the mean of the entire dataset
+-- # M2 aggregates the squared distance from the mean
+-- # count aggregates the number of samples seen so far
+-- def update(existingAggregate, newValue):
+--     (count, mean, M2) = existingAggregate
+--     count' += 1
+--     delta = newValue - mean
+--     mean' += delta / count'
+--     delta2 = newValue - mean'
+--     M2' += delta * delta2
+--     return (count', mean', M2')
+
+-- # Retrieve the mean, variance and sample variance from an aggregate
+-- def finalize(existingAggregate):
+--     (count, mean, M2) = existingAggregate
+--     if count < 2:
+--         return float("nan")
+--     else:
+--         (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))
+--         return (mean, variance, sampleVariance)
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,1 @@
+{-# OPTIONS_GHC -F -pgmF tasty-discover #-}
diff --git a/test/Statistics/Sample/WelfordOnlineTest.hs b/test/Statistics/Sample/WelfordOnlineTest.hs
new file mode 100644
--- /dev/null
+++ b/test/Statistics/Sample/WelfordOnlineTest.hs
@@ -0,0 +1,53 @@
+module Statistics.Sample.WelfordOnlineTest where
+
+import           Control.Monad
+import           Data.List                                   (foldl')
+import qualified Data.Vector                                 as VB
+import           Debug.Trace
+import           Test.Tasty.QuickCheck
+
+import           Statistics.Sample.WelfordOnlineMeanVariance
+
+main :: IO ()
+main = putStrLn "Test suite not yet implemented"
+
+
+-- | Recursively generate candles
+generateNValues :: Int -> Gen [Double]
+generateNValues n
+  | n <= 0 = return []
+  | otherwise =  do
+      x <- arbitrary
+      (x:) <$> generateNValues (n - 1)
+
+epsEqWith :: (Ord a, Fractional a, Show a) => a -> a -> a -> Property
+epsEqWith eps a b
+  | abs (a - b) <= eps = a === a
+  | otherwise = a === b
+
+
+prop_MeanAndVariance :: Gen Property
+prop_MeanAndVariance = do
+  vals <- sized $ \n -> generateNValues (10 + 5 * n)
+  let n = fromIntegral (length vals)
+      mean = sum vals / n
+      var = sum (map (\x -> (x - mean) ^ 2) vals) / (n - 1)
+      (wMean, _, wVarSample) = snd $ foldl' (nextValue . fst) (WelfordExistingAggregateEmpty, (0, 0, 0)) vals
+      eps = min 0.01 $ max 0.001 (0.001 * mean)
+  return $ epsEqWith eps wMean mean .&&. epsEqWith eps wVarSample var
+
+
+prop_MeanAndVarianceVector :: Gen Property
+prop_MeanAndVarianceVector = do
+  len <- chooseInt (1, 100)
+  vecs <- sized $ \n -> replicateM len (VB.fromList <$> generateNValues (10 + 5 * n))
+  let mkRes vals =
+        let n = fromIntegral (VB.length vals)
+            mean = VB.sum vals / n
+            var = VB.sum (VB.map (\x -> (x - mean) ^ 2) vals) / (n - 1)
+            (wMean, _, wVarSample) = snd $ foldl' (nextValue . fst) (WelfordExistingAggregateEmpty, (0, 0, 0)) (VB.toList vals)
+            eps = min 0.01 $ max 0.001 (0.001 * mean)
+        in (eps, wMean, mean, wVarSample, var)
+      ress = map mkRes vecs
+  return $
+    foldl1 (.&&.) $ map (\(eps, wMean, mean, wVarSample, var) -> epsEqWith eps wMean mean .&&. epsEqWith eps wVarSample var) ress
diff --git a/welford-online-mean-variance.cabal b/welford-online-mean-variance.cabal
new file mode 100644
--- /dev/null
+++ b/welford-online-mean-variance.cabal
@@ -0,0 +1,63 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.34.4.
+--
+-- see: https://github.com/sol/hpack
+
+name:           welford-online-mean-variance
+version:        0.1.0.0
+synopsis:       Online computation of mean and variance using the Welford algorithm.
+description:    Please see the README on GitHub at <https://github.com/githubuser/welford-online-mean-variance#readme>
+category:       Statistics
+homepage:       https://github.com/schnecki/welford-online-mean-variance#readme
+bug-reports:    https://github.com/schnecki/welford-online-mean-variance/issues
+author:         Manuel Schneckenreither
+maintainer:     manuel.schnecki@gmail.com
+copyright:      2022 Manuel Schneckenreither
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+extra-source-files:
+    README.md
+    CHANGELOG.md
+
+source-repository head
+  type: git
+  location: https://github.com/schnecki/welford-online-mean-variance
+
+library
+  exposed-modules:
+      Lib
+      Statistics.Sample.WelfordOnlineMeanVariance
+  other-modules:
+      Paths_welford_online_mean_variance
+  hs-source-dirs:
+      src
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints
+  build-depends:
+      base >=4.7 && <5
+    , cereal
+    , deepseq
+    , vector
+  default-language: Haskell2010
+
+test-suite welford-online-mean-variance-test
+  type: exitcode-stdio-1.0
+  main-is: Spec.hs
+  other-modules:
+      Statistics.Sample.WelfordOnlineTest
+      Paths_welford_online_mean_variance
+  hs-source-dirs:
+      test
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      QuickCheck
+    , base >=4.7 && <5
+    , cereal
+    , deepseq
+    , tasty
+    , tasty-discover
+    , tasty-quickcheck
+    , vector
+    , welford-online-mean-variance
+  default-language: Haskell2010