diff --git a/.gitignore b/.gitignore index 4b4562ee..d1b5bba8 100644 --- a/.gitignore +++ b/.gitignore @@ -77,3 +77,4 @@ docs/content/Contributing.md docs/content/Contributors.md docs/content/ReleaseNotes.md docs/content/ReleaseNotes-*.md +/build.fsx.lock diff --git a/src/Numerics/Distributions/BetaBinomial.cs b/src/Numerics/Distributions/BetaBinomial.cs new file mode 100644 index 00000000..13d450fc --- /dev/null +++ b/src/Numerics/Distributions/BetaBinomial.cs @@ -0,0 +1,416 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2014 Math.NET +// +// 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. +// + +// +// Andrew J. Willshire +// + + +using System; +using System.Collections.Generic; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Random; + +namespace MathNet.Numerics.Distributions +{ + /// + /// Discrete Univariate Beta-Binomial distribution. + /// The beta-binomial distribution is a family of discrete probability distributions on a finite support of non-negative integers arising + /// when the probability of success in each of a fixed or known number of Bernoulli trials is either unknown or random. + /// The beta-binomial distribution is the binomial distribution in which the probability of success at each of n trials is not fixed but randomly drawn from a beta distribution. + /// It is frequently used in Bayesian statistics, empirical Bayes methods and classical statistics to capture overdispersion in binomial type distributed data. + /// Wikipedia - Beta-Binomial distribution. + /// + public class BetaBinomial : IDiscreteDistribution + { + System.Random _random; + + readonly int _n; + readonly double _a; + readonly double _b; + + /// + /// Initializes a new instance of the class. + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + public BetaBinomial(int n, double a, double b) + { + if (!IsValidParameterSet(n, a, b)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + _random = SystemRandomSource.Default; + _n = n; + _a = a; + _b = b; + } + + /// + /// Initializes a new instance of the class. + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + /// The random number generator which is used to draw random samples. + public BetaBinomial(int n, double a, double b, System.Random randomSource) + { + if (!IsValidParameterSet(n,a,b)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + _random = randomSource ?? SystemRandomSource.Default; + _n = n; + _a = a; + _b = b; + } + + /// + /// Returns a that represents this instance. + /// + /// + /// A that represents this instance. + /// + public override string ToString() + { + return $"BetaBinomial(n = {_n}, a = {_a}, b = {_b})"; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + public static bool IsValidParameterSet(int n, double a, double b) + { + return n >= 1.0 && a > 0.0 && b > 0.0; + } + + /// + /// Tests whether the provided values are valid parameters for this distribution. + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + /// The location in the domain where we want to evaluate the probability mass function. + public static bool IsValidParameterSet(int n, double a, double b, int k) + { + return n >= 1.0 && a > 0.0 && b > 0.0 && k >=0 && k <=n; + } + + + public int N => _n; + public double A => _a; + public double B => _b; + + public System.Random RandomSource + { + get => _random; + set => _random = value ?? SystemRandomSource.Default; + } + /// + /// Gets the mean of the distribution. + /// + double IUnivariateDistribution.Mean => (_n * _a) / (_a + _b); + + /// + /// Gets the variance of the distribution. + /// + double IUnivariateDistribution.Variance => (_n*_a*_b*(_a+_b+_n))/(Math.Pow((_a+_b),2) * (_a+_b+1)); + + /// + /// Gets the standard deviation of the distribution. + /// + double IUnivariateDistribution.StdDev => Math.Sqrt((_n * _a * _b * (_a + _b + _n)) / (Math.Pow((_a + _b), 2) * (_a + _b + 1))); + + /// + /// Gets the entropy of the distribution. + /// + double IUnivariateDistribution.Entropy => throw new NotSupportedException(); + + /// + /// Gets the skewness of the distribution. + /// + double IUnivariateDistribution.Skewness => + (_a + _b + 2 * _n) * (_b - _a) / (_a + _b + 2) * Math.Sqrt((1 + _a + _b) / (_n * _a * _b * (_n + _a + _b))); + + /// + /// Gets the mode of the distribution + /// + int IDiscreteDistribution.Mode => throw new NotSupportedException(); + + /// + /// Gets the median of the distribution. + /// + double IUnivariateDistribution.Median => throw new NotSupportedException(); + + /// + /// Gets the smallest element in the domain of the distributions which can be represented by an integer. + /// + public int Minimum => 0; + + /// + /// Gets the largest element in the domain of the distributions which can be represented by an integer. + /// + public int Maximum => int.MaxValue; + + /// + /// Computes the probability mass (PMF) at k, i.e. P(X = k). + /// + /// The location in the domain where we want to evaluate the probability mass function. + /// the probability mass at location . + public double Probability(int k) + { + return PMF(_n, _a, _b, k); + } + + /// + /// Computes the log probability mass (lnPMF) at k, i.e. ln(P(X = k)). + /// + /// The location in the domain where we want to evaluate the log probability mass function. + /// the log probability mass at location . + public double ProbabilityLn(int k) + { + return PMFLn(_n, _a, _b, k); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location + public double CumulativeDistribution(double x) + { + return CDF(_n, _a, _b, (int)Math.Floor(x)); + } + + /// + /// Computes the probability mass (PMF) at k, i.e. P(X = k). + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + /// The location in the domain where we want to evaluate the probability mass function. + /// the probability mass at location . + public static double PMF(int n, double a, double b, int k) + { + if (!IsValidParameterSet(n, a, b, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + if (k > n) + { + return 0.0; + } + else + { + return Math.Exp(PMFLn(n, a, b, k)); + } + } + + /// + /// Computes the log probability mass (lnPMF) at k, i.e. ln(P(X = k)). + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + /// The location in the domain where we want to evaluate the probability mass function. + /// the log probability mass at location . + public static double PMFLn(int n, double a, double b, int k) + { + if (!IsValidParameterSet(n, a, b, k)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + return SpecialFunctions.BinomialLn((n), k) + + SpecialFunctions.BetaLn(k + a, n - k + b) + - SpecialFunctions.BetaLn(a, b); + } + + /// + /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x). + /// + /// The number of Bernoulli trials n - n is a positive integer + /// Shape parameter alpha of the Beta distribution. Range: a > 0. + /// Shape parameter beta of the Beta distribution. Range: b > 0. + /// The location at which to compute the cumulative distribution function. + /// the cumulative distribution at location . + /// + public static double CDF(int n, double a, double b, int x) + { + if (!IsValidParameterSet(n,a,b,x)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + double accumulator = 0; + + for (int i = 0; i<=x; i++) + { + accumulator += PMF(n, a, b, i); + } + + return accumulator; + } + + /// + /// Samples BetaBinomial distributed random variables by sampling a Beta distribution then passing to a Binomial distribution. + /// + /// The random number generator to use. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The number of trials (n). Range: n ≥ 0. + /// a random number from the BetaBinomial distribution. + static int SampleUnchecked(System.Random rnd, int n, double a, double b) + { + var p = Beta.SampleUnchecked(rnd, a, b); + var x = Binomial.SampleUnchecked(rnd, p, n); + return x; + } + + static void SamplesUnchecked(System.Random rnd, int[] values, int n, double a, double b) + { + for (int i = 0; i < values.Length; i++) + { + values[i] = SampleUnchecked(rnd, n, a, b); + } + } + + static IEnumerable SamplesUnchecked(System.Random rnd, int n, double a, double b) + { + while (true) + { + yield return SampleUnchecked(rnd, n, a, b); + } + } + + + /// + /// Samples a BetaBinomial distributed random variable. + /// + /// a sample from the distribution. + + public int Sample() + { + return SampleUnchecked(_random, _n, _a, _b); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + public void Samples(int[] values) + { + SamplesUnchecked(_random, values, _n, _a, _b); + } + + /// + /// Samples an array of BetaBinomial distributed random variables. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + return SamplesUnchecked(_random, _n, _a, _b); + } + + /// + /// Samples a BetaBinomial distributed random variable. + /// + /// The random number generator to use. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The number of trials (n). Range: n ≥ 0. + /// a sample from the distribution. + + public int Sample(System.Random rnd, int n, double a, double b) + { + if (!IsValidParameterSet(n,a,b)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SampleUnchecked(rnd, n, a, b); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The random number generator to use. + /// The array to fill with the samples. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The number of trials (n). Range: n ≥ 0. + public void Samples(System.Random rnd, int[] values, int n, double a, double b) + { + if (!IsValidParameterSet(n, a, b)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + SamplesUnchecked(rnd, values, n, a, b); + } + + /// + /// Samples an array of BetaBinomial distributed random variables. + /// + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The number of trials (n). Range: n ≥ 0. + /// a sequence of samples from the distribution. + public IEnumerable Samples(int n, double a, double b) + { + if (!IsValidParameterSet(n, a, b)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + return SamplesUnchecked(_random, n, a, b); + } + + /// + /// Fills an array with samples generated from the distribution. + /// + /// The array to fill with the samples. + /// The α shape parameter of the Beta distribution. Range: α ≥ 0. + /// The β shape parameter of the Beta distribution. Range: β ≥ 0. + /// The number of trials (n). Range: n ≥ 0. + public void Samples(int[] values, int n, double a, double b) + { + if (!IsValidParameterSet(n, a, b)) + { + throw new ArgumentException(Resources.InvalidDistributionParameters); + } + + SamplesUnchecked(_random, values, n, a, b); + } + } +} diff --git a/src/Numerics/Distributions/Binomial.cs b/src/Numerics/Distributions/Binomial.cs index 4a58f163..f5516089 100644 --- a/src/Numerics/Distributions/Binomial.cs +++ b/src/Numerics/Distributions/Binomial.cs @@ -359,7 +359,7 @@ namespace MathNet.Numerics.Distributions /// The success probability (p) in each trial. Range: 0 ≤ p ≤ 1. /// The number of trials (n). Range: n ≥ 0. /// The number of successful trials. - static int SampleUnchecked(System.Random rnd, double p, int n) + internal static int SampleUnchecked(System.Random rnd, double p, int n) { var k = 0; for (var i = 0; i < n; i++) diff --git a/src/Numerics/SpecialFunctions/GeneralizedHyperGeometric.cs b/src/Numerics/SpecialFunctions/GeneralizedHyperGeometric.cs new file mode 100644 index 00000000..5a7c9e9d --- /dev/null +++ b/src/Numerics/SpecialFunctions/GeneralizedHyperGeometric.cs @@ -0,0 +1,145 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2010 Math.NET +// +// 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. +// + +// +// Andrew J. Willshire +// + +using System; +using System.Linq; + +namespace MathNet.Numerics +{ + public static partial class SpecialFunctions + { + + //Rising and falling factorials - reference here: + //https://en.wikipedia.org/wiki/Falling_and_rising_factorials + + + /// + /// Computes the Rising Factorial (Pochhammer function) x -> (x)n, n>= 0. see: https://en.wikipedia.org/wiki/Falling_and_rising_factorials + /// + /// The real value of the Rising Factorial for x and n + + public static double RisingFactorial(double x, int n) + { + double accumulator = 1.0; + + for (int k = 0; k < n; k++) + { + accumulator *= (x + k); + } + return accumulator; + } + + /// + /// Computes the Falling Factorial (Pochhammer function) x -> x(n), n>= 0. see: https://en.wikipedia.org/wiki/Falling_and_rising_factorials + /// + /// The real value of the Falling Factorial for x and n + public static double FallingFactorial(double x, int n) + { + double accumulator = 1.0; + + for (int k = 0; k < n; k++) + { + accumulator *= (x - k); + } + return accumulator; + } + // + + /// + /// A generalized hypergeometric series is a power series in which the ratio of successive coefficients indexed by n is a rational function of n. + /// This is the most common pFq(a1, ..., ap; b1,...,bq; z) representation + /// see: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function + /// + /// The list of coefficients in the numerator + /// The list of coefficients in the denominator + /// The variable in the power series + /// The value of the Generalized HyperGeometric Function. + public static double GeneralizedHypergeometric(double[] a, double[] b, int z) + { + const double epsilon = 0.000000000000001; + + double cumulatives = 0.0; + double currentIncrement; + int n = 0; + + do + { + currentIncrement = HGIncrement(a, b, z, n); + cumulatives += currentIncrement; + n += 1; + } + while (Math.Abs(currentIncrement) > epsilon && Math.Abs(currentIncrement) > 0 && currentIncrement.IsFinite()); + + return cumulatives; + } + + //Calculate each iteration of the function + private static double HGIncrement(double[] a, double[] b, int z, int currentN) + { + double incrementAs = 1.0; + double incrementBs = 1.0; + + double[] incrementAArray = new double[a.Length]; + double[] incrementBArray = new double[b.Length]; + + for (int p = 0; p < a.Length; p++) + { + incrementAs *= RisingFactorial(a[p], currentN); + incrementAArray[p] = RisingFactorial(a[p], currentN); + } + + for (int q = 0; q < b.Length; q++) + { + incrementBs *= RisingFactorial(b[q], currentN); + incrementBArray[q] = RisingFactorial(b[q], currentN); + } + + double numZeros = (from x in incrementAArray where x == 0 select x).Count(); + double numPoles = (from x in incrementBArray where x == 0 select x).Count(); + + if (numZeros > 0 && numZeros >= numPoles) + { + return 0.0; + } + else if (numPoles > 0 && numPoles > numZeros) + { + return double.PositiveInfinity; + } + else + { + return incrementAs / incrementBs * Math.Pow(z, currentN) / Factorial(currentN); + } + } + + } +}