diff --git a/src/FSharp/DenseVector.fs b/src/FSharp/DenseVector.fs index 76bcd75f..84235dba 100644 --- a/src/FSharp/DenseVector.fs +++ b/src/FSharp/DenseVector.fs @@ -45,4 +45,23 @@ module DenseVector = let n = List.length fl let v = Double.DenseVector(n) fl |> List.iteri (fun i f -> v.[i] <- f) - v \ No newline at end of file + v + + /// Create a vector from a sequences. + let inline of_seq (fs: #seq) = + let n = Seq.length fs + let v = DenseVector(n) + fs |> Seq.iteri (fun i f -> v.[i] <- f) + v + + /// Create a vector with evenly spaced entries: e.g. rangef -1.0 0.5 1.0 = [-1.0 -0.5 0.0 0.5 1.0] + let inline rangef (start: float) (step: float) (stop: float) = + let n = (int ((stop - start) / step)) + 1 + let v = new DenseVector(n) + for i=0 to n-1 do + v.[i] <- (float i) * step + start + v + + /// Create a vector with integer entries in the given range. + let inline range (start: int) (stop: int) = + new DenseVector([| for i in [start .. stop] -> float i |]) \ No newline at end of file diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj index 69c3e408..5a7b7f4e 100644 --- a/src/FSharp/FSharp.fsproj +++ b/src/FSharp/FSharp.fsproj @@ -46,6 +46,7 @@ + diff --git a/src/FSharp/Vector.fs b/src/FSharp/Vector.fs new file mode 100644 index 00000000..c7e7934a --- /dev/null +++ b/src/FSharp/Vector.fs @@ -0,0 +1,169 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.LinearAlgebra.Double + +open MathNet.Numerics.LinearAlgebra + +/// A module which implements functional vector operations. +module Vector = + + /// Transform a vector into an array. + let inline to_array (v: #Vector) = + let n = v.Count + Array.init n (fun i -> v.Item(i)) + + /// Transform a vector into an array. + let inline to_list (v: #Vector) = + let n = v.Count + List.init n (fun i -> v.Item(i)) + + /// In-place mutation by applying a function to every element of the vector. + let inline mapInPlace (f: float -> float) (v: #Vector) = + for i=0 to v.Count-1 do + v.Item(i) <- f (v.Item(i)) + () + + /// In-place mutation by applying a function to every element of the vector. + let inline mapiInPlace (f: int -> float -> float) (v: #Vector) = + for i=0 to v.Count-1 do + v.Item(i) <- f i (v.Item(i)) + () + + /// In-place vector addition. + let inline addInPlace (v: #Vector) (w: #Vector) = + v.Add w + + /// In place vector subtraction. + let inline subInPlace (v: #Vector) (w: #Vector) = + v.Subtract w + + /// Functional map operator for vectors. + /// + let inline map f (v: #Vector) = + let w = v.Clone() + inplace_mapi (fun _ x -> f x) w + w + + /// Applies a function to all elements of the vector. + let inline iter (f: float -> unit) (v: #Vector) = + for i=0 to v.Count-1 do + f (v.Item i) + + /// Applies a function to all elements of the vector. + let inline iteri (f: int -> float -> unit) (v: #Vector) = + for i=0 to v.Count-1 do + f i (v.Item i) + + /// Maps a vector to a new vector by applying a function to every element. + let inline mapi (f: int -> float -> float) (v: #Vector) = + let w = v.Clone() + inplace_mapi f w + w + + /// Fold all entries of a vector. + let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (v: #Vector) = + let mutable acc = acc0 + for i=0 to v.Count-1 do + acc <- f acc (v.Item(i)) + acc + + /// Fold all entries of a vector using a position dependent folding function. + let inline foldi (f: int -> 'a -> float -> 'a) (acc0: 'a) (v: #Vector) = + let mutable acc = acc0 + for i=0 to v.Count-1 do + acc <- f i acc (v.Item(i)) + acc + + /// Checks whether a predicate is satisfied for every element in the vector. + let inline forall (p: float -> bool) (v: #Vector) = + let mutable b = true + let mutable i = 0 + while b && i < v.Count do + b <- b && (p (v.Item(i))) + i <- i+1 + b + + /// Checks whether there is an entry in the vector that satisfies a given predicate. + let inline exists (p: float -> bool) (v: #Vector) = + let mutable b = false + let mutable i = 0 + while not(b) && i < v.Count do + b <- b || (p (v.Item(i))) + i <- i+1 + b + + /// Checks whether a predicate is true for all entries in a vector. + let inline foralli (p: int -> float -> bool) (v: #Vector) = + let mutable b = true + let mutable i = 0 + while b && i < v.Count do + b <- b && (p i (v.Item(i))) + i <- i+1 + b + + /// Checks whether there is an entry in the vector that satisfies a given position dependent predicate. + let inline existsi (p: int -> float -> bool) (v: #Vector) = + let mutable b = false + let mutable i = 0 + while not(b) && i < v.Count do + b <- b || (p i (v.Item(i))) + i <- i+1 + b + + /// Scans a vector; like fold but returns the intermediate result. + let inline scan (f: float -> float -> float) (v: #Vector) = + let w = v.Clone() + let mutable p = v.Item(0) + for i=1 to v.Count-1 do + p <- f p (v.Item(i)) + w.[i] <- p + w + + /// Scans a vector; like fold but returns the intermediate result. + let inline scanBack (f: float -> float -> float) (v: #Vector) = + let w = v.Clone() + let mutable p = v.Item(v.Count-1) + for i=2 to v.Count do + p <- f (v.Item(v.Count - i)) p + w.[v.Count - i] <- p + w + + /// Reduces a vector: the result of this function will be f(...f(f(v[0],v[1]), v[2]),..., v[n]). + let inline reduce (f: float -> float -> float) (v: #Vector) = + let mutable p = v.Item(0) + for i=1 to v.Count-1 do + p <- f p (v.Item(i)) + p + + /// Reduces a vector: the result of this function will be f(v[1], ..., f(v[n-2], f(v[n-1],v[n]))...). + let inline reduceBack (f: float -> float -> float) (v: #Vector) = + let mutable p = v.Item(v.Count-1) + for i=2 to v.Count do + p <- f (v.Item(v.Count - i)) p + p \ No newline at end of file diff --git a/src/Numerics/Distributions/Discrete/Bernoulli.cs b/src/Numerics/Distributions/Discrete/Bernoulli.cs index 2188d008..cc492e8b 100644 --- a/src/Numerics/Distributions/Discrete/Bernoulli.cs +++ b/src/Numerics/Distributions/Discrete/Bernoulli.cs @@ -24,4 +24,320 @@ // 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. -// \ No newline at end of file +// + +namespace MathNet.Numerics.Distributions +{ + using System; + using System.Collections.Generic; + using Properties; + + /// + /// The Bernoulli distribution is a distribution over bits. The parameter + /// p specifies the probability that a 1 is generated. + /// + /// The distribution will use the by default. + /// Users can set the random number generator by using the property. + /// The statistics classes will check all the incoming parameters whether they are in the allowed + /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters + /// to false, all parameter checks can be turned off. + public class Bernoulli : IDiscreteDistribution + { + /// + /// The probability of generating a one. + /// + private double _p; + + /// + /// The distribution's random number generator. + /// + private Random _random; + + /// + /// Construct a new Bernoulli distribution. + /// + /// The probability of generating one. + /// If the Bernoulli parameter is not in the range [0,1]. + public Bernoulli(double p) + { + SetParameters(p); + RandomSource = new System.Random(); + } + + /// + /// A string representation of the distribution. + /// + public override string ToString() + { + return "Bernoulli(P = " + _p + ")"; + } + + /// + /// Checks whether the parameters of the distribution are valid. + /// + /// The probability of generating a one. + /// True when the parameters are valid, false otherwise. + private static bool IsValidParameterSet(double p) + { + if (p >= 0.0 && p <= 1.0) + { + return true; + } + + return false; + } + + /// + /// Sets the parameters of the distribution after checking their validity. + /// + /// The probability of generating a one. + /// When the parameters don't pass the function. + private void SetParameters(double p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + _p = p; + } + + /// + /// Gets or sets the probability of generating a one. + /// + public double P + { + get + { + return _p; + } + + set + { + SetParameters(value); + } + } + + #region IDistribution Members + + /// + /// Gets or sets the random number generator which is used to draw random samples. + /// + public Random RandomSource + { + get + { + return _random; + } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + + _random = value; + } + } + + /// + /// Gets the mean of the distribution. + /// + public double Mean + { + get { return _p; } + } + + /// + /// Gets the standard deviation of the distribution. + /// + public double StdDev + { + get { return Math.Sqrt(_p * (1.0 - _p)); } + } + + /// + /// Gets the variance of the distribution. + /// + public double Variance + { + get { return _p * (1.0 - _p); } + } + + /// + /// Gets the entropy of the distribution. + /// + public double Entropy + { + get { return -_p * Math.Log(_p) - (1.0 - _p) * Math.Log(1.0 - _p); } + } + + /// + /// Gets the skewness of the distribution. + /// + public double Skewness + { + get { return (1.0 - 2.0 * _p) / Math.Sqrt(_p * (1.0 - _p)); } + } + + /// + /// Gets the smallest element in the domain of the distributions which can be represented by an integer. + /// + public int Minimum { get { return 0; } } + + /// + /// Gets the largest element in the domain of the distributions which can be represented by an integer. + /// + public int Maximum { get { return 1; } } + + /// + /// Computes the cumulative distribution function of the Bernoulli distribution. + /// + /// The location at which to compute the cumulative density. + /// the cumulative density at . + public double CumulativeDistribution(double x) + { + if (x < 0) + { + return 0.0; + } + if (x == 0) + { + return 1.0 - _p; + } + + return 1.0; + } + + #endregion + + #region IDiscreteDistribution Members + + /// + /// The mode of the distribution. + /// + public int Mode + { + get { return _p > 0.5 ? 1 : 0; } + } + + /// + /// The median of the distribution. + /// + public int Median + { + get { throw new Exception("The median of the Bernoulli distribution is undefined."); } + } + + /// + /// Computes the probability of a specific value. + /// + public double Probability(int val) + { + if (val == 0) + { + return 1.0 - _p; + } + + if (val == 1) + { + return _p; + } + + return 0.0; + } + + /// + /// Computes the probability of a specific value. + /// + public double ProbabilityLn(int val) + { + if (val == 0) + { + return Math.Log(1.0 - _p); + } + + if (val == 1) + { + return Math.Log(_p); + } + + return Double.NegativeInfinity; + } + + /// + /// Samples a Bernoulli distributed random variable. + /// + /// A sample from the Bernoulli distribution. + public int Sample() + { + return DoSample(RandomSource, _p); + } + + /// + /// Samples an array of Bernoulli distributed random variables. + /// + /// a sequence of samples from the distribution. + public IEnumerable Samples() + { + while (true) + { + yield return DoSample(RandomSource, _p); + } + } + + #endregion + + /// + /// Samples a Bernoulli distributed random variable. + /// + /// The random number generator to use. + /// The probability of generating a 1. + /// A sample from the Bernoulli distribution. + public static int Sample(System.Random rnd, double p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + return DoSample(rnd, p); + } + + /// + /// Samples an array of Bernoulli distributed random variables. + /// + /// The random number generator to use. + /// The probability of generating a 1. + /// a sequence of samples from the distribution. + public static IEnumerable Samples(System.Random rnd, double p) + { + if (Control.CheckDistributionParameters && !IsValidParameterSet(p)) + { + throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters); + } + + while (true) + { + yield return DoSample(rnd, p); + } + } + + /// + /// Generates one sample from the Bernoulli distribution. + /// + /// The random source to use. + /// The probability of generating a one. + /// A random sample from the Bernoulli distribution. + private static int DoSample(System.Random rnd, double p) + { + if (rnd.NextDouble() < p) + { + return 1; + } + + return 0; + } + } +} \ No newline at end of file diff --git a/src/UnitTests/DistributionTests/CommonDistributionTests.cs b/src/UnitTests/DistributionTests/CommonDistributionTests.cs index ed3e3240..3d5ab6ff 100644 --- a/src/UnitTests/DistributionTests/CommonDistributionTests.cs +++ b/src/UnitTests/DistributionTests/CommonDistributionTests.cs @@ -41,12 +41,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [SetUp] public void SetupDistributions() { - dists = new IDistribution[4]; + dists = new IDistribution[5]; dists[0] = new Beta(1.0, 1.0); dists[1] = new ContinuousUniform(0.0, 1.0); dists[2] = new Gamma(1.0, 1.0); dists[3] = new Normal(0.0, 1.0); + dists[4] = new Bernoulli(0.6); } [Test] @@ -54,6 +55,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(1)] [Row(2)] [Row(3)] + [Row(4)] public void ValidateThatUnivariateDistributionsHaveRandomSource(int i) { Assert.IsNotNull(dists[i].RandomSource); @@ -64,6 +66,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(1)] [Row(2)] [Row(3)] + [Row(4)] public void CanSetRandomSource(int i) { dists[i].RandomSource = new Random(); @@ -74,6 +77,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests [Row(1)] [Row(2)] [Row(3)] + [Row(4)] [ExpectedException(typeof(ArgumentNullException))] public void FailSetRandomSourceWithNullReference(int i) { diff --git a/src/UnitTests/DistributionTests/Discrete/BernoulliTests.cs b/src/UnitTests/DistributionTests/Discrete/BernoulliTests.cs new file mode 100644 index 00000000..6d1316fe --- /dev/null +++ b/src/UnitTests/DistributionTests/Discrete/BernoulliTests.cs @@ -0,0 +1,252 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.UnitTests.DistributionTests +{ + using System; + using System.Linq; + using MbUnit.Framework; + using MathNet.Numerics.Distributions; + + [TestFixture] + public class BernoulliTests + { + [SetUp] + public void SetUp() + { + Control.CheckDistributionParameters = true; + } + + [Test] + [Row(0.0)] + [Row(0.3)] + [Row(1.0)] + public void CanCreateBernoulli(double p) + { + var bernoulli = new Bernoulli(p); + AssertEx.AreEqual(p, bernoulli.P); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + [Row(Double.NaN)] + [Row(-1.0)] + [Row(2.0)] + public void NormalCreateFailsWithBadParameters(double p) + { + var bernoulli = new Bernoulli(p); + } + + [Test] + public void ValidateToString() + { + var b = new Bernoulli(0.3); + AssertEx.AreEqual("Bernoulli(P = 0.3)", n.ToString()); + } + + [Test] + [Row(0.0)] + [Row(0.3)] + [Row(1.0)] + public void CanSetProbabilityOfOne(double p) + { + var b = new Bernoulli(0.3); + b.P = p; + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + [Row(Double.NaN)] + [Row(-1.0)] + [Row(2.0)] + public void SetProbabilityOfOneFails(double p) + { + var b = new Bernoulli(0.3); + b.P = p; + } + + [Test] + [Row(0.0)] + [Row(0.3)] + [Row(1.0)] + public void ValidateEntropy(double p) + { + var b = new Bernoulli(p); + AssertEx.AreEqual((1.0 - p) * Math.Log(1.0 - p) + p * Math.Log(p), b.Entropy); + } + + [Test] + [Row(0.0)] + [Row(0.3)] + [Row(1.0)] + public void ValidateSkewness(double p) + { + var b = new Bernoulli(p); + AssertEx.AreEqual((1.0 - 2.0 * p) / Math.Sqrt(p * (1.0 - p)), n.Skewness); + } + + [Test] + [Row(0.0, 0)] + [Row(0.3, 0)] + [Row(1.0, 1)] + public void ValidateMode(double p, double m) + { + var b = new Bernoulli(p); + AssertEx.AreEqual(mean, n.Mode); + } + + [Test] + [ExpectedException(typeof(Exception))] + public void ValidateMedian() + { + var b = new Bernoulli(0.3); + } + + [Test] + public void ValidateMinimum() + { + var b = new Bernoulli(0.3); + AssertEx.AreEqual(0.0, n.Minimum); + } + + [Test] + public void ValidateMaximum() + { + var b = new Bernoulli(0.3); + AssertEx.AreEqual(1.0, n.Maximum); + } + + [Test] + [Row(0.0, -1.0, 0.0)] + [Row(0.0, 0.0, 1.0)] + [Row(0.0, 0.5, 0.0)] + [Row(0.0, 1.0, 0.0)] + [Row(0.0, 2.0, 0.0)] + [Row(0.3, -1.0, 0.0)] + [Row(0.3, 0.0, 0.7)] + [Row(0.3, 0.5, 0.0)] + [Row(0.3, 1.0, 0.3)] + [Row(0.3, 2.0, 0.0)] + [Row(1.0, -1.0, 0.0)] + [Row(1.0, 0.0, 0.0)] + [Row(1.0, 0.5, 0.0)] + [Row(1.0, 1.0, 1.0)] + [Row(1.0, 2.0, 0.0)] + public void ValidateProbability(double p, double x, double d) + { + var b = new Bernoulli(p); + AssertEx.AreEqual(d, b.Probability(x)); + } + + [Test] + [Row(0.0, -1.0, Double.NegativeInfinity)] + [Row(0.0, 0.0, 0.0)] + [Row(0.0, 0.5, Double.NegativeInfinity)] + [Row(0.0, 1.0, Double.NegativeInfinity)] + [Row(0.0, 2.0, Double.NegativeInfinity)] + [Row(0.3, -1.0, Double.NegativeInfinity)] + [Row(0.3, 0.0, -0.35667494393873244235395440410727451457180907089949815)] + [Row(0.3, 0.5, Double.NegativeInfinity)] + [Row(0.3, 1.0, -1.2039728043259360296301803719337238685164245381839102)] + [Row(0.3, 2.0, Double.NegativeInfinity)] + [Row(1.0, -1.0, Double.NegativeInfinity)] + [Row(1.0, 0.0, Double.NegativeInfinity)] + [Row(1.0, 0.5, Double.NegativeInfinity)] + [Row(1.0, 1.0, 0.0)] + [Row(1.0, 2.0, Double.NegativeInfinity)] + public void ValidateProbabilityLn(double p, double x, double dln) + { + var b = new Bernoulli(p); + AssertEx.AreEqual(dln, b.ProbabilityLn(x)); + } + + [Test] + public void CanSampleStatic() + { + var d = Bernoulli.Sample(new Random(), 0.3); + } + + [Test] + public void CanSampleSequenceStatic() + { + var ied = Bernoulli.Samples(new Random(), 0.3); + var arr = ied.Take(5).ToArray(); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleStatic() + { + var d = Bernoulli.Sample(new Random(), -1.0); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void FailSampleSequenceStatic() + { + var ied = Bernoulli.Samples(new Random(), -1.0).First(); + } + + [Test] + public void CanSample() + { + var n = new Bernoulli(); + var d = n.Sample(); + } + + [Test] + public void CanSampleSequence() + { + var n = new Bernoulli(); + var ied = n.Samples(); + var e = ied.Take(5).ToArray(); + } + + [Test] + [Row(0.0, -1.0, 0.0)] + [Row(0.0, 0.0, 1.0)] + [Row(0.0, 0.5, 1.0)] + [Row(0.0, 1.0, 1.0)] + [Row(0.0, 2.0, 1.0)] + [Row(0.3, -1.0, 0.0)] + [Row(0.3, 0.0, 0.7)] + [Row(0.3, 0.5, 0.7)] + [Row(0.3, 1.0, 1.0)] + [Row(0.3, 2.0, 1.0)] + [Row(1.0, -1.0, 0.0)] + [Row(1.0, 0.0, 0.0)] + [Row(1.0, 0.5, 0.0)] + [Row(1.0, 1.0, 1.0)] + [Row(1.0, 2.0, 1.0)] + public void ValidateCumulativeDistribution(double p, double x, double cdf) + { + var b = new Bernoulli(p); + AssertEx.AreEqual(cdf, n.CumulativeDistribution(x)); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 79126209..11092c1f 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -68,6 +68,7 @@ + @@ -109,9 +110,6 @@ MathNet.Numerics.snk - - -