fastbayes (empty) → 0.1.0.0

raw patch · 6 files changed

+194/−0 lines, 6 filesdep +basedep +hmatrixdep +vectorsetup-changed

Dependencies added: base, hmatrix, vector

Files

• LICENSE +20/−0
• README.md +21/−0
• Setup.hs +2/−0
• fastbayes.cabal +46/−0
• src/Statistics/FastBayes.hs +3/−0
• src/Statistics/FastBayes/Linear.hs +102/−0

@@ -0,0 +1,20 @@+Copyright (c) 2014 Melinae, Inc++Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the+"Software"), to deal in the Software without restriction, including+without limitation the rights to use, copy, modify, merge, publish,+distribute, sublicense, and/or sell copies of the Software, and to+permit persons to whom the Software is furnished to do so, subject to+the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

@@ -0,0 +1,21 @@+# fastbayes++TODO: Write description here++## Installation++TODO: Write installation instructions here++## Usage++TODO: Write usage instructions here++## How to run tests++```+cabal configure --enable-tests && cabal build && cabal test+```++## Contributing++TODO: Write contribution instructions here

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,46 @@+-- Initial fastbayes.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++-- The name of the package.+name:                fastbayes+version:             0.1.0.0++-- A short (one-line) description of the package.+synopsis:            Bayesian modeling algorithms accelerated for particular model structures++-- A longer description of the package.+description:         General-purpose sampling approaches like Gibbs sampling are very useful for models that have not been studied extensively. But for some cases, specialized algorithms are available because of the model's generality (/e.g./, linear regression) or niche popularity (/e.g./, Latent Dirichlet Allocation). This package is an effort to collect such algorithms in one place.++-- URL for the project homepage or repository.+homepage:            https://github.com/cscherrer/fastbayes++license:             MIT+license-file:        LICENSE+author:              Chad Scherrer+maintainer:          chad.scherrer@gmail.com+copyright:           Copyright (c) 2014 Melinae, Inc+category:            Statistics+build-type:          Simple+extra-source-files:  README.md+cabal-version:       >=1.10+++library+  -- Modules exported by the library.+  exposed-modules:     Statistics.FastBayes, Statistics.FastBayes.Linear+  +  -- Modules included in this library but not exported.+  -- other-modules:       +  +  -- LANGUAGE extensions used by modules in this package.+  other-extensions:    UnicodeSyntax, BangPatterns+  +  -- Other library packages from which modules are imported.+  build-depends:       base >=4.7 && <4.8, vector >=0.10 && <0.11, hmatrix >=0.16 && <0.17+  +  -- Directories containing source files.+  hs-source-dirs:      src+  +  -- Base language which the package is written in.+  default-language:    Haskell2010+  

@@ -0,0 +1,3 @@+module Statistics.FastBayes where++import Statistics.FastBayes.Linear()

@@ -0,0 +1,102 @@+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE BangPatterns #-}++{-|+Module      : Statistics.FastBayes.Linear+Description : Bayesian linear regression via maximum marginal likelihood.+Copyright   : (c) Melinae, 2014+                  Chad Scherrer, 2014+License     : MIT+Maintainer  : chad.scherrer@gmail.com+Stability   : experimental+Portability : POSIX++This module gives an implementation of Bayesian linear regression, with the scale of the prior chosen by marginal likelihood.++The inputs for a Bayesian linear model are identical to those of a classical linear model, except that in addition to a design matrix and response, we must also specify a prior distribution on the weights and the noise. This leaves us with an open question of how these should be specified.++In his book /Pattern Recognition and Machine Learning/, Christopher Bishop provides details for an approach that simplifies the situation significantly, and allows for much faster inference. The structure of the linear model allows us to integrate the posterior over the weights, resulting in the /marginal likelihood/, expressed as a function of the prior precision and noise precision. This, in turn, can be easily optimized.+-}++module Statistics.FastBayes.Linear +  ( Fit+  , fit+  , design+  , response+  , priorPrecision        +  , noisePrecision        +  , numEffectiveParameters+  , logEvidence           +  , mapWeights            +  , hessian+  )+  where++import qualified Data.Vector.Storable as V+import Numeric.LinearAlgebra+++data Fit = Fit+  { design                 :: Matrix Double -- ^The design matrix used for the fit. +  , response               :: Vector Double -- ^The response vector used for the fit+  , priorPrecision         :: Double        -- ^The precision (inverse variance) of the prior distribution, determined by maximizing the marginal likelihood+  , noisePrecision         :: Double        -- ^The precision (inverse variance) of the noise+  , numEffectiveParameters :: Double        -- ^The number of effective parameters in the model+  , logEvidence            :: Double        -- ^The log of the evidence, which is useful for model comparison (different features, same response)+  , mapWeights             :: Vector Double -- ^The MAP (maximum a posteriori) values for the paramter weights+  , hessian                :: Matrix Double -- ^The Hessian (matrix of second derivatives) for the posterior distribution+  }+  deriving Show++-- |@fit lim x y@ fits a Bayesian linear model to a design matrix @x@ and response vector @y@. This is an iterative algorithm, resulting in a sequence (list) of (α,β) values. Here α is the prior precision, and β is the noise precision. The @lim@ function passed in is used to specify how the limit of this sequence should be computed.+fit :: +  ([(Double, Double)] → (Double, Double))   -- ^How to take the limit of the (α,β) sequence. A simple approach is, /e.g./, @fit (!! 1000) x y@+  → Matrix Double                           -- ^The design matrix (each column is a feature)+  → Vector Double                           -- ^The response vector+  → Fit+fit lim x y = Fit x y α β γ logEv m h+  where+  n = rows x+  p = cols x+  α0 = 1.0+  β0 = 1.0++  -- A vector of the eigenvalues of xtx+  (_,sqrtEigs,_) = compactSVD x+  eigs = V.map square sqrtEigs++  xtx = trans x <> x+  xty = trans x <> y+  getHessian a b = diag (V.replicate p a) + scale b xtx++  m = scale β $ h <\> xty+  +  go :: Double → Double → [(Double, Double)]+  go a0 b0 = (a0, b0) : go a b+    where+    h0 = getHessian a0 b0+    m0 = scale b0 $ h0 <\> xty +    c = V.sum $ V.map (\x' → x' / (a0 + x')) eigs+    a = c / (m0 <.> m0)+    b = recip $ (normSq $ y - x <> m0) / (fromIntegral n - c)+  +  γ = V.sum $ V.map (\x' → x' / (α + x')) eigs++  h = getHessian α β +  logEv = 0.5 * +    ( fromIntegral p * log α +    + fromIntegral n * log β +    - (β * normSq (y - x <> m) + α * (m <.> m))+    - logDetH +    - fromIntegral n * log (2*pi)+    )+    where+    (_,(logDetH, _)) = invlndet h++  (α, β) = lim $ go α0 β0++square :: Double → Double+square x = x * x++normSq :: Vector Double → Double+normSq x = x <.> x