diff --git a/MathNet.Numerics.sln b/MathNet.Numerics.sln index 94ee2de2..a927c196 100644 --- a/MathNet.Numerics.sln +++ b/MathNet.Numerics.sln @@ -50,7 +50,9 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Docs", "Docs", "{039229DA-A docs\content\Integration.fsx = docs\content\Integration.fsx docs\content\Interpolation.fsx = docs\content\Interpolation.fsx docs\content\LinearEquations.fsx = docs\content\LinearEquations.fsx + docs\content\MatlabFiles.fsx = docs\content\MatlabFiles.fsx docs\content\Matrix.fsx = docs\content\Matrix.fsx + docs\content\MatrixMarket.fsx = docs\content\MatrixMarket.fsx docs\content\MKL.fsx = docs\content\MKL.fsx docs\content\Probability.fsx = docs\content\Probability.fsx docs\content\Random.fsx = docs\content\Random.fsx diff --git a/src/Numerics/Settings.StyleCop b/src/Numerics/Settings.StyleCop index 34b30a8c..1bf36b60 100644 --- a/src/Numerics/Settings.StyleCop +++ b/src/Numerics/Settings.StyleCop @@ -391,5 +391,15 @@ + + + + + False + + + + + \ No newline at end of file diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs index 95b875bd..f880e504 100644 --- a/src/Numerics/Statistics/Correlation.cs +++ b/src/Numerics/Statistics/Correlation.cs @@ -69,6 +69,7 @@ namespace MathNet.Numerics.Statistics { throw new ArgumentOutOfRangeException("dataB", Resources.ArgumentArraysSameLength); } + double currentA = ieA.Current; double currentB = ieB.Current; @@ -83,8 +84,9 @@ namespace MathNet.Numerics.Statistics varA += scaleDeltaA * deltaA * (n - 1); varB += scaleDeltaB * deltaB * (n - 1); - r += ((deltaA * deltaB * (n - 1)) / n); + r += (deltaA * deltaB * (n - 1)) / n; } + if (ieB.MoveNext()) { throw new ArgumentOutOfRangeException("dataA", Resources.ArgumentArraysSameLength); @@ -103,12 +105,15 @@ namespace MathNet.Numerics.Statistics { var m = Matrix.Build.DenseIdentity(vectors.Length); for (int i = 0; i < vectors.Length; i++) + { for (int j = i + 1; j < vectors.Length; j++) { var c = Pearson(vectors[i], vectors[j]); m.At(i, j, c); m.At(j, i, c); } + } + return m; } diff --git a/src/Numerics/Statistics/DescriptiveStatistics.cs b/src/Numerics/Statistics/DescriptiveStatistics.cs index 5a6c294a..cf098eb8 100644 --- a/src/Numerics/Statistics/DescriptiveStatistics.cs +++ b/src/Numerics/Statistics/DescriptiveStatistics.cs @@ -160,9 +160,10 @@ namespace MathNet.Numerics.Statistics double variance = 0; double skewness = 0; double kurtosis = 0; - double minimum = Double.PositiveInfinity; - double maximum = Double.NegativeInfinity; + double minimum = double.PositiveInfinity; + double maximum = double.NegativeInfinity; long n = 0; + foreach (var xi in data) { double delta = xi - mean; @@ -182,11 +183,13 @@ namespace MathNet.Numerics.Statistics { minimum = xi; } + if (maximum < xi) { maximum = xi; } } + SetStatistics(mean, variance, skewness, kurtosis, minimum, maximum, n); } @@ -200,9 +203,10 @@ namespace MathNet.Numerics.Statistics double variance = 0; double skewness = 0; double kurtosis = 0; - double minimum = Double.PositiveInfinity; - double maximum = Double.NegativeInfinity; + double minimum = double.PositiveInfinity; + double maximum = double.NegativeInfinity; long n = 0; + foreach (var xi in data) { if (xi.HasValue) @@ -219,16 +223,19 @@ namespace MathNet.Numerics.Statistics skewness += tmpDelta*scaleDeltaSqr*(n - 2) - 3*scaleDelta*variance; variance += tmpDelta*scaleDelta; + if (minimum > xi) { minimum = xi.Value; } + if (maximum < xi) { maximum = xi.Value; } } } + SetStatistics(mean, variance, skewness, kurtosis, minimum, maximum, n); } @@ -242,9 +249,10 @@ namespace MathNet.Numerics.Statistics decimal variance = 0; decimal skewness = 0; decimal kurtosis = 0; - decimal minimum = Decimal.MaxValue; - decimal maximum = Decimal.MinValue; + decimal minimum = decimal.MaxValue; + decimal maximum = decimal.MinValue; long n = 0; + foreach (double x in data) { decimal xi = (decimal)x; @@ -260,15 +268,18 @@ namespace MathNet.Numerics.Statistics skewness += tmpDelta*scaleDelta2*(n - 2) - 3*scaleDelta*variance; variance += tmpDelta*scaleDelta; + if (minimum > xi) { minimum = xi; } + if (maximum < xi) { maximum = xi; } } + SetStatistics((double)mean, (double)variance, (double)skewness, (double)kurtosis, (double)minimum, (double)maximum, n); } @@ -282,9 +293,10 @@ namespace MathNet.Numerics.Statistics decimal variance = 0; decimal skewness = 0; decimal kurtosis = 0; - decimal minimum = Decimal.MaxValue; - decimal maximum = Decimal.MinValue; + decimal minimum = decimal.MaxValue; + decimal maximum = decimal.MinValue; long n = 0; + foreach (double? x in data) { if (x.HasValue) @@ -302,16 +314,19 @@ namespace MathNet.Numerics.Statistics skewness += tmpDelta*scaleDeltaSQR*(n - 2) - 3*scaleDelta*variance; variance += tmpDelta*scaleDelta; + if (minimum > xi) { minimum = xi; } + if (maximum < xi) { maximum = xi; } } } + SetStatistics((double)mean, (double)variance, (double)skewness, (double)kurtosis, (double)minimum, (double)maximum, n); } @@ -341,11 +356,13 @@ namespace MathNet.Numerics.Statistics { Minimum = minimum; Maximum = maximum; + if (n > 1) { Variance = variance/(n - 1); StandardDeviation = Math.Sqrt(Variance); } + if (Variance != 0) { if (n > 2) diff --git a/src/Numerics/Statistics/Histogram.cs b/src/Numerics/Statistics/Histogram.cs index 78dbb130..f8c956ed 100644 --- a/src/Numerics/Statistics/Histogram.cs +++ b/src/Numerics/Statistics/Histogram.cs @@ -36,7 +36,7 @@ namespace MathNet.Numerics.Statistics using Properties; /// - /// A consists of a series of s, + /// A consists of a series of s, /// each representing a region limited by a lower bound (exclusive) and an upper bound (inclusive). /// [Serializable] @@ -45,7 +45,7 @@ namespace MathNet.Numerics.Statistics IComparable #else IComparable, ICloneable -#endif +#endif { /// /// This IComparer performs comparisons between a point and a bucket. @@ -159,7 +159,7 @@ namespace MathNet.Numerics.Statistics throw new ArgumentException(Resources.PartialOrderException); } - if (UpperBound.AlmostEqual(bucket.UpperBound) + if (UpperBound.AlmostEqual(bucket.UpperBound) && LowerBound.AlmostEqual(bucket.LowerBound)) { return 0; @@ -175,7 +175,7 @@ namespace MathNet.Numerics.Statistics /// /// Checks whether two Buckets are equal; this method tolerates a difference in lowerbound, upperbound - /// and count given by . + /// and count given by . /// public override bool Equals(object obj) { @@ -206,7 +206,7 @@ namespace MathNet.Numerics.Statistics return "(" + LowerBound + ";" + UpperBound + "] = " + Count; } } - + /// /// A class which computes histograms of data. /// @@ -359,7 +359,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Returns the Bucket that contains the value v. + /// Returns the Bucket that contains the value v. /// /// The point to search the bucket for. /// A copy of the bucket containing point . @@ -427,7 +427,7 @@ namespace MathNet.Numerics.Statistics return (Bucket) _buckets[n].Clone(); } } - + /// /// Gets the number of buckets. /// @@ -444,7 +444,7 @@ namespace MathNet.Numerics.Statistics get { double totalCount = 0; - + for (int i = 0; i < BucketCount; i++) { totalCount += this[i].Count; @@ -464,6 +464,7 @@ namespace MathNet.Numerics.Statistics { sb.Append(b); } + return sb.ToString(); } } diff --git a/src/Numerics/Statistics/RunningStatistics.cs b/src/Numerics/Statistics/RunningStatistics.cs index 1458e52b..2a2e8836 100644 --- a/src/Numerics/Statistics/RunningStatistics.cs +++ b/src/Numerics/Statistics/RunningStatistics.cs @@ -47,8 +47,8 @@ namespace MathNet.Numerics.Statistics double _m2; double _m3; double _m4; - double _min = Double.PositiveInfinity; - double _max = Double.NegativeInfinity; + double _min = double.PositiveInfinity; + double _max = double.NegativeInfinity; public RunningStatistics() { @@ -194,6 +194,7 @@ namespace MathNet.Numerics.Statistics { _min = value; } + if (_max < value) { _max = value; diff --git a/src/Numerics/Statistics/SortedArrayStatistics.cs b/src/Numerics/Statistics/SortedArrayStatistics.cs index 07ae7b26..c29b17e2 100644 --- a/src/Numerics/Statistics/SortedArrayStatistics.cs +++ b/src/Numerics/Statistics/SortedArrayStatistics.cs @@ -46,7 +46,10 @@ namespace MathNet.Numerics.Statistics /// Sample array, must be sorted ascendingly. public static double Minimum(double[] data) { - if (data.Length == 0) return double.NaN; + if (data.Length == 0) + { + return double.NaN; + } return data[0]; } @@ -57,7 +60,10 @@ namespace MathNet.Numerics.Statistics /// Sample array, must be sorted ascendingly. public static double Maximum(double[] data) { - if (data.Length == 0) return double.NaN; + if (data.Length == 0) + { + return double.NaN; + } return data[data.Length - 1]; } @@ -69,7 +75,10 @@ namespace MathNet.Numerics.Statistics /// One-based order of the statistic, must be between 1 and N (inclusive). public static double OrderStatistic(double[] data, int order) { - if (order < 1 || order > data.Length) return double.NaN; + if (order < 1 || order > data.Length) + { + return double.NaN; + } return data[order - 1]; } @@ -81,7 +90,10 @@ namespace MathNet.Numerics.Statistics /// Sample array, must be sorted ascendingly. public static double Median(double[] data) { - if (data.Length == 0) return double.NaN; + if (data.Length == 0) + { + return double.NaN; + } var k = data.Length/2; return data.Length.IsOdd() @@ -138,7 +150,11 @@ namespace MathNet.Numerics.Statistics /// Sample array, must be sorted ascendingly. public static double[] FiveNumberSummary(double[] data) { - if (data.Length == 0) return new[] { double.NaN, double.NaN, double.NaN, double.NaN, double.NaN }; + if (data.Length == 0) + { + return new[] { double.NaN, double.NaN, double.NaN, double.NaN, double.NaN }; + } + return new[] { data[0], Quantile(data, 0.25), Quantile(data, 0.50), Quantile(data, 0.75), data[data.Length - 1] }; } @@ -157,9 +173,20 @@ namespace MathNet.Numerics.Statistics /// public static double Quantile(double[] data, double tau) { - if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; - if (tau == 0d || data.Length == 1) return data[0]; - if (tau == 1d) return data[data.Length - 1]; + if (tau < 0d || tau > 1d || data.Length == 0) + { + return double.NaN; + } + + if (tau == 0d || data.Length == 1) + { + return data[0]; + } + + if (tau == 1d) + { + return data[data.Length - 1]; + } double h = (data.Length + 1/3d)*tau + 1/3d; var hf = (int)h; @@ -182,7 +209,10 @@ namespace MathNet.Numerics.Statistics /// d-parameter public static double QuantileCustom(double[] data, double tau, double a, double b, double c, double d) { - if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; + if (tau < 0d || tau > 1d || data.Length == 0) + { + return double.NaN; + } var x = a + (data.Length + b)*tau - 1; #if PORTABLE @@ -213,9 +243,20 @@ namespace MathNet.Numerics.Statistics /// Quantile definition, to choose what product/definition it should be consistent with public static double QuantileCustom(double[] data, double tau, QuantileDefinition definition) { - if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; - if (tau == 0d || data.Length == 1) return data[0]; - if (tau == 1d) return data[data.Length - 1]; + if (tau < 0d || tau > 1d || data.Length == 0) + { + return double.NaN; + } + + if (tau == 0d || data.Length == 1) + { + return data[0]; + } + + if (tau == 1d) + { + return data[data.Length - 1]; + } switch (definition) { @@ -224,16 +265,19 @@ namespace MathNet.Numerics.Statistics double h = data.Length*tau + 0.5d; return data[(int)Math.Ceiling(h - 0.5d) - 1]; } + case QuantileDefinition.R2: { double h = data.Length*tau + 0.5d; return (data[(int)Math.Ceiling(h - 0.5d) - 1] + data[(int)(h + 0.5d) - 1])*0.5d; } + case QuantileDefinition.R3: { double h = data.Length*tau; return data[Math.Max((int)Math.Round(h) - 1, 0)]; } + case QuantileDefinition.R4: { double h = data.Length*tau; @@ -242,6 +286,7 @@ namespace MathNet.Numerics.Statistics var upper = data[Math.Min(hf, data.Length - 1)]; return lower + (h - hf)*(upper - lower); } + case QuantileDefinition.R5: { double h = data.Length*tau + 0.5d; @@ -250,6 +295,7 @@ namespace MathNet.Numerics.Statistics var upper = data[Math.Min(hf, data.Length - 1)]; return lower + (h - hf)*(upper - lower); } + case QuantileDefinition.R6: { double h = (data.Length + 1)*tau; @@ -258,6 +304,7 @@ namespace MathNet.Numerics.Statistics var upper = data[Math.Min(hf, data.Length - 1)]; return lower + (h - hf)*(upper - lower); } + case QuantileDefinition.R7: { double h = (data.Length - 1)*tau + 1d; @@ -266,6 +313,7 @@ namespace MathNet.Numerics.Statistics var upper = data[Math.Min(hf, data.Length - 1)]; return lower + (h - hf)*(upper - lower); } + case QuantileDefinition.R8: { double h = (data.Length + 1/3d)*tau + 1/3d; @@ -274,6 +322,7 @@ namespace MathNet.Numerics.Statistics var upper = data[Math.Min(hf, data.Length - 1)]; return lower + (h - hf)*(upper - lower); } + case QuantileDefinition.R9: { double h = (data.Length + 0.25d)*tau + 0.375d; @@ -282,6 +331,7 @@ namespace MathNet.Numerics.Statistics var upper = data[Math.Min(hf, data.Length - 1)]; return lower + (h - hf)*(upper - lower); } + default: throw new NotSupportedException(); } @@ -294,8 +344,15 @@ namespace MathNet.Numerics.Statistics /// The value where to estimate the CDF at. public static double EmpiricalCDF(double[] data, double x) { - if (x < data[0]) return 0.0; - if (x >= data[data.Length - 1]) return 1.0; + if (x < data[0]) + { + return 0.0; + } + + if (x >= data[data.Length - 1]) + { + return 1.0; + } int right = Array.BinarySearch(data, x); if (right >= 0) @@ -304,10 +361,11 @@ namespace MathNet.Numerics.Statistics { right++; } - return (right + 1)/(double)(data.Length); + + return (right + 1)/(double)data.Length; } - return (~right)/(double)(data.Length); + return (~right)/(double)data.Length; } /// @@ -321,17 +379,26 @@ namespace MathNet.Numerics.Statistics /// Rank definition, to choose how ties should be handled and what product/definition it should be consistent with public static double QuantileRank(double[] data, double x, RankDefinition definition = RankDefinition.Default) { - if (x < data[0]) return 0.0; - if (x >= data[data.Length - 1]) return 1.0; + if (x < data[0]) + { + return 0.0; + } + + if (x >= data[data.Length - 1]) + { + return 1.0; + } int right = Array.BinarySearch(data, x); if (right >= 0) { int left = right; + while (left > 0 && data[left - 1] == data[left]) { left--; } + while (right < data.Length - 1 && data[right + 1] == data[right]) { right++; @@ -340,13 +407,17 @@ namespace MathNet.Numerics.Statistics switch (definition) { case RankDefinition.EmpiricalCDF: - return (right + 1)/(double)(data.Length); + return (right + 1)/(double)data.Length; + case RankDefinition.Max: return right/(double)(data.Length - 1); + case RankDefinition.Min: return left/(double)(data.Length - 1); + case RankDefinition.Average: return (left/(double)(data.Length - 1) + right/(double)(data.Length - 1))/2; + default: throw new NotSupportedException(); } @@ -359,7 +430,8 @@ namespace MathNet.Numerics.Statistics switch (definition) { case RankDefinition.EmpiricalCDF: - return (left + 1)/(double)(data.Length); + return (left + 1)/(double)data.Length; + default: { var a = left/(double)(data.Length - 1); @@ -385,6 +457,7 @@ namespace MathNet.Numerics.Statistics { ranks[i] = i + 1; } + return ranks; } @@ -424,16 +497,19 @@ namespace MathNet.Numerics.Statistics rank = (b + a - 1)/2d + 1; break; } + case RankDefinition.Min: { rank = a + 1; break; } + case RankDefinition.Max: { rank = b; break; } + default: throw new NotSupportedException(); } diff --git a/src/Numerics/Statistics/Statistics.cs b/src/Numerics/Statistics/Statistics.cs index ce4122da..ee2a1393 100644 --- a/src/Numerics/Statistics/Statistics.cs +++ b/src/Numerics/Statistics/Statistics.cs @@ -721,7 +721,6 @@ namespace MathNet.Numerics.Statistics return order => SortedArrayStatistics.OrderStatistic(array, order); } - /// /// Evaluates the rank of each entry of the provided samples. /// The rank definition can be specificed to be compatible @@ -747,7 +746,6 @@ namespace MathNet.Numerics.Statistics return Ranks(data.Where(d => d.HasValue).Select(d => d.Value), definition); } - /// /// Estimates the quantile tau from the provided samples. /// The tau-th quantile is the data value where the cumulative distribution @@ -806,7 +804,6 @@ namespace MathNet.Numerics.Statistics return QuantileRankFunc(data.Where(d => d.HasValue).Select(d => d.Value), definition); } - /// /// Estimates the empirical cummulative distribution function (CDF) at x from the provided samples. /// @@ -849,7 +846,6 @@ namespace MathNet.Numerics.Statistics return EmpiricalCDFFunc(data.Where(d => d.HasValue).Select(d => d.Value)); } - /// /// Estimates the empirical inverse CDF at tau from the provided samples. /// diff --git a/src/Numerics/Statistics/StreamingStatistics.cs b/src/Numerics/Statistics/StreamingStatistics.cs index 0c09e65c..4c593740 100644 --- a/src/Numerics/Statistics/StreamingStatistics.cs +++ b/src/Numerics/Statistics/StreamingStatistics.cs @@ -52,14 +52,17 @@ namespace MathNet.Numerics.Statistics { var min = double.PositiveInfinity; bool any = false; + foreach (var d in stream) { if (d < min || double.IsNaN(d)) { min = d; } + any = true; } + return any ? min : double.NaN; } @@ -72,14 +75,17 @@ namespace MathNet.Numerics.Statistics { var max = double.NegativeInfinity; bool any = false; + foreach (var d in stream) { if (d > max || double.IsNaN(d)) { max = d; } + any = true; } + return any ? max : double.NaN; } @@ -93,11 +99,13 @@ namespace MathNet.Numerics.Statistics double mean = 0; ulong m = 0; bool any = false; + foreach (var d in stream) { mean += (d - mean)/++m; any = true; } + return any ? mean : double.NaN; } @@ -112,6 +120,7 @@ namespace MathNet.Numerics.Statistics double variance = 0; double sum = 0; ulong count = 0; + using (var iterator = samples.GetEnumerator()) { if (iterator.MoveNext()) @@ -129,6 +138,7 @@ namespace MathNet.Numerics.Statistics variance += (diff*diff)/(count*(count - 1)); } } + return count > 1 ? variance/(count - 1) : double.NaN; } @@ -143,6 +153,7 @@ namespace MathNet.Numerics.Statistics double variance = 0; double sum = 0; ulong count = 0; + using (var iterator = population.GetEnumerator()) { if (iterator.MoveNext()) @@ -160,6 +171,7 @@ namespace MathNet.Numerics.Statistics variance += (diff*diff)/(count*(count - 1)); } } + return variance/count; } @@ -197,6 +209,7 @@ namespace MathNet.Numerics.Statistics double variance = 0; double sum = 0; ulong count = 0; + using (var iterator = samples.GetEnumerator()) { if (iterator.MoveNext()) @@ -215,6 +228,7 @@ namespace MathNet.Numerics.Statistics mean += (xi - mean) / count; } } + return new Tuple( count > 0 ? mean : double.NaN, count > 1 ? variance/(count - 1) : double.NaN); @@ -246,6 +260,7 @@ namespace MathNet.Numerics.Statistics var mean1 = 0.0; var mean2 = 0.0; var comoment = 0.0; + using (var s1 = samples1.GetEnumerator()) using (var s2 = samples2.GetEnumerator()) { @@ -261,13 +276,14 @@ namespace MathNet.Numerics.Statistics mean1 += (s1.Current - mean1)/n; mean2 += (s2.Current - mean2)/n; comoment += (s1.Current - mean1)*(s2.Current - mean2Prev); - } + if (s2.MoveNext()) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } } + return n > 1 ? comoment/(n - 1) : double.NaN; } @@ -285,6 +301,7 @@ namespace MathNet.Numerics.Statistics var mean1 = 0.0; var mean2 = 0.0; var comoment = 0.0; + using (var p1 = population1.GetEnumerator()) using (var p2 = population2.GetEnumerator()) { @@ -300,18 +317,19 @@ namespace MathNet.Numerics.Statistics mean1 += (p1.Current - mean1) / n; mean2 += (p2.Current - mean2) / n; comoment += (p1.Current - mean1) * (p2.Current - mean2Prev); - } + if (p2.MoveNext()) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } } + return comoment/n; } /// - /// Calculates the entropy of a stream of double values. + /// Calculates the entropy of a stream of double values. /// Returns NaN if any of the values in the stream are NaN. /// /// The input stream to evaluate. @@ -326,20 +344,26 @@ namespace MathNet.Numerics.Statistics int totalCount = 0; foreach (double value in stream) { - if (double.IsNaN(value)) return double.NaN; + if (double.IsNaN(value)) + { + return double.NaN; + } - double currentValueCount = 0; - + double currentValueCount; if (index.TryGetValue(value, out currentValueCount)) + { index[value] = ++currentValueCount; + } else + { index.Add(value, 1); - + } + ++totalCount; } // calculate the entropy of the stream - double entropy = 0; + double entropy = 0; foreach (var item in index) { double p = item.Value / totalCount;