probability-dist (empty) → 0.1.0.0
raw patch · 9 files changed
+1750/−0 lines, 9 filesdep +basedep +probability-dist
Dependencies added: base, probability-dist
Files
- CHANGELOG.md +23/−0
- LICENSE +29/−0
- README.md +351/−0
- probability-dist.cabal +47/−0
- src/Probability/Continuous.hs +340/−0
- src/Probability/Discrete.hs +420/−0
- src/Probability/Error.hs +19/−0
- src/Probability/Math.hs +128/−0
- test/Main.hs +393/−0
+ CHANGELOG.md view
@@ -0,0 +1,23 @@+# Changelog++All notable changes to this project will be documented in this file.++## 0.1.0.0 — Initial release++* Initial public release.+* Added discrete probability distributions:++ * Binomial+ * Negative Binomial+ * Geometric+ * Multinomial+* Added continuous probability distributions:++ * Normal+ * Exponential+ * Gamma+ * Uniform+* Added probability mass, density, and cumulative distribution functions where applicable.+* Added explicit error handling through `Either`.+* Added numerical support for logarithmic factorial, combination, Gamma, Beta, and error functions.+* Added Cabal test suite covering numerical results, boundary conditions, and invalid inputs.
+ LICENSE view
@@ -0,0 +1,29 @@+BSD 3-Clause License++Copyright (c) 2026, Barış Barış++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice,+ this list of conditions and the following disclaimer.++2. 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.++3. Neither the name of the copyright holder nor the names of its+ 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 HOLDER 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.
+ README.md view
@@ -0,0 +1,351 @@+# probability-dist++A Haskell library providing probability distributions and related probability functions.++The library is designed around a small, explicit public API. Discrete and continuous probability distributions are exposed through separate modules, allowing applications to depend only on the functionality they need.++## Features++* Discrete probability distributions+* Continuous probability distributions+* Probability mass functions (PMF)+* Probability density functions (PDF)+* Cumulative distribution functions (CDF)+* Distribution-specific moments where implemented+* Numerically stable logarithmic calculations for combinatorial and special functions+* Explicit error handling through `Either`++## Requirements++* GHC 9.10.3 or compatible GHC version+* Cabal 3.16.1.0 or compatible Cabal version+* `base >= 4.20 && < 5`++## Installation++Clone the repository and build it with Cabal:++```bash+git clone <repository-url>+cd probability-dist+cabal build+```++Run the test suite with:++```bash+cabal test+```++The package can also be built as a source distribution:++```bash+cabal sdist+```++## Public API++The library exposes two public modules.++### Discrete distributions++```haskell+import Probability.Discrete+```++This module contains the discrete probability distributions implemented by the package.++Current distributions include:++* Binomial+* Negative Binomial+* Geometric+* Multinomial++### Continuous distributions++```haskell+import Probability.Continuous+```++This module contains the continuous probability distributions implemented by the package.++Current distributions include:++* Normal+* Exponential+* Gamma+* Uniform++The modules are intentionally separated so that an application requiring only discrete or only continuous distributions does not need to import both APIs.++## Basic Usage++### Binomial distribution++The binomial PMF is exposed through `binomialPMF`.++```haskell+import Probability.Discrete++main :: IO ()+main = do+ print (binomialPMF 3 10 0.5)+```++The parameters are:++```text+k n p+```++where:++* `k` is the number of successes+* `n` is the number of trials+* `p` is the probability of success++The result is returned as an `Either` value, allowing invalid parameters to be handled explicitly.++For example:++```haskell+binomialPMF 3 10 0.5+```++returns a `Right` value containing the probability.++An invalid success count produces an error:++```haskell+binomialPMF 11 10 0.5+```++which returns a `Left` value.++## Negative Binomial and Geometric Distributions++The geometric distribution is implemented as the special case of the negative binomial distribution where:++```text+r = 1+```++For example:++```haskell+import Probability.Discrete++main :: IO ()+main = do+ print (negativeBinomialPMF 2 1 0.5)+ print (geometricPMF 2 0.5)+```++The geometric distribution is therefore provided as a convenience function rather than requiring users to manually express it as a negative binomial distribution with `r = 1`.++## Multinomial Distribution++The multinomial PMF accepts the total number of trials, a vector of category counts, and a corresponding probability vector.++```haskell+import Probability.Discrete++main :: IO ()+main = do+ print $+ multinomialPMF+ 10+ [4, 3, 3]+ [0.4, 0.3, 0.3]+```++The count vector and probability vector must be dimensionally compatible.++Invalid input is represented through the package's error type rather than silently producing an invalid probability.++## Normal Distribution++The normal distribution provides PDF and CDF functionality.++```haskell+import Probability.Continuous++main :: IO ()+main = do+ print (normalPDF 0.0 0.0 1.0)+ print (normalCDF 0.0 0.0 1.0)+```++For the standard normal distribution:++```text+μ = 0+σ = 1+```++the CDF at zero is:++```text+0.5+```++The standard deviation is validated explicitly; non-positive values result in an error.++## Exponential Distribution++The exponential distribution is parameterized by its rate `λ`.++```haskell+import Probability.Continuous++main :: IO ()+main = do+ print (exponentialPDF 1.0 2.0)+ print (exponentialCDF 1.0 2.0)+```++The implementation uses numerically appropriate calculations for expressions such as:++```text+1 - exp(-λx)+```++in order to improve numerical behavior for small values.++## Gamma Distribution++The Gamma distribution is parameterized by shape and rate parameters.++```haskell+import Probability.Continuous++main :: IO ()+main = do+ print (gammaPDF 1.0 2.0 1.0)+ print (gammaMean 2.0 1.0)+ print (gammaVar 2.0 1.0)+```++The package also provides the corresponding mean and variance functions.++The exponential distribution can be viewed as a special case of the Gamma distribution with shape parameter equal to one.++## Uniform Distribution++The continuous uniform distribution is parameterized by its lower and upper bounds.++```haskell+import Probability.Continuous++main :: IO ()+main = do+ print (uniformPDF 0.5 0.0 1.0)+ print (uniformCDF 0.5 0.0 1.0)+ print (uniformMean 0.0 1.0)+ print (uniformVar 0.0 1.0)+```++Invalid bounds are reported through the package's error handling mechanism.++## Error Handling++Public distribution functions return results using `Either`.++This makes invalid statistical parameters explicit instead of relying on exceptions or silently returning meaningless numerical values.++For example:++```haskell+case binomialPMF 11 10 0.5 of+ Right probability ->+ print probability++ Left err ->+ print err+```++The package defines distribution-related errors in its internal error module and exposes the resulting error values through the public functions.++## Numerical Design++Several calculations involved in probability distributions can become numerically unstable when performed directly.++The package therefore uses logarithmic forms where appropriate, particularly for:++* factorial-related calculations+* combinations+* Gamma functions+* Beta functions+* probability expressions involving products of many terms++Internally, the package provides mathematical support functions such as:++* `logFactorial`+* `logCombination`+* `logGamma`+* `logBeta`+* `erf`++These functions are implementation details and are kept outside the public API.++This separation allows the public distribution modules to remain focused on probability distributions while the numerical machinery remains internal to the package.++## Testing++The project includes a Cabal test suite covering:++* ordinary distribution values+* boundary conditions+* invalid parameters+* probability constraints+* dimensional consistency+* numerical results+* degenerate cases+* error handling++Run all tests with:++```bash+cabal test+```++The package is also checked with:++```bash+cabal check+```++and can be packaged with:++```bash+cabal sdist+```++## Project Structure++```text+probability-dist/+├── src/+│ └── Probability/+│ ├── Continuous.hs+│ ├── Discrete.hs+│ ├── Error.hs+│ └── Math.hs+├── test/+│ └── Main.hs+├── LICENSE+├── README.md+├── probability-dist.cabal+└── CHANGELOG.md+```++`Probability.Discrete` and `Probability.Continuous` form the public API.++`Probability.Math` and `Probability.Error` are internal implementation modules.++## License++This project is licensed under the BSD 3-Clause License.++See the `LICENSE` file for the complete license text.
+ probability-dist.cabal view
@@ -0,0 +1,47 @@+cabal-version: 3.0++name: probability-dist+version: 0.1.0.0+synopsis: Probability distributions in Haskell+description:+ A Haskell library providing discrete and continuous probability+ distributions.+category: Mathematics+license: BSD-3-Clause+license-file: LICENSE+author: BARIŞ BARIŞ+maintainer: barisbaris2005@gmail.com+build-type: Simple+extra-doc-files:+ README.md+ CHANGELOG.md++library+ hs-source-dirs: src++ exposed-modules:+ Probability.Discrete+ Probability.Continuous++ other-modules:+ Probability.Error+ Probability.Math++ build-depends:+ base >= 4.18 && < 5+ default-language: Haskell2010++test-suite probability-dist-test+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Main.hs++ build-depends:+ base,+ probability-dist++ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/barisbarisgithub/probability-dist
+ src/Probability/Continuous.hs view
@@ -0,0 +1,340 @@+module Probability.Continuous+(+ normalPDF+ ,normalLogPDF+ ,normalCDF+ ,normalMean+ ,normalVar+ ,exponentialLogPDF+ ,exponentialPDF+ ,exponentialCDF+ ,exponentialMean+ ,exponentialVar+ ,gammaLogPDF+ ,gammaPDF+ ,gammaMean+ ,gammaVar+ ,uniformLogPDF+ ,uniformPDF+ ,uniformCDF+ ,uniformMean+ ,uniformVar+)where++import Probability.Error (ProbabilityError(..))+import Probability.Math+ (logFactorial+ , logCombination+ , logGamma+ , logBeta+ ,erf+ ,expm1)+normalLogPDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+normalLogPDF x mu sigma+ | sigma <= 0 =+ Left (InvalidStandardDeviation sigma)++ | otherwise =+ let z = (x - mu) / sigma++ logP =+ - log sigma+ - 0.5 * log (2.0 * pi)+ - 0.5 * z * z++ in Right logP+++normalPDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+normalPDF x mu sigma+ | sigma <= 0 =+ Left (InvalidStandardDeviation sigma)++ | otherwise = do+ logP <- normalLogPDF x mu sigma+ pure (exp logP)+++normalCDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+normalCDF x mu sigma+ | sigma <= 0 =+ Left (InvalidStandardDeviation sigma)++ | otherwise =+ let z = (x - mu) / (sigma * sqrt 2.0)++ in Right+ (0.5 * (1.0 + erf z))+++normalMean+ :: Double+ -> Double+ -> Either ProbabilityError Double+normalMean mu sigma+ | sigma <= 0 =+ Left (InvalidStandardDeviation sigma)++ | otherwise =+ Right mu+++normalVar+ :: Double+ -> Double+ -> Either ProbabilityError Double+normalVar mu sigma+ | sigma <= 0 =+ Left (InvalidStandardDeviation sigma)++ | otherwise =+ Right (sigma * sigma)+++exponentialLogPDF+ :: Double+ -> Double+ -> Either ProbabilityError Double+exponentialLogPDF x lambda+ | lambda <= 0 =+ Left (InvalidRate lambda)++ | x < 0 =+ Right (-1.0 / 0.0)++ | otherwise =+ Right (log lambda - lambda * x)+++exponentialPDF+ :: Double+ -> Double+ -> Either ProbabilityError Double+exponentialPDF x lambda+ | lambda <= 0 =+ Left (InvalidRate lambda)++ | x < 0 =+ Right 0.0++ | otherwise = do+ logP <- exponentialLogPDF x lambda+ pure (exp logP)+++exponentialCDF+ :: Double+ -> Double+ -> Either ProbabilityError Double+exponentialCDF x lambda+ | lambda <= 0 =+ Left (InvalidRate lambda)++ | x < 0 =+ Right 0.0++ | otherwise =+ Right (- expm1 (-lambda * x))+++exponentialMean+ :: Double+ -> Either ProbabilityError Double+exponentialMean lambda+ | lambda <= 0 =+ Left (InvalidRate lambda)++ | otherwise =+ Right (1.0 / lambda)+++exponentialVar+ :: Double+ -> Either ProbabilityError Double+exponentialVar lambda+ | lambda <= 0 =+ Left (InvalidRate lambda)++ | otherwise =+ Right (1.0 / (lambda * lambda))++gammaLogPDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+gammaLogPDF x alpha lambda+ | alpha <= 0 =+ Left (InvalidShape alpha)++ | lambda <= 0 =+ Left (InvalidRate lambda)++ | x < 0 =+ Right (-1.0 / 0.0)++ | x == 0 && alpha < 1 =+ Right (1.0 / 0.0)++ | x == 0 && alpha == 1 =+ Right (log lambda)++ | x == 0 =+ Right (-1.0 / 0.0)++ | otherwise = do+ logG <- logGamma alpha++ let logP =+ alpha * log lambda+ - logG+ + (alpha - 1.0) * log x+ - lambda * x++ pure logP+++gammaPDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+gammaPDF x alpha lambda+ | alpha <= 0 =+ Left (InvalidShape alpha)++ | lambda <= 0 =+ Left (InvalidRate lambda)++ | x < 0 =+ Right 0.0++ | x == 0 && alpha < 1 =+ Right (1.0 / 0.0)++ | x == 0 && alpha == 1 =+ Right lambda++ | x == 0 =+ Right 0.0++ | otherwise = do+ logP <- gammaLogPDF x alpha lambda+ pure (exp logP)+++gammaMean+ :: Double+ -> Double+ -> Either ProbabilityError Double+gammaMean alpha lambda+ | alpha <= 0 =+ Left (InvalidShape alpha)++ | lambda <= 0 =+ Left (InvalidRate lambda)++ | otherwise =+ Right (alpha / lambda)+++gammaVar+ :: Double+ -> Double+ -> Either ProbabilityError Double+gammaVar alpha lambda+ | alpha <= 0 =+ Left (InvalidShape alpha)++ | lambda <= 0 =+ Left (InvalidRate lambda)++ | otherwise =+ Right (alpha / (lambda * lambda))++uniformLogPDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+uniformLogPDF x a b+ | b <= a =+ Left (InvalidBounds a b)++ | x < a || x > b =+ Right (-1.0 / 0.0)++ | otherwise =+ Right (- log (b - a))+++uniformPDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+uniformPDF x a b+ | b <= a =+ Left (InvalidBounds a b)++ | x < a || x > b =+ Right 0.0++ | otherwise =+ Right (1.0 / (b - a))+++uniformCDF+ :: Double+ -> Double+ -> Double+ -> Either ProbabilityError Double+uniformCDF x a b+ | b <= a =+ Left (InvalidBounds a b)++ | x < a =+ Right 0.0++ | x > b =+ Right 1.0++ | otherwise =+ Right ((x - a) / (b - a))+++uniformMean+ :: Double+ -> Double+ -> Either ProbabilityError Double+uniformMean a b+ | b <= a =+ Left (InvalidBounds a b)++ | otherwise =+ Right ((a + b) / 2.0)+++uniformVar+ :: Double+ -> Double+ -> Either ProbabilityError Double+uniformVar a b+ | b <= a =+ Left (InvalidBounds a b)++ | otherwise =+ let width = b - a+ in Right (width * width / 12.0)
+ src/Probability/Discrete.hs view
@@ -0,0 +1,420 @@+module Probability.Discrete+ ( bernoulliPMF+ , bernoulliMean+ , bernoulliVar+ ,binomialPMF+ ,binomialLogPMF+ ,binomialMean+ ,binomialVar+ ,poissonPMF+ ,poissonLogPMF+ ,poissonMean+ ,poissonVar+ ,negativeBinomialPMF+ ,negativeBinomialLogPMF+ ,negativeBinomialMean+ ,negativeBinomialVar+ ,geometricLogPMF+ ,geometricPMF+ ,geometricMean+ ,geometricVar+ ,hypergeometricPMF+ ,hypergeometricLogPMF+ ,hypergeometricMean+ ,hypergeometricVar+ ,multinomialLogPMF+ ,multinomialPMF+ ) where++import Probability.Error (ProbabilityError(..))+import Probability.Math+ (logCombination+ ,logFactorial)++bernoulliPMF :: Int -> Double -> Either ProbabilityError Double+bernoulliPMF x p+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | x == 0 = Right (1 - p)+ | x == 1 = Right p+ | otherwise = Right 0.0+++bernoulliMean :: Double -> Either ProbabilityError Double+bernoulliMean p+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | otherwise = Right p+++bernoulliVar :: Double -> Either ProbabilityError Double+bernoulliVar p+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | otherwise = Right (p * (1 - p))++++binomialLogPMF+ :: Int+ -> Int+ -> Double+ -> Either ProbabilityError Double+binomialLogPMF k n p+ | n < 0 = Left (InvalidSampleSize n)+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | k < 0 || k > n = Right 0.0+ | otherwise = do+ logC <- logCombination n k++ let termK =+ if k == 0+ then 0.0+ else fromIntegral k * log p++ termNMinusK =+ if k == n+ then 0.0+ else fromIntegral (n - k) * log (1 - p)++ return (logC + termK + termNMinusK)+++binomialPMF+ :: Int+ -> Int+ -> Double+ -> Either ProbabilityError Double+binomialPMF k n p+ | n < 0 = Left (InvalidSampleSize n)+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | k < 0 || k > n = Left(InvalidSuccessCount k)+ | otherwise = do+ logP <- binomialLogPMF k n p+ pure (exp logP)++binomialMean+ :: Int+ -> Double+ -> Either ProbabilityError Double+binomialMean n p+ | n < 0 = Left (InvalidSampleSize n)+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | otherwise = Right (fromIntegral n * p)++binomialVar+ :: Int+ -> Double+ -> Either ProbabilityError Double+binomialVar n p+ | n < 0 = Left (InvalidSampleSize n)+ | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | otherwise = Right (fromIntegral n * p * (1 - p))+++poissonLogPMF+ :: Int+ -> Double+ -> Either ProbabilityError Double+poissonLogPMF k lambda+ | lambda <= 0 = Left (InvalidRate lambda)+ | k < 0 = Right 0.0+ | otherwise = do+ logFact <- logFactorial k++ let logP =+ -lambda+ + if k == 0+ then 0.0+ else fromIntegral k * log lambda+ - logFact++ pure logP+++poissonPMF+ :: Int+ -> Double+ -> Either ProbabilityError Double+poissonPMF k lambda+ | lambda <= 0 = Left (InvalidRate lambda)+ | k < 0 = Right 0.0+ | otherwise = do+ logP <- poissonLogPMF k lambda+ pure (exp logP)+++poissonMean+ :: Double+ -> Either ProbabilityError Double+poissonMean lambda+ | lambda <= 0 = Left (InvalidRate lambda)+ | otherwise = Right lambda+++poissonVar+ :: Double+ -> Either ProbabilityError Double+poissonVar lambda+ | lambda <= 0 = Left (InvalidRate lambda)+ | otherwise = Right lambda+++negativeBinomialLogPMF+ :: Int+ -> Int+ -> Double+ -> Either ProbabilityError Double+negativeBinomialLogPMF k r p+ | r <= 0 = Left (InvalidSuccessCount r)+ | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | k < 0 = Left (InvalidFailureCount k)+ | otherwise = do+ logC <- logCombination (r + k - 1) k++ let logP =+ logC+ + fromIntegral r * log p+ + if k == 0+ then 0.0+ else fromIntegral k * log (1 - p)++ pure logP+++negativeBinomialPMF+ :: Int+ -> Int+ -> Double+ -> Either ProbabilityError Double+negativeBinomialPMF k r p+ | r <= 0 = Left (InvalidSuccessCount r)+ | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | k < 0 = Left (InvalidFailureCount k)+ | otherwise = do+ logP <- negativeBinomialLogPMF k r p+ pure (exp logP)+++negativeBinomialMean+ :: Int+ -> Double+ -> Either ProbabilityError Double+negativeBinomialMean r p+ | r <= 0 = Left (InvalidSuccessCount r)+ | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | otherwise =+ Right (fromIntegral r * (1 - p) / p)+++negativeBinomialVar+ :: Int+ -> Double+ -> Either ProbabilityError Double+negativeBinomialVar r p+ | r <= 0 = Left (InvalidSuccessCount r)+ | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)+ | otherwise =+ Right (fromIntegral r * (1 - p) / (p * p))+++geometricLogPMF+ :: Int+ -> Double+ -> Either ProbabilityError Double+geometricLogPMF k p =+ negativeBinomialLogPMF k 1 p+++geometricPMF+ :: Int+ -> Double+ -> Either ProbabilityError Double+geometricPMF k p =+ negativeBinomialPMF k 1 p+++geometricMean+ :: Double+ -> Either ProbabilityError Double+geometricMean p =+ negativeBinomialMean 1 p+++geometricVar+ :: Double+ -> Either ProbabilityError Double+geometricVar p =+ negativeBinomialVar 1 p+++hypergeometricLogPMF+ :: Int+ -> Int+ -> Int+ -> Int+ -> Either ProbabilityError Double+hypergeometricLogPMF k sampleSize successes populationSize+ | populationSize <= 0 =+ Left (InvalidPopulation populationSize) --size++ | successes < 0 || successes > populationSize =+ Left (InvalidSuccessCount successes)++ | sampleSize < 0 || sampleSize > populationSize =+ Left (InvalidSampleSize sampleSize)++ | k < 0+ || k > sampleSize+ || k > successes+ || sampleSize - k > populationSize - successes =+ Right 0.0++ | otherwise = do+ logC1 <- logCombination successes k+ logC2 <- logCombination+ (populationSize - successes)+ (sampleSize - k)+ logC3 <- logCombination populationSize sampleSize++ pure (logC1 + logC2 - logC3)+++hypergeometricPMF+ :: Int+ -> Int+ -> Int+ -> Int+ -> Either ProbabilityError Double+hypergeometricPMF k sampleSize successes populationSize+ | populationSize <= 0 =+ Left (InvalidPopulation populationSize) --size++ | successes < 0 || successes > populationSize =+ Left (InvalidSuccessCount successes)++ | sampleSize < 0 || sampleSize > populationSize =+ Left (InvalidSampleSize sampleSize)++ | k < 0+ || k > sampleSize+ || k > successes+ || sampleSize - k > populationSize - successes =+ Right 0.0++ | otherwise = do+ logP <- hypergeometricLogPMF+ k+ sampleSize+ successes+ populationSize++ pure (exp logP)+++hypergeometricMean+ :: Int+ -> Int+ -> Int+ -> Either ProbabilityError Double+hypergeometricMean sampleSize successes populationSize+ | populationSize <= 0 =+ Left (InvalidPopulation populationSize) -- size kısımını düzelttik++ | successes < 0 || successes > populationSize =+ Left (InvalidSuccessCount successes)++ | sampleSize < 0 || sampleSize > populationSize =+ Left (InvalidSampleSize sampleSize)++ | otherwise =+ Right+ ( fromIntegral sampleSize+ * fromIntegral successes+ / fromIntegral populationSize+ )+++hypergeometricVar+ :: Int+ -> Int+ -> Int+ -> Either ProbabilityError Double+hypergeometricVar sampleSize successes populationSize+ | populationSize <= 1 =+ Left (InvalidPopulation populationSize) -- düzeltme yaptık++ | successes < 0 || successes > populationSize =+ Left (InvalidSuccessCount successes)++ | sampleSize < 0 || sampleSize > populationSize =+ Left (InvalidSampleSize sampleSize)++ | otherwise =+ let n = fromIntegral sampleSize+ k = fromIntegral successes+ nPop = fromIntegral populationSize++ p = k / nPop++ in Right+ ( n+ * p+ * (1 - p)+ * ((nPop - n) / (nPop - 1))+ )++multinomialLogPMF+ :: Int+ -> [Int]+ -> [Double]+ -> Either ProbabilityError Double+multinomialLogPMF n counts probabilities+ | n < 0 =+ Left (InvalidSampleSize n)++ | null probabilities =+ Left (InvalidProbabilityVector probabilities)++ | length counts /= length probabilities =+ Left (DimensionMismatch+ (length counts)+ (length probabilities))++ | any (< 0) counts =+ Left (InvalidCountVector counts)++ | sum counts /= n =+ Left (InvalidSampleSize n)++ | any (< 0) probabilities =+ Left (InvalidProbabilityVector probabilities)++ | abs (sum probabilities - 1.0) > 1.0e-12 =+ Left (InvalidProbabilityVector probabilities)++ | otherwise = do+ logNFact <- logFactorial n++ logCountsFact <- mapM logFactorial counts++ let logCountTerm =+ logNFact - sum logCountsFact++ logProbabilityTerm =+ sum+ [ if k == 0+ then 0.0+ else fromIntegral k * log p+ | (k, p) <- zip counts probabilities+ ]++ pure (logCountTerm + logProbabilityTerm)+++multinomialPMF+ :: Int+ -> [Int]+ -> [Double]+ -> Either ProbabilityError Double+multinomialPMF n counts probabilities+ | otherwise = do+ logP <- multinomialLogPMF n counts probabilities+ pure (exp logP)
+ src/Probability/Error.hs view
@@ -0,0 +1,19 @@+module Probability.Error+(ProbabilityError(..)) where++data ProbabilityError+ = ProbabilityOutOfRange Double+ | NegativeValue Double+ | InvalidBounds Double Double+ | InvalidSampleSize Int+ | InvalidSuccessCount Int+ | InvalidFailureCount Int+ | InvalidRate Double+ | InvalidScale Double+ | InvalidShape Double+ | InvalidPopulation Int+ |InvalidProbabilityVector [Double]+ |InvalidCountVector [Int]+ |DimensionMismatch Int Int+ |InvalidStandardDeviation Double+ deriving(Show, Eq)
+ src/Probability/Math.hs view
@@ -0,0 +1,128 @@+module Probability.Math+ ( logFactorial+ , logCombination+ , logGamma+ , logBeta+ ,erf+ ,expm1+ ) where++import Probability.Error (ProbabilityError(..))+++logFactorial :: Int -> Either ProbabilityError Double+logFactorial n+ | n < 0 = Left (NegativeValue (fromIntegral n))+ | n < 35 = Right (sum [log (fromIntegral i) | i <- [1 .. n]])+ | otherwise = logGamma (fromIntegral n + 1)+++logCombination :: Int -> Int -> Either ProbabilityError Double+logCombination n r+ | n < 0 = Left (NegativeValue (fromIntegral n))+ | r < 0 = Left (NegativeValue (fromIntegral r))+ | r > n = Left (InvalidSuccessCount r)+ | otherwise = do+ lnN <- logFactorial n+ lnR <- logFactorial r+ lnNR <- logFactorial (n - r)++ pure (lnN - lnR - lnNR)+++logGamma :: Double -> Either ProbabilityError Double+logGamma z+ | z <= 0 = Left (NegativeValue z)+ | otherwise = Right result+ where+ g :: Double+ g = 7.0++ coefficients :: [Double]+ coefficients =+ [ 0.99999999999980993+ , 676.5203681218851+ , -1259.1392167224028+ , 771.32342877765313+ , -176.61502916214059+ , 12.507343278686905+ , -0.13857109526572012+ , 9.9843695780195716e-6+ , 1.5056327351493116e-7+ ]++ x :: Double+ x = z - 1.0++ a :: Double+ a =+ foldl+ (\acc (i, c) ->+ acc + c / (x + fromIntegral i))+ (head coefficients)+ (zip [1 ..] (tail coefficients))++ t :: Double+ t = x + g + 0.5++ result :: Double+ result =+ 0.5 * log (2.0 * pi)+ + (x + 0.5) * log t+ - t+ + log a+++logBeta :: Double -> Double -> Either ProbabilityError Double+logBeta a b+ | a <= 0 = Left (InvalidShape a)+ | b <= 0 = Left (InvalidShape b)+ | otherwise = do+ logA <- logGamma a+ logB <- logGamma b+ logAB <- logGamma (a + b)++ pure (logA + logB - logAB)++erf+ :: Double+ -> Double+erf x+ | x == 0.0 = 0.0+ | otherwise =+ let sign = if x < 0.0 then -1.0 else 1.0+ ax = abs x++ p = 0.3275911+ t = 1.0 / (1.0 + p * ax)++ poly =+ t * exp+ ( -ax * ax+ -1.26551223+ + t * ( 1.00002368+ + t * ( 0.37409196+ + t * ( 0.09678418+ + t * (-0.18628806+ + t * (0.27886807+ + t * (-1.13520398+ + t * (1.48851587+ + t * (-0.82215223+ + t * 0.17087277)))))))))++ in sign * (1.0 - poly)++expm1+ :: Double+ -> Double+expm1 x+ | abs x < 1.0e-5 =+ x+ + x^2 / 2.0+ + x^3 / 6.0+ + x^4 / 24.0+ + x^5 / 120.0+ + x^6 / 720.0++ | otherwise =+ exp x - 1.0
+ test/Main.hs view
@@ -0,0 +1,393 @@+module Main where++import Control.Monad (unless)+import Data.Either (isLeft)+import System.Exit (exitFailure)++--import Probability.Math+import Probability.Discrete+import Probability.Continuous+++-- ============================================================+-- Test helpers+-- ============================================================++epsilon :: Double+epsilon = 1.0e-12+++approxEqual :: Double -> Double -> Bool+approxEqual x y =+ abs (x - y) < epsilon+++assert :: Bool -> String -> IO ()+assert condition message =+ unless condition $ do+ putStrLn ("FAIL: " ++ message)+ exitFailure+++assertApprox :: Double -> Double -> String -> IO ()+assertApprox actual expected message =+ assert (approxEqual actual expected) message+++assertRightApprox+ :: Either e Double+ -> Double+ -> String+ -> IO ()+assertRightApprox result expected message =+ case result of+ Right value ->+ assertApprox value expected message+ Left _ ->+ assert False (message ++ " returned Left")+++assertLeft+ :: Either e a+ -> String+ -> IO ()+assertLeft result message =+ assert (isLeft result) message++{-+-- ============================================================+-- Math+-- ============================================================++testMath :: IO ()+testMath = do+ putStrLn "Testing Math..."++ -- logGamma+ assertApprox+ (logGamma 4)+ 1.791759469228055+ "logGamma 4"++ assertApprox+ (logGamma 5)+ 3.178053830347944+ "logGamma 5"++ assertApprox+ (logGamma 10)+ 12.801827480081474+ "logGamma 10"++ -- logFactorial+ assertApprox+ (logFactorial 0)+ 0.0+ "logFactorial 0"++ assertApprox+ (logFactorial 10)+ 15.104412573075516+ "logFactorial 10"++ assertApprox+ (logFactorial 40)+ 110.32063971475739+ "logFactorial 40"++ -- logCombination+ assertRightApprox+ (logCombination 10 0)+ 0.0+ "logCombination 10 0"++ assertRightApprox+ (logCombination 10 10)+ 0.0+ "logCombination 10 10"++ assertLeft+ (logCombination 10 11)+ "logCombination 10 11"++ assertLeft+ (logCombination (-1) 2)+ "logCombination -1 2"++ -- logBeta+ assertApprox+ (logBeta 2 2)+ (-1.791759469228055)+ "logBeta 2 2"++ assertApprox+ (logBeta 1 2)+ (-0.6931471805599472)+ "logBeta 1 2"++ -- erf+ assertApprox+ (erf 0.0)+ 0.0+ "erf 0"++ -- expm1+ assertApprox+ (expm1 0.0)+ 0.0+ "expm1 0"++ assertApprox+ (expm1 1.0)+ (exp 1.0 - 1.0)+ "expm1 1"++ putStrLn "Math: OK"++-}+-- ============================================================+-- Discrete+-- ============================================================++testDiscrete :: IO ()+testDiscrete = do+ putStrLn "Testing Discrete..."++ -- --------------------------------------------------------+ -- Binomial+ -- --------------------------------------------------------++ assertRightApprox+ (binomialPMF 0 10 0.5)+ 0.0009765625+ "binomialPMF 10 0 0.5"++ assertRightApprox+ (binomialPMF 10 10 0.5)+ 0.0009765625+ "binomialPMF 10 10 0.5"++ assertLeft+ (binomialPMF 11 10 0.5)+ "binomialPMF invalid success count"++ assertLeft+ (binomialPMF (-1) 2 0.5)+ "binomialPMF negative population"++ -- --------------------------------------------------------+ -- Geometric+ -- --------------------------------------------------------++ assertRightApprox+ (geometricPMF 1 0.5)+ 0.25+ "geometricPMF 1 0.5"++ -- --------------------------------------------------------+ -- Negative Binomial+ -- --------------------------------------------------------++ assertRightApprox+ (negativeBinomialPMF 0 1 0.5)+ 0.5+ "negativeBinomialPMF r=1"++ -- --------------------------------------------------------+ -- Multinomial+ -- --------------------------------------------------------++ assertRightApprox+ (multinomialPMF+ 10+ [10, 0]+ [1.0, 0.0])+ 1.0+ "multinomialPMF degenerate case"++ -- --------------------------------------------------------+ -- Poisson+ -- --------------------------------------------------------++ assertRightApprox+ (poissonPMF 0 2.0)+ (exp (-2.0))+ "poissonPMF 0 2"++ -- --------------------------------------------------------+ -- Bernoulli+ -- --------------------------------------------------------++ assertRightApprox+ (bernoulliPMF 1 0.7)+ 0.7+ "bernoulliPMF success"++ assertRightApprox+ (bernoulliPMF 0 0.7)+ 0.3+ "bernoulliPMF failure"++ putStrLn "Discrete: OK"+++-- ============================================================+-- Continuous+-- ============================================================++testContinuous :: IO ()+testContinuous = do+ putStrLn "Testing Continuous..."++ -- --------------------------------------------------------+ -- Normal+ -- --------------------------------------------------------++ assertRightApprox+ (normalPDF 0.0 0.0 1.0)+ 0.3989422804014327+ "normalPDF standard normal at 0"++ assertRightApprox+ (normalCDF 0.0 0.0 1.0)+ 0.5+ "normalCDF standard normal at 0"++ -- Symmetry+ let leftCDF = normalCDF (-1.0) 0.0 1.0+ rightCDF = normalCDF 1.0 0.0 1.0++ case (leftCDF, rightCDF) of+ (Right l, Right r) ->+ assertApprox+ (l + r)+ 1.0+ "normalCDF symmetry"+ _ ->+ assert False "normalCDF symmetry returned Left"++ -- Invalid standard deviation+ assertLeft+ (normalPDF 0.0 0.0 (-1.0))+ "normalPDF invalid standard deviation"++ -- --------------------------------------------------------+ -- Exponential+ -- --------------------------------------------------------++ assertRightApprox+ (exponentialPDF 0.0 2.0)+ 2.0+ "exponentialPDF x=0"++ assertRightApprox+ (exponentialCDF 0.0 2.0)+ 0.0+ "exponentialCDF x=0"++ assertRightApprox+ (exponentialCDF 1.0 2.0)+ (1.0 - exp (-2.0))+ "exponentialCDF x=1 lambda=2"++ assertLeft+ (exponentialPDF 1.0 (-1.0))+ "exponentialPDF invalid lambda"++ -- --------------------------------------------------------+ -- Gamma+ -- --------------------------------------------------------++ -- Gamma(1, lambda) = Exponential(lambda)+ assertRightApprox+ (gammaPDF 0.0 1.0 2.0)+ 2.0+ "gammaPDF alpha=1 x=0"++ assertRightApprox+ (gammaPDF 1.0 2.0 1.0)+ 0.36787944117144233+ "gammaPDF alpha=2 lambda=1 x=1"++ assertRightApprox+ (gammaMean 2.0 1.0)+ 2.0+ "gammaMean"++ assertRightApprox+ (gammaVar 2.0 1.0)+ 2.0+ "gammaVar"++ assertLeft+ (gammaPDF 1.0 (-1.0) 1.0)+ "gammaPDF invalid shape"++ assertLeft+ (gammaPDF 1.0 1.0 (-1.0))+ "gammaPDF invalid rate"++ -- --------------------------------------------------------+ -- Uniform+ -- --------------------------------------------------------++ assertRightApprox+ (uniformPDF 0.5 0.0 1.0)+ 1.0+ "uniformPDF U(0,1)"++ assertRightApprox+ (uniformCDF 0.5 0.0 1.0)+ 0.5+ "uniformCDF U(0,1)"++ assertRightApprox+ (uniformMean 0.0 1.0)+ 0.5+ "uniformMean U(0,1)"++ assertRightApprox+ (uniformVar 0.0 1.0)+ (1.0 / 12.0)+ "uniformVar U(0,1)"++ -- CDF boundaries+ assertRightApprox+ (uniformCDF (-1.0) 0.0 1.0)+ 0.0+ "uniformCDF below lower bound"++ assertRightApprox+ (uniformCDF 2.0 0.0 1.0)+ 1.0+ "uniformCDF above upper bound"++ -- Invalid bounds+ assertLeft+ (uniformPDF 0.0 1.0 1.0)+ "uniformPDF invalid bounds"++ assertLeft+ (uniformCDF 0.0 2.0 1.0)+ "uniformCDF invalid bounds"++ putStrLn "Continuous: OK"+++-- ============================================================+-- Main+-- ============================================================++main :: IO ()+main = do+ putStrLn "========================================"+ putStrLn " probability-dist test suite"+ putStrLn "========================================"++ --testMath+ testDiscrete+ testContinuous++ putStrLn "========================================"+ putStrLn " All tests passed."+ putStrLn "========================================"