diff --git a/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs b/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs index 6354391..e420497 100644 --- a/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs +++ b/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs @@ -141,20 +141,22 @@ private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyLi double value = f.Evaluate(point); double delta = model.MinimumValue - value; - double tol = settings.ComputePrecision(value); + double tol = settings.ComputePrecision(Math.Min(model.MinimumValue, value)); + // Note value can be way off, so use better of old best and new value to compute tol. + // When we didn't do this before, we got value = infinity, so tol = infinity, and thus terminated! - if (delta > 0 && settings.Listener != null) { + if (delta > 0.0 && settings.Listener != null) { MultiExtremum report = new MultiExtremum(f.EvaluationCount, settings, point, value, Math.Max(Math.Abs(delta), 0.75 * tol), model.GetHessian()); settings.Listener(report); } // To terminate, we demand: a reduction, that the reduction be small, that the reduction be in line with - // its expected value, that we have not run up against trust boundary, and that the gradient is small. + // its expected value, that we have not run up against the trust boundary, and that the gradient is small. // I had wanted to demand delta > 0, but we run into some cases where delta keeps being very slightly // negative, typically orders of magnitude less than tol, causing the trust radius to shrink in an // endless cycle that causes our approximation to ultimately go sour, even though terminating on the original // very slightly negative delta would have produced an accurate estimate. So we tolerate this case for now. - if ((-tol / 4.0 <= delta) && (delta <= tol)) { + if ((delta <= tol) && (-0.25 * tol <= delta)) { // We demand that the model be decent, i.e. that the expected delta was within tol of the measured delta. if (Math.Abs(delta - deltaExpected) <= tol) { // We demand that the step not just be small because it ran up against the trust radius. @@ -181,7 +183,7 @@ private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyLi // If the actual change was very far from the expected change, reduce the trust radius. // If the expected change did a good job of predicting the actual change, increase the trust radius. if ((delta < 0.25 * deltaExpected) /*|| (8.0 * deltaExpected < delta)*/) { - trustRadius = trustRadius / 2.0; + trustRadius = 0.5 * trustRadius; } else if ((0.75 * deltaExpected <= delta) /*&& (delta <= 2.0 * deltaExpected)*/) { trustRadius = 2.0 * trustRadius; } @@ -195,8 +197,9 @@ private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyLi double bad = model.ComputeBadness(i, z, point, value); if (bad > fBad) { iBad = i; fBad = bad; } } + // Use the new point as long as it is better than our worst existing point. if (value < fMax) { - //Debug.WriteLine("iMax={0}, iBad={1}", iMax, iBad); + Debug.Assert(!Double.IsPositiveInfinity(value) && !Double.IsNaN(value)); model.ReplacePoint(iBad, point, z, value); } // There is some question about how best to choose which point to replace. @@ -710,21 +713,28 @@ public double[] ConvertPoint (double[] z) { return (x); } + // Badness is used as a criteria for which point to throw out of collection. + // Basic idea is to compute a number which (i) increases with value above + // minimum and also (ii) increases with distance from minimum. + public double ComputeBadness (int index, double[] z, double[] point, double value) { double s = 0.0; if (value < MinimumValue) { + // If new candidate is best, compute badness relative to it. for (int i = 0; i < point.Length; i++) { s += MoreMath.Sqr(points[index][i] - point[i]); } s = Math.Pow(s, 3.0 / 2.0) * (values[index] - value); } else { + // Otherwise compute badness relative to best existing point. if (index == minValueIndex) return (0.0); for (int i = 0; i < point.Length; i++) { s += MoreMath.Sqr(points[index][i] - origin[i]); } s = Math.Pow(s, 3.0 / 2.0) * (values[index] - values[minValueIndex]); } - return (s); + Debug.Assert(s >= 0.0); + return s; } public double FindMinimumSeperation () { diff --git a/Numerics/Analysis/NamespaceDoc.cs b/Numerics/Analysis/NamespaceDoc.cs index d3d5c38..e2f85a5 100644 --- a/Numerics/Analysis/NamespaceDoc.cs +++ b/Numerics/Analysis/NamespaceDoc.cs @@ -3,7 +3,7 @@ namespace Meta.Numerics.Analysis { /// - /// Contains types used to analyze user-supplied functions. + /// Contains types used to solve, maximize, integrate, and otherwise perform analysis on user-supplied functions. /// /// /// This namespace contains types for the numerical analysis of user-supplied functions. Examples diff --git a/Numerics/Core/DiscreteInterval.cs b/Numerics/Core/DiscreteInterval.cs index 844b774..55b5326 100644 --- a/Numerics/Core/DiscreteInterval.cs +++ b/Numerics/Core/DiscreteInterval.cs @@ -1,8 +1,6 @@ using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Threading.Tasks; +using System.Diagnostics; +using System.Globalization; namespace Meta.Numerics { @@ -12,6 +10,7 @@ namespace Meta.Numerics { public struct DiscreteInterval : IEquatable { internal DiscreteInterval (int min, int max) { + Debug.Assert(max >= min); this.min = min; this.max = max; } @@ -63,31 +62,68 @@ internal int Width { // Equality - public bool Equals (DiscreteInterval other) { - return ((this.min == other.min) && (this.max == other.max)); + private static bool Equals (DiscreteInterval u, DiscreteInterval v) { + return (u.min == v.min) && (u.max == v.max); } + + /// + /// Tests whether two discrete intervals are equal. + /// + /// The first discrete interval. + /// The second discrete interval. + /// if and are equal, otherwise . public static bool operator == (DiscreteInterval a, DiscreteInterval b) { - return (a.Equals(b)); + return Equals(a, b); } + /// + /// Tests whether two intervals are not equal. + /// + /// The first interval. + /// The second interval. + /// if and are equal, otherwise . public static bool operator != (DiscreteInterval a, DiscreteInterval b) { - return (!a.Equals(b)); + return !Equals(a, b); } + /// + /// Tests wither the current instance is equal to another discrete interval. + /// + /// Another discrete interval. + /// if the current instance equals , otherwise . + public bool Equals (DiscreteInterval other) { + return Equals(this, other); + } + + /// + /// Tests whether a given object is equal to the current interval. + /// + /// An object. + /// if is an equal , otherwise . public override bool Equals (object obj) { if (obj is DiscreteInterval other) { - return (this.Equals(other)); + return Equals(this, other); } else { - return (false); + return false; } } + /// public override int GetHashCode () { unchecked { return (min.GetHashCode() + 31 * max.GetHashCode()); } } + /// + public override string ToString () { + return ToString(CultureInfo.CurrentCulture); + } + + private string ToString (IFormatProvider format) { + return String.Format(format, "{{{0}..{1}}}", min, max); + } + } } diff --git a/Numerics/Core/Interval.cs b/Numerics/Core/Interval.cs index 8f0303d..5c11232 100644 --- a/Numerics/Core/Interval.cs +++ b/Numerics/Core/Interval.cs @@ -13,7 +13,7 @@ namespace Meta.Numerics { /// public struct Interval : IEquatable { - private double a, b, w; + private readonly double a, b, w; private Interval (double a, double b, double w) { this.a = a; @@ -120,55 +120,67 @@ public static Interval FromMidpointAndWidth (double midpoint, double width) { // equality + private static bool Equals (Interval u, Interval v) { + return ((u.a == v.a) && (u.w == v.w)); + } + + /// + /// Tests wither the current instance is equal to another interval. + /// + /// Another iteterval. + /// if the current instance equals , otherwise . public bool Equals (Interval other) { - return ((this.a == other.a) && (this.w == other.w)); + return Equals(this, other); } /// - /// Determines whether two intervals are equal. + /// Tests whether two intervals are equal. /// - /// The first interval. - /// The second interval. - /// True if and are equal, otherwise false. - public static bool operator == (Interval i1, Interval i2) { - return (i1.Equals(i2)); + /// The first interval. + /// The second interval. + /// if and are equal, otherwise . + public static bool operator == (Interval u, Interval v) { + return Equals(u, v); } /// - /// Determines whether two intervals are not equal. + /// Tests whether two intervals are not equal. /// - /// The first interval. - /// The second interval. - /// True if and are not equal, otherwise false. - public static bool operator != (Interval i1, Interval i2) { - return (!i1.Equals(i2)); + /// The first interval. + /// The second interval. + /// if and are equal, otherwise . + public static bool operator != (Interval u, Interval v) { + return !Equals(u, v); } /// - /// Determines whether a given object is an equal interval. + /// Tests whether a given object is equal to the current interval. /// /// An object. - /// True if is an equal , otherwise false. + /// if is an equal , otherwise . public override bool Equals (object obj) { if (obj is Interval other) { - return(this.Equals(other)); + return Equals(this, other); } else { - return (false); + return false; } } /// public override int GetHashCode () { - return (a.GetHashCode() ^ w.GetHashCode()); + return a.GetHashCode() ^ w.GetHashCode(); } - // text representation /// /// Produces a string representation of the interval. /// /// A string representation of the interval. public override string ToString () { - return (String.Format(CultureInfo.CurrentCulture, "[{0},{1}]", a, b)); + return ToString(CultureInfo.CurrentCulture); + } + + private string ToString (IFormatProvider format) { + return String.Format(CultureInfo.CurrentCulture, "[{0},{1}]", a, b); } #if SHO diff --git a/Numerics/Core/Polynomial.cs b/Numerics/Core/Polynomial.cs index 733735c..02b814f 100644 --- a/Numerics/Core/Polynomial.cs +++ b/Numerics/Core/Polynomial.cs @@ -17,7 +17,7 @@ public class Polynomial { /// The specified polynomial. /// /// Coefficients should be arranged from low to high order, so that the kth entry is the coefficient of xk. For example, - /// to specify the polynomial 5 - 6 x + 7 x2, give the values 5, -6, 7. + /// to specify the polynomial 5 - 6 x + 7 x3, give the values 5, -6, , 0, 7. /// public static Polynomial FromCoefficients (params double[] coefficients) { if (coefficients == null) throw new ArgumentNullException(nameof(coefficients)); diff --git a/Numerics/Data/FrameColumn.cs b/Numerics/Data/FrameColumn.cs index 2e8dbb7..471e054 100644 --- a/Numerics/Data/FrameColumn.cs +++ b/Numerics/Data/FrameColumn.cs @@ -17,6 +17,8 @@ namespace Meta.Numerics.Data public sealed class FrameColumn : IEnumerable { internal FrameColumn(FrameView frame, int c) { + Debug.Assert(frame != null); + Debug.Assert(c >= 0 && c < frame.Columns.Count); this.column = frame.columns[c]; this.map = frame.map; } @@ -77,8 +79,7 @@ public object this[int r] { /// public IReadOnlyList As () { // If the requested column is of the requested type, expose it directly. - IReadOnlyList typedColumn = column as IReadOnlyList; - if (typedColumn != null) { + if (column is IReadOnlyList typedColumn) { return (new TypedFrameColumn(column.Name, typedColumn, map)); } else { // If the requested column is cast-able to the requested type, cast it. diff --git a/Numerics/Data/FrameColumnCollection.cs b/Numerics/Data/FrameColumnCollection.cs index 828c7a2..63b4414 100644 --- a/Numerics/Data/FrameColumnCollection.cs +++ b/Numerics/Data/FrameColumnCollection.cs @@ -33,10 +33,11 @@ public int Count { /// /// The index of the column. /// The column with index . + /// is less than zero or greater than or equal to the number of rows . public FrameColumn this [int index] { get { - if ((index < 0) || (index >= frame.columns.Count)) throw new ArgumentOutOfRangeException(nameof(index)); - return (new FrameColumn(frame, index)); + if ((index < 0) || (index >= frame.columns.Count)) throw new IndexOutOfRangeException(); + return new FrameColumn(frame, index); } } @@ -45,10 +46,12 @@ public FrameColumn this [int index] { /// /// The name of the column. /// The column with the given name. + /// There is no column with the name in the collection. public FrameColumn this [string name] { get { int index = frame.GetColumnIndex(name); - return (new FrameColumn(frame, index)); + if (index < 0) throw new IndexOutOfRangeException(); + return new FrameColumn(frame, index); } } diff --git a/Numerics/Data/FrameRow.cs b/Numerics/Data/FrameRow.cs index df742c1..0f22a1c 100644 --- a/Numerics/Data/FrameRow.cs +++ b/Numerics/Data/FrameRow.cs @@ -52,9 +52,10 @@ IEnumerable IReadOnlyDictionary.Values { /// /// The index of the column. /// The value of the th column. + /// is less than zero or greater than or equal to the number of columns . public object this[int c] { get { - if ((c < 0) || (c >= Count)) throw new ArgumentOutOfRangeException(nameof(c)); + if ((c < 0) || (c >= Count)) throw new IndexOutOfRangeException(); return (frame.columns[c].GetItem(frame.map[r])); } } @@ -64,6 +65,7 @@ public object this[int c] { /// /// The name of the column. /// The value. + /// There is no column in the row with the given . public object this[string name] { get { int c = frame.GetColumnIndex(name); diff --git a/Numerics/Extended/AdvancedDoubleDoubleMath_Gamma.cs b/Numerics/Extended/AdvancedDoubleDoubleMath_Gamma.cs index d2cf2f5..01182b1 100644 --- a/Numerics/Extended/AdvancedDoubleDoubleMath_Gamma.cs +++ b/Numerics/Extended/AdvancedDoubleDoubleMath_Gamma.cs @@ -132,7 +132,7 @@ private static DoubleDouble LogGamma_ShiftToZetaSeries (DoubleDouble x) { -((DoubleDouble) 74611) / 330, ((DoubleDouble) 854513) / 138, -((DoubleDouble) 236364091) / 2730, ((DoubleDouble) 8553103) / 6, -((DoubleDouble) 23749461029) / 870 }; - public static DoubleDouble BernoulliSum (DoubleDouble x) { + private static DoubleDouble BernoulliSum (DoubleDouble x) { DoubleDouble rxPower = 1.0 / x; DoubleDouble rxSquared = rxPower * rxPower; DoubleDouble f = 0.5 * Bernoulli[1] * rxPower; diff --git a/Numerics/Extended/Int128.cs b/Numerics/Extended/Int128.cs index e3add1c..49f19e8 100644 --- a/Numerics/Extended/Int128.cs +++ b/Numerics/Extended/Int128.cs @@ -52,10 +52,19 @@ private Int128 (UInt128 u) { /// public static readonly Int128 Zero = new Int128(0L, 0UL); + /// + /// The unit value of the signed 128-bit integer. + /// public static readonly Int128 One = new Int128(0L, 1UL); + /// + /// The maximum value of the signed 128-bit integer. + /// public static readonly Int128 MaxValue = new Int128(ulong.MaxValue >> 1, ulong.MaxValue); + /// + /// The minimum value of the signed 128-bit integer. + /// public static readonly Int128 MinValue = new Int128(1UL << 63, 0UL); // Equality @@ -103,7 +112,7 @@ public bool Equals (Int128 other) { /// Tests whether the current instance equals another object. /// /// The other object. - /// if the current instance equals , otherwise. + /// if the current instance equals , otherwise. public override bool Equals (object obj) { if (obj is Int128) { return Equals(this, (Int128) obj); @@ -112,6 +121,10 @@ public override bool Equals (object obj) { } } + /// + /// Return a hash code for the 128-bit integer. + /// + /// A hash code for the value. public override int GetHashCode () { return u.GetHashCode() + 31; } @@ -178,15 +191,20 @@ public static int Compare (Int128 x, Int128 y) { /// /// Compares the current instance to another 128-bit integer. /// - /// The other integer. - /// 1 if the current instance is greater, -1 if is greater, - /// and 0 if the current instance is equal to . + /// The other integer. + /// 1 if the current instance is greater, -1 if is greater, + /// and 0 if the current instance is equal to . public int CompareTo (Int128 x) { return Compare(this, x); } // Implicit conversions + /// + /// Converts a 64-bit integer into a 128-bit integer. + /// + /// The 64-bit integer to convert. + /// THe corresponding 128-bit integer. public static implicit operator Int128 (long s) { if (s < 0L) { return -(new Int128(0UL, (ulong) (-s))); @@ -196,6 +214,11 @@ public static implicit operator Int128 (long s) { } + /// + /// Converts a 128-bit integer into a floating point value. + /// + /// The 128-bit integer to convert. + /// The correspoding floating point value. public static implicit operator double (Int128 s) { if (s.u.IsNegative) { return -((double) s.u.Negate()); @@ -204,6 +227,11 @@ public static implicit operator double (Int128 s) { } } + /// + /// Converts a 128-bit integer into an arbitrary-size big integer. + /// + /// The 128-bit integer to convert. + /// The corresponding big integer. public static implicit operator BigInteger (Int128 s) { byte[] bytes = s.u.GetBytesForBigInteger(true); return new BigInteger(bytes); @@ -211,10 +239,20 @@ public static implicit operator BigInteger (Int128 s) { // Explicit conversions + /// + /// Converts a 128-bit integer into a 64-bit integer. + /// + /// The 128-bit integer to convert. + /// The 64-bit integer with the same lower 64 bits. public static explicit operator long (Int128 s) { return unchecked((long) (ulong) s.u); } + /// + /// Converts an aribitrary-sized big integer into a 128-bit integer. + /// + /// Te big integer to be converted. + /// The 128-bit integer with the same lower 128 bits. public static explicit operator Int128 (BigInteger b) { UInt128 u = (UInt128) b; return new Int128(u); @@ -222,6 +260,10 @@ public static explicit operator Int128 (BigInteger b) { // Serialization and deserialization + /// + /// Produces a string representation of the 128-bit integer. + /// + /// A string containing the base-10 representation of the 128-bit integer value. public override string ToString () { if (u.IsNegative) { UInt128 t = u.Negate(); @@ -231,6 +273,14 @@ public override string ToString () { } } + /// + /// Produces a 128-bit integer from its string representation. + /// + /// A string containing the base-10 representation of a 128-bit integer value. + /// The 128-bit integer represented by . + /// was . + /// did not contain a valid representation of an integer. + /// represented an integer outside the range of . public static Int128 Parse (string text) { ParseResult result = UInt128.TryParse(text, true, out UInt128 u); if (result == ParseResult.Null) throw new ArgumentNullException(nameof(text)); @@ -239,6 +289,12 @@ public static Int128 Parse (string text) { return new Int128(u); } + /// + /// Attempts to produce a 128-bit integer from a string representation. + /// + /// A string containing the base-10 representation of a 128-bit integer value. + /// If the parse was successful, the 128-bit integer represented by . + /// if the parse was successful, otherwise . public static bool TryParse (string text, out Int128 value) { ParseResult result = UInt128.TryParse(text, true, out UInt128 u); if (result == ParseResult.Success) { @@ -368,8 +424,9 @@ public static Int128 DivRem (Int128 x, Int128 y, out Int128 r) { /// The argument. /// The absolute value of the argument. /// Because is one unit smaller in absolute value than , this - /// method throws a when passed . Built in types such as + /// method throws an when passed . Built in types such as /// have analogous behavior. All other values are supported. + /// was . public static Int128 Abs (Int128 x) { UInt128 xu = x.u; if (xu.IsNegative) { diff --git a/Numerics/Extended/NamespaceDoc.cs b/Numerics/Extended/NamespaceDoc.cs index 1a2a41a..345da7e 100644 --- a/Numerics/Extended/NamespaceDoc.cs +++ b/Numerics/Extended/NamespaceDoc.cs @@ -5,7 +5,7 @@ namespace Meta.Numerics.Extended { /// - /// Contains types that enable extended-precision floating point calculations. + /// Contains types that enable extended-precision integer and floating point calculations. /// /// /// The structure supports floating-point arithmetic with diff --git a/Numerics/Extended/UInt128.cs b/Numerics/Extended/UInt128.cs index 2ed704e..f0c4b41 100644 --- a/Numerics/Extended/UInt128.cs +++ b/Numerics/Extended/UInt128.cs @@ -77,19 +77,19 @@ internal UInt128 (ulong s1, ulong s0) { // Implicit (widening) casts /// - /// Creates an unsigned 128-bit integer from an unsigned 64-bit integer. + /// Converts an unsigned 64-bit integer to an unsigned 128-bit integer. /// - /// The unsigned 64-bit integer. - /// This is x widening cast, which preserves the value and is thus safe to perform without an explicit cast. + /// The unsigned 64-bit integer. + /// The equivilent 128-bit unsigned integer. public static implicit operator UInt128 (ulong u) { return (new UInt128(0UL, u)); } /// - /// Creates a BigInteger from an unsigned 128-bit integer. + /// Converts an unsigned 128-bit integer into an arbitrary-size big integer. /// - /// The unishged 128-bit integer. - /// This is x widening cast, which preserves the value and is thus safe to perform without an explicit cast. + /// The unisgned 128-bit integer. + /// The equivilent . public static implicit operator BigInteger (UInt128 u) { byte[] bytes = u.GetBytesForBigInteger(false); return new BigInteger(bytes); @@ -113,9 +113,10 @@ internal byte[] GetBytesForBigInteger (bool isSigned) { } /// - /// Creates an double from an unsigned 128-bit integer. + /// Converts an unsigned 128-bit integer into a floating point value. /// /// The unsigned 128-bit integer. + /// The nearest floating point value. /// This cast will never fail, but it will loose precision /// for values larger than about 1016, because /// maintains only about 52 bits of precision. @@ -127,9 +128,10 @@ public static implicit operator double (UInt128 u) { // Explicit (narrowing) casts /// - /// Creates an unsigned 64-bit integer from an unsigned 128-bit integer. + /// Converts an unsigned 128-bit integer into an unsigned 64-bit integer. /// - /// The unsigned 128-bit integer. + /// The unsigned 128-bit integer. + /// The unsigned 64-bit integer with the same lower 64 bits. /// This is x narrowing cast, which discards high-order bits if /// the 128-bit integer is greater than . public static explicit operator ulong (UInt128 u) { @@ -137,11 +139,10 @@ public static explicit operator ulong (UInt128 u) { } /// - /// Creates an unsigned 128-bit integer from a BigInteger. + /// Converts an arbitrary-size big integer into an unsigned 128-bit integer. /// - /// The BigInteger - /// This is a narrowing cast, because the range of is larger than - /// the range of . + /// A value. + /// The unsigned 128-bit integer with the same lower 128 bits. public static explicit operator UInt128 (BigInteger b) { byte[] bytes = b.ToByteArray(); // BigInteger produces the minimum length byte array to represent the number. @@ -162,6 +163,12 @@ public static explicit operator UInt128 (BigInteger b) { return new UInt128(s1, s0); } + /// + /// Converts a floating-point value into an unsigned 128-bit integer. + /// + /// A floating-point number. + /// The lower 128 bits of the integer part of the floating-point number. + /// is negative, or NaN, or infinite. public static explicit operator UInt128 (double x) { DoubleInfo s = new DoubleInfo(Math.Floor(x)); if (s.IsNegative || !s.IsFinite) throw new InvalidCastException(); @@ -179,7 +186,7 @@ public static explicit operator UInt128 (double x) { /// /// The first integer. /// The second integer. - /// True if equals , otherwise false. + /// if equals , otherwise . public static bool Equals (UInt128 x, UInt128 y) { return (x.s0 == y.s0) && (x.s1 == y.s1); } @@ -189,7 +196,7 @@ public static bool Equals (UInt128 x, UInt128 y) { /// /// The first integer. /// The second integer. - /// True if equals , otherwise false. + /// if equals , otherwise . public static bool operator == (UInt128 x, UInt128 y) { return Equals(x, y); } @@ -199,7 +206,7 @@ public static bool Equals (UInt128 x, UInt128 y) { /// /// The first integer. /// The second integer. - /// True if does not equal , otherwise false. + /// if equals , otherwise . public static bool operator != (UInt128 x, UInt128 y) { return !Equals(x, y); } @@ -208,7 +215,7 @@ public static bool Equals (UInt128 x, UInt128 y) { /// Tests whether the current instance is equal to another unsigned 128-bit integer. /// /// The other integer. - /// True if this equals , otherwise false. + /// if the current instance equals , otherwise . public bool Equals (UInt128 x) { return Equals(this, x); } @@ -217,7 +224,7 @@ public bool Equals (UInt128 x) { /// Tests whether the current instance is equal to another object. /// /// The other object. - /// True if the other object is an equal , otherwise false. + /// if the other object is an equal , otherwise . public override bool Equals (object obj) { if (obj is UInt128) { return (Equals(this, (UInt128) obj)); @@ -226,6 +233,10 @@ public override bool Equals (object obj) { } } + /// + /// Gets a hash code for the current instance. + /// + /// A hash code for the current instance. public override int GetHashCode () { unchecked { return (s0.GetHashCode() + 13 * s1.GetHashCode()); @@ -366,7 +377,7 @@ public int CompareTo (UInt128 u) { /// /// The dividend. /// The divisor. - /// The integer quotient / . + /// The integer quotient / . public static UInt128 operator / (UInt128 x, UInt128 y) { Int128Calculator.Divide128By128(x.s1, x.s0, y.s1, y.s0, out ulong s1, out ulong s0, out _, out _); return new UInt128(s1, s0); @@ -650,6 +661,14 @@ public static UInt128 DivRem (UInt128 x, uint y, out uint r) { // Serialization and deserialization + /// + /// Creates an unsigned 128-bit integer from its string representation. + /// + /// The base-10 string representation of an unsigned 128-bit integer. + /// The integer value it represents. + /// is . + /// is not a valid base-10 representation of an unsigned integer. + /// represents an unsigned integer outside the range of . public static UInt128 Parse (string text) { ParseResult result = TryParse(text, false, out UInt128 value); if (result == ParseResult.Null) throw new ArgumentNullException(nameof(text)); @@ -663,8 +682,8 @@ public static UInt128 Parse (string text) { /// /// The string to parse. /// The integer value it represents. - /// True if is x valid representation of x - /// x non-negative integer, otherwise false. + /// if was successfully + /// parsed as an unsigned 128-bit integer, otherwise . public static bool TryParse (string text, out UInt128 value) { ParseResult result = TryParse(text, false, out value); return (result == ParseResult.Success); @@ -711,6 +730,10 @@ internal static ParseResult TryParse (string text, bool allowNegative, out UInt1 } + /// + /// Produces a string representation of the unsigned 128-bit integer. + /// + /// A string containing the base-10 representation of the integer value. public override string ToString () { if (this == UInt128.Zero) return "0"; diff --git a/Numerics/Functions/AdvancedMath_Exponential.cs b/Numerics/Functions/AdvancedMath_Exponential.cs index adbaec8..ff46798 100644 --- a/Numerics/Functions/AdvancedMath_Exponential.cs +++ b/Numerics/Functions/AdvancedMath_Exponential.cs @@ -62,11 +62,11 @@ private static double IntegralEi_Asymptotic (double x) { /// /// is negative. /// - /// + /// /// /// public static double IntegralEi (double x) { - if (x < 0) { + if (x < 0.0) { throw new ArgumentOutOfRangeException(nameof(x)); } else if (x < 40.0) { return (EulerGamma + Math.Log(x) + IntegralEi_Series(x)); @@ -125,7 +125,7 @@ public static double IntegralE (int n, double x) { return (Math.Exp(-x) / x); } - // Other x > 0 and n > 0. + // Now we are sure x > 0 and n > 0. if (x < 2.0) { return (IntegralE_Series(n, x)); } else if (x < expLimit) { @@ -208,7 +208,7 @@ private static double IntegralE_ContinuedFraction (int n, double x) { /// is negative. /// /// - /// + /// public static double IntegralCi (double x) { if (x < 0.0) { throw new ArgumentOutOfRangeException(nameof(x)); @@ -257,7 +257,7 @@ public static double IntegralCin (double x) { /// /// /// - /// + /// public static double IntegralSi (double x) { if (x < 0.0) { return (-IntegralSi(-x)); diff --git a/Numerics/Functions/AdvancedMath_Jacobi.cs b/Numerics/Functions/AdvancedMath_Jacobi.cs index 0749125..16ece5b 100644 --- a/Numerics/Functions/AdvancedMath_Jacobi.cs +++ b/Numerics/Functions/AdvancedMath_Jacobi.cs @@ -34,8 +34,7 @@ public static double JacobiSn (double u, double k) { if (k == 1.0) return (Math.Tanh(u)); // Decompose u = u_0 K + u_1, where -K/2 < u_1 < K/2, and compute sn(u_1) - long u0; double u1; - double sn = JacobiSn_ViaRangeReduction(u, k, out u0, out u1); + double sn = JacobiSn_ViaRangeReduction(u, k, out long u0, out double u1); // Transform to the appropriate quadrant switch (MoreMath.Mod(u0, 4)) { @@ -74,9 +73,7 @@ public static double JacobiCn (double u, double k) { if (k == 1.0) return (1.0 / Math.Cosh(u)); // Decompose u = u_0 K + u_1, where -K/2 < u_1 < K/2, and compute sn(u_1) - long u0; - double u1; - double sn = JacobiSn_ViaRangeReduction(u, k, out u0, out u1); + double sn = JacobiSn_ViaRangeReduction(u, k, out long u0, out double u1); switch (MoreMath.Mod(u0, 4)) { case 0: @@ -114,8 +111,7 @@ public static double JacobiDn (double u, double k) { if (k == 1.0) return (1.0 / Math.Cosh(u)); // Decompose u = u_0 (K / 2) + u_1, where -K/4 < u_1 < K/4, and compute sn(u_1) - long u0; double u1; - double sn = JacobiSn_ViaRangeReduction(u, k, out u0, out u1); + double sn = JacobiSn_ViaRangeReduction(u, k, out long u0, out double u1); double t = k * sn; /* diff --git a/Numerics/Functions/AdvancedMath_ModifiedBessel.cs b/Numerics/Functions/AdvancedMath_ModifiedBessel.cs index f852741..8ca3d82 100644 --- a/Numerics/Functions/AdvancedMath_ModifiedBessel.cs +++ b/Numerics/Functions/AdvancedMath_ModifiedBessel.cs @@ -112,7 +112,8 @@ private static void ModifiedBesselK_RecurrUpward (double mu, double x, ref doubl /// regeime, you can still obtain the value of the scaled modified e-x Iν(x) /// by calling . /// - /// + /// + /// public static double ModifiedBesselI (double nu, double x) { return ModifiedBesselI(nu, x, false); } @@ -123,6 +124,8 @@ public static double ModifiedBesselI (double nu, double x) { /// The order parameter. /// The argument, which must be non-negative. /// The value of e-x Iν(x). + /// + /// public static double ScaledModifiedBesselI (double nu, double x) { return ModifiedBesselI(nu, x, true); } diff --git a/Numerics/Functions/SolutionPair.cs b/Numerics/Functions/SolutionPair.cs index a046025..d22238b 100644 --- a/Numerics/Functions/SolutionPair.cs +++ b/Numerics/Functions/SolutionPair.cs @@ -19,7 +19,7 @@ namespace Meta.Numerics.Functions { /// a matter of convention. When one solution is regular (finite) at the origin and the other is not, we take the regular solution /// to be the first. /// - public struct SolutionPair : IEquatable { + public struct SolutionPair { private readonly double j, jPrime, y, yPrime; @@ -90,7 +90,8 @@ public SolutionPair (double firstSolutionValue, double firstSolutionDerivative, } // Equality - + // These are not yet public and of very questionable necessity, so commenting them out for now. + /* public override bool Equals (object obj) { if (obj is SolutionPair other) { return (Equals(this, other)); @@ -120,6 +121,7 @@ public override int GetHashCode () { return (j.GetHashCode() + 5 * jPrime.GetHashCode() + 7 * y.GetHashCode() + 13 * yPrime.GetHashCode()); } } + */ // Fitting /* diff --git a/Numerics/Functions/SpinState.cs b/Numerics/Functions/SpinState.cs index 5458778..881ee96 100644 --- a/Numerics/Functions/SpinState.cs +++ b/Numerics/Functions/SpinState.cs @@ -8,7 +8,7 @@ namespace Meta.Numerics.Functions { /// /// Represents the state of a spinor. /// - public struct SpinState { + public struct SpinState : IEquatable { /// /// Instantiates a new SpinState with the given spin and magnetic quantum numbers. @@ -105,7 +105,16 @@ public Spin Representation { // equality private static bool Equals (SpinState s1, SpinState s2) { - return ((s1.spin == s2.spin) && (s1.twoM == s2.twoM)); + return (s1.spin == s2.spin) && (s1.twoM == s2.twoM); + } + + /// + /// Tests whether another instance equals the current instance. + /// + /// The other instance. + /// if the current instance is equal to , otherwise . + public bool Equals (SpinState other) { + return Equals(this, other); } /// @@ -113,9 +122,9 @@ private static bool Equals (SpinState s1, SpinState s2) { /// /// The first spin state. /// The second spin state. - /// True if and are equal, otherwise false. + /// if and are equal, otherwise . public static bool operator == (SpinState s1, SpinState s2) { - return (Equals(s1, s2)); + return Equals(s1, s2); } /// @@ -123,18 +132,18 @@ private static bool Equals (SpinState s1, SpinState s2) { /// /// The first spin state. /// The second spin state. - /// False if and are equal, otherwise true. + /// if and are equal, otherwise . public static bool operator != (SpinState s1, SpinState s2) { - return (!Equals(s1, s2)); + return !Equals(s1, s2); } /// /// Determines whether the given object represents the same spin state. /// /// The object to compare. - /// True if is an equal spin state, otherwise false. + /// if is an equal spin state, otherwise . public override bool Equals (object obj) { - return ((obj is SpinState) && Equals(this, (SpinState)obj)); + return (obj is SpinState) && Equals(this, (SpinState) obj); } /// @@ -142,8 +151,8 @@ public override bool Equals (object obj) { /// /// An integer which is guaranteed equal for equal spin states, an unlikely to be equal for unequal spin states. public override int GetHashCode () { - // left-shift 2j by about half the with of an int (which is 32 bits) so 2j and 2m values are unlikely to interfere - return ((TwoJ << 14) ^ (TwoM)); + // left-shift 2j by about half the width of an int (which is 32 bits) so 2j and 2m values are unlikely to interfere + return (TwoJ << 14) ^ (TwoM); } /// @@ -151,15 +160,20 @@ public override int GetHashCode () { /// /// A string representation of the spin state. public override string ToString () { + return ToString(CultureInfo.CurrentCulture); + } + + private string ToString (IFormatProvider provider) { + NumberFormatInfo numberFormat = (NumberFormatInfo) provider.GetFormat(typeof(NumberFormatInfo)); StringBuilder text = new StringBuilder(spin.ToString()); - text.AppendFormat(CultureInfo.CurrentCulture, ","); + text.Append(","); if (twoM > 0) { - text.Append(CultureInfo.CurrentCulture.NumberFormat.PositiveSign); + text.Append(numberFormat.PositiveSign); } else if (twoM < 0) { - text.Append(CultureInfo.CurrentCulture.NumberFormat.NegativeSign); + text.Append(numberFormat.NegativeSign); } text.Append(SpinString(Math.Abs(twoM))); - return (text.ToString()); + return text.ToString(); } private static string SpinString (int twoJ) { diff --git a/Numerics/Numerics.csproj b/Numerics/Numerics.csproj index 8de7f42..a14e03e 100644 --- a/Numerics/Numerics.csproj +++ b/Numerics/Numerics.csproj @@ -1,28 +1,33 @@ - + netstandard1.1 Meta.Numerics Meta.Numerics - 4.0.8 - Copyright © David Wright 2008-2018 + 4.1.0-alpha + Copyright © David Wright 2008-2020 David Wright Meta Numerics A library for numeric computing with support for data manipulation, statistics, matrix algebra, and advanced functions. http://www.meta-numerics.net - numeric scientific technical math statistics matrix fft + numeric scientific technical math statistics matrix fft data optimization http://www.meta-numerics.net/Images/ComplexPsiPlot.png https://github.com/dcwuser/metanumerics - Full + portable + true + https://opensource.org/licenses/MS-PL + false TRACE;RELEASE;NETSTANDARD1_1 Numerics1.ruleset + true Numerics1.ruleset + diff --git a/Numerics/Statistics/AR1FitResult.cs b/Numerics/Statistics/AR1FitResult.cs index 4c90ffb..47c814b 100644 --- a/Numerics/Statistics/AR1FitResult.cs +++ b/Numerics/Statistics/AR1FitResult.cs @@ -29,7 +29,7 @@ internal AR1FitResult(UncertainValue mu, UncertainValue alpha, UncertainValue si /// mean of the underlying population from which the series was drawn. public UncertainValue Mu { get { - return (mu); + return mu; } } @@ -40,7 +40,7 @@ public UncertainValue Mu { /// lag-1 auto-correlation of the underlying population from which the series was drawn. public UncertainValue Alpha { get { - return (alpha); + return alpha; } } @@ -51,16 +51,17 @@ public UncertainValue Alpha { /// per-step Gaussian noise in the underlying population from which the series was drawn. public UncertainValue Sigma { get { - return (sigma); + return sigma; } } /// - /// Gets the residuals of the time series fit. + /// Gets a list of the residuals of the fit. /// + /// A read-only list, in the same order as the original data, of the difference between each measured and predicted value. public IReadOnlyList Residuals { get { - return (residuals); + return residuals; } } diff --git a/Numerics/Statistics/Distributions/BetaDistribution.cs b/Numerics/Statistics/Distributions/BetaDistribution.cs index 1ba670b..1268f20 100644 --- a/Numerics/Statistics/Distributions/BetaDistribution.cs +++ b/Numerics/Statistics/Distributions/BetaDistribution.cs @@ -240,8 +240,7 @@ private IEnumerator CentralMoments () { public override double InverseLeftProbability (double P) { if ((P < 0.0) || (P > 1.0)) throw new ArgumentOutOfRangeException(nameof(P)); - double x, y; - betaInverter.InverseRegularizedBeta(P, 1.0 - P, out x, out y); + betaInverter.InverseRegularizedBeta(P, 1.0 - P, out double x, out double y); return (x); } @@ -249,8 +248,7 @@ public override double InverseLeftProbability (double P) { public override double InverseRightProbability (double Q) { if ((Q < 0.0) || (Q > 1.0)) throw new ArgumentOutOfRangeException(nameof(Q)); - double x, y; - betaInverter.InverseRegularizedBeta(1.0 - Q, Q, out x, out y); + betaInverter.InverseRegularizedBeta(1.0 - Q, Q, out double x, out double y); return (x); } @@ -462,8 +460,7 @@ private void RefineInverseRegularizedBeta (double P0, double Q0, ref double x, r //double y_old = y; // Evaluate I and I' at the coordinates - double P, Q, D; - IncompleteBetaAndDerivative(x, y, out P, out Q, out D); + IncompleteBetaAndDerivative(x, y, out double P, out double Q, out double D); Debug.Assert(P + Q == 1.0); // Compute the deficit from the smaller of P and Q, to reduce inaccuracies when P or Q is near 1 @@ -529,7 +526,7 @@ public void InverseRegularizedBeta (double P, double Q, out double x, out double // R. C. H. Cheng, Generating Beta Variates with Nonintegral Shape Parameters, // Communications of the ACM 21 (1978) 317 - public class ChengBetaGenerator : IDeviateGenerator { + internal class ChengBetaGenerator : IDeviateGenerator { public ChengBetaGenerator (double a, double b) { Debug.Assert(a > 1.0); diff --git a/Numerics/Statistics/Distributions/ChiDistribution.cs b/Numerics/Statistics/Distributions/ChiDistribution.cs index ce544dc..5968d9b 100644 --- a/Numerics/Statistics/Distributions/ChiDistribution.cs +++ b/Numerics/Statistics/Distributions/ChiDistribution.cs @@ -28,7 +28,7 @@ public sealed class ChiDistribution : ContinuousDistribution { public ChiDistribution (int nu) { if (nu < 1) throw new ArgumentOutOfRangeException(nameof(nu)); this.nu = nu; - this.gamma = new GammaDistribution(nu / 2.0); + this.gamma = new GammaDistribution(0.5 * nu); } private readonly int nu; @@ -55,7 +55,7 @@ public override double ProbabilityDensity (double x) { if (x < 0.0) { return (0.0); } else { - return (x * gamma.ProbabilityDensity(x * x / 2.0)); + return (x * gamma.ProbabilityDensity(0.5 * x * x)); } } @@ -64,7 +64,7 @@ public override double LeftProbability (double x) { if (x < 0.0) { return (0.0); } else { - return (gamma.LeftProbability(x * x / 2.0)); + return (gamma.LeftProbability(0.5 * x * x)); } } @@ -73,7 +73,7 @@ public override double RightProbability (double x) { if (x <= 0.0) { return (1.0); } else { - return (gamma.RightProbability(x * x / 2.0)); + return (gamma.RightProbability(0.5 * x * x)); } } @@ -94,7 +94,7 @@ public override double RawMoment (int r) { if (r < 0) { throw new ArgumentOutOfRangeException(nameof(r)); } else { - return (Math.Pow(2.0, r / 2.0) * AdvancedMath.Pochhammer(nu / 2.0, r / 2.0)); + return (Math.Pow(2.0, 0.5 * r) * AdvancedMath.Pochhammer(0.5 * nu, 0.5 * r)); } } @@ -108,7 +108,7 @@ public override double Variance { // C_2 = \frac{1}{2} - \frac{1}{8\nu} - \frac{1}{16\nu^2} + \frac{5}{128\nu^3} + \cdots // which we use for large \nu. if (nu < 32) { - return (nu - 2.0 * MoreMath.Sqr(AdvancedMath.Pochhammer(nu / 2.0, 1.0 / 2.0))); + return (nu - 2.0 * MoreMath.Sqr(AdvancedMath.Pochhammer(0.5 * nu, 0.5))); } else { double C = varianceCoefficients[0]; double nuk = 1.0; diff --git a/Numerics/Statistics/Distributions/DiscreteAsContinuousDistribution.cs b/Numerics/Statistics/Distributions/DiscreteAsContinuousDistribution.cs index c2cae08..57a3150 100644 --- a/Numerics/Statistics/Distributions/DiscreteAsContinuousDistribution.cs +++ b/Numerics/Statistics/Distributions/DiscreteAsContinuousDistribution.cs @@ -10,8 +10,8 @@ namespace Meta.Numerics.Statistics.Distributions { // Here is a way that we can get everything we want and have P(x) increase continuously rather than in jumps. Define // z = (x - x_min) / (x_max - x_min) - // to be the fraction of the way through the continuos interval and define - // k = k_min + z (k_max - k_min + 1) + // to be the fraction of the way through the continuous interval and define + // k = (1 - z) * k_min + z * (k_max + 1) // to be the corresponding "effective k". Overall this means // k = k_min + m (x - x_min) m = (k_max - k_min + 1)/(x_max - x_min) // Split k = [k] + {k} into into integer [k] and fractional {k} parts. Then let @@ -19,8 +19,8 @@ namespace Meta.Numerics.Statistics.Distributions { // where P_L([k]) is the exclusive left probability for bin [k] and P([k]) is the probability for bin [k]. // Then - // x = a + e => z = 0 + e => k = k_min + e => P = P_L(k_min) + e P(k_min) = 0 + e - // x = b - e => z = 1 - e => k = k_max + 1 - e => P = P_L(k_max) + (1 - e) P(k_max) = 1 - e + // x = a + e => z = 0 + e => k = k_min + e => P = P_L(k_min) + e P(k_min) = e P(k_min) + // x = b - e => z = 1 - e => k = (k_max + 1) - e => P = P_L(k_max) + (1 - e) P(k_max) = P_R(k_max) - e P(k_max) // x = (a+b)/2 => z = 1/2 => k = (k_max + k_min + 1) / 2 => P = P_L((k_max + k_min)/2) + 1/2 P((k_max + k_min)/2) = half of mid-bin // which is what we wanted. @@ -55,11 +55,12 @@ public DiscreteAsContinuousDistribution (DiscreteDistribution distribution, Inte private DiscreteDistribution d; private Interval xSupport; - private void ComputeEffectiveBin (double x, out int ki, out double kf) { + private void ComputeEffectiveBin (double x, out int k0, out double k1) { double z = (x - xSupport.LeftEndpoint) / xSupport.Width; - double k = d.Support.LeftEndpoint + z * d.Support.Width; - ki = (int) Math.Floor(k); - kf = k - ki; + double k = (1.0 - z) * d.Support.LeftEndpoint + z * (d.Support.RightEndpoint + 1); + k0 = (int) Math.Floor(k); + k1 = k - k0; + Debug.Assert(0.0 <= k1 && k1 < 1.0); } private double ComputeEffectivePoint (double k) { @@ -77,9 +78,8 @@ public override double LeftProbability (double x) { } else if (x > xSupport.RightEndpoint) { return (1.0); } else { - int ki; double kf; - ComputeEffectiveBin(x, out ki, out kf); - return (d.LeftExclusiveProbability(ki) + kf * d.ProbabilityMass(ki)); + ComputeEffectiveBin(x, out int k0, out double k1); + return (d.LeftExclusiveProbability(k0) + k1 * d.ProbabilityMass(k0)); } } @@ -90,9 +90,8 @@ public override double RightProbability (double x) { } else if (x > xSupport.RightEndpoint) { return (0.0); } else { - int ki; double kf; - ComputeEffectiveBin(x, out ki, out kf); - return ((1.0 - kf) * d.ProbabilityMass(ki) + d.RightExclusiveProbability(ki)); + ComputeEffectiveBin(x, out int k0, out double k1); + return ((1.0 - k1) * d.ProbabilityMass(k0) + d.RightExclusiveProbability(k0)); } } diff --git a/Numerics/Statistics/Distributions/DiscreteDistribution.cs b/Numerics/Statistics/Distributions/DiscreteDistribution.cs index e606c19..a29b25d 100644 --- a/Numerics/Statistics/Distributions/DiscreteDistribution.cs +++ b/Numerics/Statistics/Distributions/DiscreteDistribution.cs @@ -13,8 +13,6 @@ namespace Meta.Numerics.Statistics.Distributions { /// public abstract class DiscreteDistribution : UnivariateDistribution { - // replace ProbabilityDensity with ProbabilityMass - /// /// Returns the probability of the obtaining the given value. /// diff --git a/Numerics/Statistics/Distributions/GeometricDistribution.cs b/Numerics/Statistics/Distributions/GeometricDistribution.cs index 8523807..feb274e 100644 --- a/Numerics/Statistics/Distributions/GeometricDistribution.cs +++ b/Numerics/Statistics/Distributions/GeometricDistribution.cs @@ -168,11 +168,9 @@ private static double EulerianPolynomial (int n, double x) { } - /// /// - /// public override int GetRandomValue (Random rng) { - if (rng == null) throw new ArgumentNullException(nameof(rng)); + if (rng is null) throw new ArgumentNullException(nameof(rng)); double u = rng.NextDouble(); return ((int) Math.Min(Math.Floor(Math.Log(u) / lnq), Int32.MaxValue)); } diff --git a/Numerics/Statistics/Distributions/IDeviateGenerator.cs b/Numerics/Statistics/Distributions/IDeviateGenerator.cs index d4cb8fc..aadd745 100644 --- a/Numerics/Statistics/Distributions/IDeviateGenerator.cs +++ b/Numerics/Statistics/Distributions/IDeviateGenerator.cs @@ -202,14 +202,14 @@ internal class CauchyGenerator : IDeviateGenerator { public double GetNext (Random rng) { - // pick a random point within the unit semicircle using rejection + // Pick a random point within the unit semicircle using rejection. double x, y; do { x = 2.0 * rng.NextDouble() - 1.0; y = rng.NextDouble(); } while ((x * x + y * y > 1.0) || (y == 0.0)); - // its tangent is the tangent of a random angle and is thus Cauchy distributed + // Its tangent is the tangent of a random angle and is thus Cauchy distributed. return (x / y); } @@ -229,7 +229,7 @@ public AhrensDieterLowAlphaGammaGenerator (double alpha) { b = (Math.E + a) / Math.E; } - double a, b; + private readonly double a, b; public double GetNext (Random rng) { diff --git a/Numerics/Statistics/Distributions/MannWhitneyDistribution.cs b/Numerics/Statistics/Distributions/MannWhitneyDistribution.cs index 05b8082..0df762e 100644 --- a/Numerics/Statistics/Distributions/MannWhitneyDistribution.cs +++ b/Numerics/Statistics/Distributions/MannWhitneyDistribution.cs @@ -1,6 +1,8 @@ using System; using System.Diagnostics; +using Meta.Numerics.Extended; + namespace Meta.Numerics.Statistics.Distributions { /// @@ -29,28 +31,26 @@ public MannWhitneyExactDistribution (int m, int n) { } - private int m, n; - private double total; - private decimal[] counts; + private readonly int m, n; + private readonly Int128[] counts; + private readonly double total; - // this routine is based on the recursion - // [ m n ] = ( 1 - q^m ) / ( 1 - q^(m-n) ) [ m-1 n ] - // and the starting point [ n n ] = 1 + // This routine is based on the recursion + // [ m n ] = ( 1 - q^m ) / ( 1 - q^(m-n) ) [ m-1 n ] + // and the starting point [ n n ] = 1. - // the coefficients are integers and get large quickly as m and n increase - // we use decimal because it handles larger integers than long - // we can't use double because the calculation requires delicate cancelations - // among large intermediate values, thus necessitating exact integer arithmetic - // look into using an arbitrary-sized integer structure in the future + // The coefficients are integers and get large quickly as m and n increase. + // First we used decimal because it was the has the largest precise integer + // values of any built-in, but now that we have Int128, we use that instead. - private static decimal[] GaussianBinomialCoefficients (int m, int n) { + private static Int128[] GaussianBinomialCoefficients (int m, int n) { Debug.Assert(m >= 0); Debug.Assert(n >= 0); Debug.Assert(m >= n); // create an array to hold our coefficients - decimal[] c = new decimal[(m - n) * n + 1]; + Int128[] c = new Int128[(m - n) * n + 1]; // start with [n n] = 1 * q^0 c[0] = 1; @@ -60,7 +60,7 @@ private static decimal[] GaussianBinomialCoefficients (int m, int n) { // create a scratch array for intermediate use // it needs to be larger than the previous array by (m-n) to hold intermediate polynomials - decimal[] b = new decimal[c.Length + (m - n)]; + Int128[] b = new Int128[c.Length + (m - n)]; // iterate from [n n] up to [m n] for (int k = n + 1; k <= m; k++) { diff --git a/Numerics/Statistics/Distributions/MomentMath.cs b/Numerics/Statistics/Distributions/MomentMath.cs index cd0c577..f334a61 100644 --- a/Numerics/Statistics/Distributions/MomentMath.cs +++ b/Numerics/Statistics/Distributions/MomentMath.cs @@ -34,8 +34,6 @@ public static double[] RawToCentral (double[] M) { C[1] = 0.0; - if (C.Length == 2) return (C); - for (int r = 2; r < C.Length; r++) { C[r] = RawToCentral(M, r); } @@ -48,6 +46,10 @@ public static double[] RawToCentral (double[] M) { internal static double RawToCentral (double[] M, int r) { + Debug.Assert(M.Length > 1); + Debug.Assert(0 <= r && r < M.Length); + Debug.Assert(M[0] == 1.0); + double mmu = -M[1]; IEnumerator B = AdvancedIntegerMath.BinomialCoefficients(r).GetEnumerator(); @@ -95,8 +97,6 @@ public static double[] CentralToRaw (double mu, double[] C) { M[1] = mu; - if (M.Length == 2) return (M); - for (int r = 2; r < M.Length; r++) { M[r] = CentralToRaw(mu, C, r); } @@ -109,6 +109,11 @@ public static double[] CentralToRaw (double mu, double[] C) { internal static double CentralToRaw (double mu, double[] C, int r) { + Debug.Assert(C.Length > 0); + Debug.Assert(0 <= r && r < C.Length); + Debug.Assert(C[0] == 1.0); + //Debug.Assert(C[1] == 0.0); + IEnumerator B = AdvancedIntegerMath.BinomialCoefficients(r).GetEnumerator(); // k = 0 term @@ -149,8 +154,6 @@ public static double[] CumulantToRaw (double[] K) { M[0] = 1.0; - if (M.Length == 1) return (M); - for (int r = 1; r < M.Length; r++) { M[r] = CumulantToRaw(K, r); } @@ -175,14 +178,16 @@ public static double[] CumulantToCentral (double[] K) { if (K == null) throw new ArgumentNullException(nameof(K)); double[] C = new double[K.Length]; + if (C.Length == 0) return (C); if (K[0] != 0.0) throw new ArgumentOutOfRangeException(nameof(K)); + C[0] = 1.0; + if (C.Length == 1) return (C); C[1] = 0.0; - if (C.Length == 2) return (C); for (int r = 2; r < C.Length; r++) { C[r] = CumulantToCentral(K, r); diff --git a/Numerics/Statistics/Distributions/NegativeBinomialDistribution.cs b/Numerics/Statistics/Distributions/NegativeBinomialDistribution.cs index 9b6bc9e..8dd5a3e 100644 --- a/Numerics/Statistics/Distributions/NegativeBinomialDistribution.cs +++ b/Numerics/Statistics/Distributions/NegativeBinomialDistribution.cs @@ -68,10 +68,19 @@ public override double Mean { /// public override double LeftInclusiveProbability (int k) { - if (k < 0) { - return (0.0); + if (k == Int32.MaxValue) { + return 1.0; + } else { + return LeftExclusiveProbability(k + 1); + } + } + + /// + public override double LeftExclusiveProbability (int k) { + if (k < 1) { + return 0.0; } else { - return (AdvancedMath.LeftRegularizedBeta(r, k + 1, q)); + return AdvancedMath.LeftRegularizedBeta(r, k, q); } } @@ -106,6 +115,7 @@ public override double Skewness { } /// + /// The order of moment to compute. public override double RawMoment (int k) { if (k < 0) { diff --git a/Numerics/Statistics/Distributions/PoissonDistribution.cs b/Numerics/Statistics/Distributions/PoissonDistribution.cs index d0a7207..87d89c2 100644 --- a/Numerics/Statistics/Distributions/PoissonDistribution.cs +++ b/Numerics/Statistics/Distributions/PoissonDistribution.cs @@ -20,7 +20,7 @@ public PoissonDistribution (double mu) { this.mu = mu; if (mu < 4.0) { - poissonGenerator = new PoissonGeneratorAdditive(mu); + poissonGenerator = new PoissonGeneratorMultiplicative(mu); } else if (mu < 16.0) { poissonGenerator = new PoissonGeneratorTabulated(mu); } else { @@ -84,27 +84,30 @@ public override double ExcessKurtosis { /// public override double LeftInclusiveProbability (int k) { - if (k < 0) { - return (0.0); + if (k == Int32.MaxValue) { + return 1.0; } else { - // these are equivalent expressions, but the explicit sum is faster for if the number of terms is small - if (k < 16) { - double ds = 1.0; - double s = ds; - for (int i = 1; i <= k; i++) { - ds *= mu / i; - s += ds; - } - return (Math.Exp(-mu) * s); - } else { - return(AdvancedMath.RightRegularizedGamma(k + 1, mu)); - } + return LeftExclusiveProbability(k + 1); } } /// public override double LeftExclusiveProbability (int k) { - return(LeftInclusiveProbability(k-1)); + if (k < 1) { + return 0.0; + } else if (k < 16) { + // For small k, it's fastest to sum probabilities directly. + double ds = 1.0; + double s = ds; + for (int i = 1; i < k; i++) { + ds *= mu / i; + s += ds; + } + return Math.Exp(-mu) * s; + } else { + // Otherwise, fall back to regularized incomplete Gamma. + return AdvancedMath.RightRegularizedGamma(k, mu); + } } /// @@ -317,30 +320,32 @@ public override double Cumulant (int r) { // form of simple inversion: binary search, starting from the median, using a pre-computed table // of CDF values that covers most of the space. + /// public override int GetRandomValue (Random rng) { - if (rng == null) throw new ArgumentNullException(nameof(rng)); + if (rng is null) throw new ArgumentNullException(nameof(rng)); return (poissonGenerator.GetNext(rng)); } } - internal class PoissonGeneratorAdditive : IDeviateGenerator { + internal class PoissonGeneratorMultiplicative : IDeviateGenerator { - public PoissonGeneratorAdditive (double lambda) { - this.lambda = lambda; + public PoissonGeneratorMultiplicative (double lambda) { + expMinusLambda = Math.Exp(-lambda); } - private readonly double lambda; + private readonly double expMinusLambda; public int GetNext (Random rng) { int k = 0; double t = rng.NextDouble(); - while (t < lambda) { + while (t > expMinusLambda) { k++; - t += rng.NextDouble(); + t *= rng.NextDouble(); } return (k); } + } internal class PoissonGeneratorTabulated : IDeviateGenerator { diff --git a/Numerics/Statistics/GeneralLinearRegressionResult.cs b/Numerics/Statistics/GeneralLinearRegressionResult.cs index 47b1403..965046f 100644 --- a/Numerics/Statistics/GeneralLinearRegressionResult.cs +++ b/Numerics/Statistics/GeneralLinearRegressionResult.cs @@ -7,7 +7,7 @@ namespace Meta.Numerics.Statistics { /// /// Describes the result of any generalized linear regression. /// - public abstract class GeneralLinearRegressionResult : FitResult { + public abstract class GeneralLinearRegressionResult : ResidualsResult { internal GeneralLinearRegressionResult () : base() { this.anova = new Lazy(CreateAnova); @@ -25,7 +25,7 @@ internal GeneralLinearRegressionResult () : base() { /// public virtual double RSquared { get { - return (this.Anova.Factor.SumOfSquares / this.Anova.Total.SumOfSquares); + return anova.Value.Factor.SumOfSquares / anova.Value.Total.SumOfSquares; } } @@ -35,7 +35,7 @@ public virtual double RSquared { /// public virtual TestResult F { get { - return (this.Anova.Factor.Result); + return anova.Value.Factor.Result; } } @@ -44,7 +44,14 @@ public virtual TestResult F { /// public virtual OneWayAnovaResult Anova { get { - return (anova.Value); + return anova.Value; + } + } + + /// + public override double SumOfSquaredResiduals { + get { + return anova.Value.Residual.SumOfSquares; } } diff --git a/Numerics/Statistics/LinearRegressionResult.cs b/Numerics/Statistics/LinearRegressionResult.cs index e885916..5b67c7b 100644 --- a/Numerics/Statistics/LinearRegressionResult.cs +++ b/Numerics/Statistics/LinearRegressionResult.cs @@ -65,7 +65,7 @@ internal LinearRegressionResult (IReadOnlyList x, IReadOnlyList /// public override UncertainValue Intercept { get { - return (new UncertainValue(a, Math.Sqrt(caa))); + return new UncertainValue(a, Math.Sqrt(caa)); } } @@ -74,7 +74,7 @@ public override UncertainValue Intercept { /// public UncertainValue Slope { get { - return (new UncertainValue(b, Math.Sqrt(cbb))); + return new UncertainValue(b, Math.Sqrt(cbb)); } } @@ -86,16 +86,20 @@ public UncertainValue Slope { public UncertainValue Predict (double x) { double y = a + b * x; double yVariance = sigmaSquared * (1.0 + (MoreMath.Sqr(x - xMean) / xVariance + 1.0) / n); - return (new UncertainValue(y, Math.Sqrt(yVariance))); + return new UncertainValue(y, Math.Sqrt(yVariance)); } - /// - /// Gets the residuals. - /// - /// A list of the difference between the actual and predicted value for each point. - public IReadOnlyList Residuals { + /// + public override IReadOnlyList Residuals { + get { + return residuals; + } + } + + /// + public override double SumOfSquaredResiduals { get { - return (residuals); + return SSR; } } @@ -104,7 +108,7 @@ public IReadOnlyList Residuals { /// public TestResult R { get { - return (rTest.Value); + return rTest.Value; } } @@ -112,7 +116,7 @@ internal override ParameterCollection CreateParameters () { ParameterCollection parameters = new ParameterCollection( nameof(Intercept), a, caa, nameof(Slope), b, cbb, cab ); - return (parameters); + return parameters; } internal override OneWayAnovaResult CreateAnova () { @@ -120,7 +124,7 @@ internal override OneWayAnovaResult CreateAnova () { AnovaRow residual = new AnovaRow(SSR, n - 2); AnovaRow total = new AnovaRow(SST, n - 1); OneWayAnovaResult anova = new OneWayAnovaResult(fit, residual, total); - return (anova); + return anova; } } diff --git a/Numerics/Statistics/MA1FitResult.cs b/Numerics/Statistics/MA1FitResult.cs index d005a11..b86468b 100644 --- a/Numerics/Statistics/MA1FitResult.cs +++ b/Numerics/Statistics/MA1FitResult.cs @@ -29,7 +29,7 @@ internal MA1FitResult (UncertainValue mu, UncertainValue beta, UncertainValue si /// mean of the underlying population from which the series was drawn. public UncertainValue Mu { get { - return (mu); + return mu; } } @@ -40,7 +40,7 @@ public UncertainValue Mu { /// lag-1 auto-correlation of the underlying population from which the series was drawn. public UncertainValue Beta { get { - return (beta); + return beta; } } @@ -51,16 +51,17 @@ public UncertainValue Beta { /// per-step Gaussian noise in the underlying population from which the series was drawn. public UncertainValue Sigma { get { - return (sigma); + return sigma; } } /// - /// Gets the residuals of the time series fit. + /// Gets a list of the residuals of the fit. /// + /// A read-only list, in the same order as the original data, of the difference between each measured and predicted value. public IReadOnlyList Residuals { get { - return (residuals); + return residuals; } } @@ -70,7 +71,7 @@ public IReadOnlyList Residuals { /// A Ljung-Box test for non-correlation of the residuals. public TestResult GoodnessOfFit { get { - return (goodnessOfFit.Value); + return goodnessOfFit.Value; } } diff --git a/Numerics/Statistics/MultiLinearRegressionResult.cs b/Numerics/Statistics/MultiLinearRegressionResult.cs index 69960a8..927f7a8 100644 --- a/Numerics/Statistics/MultiLinearRegressionResult.cs +++ b/Numerics/Statistics/MultiLinearRegressionResult.cs @@ -154,11 +154,8 @@ public UncertainValue Predict (IReadOnlyList x) { } - /// - /// Gets the residuals - /// - /// A list of the differences between each measured and predicted value. - public IReadOnlyList Residuals { + /// + public override IReadOnlyList Residuals { get { return (residuals); } diff --git a/Numerics/Statistics/NonlinearRegressionResult.cs b/Numerics/Statistics/NonlinearRegressionResult.cs index a9240d4..97bff69 100644 --- a/Numerics/Statistics/NonlinearRegressionResult.cs +++ b/Numerics/Statistics/NonlinearRegressionResult.cs @@ -11,7 +11,7 @@ namespace Meta.Numerics.Statistics { /// /// Describes the result of a fit to a non-linear function. /// - public sealed class NonlinearRegressionResult : FitResult { + public sealed class NonlinearRegressionResult : ResidualsResult { internal NonlinearRegressionResult( IReadOnlyList x, IReadOnlyList y, @@ -47,10 +47,12 @@ internal NonlinearRegressionResult( C = cholesky.Inverse(); C = (2.0 * min.Value / (n - d)) * C; + sumOfSquaredResiduals = 0.0; residuals = new List(n); for (int i = 0; i < n; i++) { - double r = y[i] - function(b, x[i]); - residuals.Add(r); + double z = y[i] - function(b, x[i]); + sumOfSquaredResiduals += z * z; + residuals.Add(z); } this.names = names; @@ -62,13 +64,19 @@ internal NonlinearRegressionResult( private readonly SymmetricMatrix C; private readonly Func, double, double> function; private readonly List residuals; + private readonly double sumOfSquaredResiduals; - /// - /// Gets the residuals from the fit. - /// - public IReadOnlyList Residuals { + /// + public override IReadOnlyList Residuals { + get { + return residuals; + } + } + + /// + public override double SumOfSquaredResiduals { get { - return (residuals); + return sumOfSquaredResiduals; } } @@ -78,11 +86,11 @@ public IReadOnlyList Residuals { /// The value of the independent variable. /// The predicted value of the dependent variable. public double Predict (double x) { - return (function(b, x)); + return function(b, x); } internal override ParameterCollection CreateParameters () { - return (new ParameterCollection(names, b, C)); + return new ParameterCollection(names, b, C); } } diff --git a/Numerics/Statistics/PolynomialRegressionResult.cs b/Numerics/Statistics/PolynomialRegressionResult.cs index fb1e2d5..8ff3869 100644 --- a/Numerics/Statistics/PolynomialRegressionResult.cs +++ b/Numerics/Statistics/PolynomialRegressionResult.cs @@ -123,11 +123,8 @@ public UncertainValue Predict (double x) { } - /// - /// Gets the residuals - /// - /// A list of the differences between each measured and predicted value. - public IReadOnlyList Residuals { + /// + public override IReadOnlyList Residuals { get { return (residuals); } diff --git a/Numerics/Statistics/ResidualsResult.cs b/Numerics/Statistics/ResidualsResult.cs new file mode 100644 index 0000000..ac57e5f --- /dev/null +++ b/Numerics/Statistics/ResidualsResult.cs @@ -0,0 +1,29 @@ +using System; +using System.Collections.Generic; +using System.Text; + +namespace Meta.Numerics.Statistics { + + /// + /// Describes the result of a fit with residuals. + /// + public abstract class ResidualsResult : FitResult { + + internal ResidualsResult () : base() { } + + /// + /// Gets a list of the residuals of the fit. + /// + /// A read-only list, in the same order as the original data, of the difference between each measured and predicted value. + public abstract IReadOnlyList Residuals { get; } + + /// + /// Gets the sum of the squares of all residuals. + /// + /// The sum of the squares of all . + public abstract double SumOfSquaredResiduals { get; } + + // We could implement this, but everyone who inherits would override, so we don't. + + } +} diff --git a/Numerics/Statistics/Univariate.cs b/Numerics/Statistics/Univariate.cs index 674de2b..928b8fa 100644 --- a/Numerics/Statistics/Univariate.cs +++ b/Numerics/Statistics/Univariate.cs @@ -58,9 +58,7 @@ public static double Mean (this IReadOnlyCollection sample) { if (sample == null) throw new ArgumentNullException(nameof(sample)); if (sample.Count < 1) return (Double.NaN); - int n; - double mean; - ComputeMomentsUpToFirst(sample, out n, out mean); + ComputeMomentsUpToFirst(sample, out int n, out double mean); return (mean); } @@ -82,9 +80,7 @@ public static double Variance (this IReadOnlyCollection sample) { if (sample == null) throw new ArgumentNullException(nameof(sample)); if (sample.Count < 1) return (Double.NaN); - int n; - double mean, sumOfSquares; - ComputeMomentsUpToSecond(sample, out n, out mean, out sumOfSquares); + ComputeMomentsUpToSecond(sample, out int n, out double mean, out double sumOfSquares); Debug.Assert(sample.Count == n); return (sumOfSquares / n); @@ -101,9 +97,7 @@ public static double Skewness (this IReadOnlyCollection sample) { if (sample == null) throw new ArgumentNullException(nameof(sample)); if (sample.Count < 1) return (Double.NaN); - int n; - double m1, c2, c3; - ComputeMomentsUpToThird(sample, out n, out m1, out c2, out c3); + ComputeMomentsUpToThird(sample, out int n, out double m1, out double c2, out double c3); Debug.Assert(sample.Count == n); Debug.Assert(c2 >= 0.0); @@ -231,9 +225,12 @@ private static void ComputeMomentsUpToFourth(IEnumerable values, out int /// The standard deviation of the sample data. /// /// Note that a few authors use the term "sample standard deviation" to mean "the square root of the population variance as estimated by the sample". - /// We do not adopt this aberrant convention, which is contrary to all other conventional meanings of sample moments. In other words, what is returned - /// by this methods is the square root of the sum of squared deviations divided by n, not divided by (n-1). If you want an estimate of the population - /// standard deviation, use the method. + /// We do not adopt this aberrant convention, which is contrary to the conventional meaning of "sample moment" for all other moments. + /// In other words, what is returned by this methods is the square root of the sum of the squared deviations divided by n, not divided by (n-1). + /// If you want the square root of the sum of the squared deviations divided by n, which is also called the Bessel-corrected standard deviation, + /// use the method. Note, however, that, unlike the variance, the standard deviation + /// corrected in this way is not unbiased. + /// If you want a bias-corrected estimate of the standard deviation of the populace that was sampled, use the method. /// /// is . public static double StandardDeviation (this IReadOnlyCollection sample) { @@ -247,8 +244,8 @@ public static double StandardDeviation (this IReadOnlyCollection sample) /// The Bessel-corrected standard deviation. /// /// This probably isn't the quantity you want, even though it's what many software - /// packages refer to as the "standard deviation". It is the square root of the sum of - /// squared deviations from the mean, divided by n-1 instead of n. + /// packages refer to as the "standard deviation", and some authors even misleadingly call the "sample standard deviation". + /// It is the square root of the sum of the squared deviations from the mean, divided by n-1 instead of n. /// Using n-1 instead of n in the formula for variance produces an /// unbiased estimator of the . But using n-1 instead of n /// in the formula for standard deviation, which is how this property is @@ -268,9 +265,7 @@ public static double CorrectedStandardDeviation (this IReadOnlyCollection x = new List(); + List y = new List(); + for (int i = 0; i < e1[0, 0]; i++) { x.Add(0); y.Add(0); } + for (int i = 0; i < e1[0, 1]; i++) { x.Add(0); y.Add(1); } + for (int i = 0; i < e1[1, 0]; i++) { x.Add(1); y.Add(0); } + for (int i = 0; i < e1[1, 1]; i++) { x.Add(1); y.Add(1); } + double s = x.CorrelationCoefficient(y); + Assert.IsTrue(TestUtilities.IsNearlyEqual(s, e1.Binary.Phi)); } } diff --git a/Test/BivariateSampleTest.cs b/Test/BivariateSampleTest.cs index 65775fc..e0e9edf 100644 --- a/Test/BivariateSampleTest.cs +++ b/Test/BivariateSampleTest.cs @@ -208,7 +208,7 @@ public void LinearLogisticRegression () { // record estimated covariances caa += r.Parameters.CovarianceMatrix[0, 0]; cbb += r.Parameters.CovarianceMatrix[1, 1]; - cab += r.Parameters.CovarianceMatrix[0, 1]; + cab += r.Parameters.CovarianceMatrix[0, 1]; } @@ -277,6 +277,7 @@ public void LinearRegressionSimple () { double SSR = 0.0; foreach (double z in result.Residuals) SSR += z * z; Assert.IsTrue(TestUtilities.IsNearlyEqual(SSR, result.Anova.Residual.SumOfSquares)); + Assert.IsTrue(TestUtilities.IsNearlyEqual(SSR, result.SumOfSquaredResiduals)); // R is same as correlation coefficient Assert.IsTrue(TestUtilities.IsNearlyEqual(x.CorrelationCoefficient(y), result.R.Statistic.Value)); @@ -592,8 +593,7 @@ public void BivariateNonlinearFitVariances () { } [TestMethod] - public void BivariateAssociationDiscreteNullDistribution () - { + public void BivariateAssociationDiscreteNullDistribution () { Random rng = new Random(1); // Pick very non-normal distributions for our non-parameteric tests @@ -601,8 +601,7 @@ public void BivariateAssociationDiscreteNullDistribution () ContinuousDistribution yd = new CauchyDistribution(); // Pick small sample sizes to get exact distributions - foreach (int n in TestUtilities.GenerateIntegerValues(4, 24, 4)) - { + foreach (int n in TestUtilities.GenerateIntegerValues(4, 24, 4)) { // Do a bunch of test runs, recording reported statistic for each. List spearmanStatistics = new List(); @@ -610,38 +609,32 @@ public void BivariateAssociationDiscreteNullDistribution () DiscreteDistribution spearmanDistribution = null; DiscreteDistribution kendallDistribution = null; - for (int i = 0; i < 512; i++) - { + for (int i = 0; i < 512; i++) { List x = new List(); List y = new List(); - for (int j = 0; j < n; j++) - { + for (int j = 0; j < n; j++) { x.Add(xd.GetRandomValue(rng)); y.Add(yd.GetRandomValue(rng)); } DiscreteTestStatistic spearman = Bivariate.SpearmanRhoTest(x, y).UnderlyingStatistic; - if (spearman != null) - { + if (spearman != null) { spearmanStatistics.Add(spearman.Value); spearmanDistribution = spearman.Distribution; } DiscreteTestStatistic kendall = Bivariate.KendallTauTest(x, y).UnderlyingStatistic; - if (kendall != null) - { + if (kendall != null) { kendallStatistics.Add(kendall.Value); kendallDistribution = kendall.Distribution; } } // Test whether statistics are actually distributed as claimed - if (spearmanDistribution != null) - { + if (spearmanDistribution != null) { TestResult spearmanChiSquared = spearmanStatistics.ChiSquaredTest(spearmanDistribution); Assert.IsTrue(spearmanChiSquared.Probability > 0.01); } - if (kendallDistribution != null) - { + if (kendallDistribution != null) { TestResult kendallChiSquared = kendallStatistics.ChiSquaredTest(kendallDistribution); Assert.IsTrue(kendallChiSquared.Probability > 0.01); } @@ -649,5 +642,17 @@ public void BivariateAssociationDiscreteNullDistribution () } } + [TestMethod] + public void BivariateDimensionMismatch () { + + double[] x = new double[3]; + double[] y = new double[4]; + + Assert.ThrowsException(() => x.CorrelationCoefficient(y)); + Assert.ThrowsException(() => x.Covariance(y)); + Assert.ThrowsException(() => y.LinearRegression(x)); + + } + } } diff --git a/Test/Bug53.txt b/Test/Bug53.txt new file mode 100644 index 0000000..73da64a --- /dev/null +++ b/Test/Bug53.txt @@ -0,0 +1,528 @@ +Length;Age +367;3 +398;3 +365;3 +470;4 +375;3 +393;3 +423;3 +365;3 +373;3 +393;4 +425;3 +460;4 +350;3 +365;3 +363;3 +465;5 +380;3 +445;4 +375;3 +408;4 +425;4 +390;4 +700;7 +633;6 +710;7 +680;7 +475;5 +375;3 +360;3 +363;3 +315;2 +545;6 +360;3 +450;4 +475;4 +430;4 +415;4 +442;4 +415;3 +427;4 +440;4 +410;4 +417;3 +397;4 +510;4 +630;6 +700;7 +435;3 +350;2 +625;6 +520;5 +415;4 +660;6 +575;6 +475;5 +595;6 +353;3 +387;3 +370;3 +495;4 +535;5 +415;4 +380;3 +365;3 +475;4 +382;3 +430;4 +420;3 +595;6 +475;4 +622;6 +455;5 +475;4 +385;3 +345;3 +580;6 +600;6 +530;5 +475;5 +590;6 +520;5 +565;5 +560;5 +455;4 +425;4 +495;5 +570;5 +555;5 +495;4 +550;5 +630;7 +510;4 +500;4 +580;6 +510;4 +660;7 +635;7 +650;7 +730;8 +480;4 +644;6 +564;5 +615;6 +650;6 +535;5 +380;4 +342;3 +350;3 +490;5 +333;3 +395;4 +420;4 +520;5 +550;5 +435;4 +410;3 +387;3 +445;4 +475;4 +355;3 +370;3 +363;3 +455;4 +405;3 +545;5 +440;4 +480;4 +430;4 +450;4 +480;4 +482;4 +507;4 +483;4 +510;4 +570;5 +410;3 +420;3 +520;5 +575;6 +590;6 +540;5 +590;6 +580;5 +580;6 +635;6 +633;6 +550;5 +430;3 +700;7 +365;3 +730;8 +750;8 +580;5 +380;3 +640;6 +515;5 +655;6 +620;7 +685;8 +680;8 +490;4 +615;7 +356;2 +657;7 +570;5 +600;7 +620;7 +590;6 +370;3 +380;3 +690;7 +335;3 +665;7 +630;8 +615;6 +615;6 +630;6 +625;6 +540;5 +510;4 +510;4 +640;6 +440;4 +440;4 +495;4 +415;3 +615;6 +515;4 +475;4 +600;6 +395;3 +490;4 +528;5 +545;5 +473;4 +450;4 +429;3 +460;4 +425;3 +351;2 +568;5 +530;5 +475;4 +385;3 +370;3 +425;3 +430;3 +477;4 +433;3 +485;4 +560;5 +540;5 +520;5 +475;5 +420;3 +530;5 +525;5 +359;3 +480;4 +350;2 +404;3 +435;3 +430;4 +310;3 +410;4 +370;3 +455;4 +420;4 +500;4 +490;4 +580;5 +475;4 +560;5 +470;4 +360;3 +445;4 +365;4 +375;3 +421;3 +375;3 +565;5 +478;4 +610;6 +620;6 +620;6 +465;5 +481;4 +570;5 +400;3 +550;6 +462;4 +520;5 +465;4 +580;5 +510;5 +640;6 +600;6 +670;7 +625;8 +600;6 +570;5 +635;6 +590;6 +565;5 +470;4 +402;4 +495;5 +490;5 +390;4 +430;3 +450;4 +620;6 +667;7 +450;4 +465;4 +570;5 +510;5 +505;5 +528;6 +475;4 +520;5 +485;4 +504;5 +460;4 +456;4 +510;5 +380;3 +510;4 +500;4 +460;4 +500;4 +594;6 +446;4 +457;4 +536;5 +380;3 +470;4 +460;4 +635;6 +460;5 +447;4 +500;4 +665;7 +445;5 +475;4 +510;4 +515;4 +545;5 +540;5 +540;5 +535;5 +371;3 +510;4 +500;4 +500;4 +630;6 +770;8 +475;4 +485;4 +500;4 +630;6 +536;5 +455;4 +610;6 +550;5 +440;4 +485;4 +510;4 +485;4 +438;3 +500;5 +480;4 +567;5 +730;8 +650;7 +490;4 +440;4 +655;6 +489;4 +630;6 +460;4 +845;8 +400;3 +540;5 +526;5 +490;5 +450;4 +550;5 +426;3 +430;4 +575;5 +662;7 +395;4 +474;4 +400;3 +450;4 +397;3 +527;5 +621;6 +530;5 +340;3 +650;6 +437;4 +405;3 +510;4 +420;4 +645;7 +518;5 +500;4 +372;4 +710;8 +523;5 +495;4 +495;5 +450;4 +485;4 +375;3 +493;5 +443;4 +1065;16 +580;6 +415;4 +820;9 +485;4 +618;6 +650;6 +540;5 +795;7 +530;5 +540;5 +452;4 +699;6 +550;5 +462;4 +480;4 +424;3 +650;6 +412;4 +420;3 +525;6 +465;5 +550;5 +680;7 +590;6 +455;4 +480;4 +700;7 +705;7 +660;7 +515;5 +570;5 +540;6 +370;4 +455;4 +400;3 +421;3 +480;4 +487;4 +700;7 +663;7 +470;4 +490;4 +475;4 +573;5 +470;4 +565;6 +665;7 +440;4 +560;5 +495;5 +495;4 +445;4 +487;4 +457;4 +570;6 +650;7 +360;3 +420;4 +495;5 +560;6 +517;5 +435;4 +611;6 +570;5 +425;3 +400;3 +775;9 +760;8 +775;8 +357;3 +443;4 +465;5 +425;4 +570;6 +611;6 +465;4 +780;8 +540;5 +493;4 +515;5 +483;5 +361;3 +363;3 +890;11 +680;7 +390;3 +373;5 +970;12 +1005;11 +1010;12 +415;4 +420;4 +600;6 +375;4 +490;4 +525;5 +550;5 +675;7 +675;7 +640;6 +670;8 +1010;15 +670;7 +510;4 +595;6 +870;10 +530;5 +491;4 +470;5 +490;4 +775;8 +360;3 +480;5 +379;3 +327;2 +370;3 +465;4 +610;6 +755;8 +735;8 +730;8 +500;4 +560;6 +560;6 +450;4 +540;5 +380;3 +430;4 +460;4 +485;4 +618;6 +550;5 +415;4 +590;6 +540;5 +465;4 +355;3 +385;3 +585;5 +380;3 +550;5 +475;4 +470;4 +684;6 +380;3 +655;6 +623;6 diff --git a/Test/Bugs.cs b/Test/Bugs.cs index 1e76026..855b778 100644 --- a/Test/Bugs.cs +++ b/Test/Bugs.cs @@ -15,6 +15,45 @@ namespace Test { [TestClass] public class BugTests { + [TestMethod] + public void Bug56 () { + + RectangularMatrix M = new RectangularMatrix(new double[,] { + { 44.6667, -392.0000, -66.0000 }, + { -392.0000, 3488.0000, 504.0001 }, + { -66.0000, 504.0001, 216.0001 } + }); + //M = M.Transpose; + SingularValueDecomposition S = M.SingularValueDecomposition(); + + } + + //[TestMethod] + [DeploymentItem("Bug53.txt")] + public void Bug53 () { + // User sent data for which non-linear fit failed with Nonconvergence exception. + List x = new List(); + List y = new List(); + using (System.IO.StreamReader reader = System.IO.File.OpenText(@"Bug53.txt")) { + reader.ReadLine(); + while (!reader.EndOfStream) { + string line = reader.ReadLine(); + string[] values = line.Split(';'); + y.Add(Double.Parse(values[0])); + x.Add(Double.Parse(values[1])); + } + } + NonlinearRegressionResult r = y.NonlinearRegression( + x, + (IReadOnlyDictionary p, double u) => { + return p["0"] * (1.0 - Math.Exp(-p["1"] * (u - p["2"]))); + }, + //new Dictionary() { { "0", 100.0 }, { "1", 1.0 }, { "2", 1.0 } } + new Dictionary() { {"0", 1.0 }, {"1", 1.0}, {"2", 0.0} } + ); + + } + [TestMethod] [DeploymentItem("Bug8021.csv")] public void Bug8021 () { diff --git a/Test/DiscreteDistributionTest.cs b/Test/DiscreteDistributionTest.cs index c12e454..ff9a65c 100644 --- a/Test/DiscreteDistributionTest.cs +++ b/Test/DiscreteDistributionTest.cs @@ -22,11 +22,12 @@ public static DiscreteDistribution[] GetDistributions () { return (new DiscreteDistribution[] { new BernoulliDistribution(0.1), new BinomialDistribution(0.2, 30), new BinomialDistribution(0.4, 5), - new PoissonDistribution(4.5), new PoissonDistribution(400.0), + new PoissonDistribution(0.54), new PoissonDistribution(5.4), new PoissonDistribution(540.0), new DiscreteUniformDistribution(5, 11), new GeometricDistribution(0.6), new NegativeBinomialDistribution(7.8, 0.4), - new HypergeometricDistribution(9, 3, 5) + new HypergeometricDistribution(9, 3, 5), + new DiscreteTestDistribution() }); } @@ -58,6 +59,42 @@ public void DiscreteDistributionVariance () { } } + [TestMethod] + public void DiscreteDistributionSupport () { + + foreach (DiscreteDistribution distribution in distributions) { + + // To left of support + if (distribution.Support.LeftEndpoint > Int32.MinValue) { + int k0 = distribution.Support.LeftEndpoint - 1; + Assert.IsTrue(distribution.LeftInclusiveProbability(k0) == 0.0); + Assert.IsTrue(distribution.LeftExclusiveProbability(k0) == 0.0); + Assert.IsTrue(distribution.RightExclusiveProbability(k0) == 1.0); + Assert.IsTrue(distribution.ProbabilityMass(k0) == 0.0); + } + + // At left end of support + int k1 = distribution.Support.LeftEndpoint; + Assert.IsTrue(distribution.LeftInclusiveProbability(k1) == distribution.ProbabilityMass(k1)); + Assert.IsTrue(distribution.LeftExclusiveProbability(k1) == 0.0); + Assert.IsTrue(distribution.InverseLeftProbability(0.0) == k1); + + // At right end of support + int k2 = distribution.Support.RightEndpoint; + Assert.IsTrue(distribution.LeftInclusiveProbability(k2) == 1.0); + //Assert.IsTrue(distribution.RightExclusiveProbability(k2) == 0.0); + + // To right of support + if (distribution.Support.RightEndpoint < Int32.MaxValue) { + int k3 = distribution.Support.RightEndpoint + 1; + Assert.IsTrue(distribution.LeftInclusiveProbability(k3) == 1.0); + Assert.IsTrue(distribution.LeftExclusiveProbability(k3) == 1.0); + Assert.IsTrue(distribution.RightExclusiveProbability(k3) == 0.0); + } + + } + } + [TestMethod] public void DiscreteDistributionProbabilityAxioms () { @@ -66,19 +103,16 @@ public void DiscreteDistributionProbabilityAxioms () { // some of these values will be outside the support, but that's fine, our results should still be consistent with probability axioms foreach (int k in TestUtilities.GenerateUniformIntegerValues(-10, +100, 8)) { - Console.WriteLine("{0} {1}", distribution.GetType().Name, k); - double DP = distribution.ProbabilityMass(k); - Assert.IsTrue(DP >= 0.0); Assert.IsTrue(DP <= 1.0); + Assert.IsTrue(0.0 <= DP && DP <= 1.0); double P = distribution.LeftInclusiveProbability(k); + Assert.IsTrue(0.0 <= P && P <= 1.0); + double Q = distribution.RightExclusiveProbability(k); - Console.WriteLine("{0} {1} {2}", P, Q, P + Q); + Assert.IsTrue(0.0 <= Q && Q <= 1.0); - Assert.IsTrue(P >= 0.0); Assert.IsTrue(P <= 1.0); - Assert.IsTrue(Q >= 0.0); Assert.IsTrue(Q <= 1.0); Assert.IsTrue(TestUtilities.IsNearlyEqual(P + Q, 1.0)); - } } @@ -95,8 +129,6 @@ public void DiscreteDistributionInverseCDF () { foreach (DiscreteDistribution distribution in distributions) { int k = distribution.InverseLeftProbability(P); - Console.WriteLine("{0} P={1} K={2} P(k sample = new List(); @@ -260,6 +293,9 @@ public void BernoulliDistributionFit () BernoulliFitResult fit = BernoulliDistribution.FitToSample(n0, n1); Assert.IsTrue(fit.P.ConfidenceInterval(0.95).ClosedContains(p)); Assert.IsTrue(fit.GoodnessOfFit.Probability > 0.05); + + Assert.IsTrue(fit.Parameters.Count == 1); + Assert.IsTrue(fit.Parameters[0].Estimate == fit.P); } [TestMethod] diff --git a/Test/DistributionTest.cs b/Test/DistributionTest.cs index 8cf9a31..36a6ad8 100644 --- a/Test/DistributionTest.cs +++ b/Test/DistributionTest.cs @@ -13,7 +13,7 @@ using Meta.Numerics.Statistics.Distributions; namespace Test { - + [TestClass] public class DistributionTest { @@ -23,7 +23,7 @@ public class DistributionTest { private static List CreateDistributions () { - List distributions = new List( new ContinuousDistribution[] { + List distributions = new List(new ContinuousDistribution[] { new NoncentralChiSquaredDistribution(2, 3.0), new CauchyDistribution(1.0, 2.0), new UniformDistribution(Interval.FromEndpoints(-2.0,1.0)), new UniformDistribution(Interval.FromEndpoints(7.0, 9.0)), @@ -131,7 +131,7 @@ public void DistributionCentralInequality () { public void DistributionMonotonicity () { foreach (ContinuousDistribution distribution in distributions) { for (int i = 0; i < (probabilities.Length - 1); i++) { - Assert.IsTrue(distribution.InverseLeftProbability(probabilities[i]) < distribution.InverseLeftProbability(probabilities[i+1])); + Assert.IsTrue(distribution.InverseLeftProbability(probabilities[i]) < distribution.InverseLeftProbability(probabilities[i + 1])); Assert.IsTrue(distribution.InverseRightProbability(probabilities[i]) > distribution.InverseRightProbability(probabilities[i + 1])); } } @@ -178,7 +178,7 @@ public void DistributionMeanIntegral () { double mean = distribution.Mean; if (Double.IsNaN(mean) || Double.IsInfinity(mean)) continue; if (distribution is BetaDistribution b && ((b.Alpha < 1.0) || (b.Beta < 1.0))) continue; - Func f = delegate(double x) { + Func f = delegate (double x) { return (distribution.ProbabilityDensity(x) * x); }; double M1 = FunctionMath.Integrate(f, distribution.Support); @@ -203,7 +203,7 @@ public void DistributionVarianceIntegral () { double e = TestUtilities.TargetPrecision; GammaDistribution gammaDistribution = distribution as GammaDistribution; if ((gammaDistribution != null) && (gammaDistribution.Shape < 1.0)) e = Math.Sqrt(e); - Func f = delegate(double x) { + Func f = delegate (double x) { double z = x - distribution.Mean; return (distribution.ProbabilityDensity(x) * z * z); }; @@ -220,7 +220,7 @@ public void DistributionRawMomentIntegral () { foreach (int r in TestUtilities.GenerateIntegerValues(2, 32, 8)) { double M = distribution.RawMoment(r); if (Double.IsInfinity(M) || Double.IsNaN(M)) continue; // don't try to do a non-convergent integral - Func f = delegate(double x) { + Func f = delegate (double x) { return (distribution.ProbabilityDensity(x) * Math.Pow(x, r)); }; try { @@ -268,7 +268,7 @@ public void DistributionCentralMomentIntegral () { // do the integral double m = distribution.Mean; - Func f = delegate(double x) { + Func f = delegate (double x) { return (distribution.ProbabilityDensity(x) * MoreMath.Pow(x - m, n)); }; try { @@ -313,7 +313,7 @@ public void DistributionProbabilityIntegral () { if (Double.IsNegativeInfinity(distribution.Support.LeftEndpoint) && Double.IsPositiveInfinity(distribution.Support.RightEndpoint)) { // pick an exponentially distributed random point with a random sign double y = rng.NextDouble(); - x = - Math.Log(y); + x = -Math.Log(y); if (rng.NextDouble() < 0.5) x = -x; } else if (Double.IsPositiveInfinity(distribution.Support.RightEndpoint)) { // pick an exponentially distributed random point @@ -418,7 +418,7 @@ public void FisherFromChiSquared () { for (int i = 0; i < 250; i++) { double x1 = d1.GetRandomValue(rng); double x2 = d2.GetRandomValue(rng); - double x = (x1/n1) / (x2/n2); + double x = (x1 / n1) / (x2 / n2); s.Add(x); } @@ -525,7 +525,7 @@ public void GammaFromExponential () { public void InverseGaussianSummation () { // X_i ~ IG(\mu,\lambda) \rightarrow \sum_{i=0}^{n} X_i ~ IG(n \mu, n^2 \lambda) - + Random rng = new Random(1); WaldDistribution d0 = new WaldDistribution(1.0, 2.0); List s = new List(); @@ -552,7 +552,7 @@ public void CauchyStudentAgreement () { Assert.IsTrue(TestUtilities.IsNearlyEqual(xS, xC)); Assert.IsTrue(TestUtilities.IsNearlyEqual(S.ProbabilityDensity(xS), C.ProbabilityDensity(xC))); } - + } [TestMethod] @@ -564,6 +564,20 @@ public void CauchyFWHM () { Assert.IsTrue(TestUtilities.IsNearlyEqual(D.ProbabilityDensity(D.Median + D.FullWithAtHalfMaximum / 2.0), p / 2.0)); } + [TestMethod] + public void ChiAndChiSquared () { + + ChiDistribution chi = new ChiDistribution(3); + ChiSquaredDistribution chi2 = new ChiSquaredDistribution(3.0); + + Assert.IsTrue(chi.DegreesOfFreedom == chi2.DegreesOfFreedom); + + foreach (double x in TestUtilities.GenerateRealValues(1.0E-1, 1.0E1, 4)) { + Assert.IsTrue(TestUtilities.IsNearlyEqual(chi.LeftProbability(x), chi2.LeftProbability(x * x))); + } + + } + [TestMethod] public void DistributionRandomDeviates () { // Check that random deviates generated match distribution diff --git a/Test/FutureTest.cs b/Test/FutureTest.cs index 08b04dd..3428153 100644 --- a/Test/FutureTest.cs +++ b/Test/FutureTest.cs @@ -24,267 +24,7 @@ internal interface IDeviateGenerator { } - public class PoissonGeneratorMultiplicative : IDeviateGenerator { - - public PoissonGeneratorMultiplicative (double lambda) { - expMinusLambda = Math.Exp(-lambda); - } - - private readonly double expMinusLambda; - - public int GetNext (Random rng) { - int k = 0; - double t = rng.NextDouble(); - while (t > expMinusLambda) { - k++; - t *= rng.NextDouble(); - } - return (k); - } - - } - - public class PoissonGeneratorAdditive : IDeviateGenerator { - - public PoissonGeneratorAdditive (double lambda) { - this.lambda = lambda; - } - - private readonly double lambda; - - public int GetNext (Random rng) { - int k = 0; - double t = rng.NextDouble(); - while (t < lambda) { - k++; - t += rng.NextDouble(); - } - return (k); - } - } - - public class PoissonGeneratorTabulated : IDeviateGenerator { - - public PoissonGeneratorTabulated (double lambda) { - cdf = new double[(int) Math.Ceiling(4.0 * lambda)]; - - double P = Math.Exp(-lambda); - cdf[0] = P; - for (int k = 1; k < cdf.Length; k++) { - P *= lambda / k; - cdf[k] = cdf[k - 1] + P; - } - } - - private readonly double[] cdf; - - public int GetNext (Random rng) { - - double u = rng.NextDouble(); - - int k = Array.BinarySearch(cdf, u); - if (k < 0) k = ~k; - - return (k); - } - - } - - public class PoissonGeneratorNR : IDeviateGenerator { - - public PoissonGeneratorNR (double lambda) { - this.lambda = lambda; - this.sqlam = Math.Sqrt(lambda); - this.loglam = Math.Log(lambda); - } - - private readonly double lambda, sqlam, loglam; - - public int GetNext (Random rng) { - while (true) { - double u = 0.64 * rng.NextDouble(); - double v = -0.68 + 1.28 * rng.NextDouble(); - int k = (int) Math.Floor(sqlam * (v / u) + lambda + 0.5); - if (k < 0) continue; - double u2 = u * u; - double lfac = AdvancedIntegerMath.LogFactorial(k); - double p = sqlam * Math.Exp(-lambda + k * loglam - lfac); - if (u2 < p) return (k); - } - } - } - - public class PoissonGeneratorPTRS : IDeviateGenerator { - - public PoissonGeneratorPTRS (double lambda) { - this.lambda = lambda; - this.lnmu = Math.Log(lambda); - this.b = 0.931 + 2.53 * Math.Sqrt(lambda); - this.a = -0.059 + 0.02483 * b; - this.vr = 0.9277 - 3.6224 / (b - 2.0); - } - - private readonly double lambda, lnmu, b, a, vr; - - public int GetNext (Random rng) { - - while (true) { - double u = rng.NextDouble() - 0.5; - double v = rng.NextDouble(); - - double us = 0.5 - Math.Abs(u); - int k = (int) Math.Floor((2.0 * a / us + b) * u + lambda + 0.43); - if (us > 0.07 && v < vr) return (k); - if (k < 0) continue; - if (us < 0.013 && v > us) continue; - double ai = 1.1239 + 1.1328 / (b - 3.4); - if (Math.Log(v * ai / (a / (us * us) + b)) <= -lambda + k * lnmu - AdvancedIntegerMath.LogFactorial(k)) { - return (k); - } - } - - } - } - - public static class DoubleDoubleGamma { - - public static readonly DoubleDouble[] Bernoulli = new DoubleDouble[] { - DoubleDouble.One, DoubleDouble.One / 6, -DoubleDouble.One / 30, DoubleDouble.One / 42, -DoubleDouble.One / 30, - ((DoubleDouble) 5) / 66, -((DoubleDouble) 691) / 2730, ((DoubleDouble) 7) / 6, -((DoubleDouble) 3617) / 510, ((DoubleDouble) 43867) / 798, - -((DoubleDouble) 74611) / 330, ((DoubleDouble) 854513) / 138, -((DoubleDouble) 236364091) / 2730, ((DoubleDouble) 8553103) / 6, -((DoubleDouble) 23749461029) / 870 - }; - - public static DoubleDouble Sum (DoubleDouble x) { - DoubleDouble rxPower = 1.0 / x; - DoubleDouble rxSquared = rxPower * rxPower; - DoubleDouble f = 0.5 * Bernoulli[1] * rxPower; - for (int k = 2; k < Bernoulli.Length; k++) { - DoubleDouble f_old = f; - rxPower *= rxSquared; - f += Bernoulli[k] / ((2 * k) * (2 * k - 1)) * rxPower; - if (f == f_old) { - return (f); - } - } - throw new NonconvergenceException(); - } - - public static DoubleDouble LogGamma_Asymptotic (DoubleDouble x) { - // Sum from smallest to largest terms to minimize error. - return (Sum(x) + halfLogTwoPi - x + (x - 0.5) * DoubleDouble.Log(x)); - } - - public static DoubleDouble LogGammaFromAsymptotic (DoubleDouble x) { - - Debug.Assert(x > 0.0); - - DoubleDouble s = DoubleDouble.Zero; - while (x < 34.0) { - s += DoubleDouble.Log(x); - x += DoubleDouble.One; - } - - return (LogGamma_Asymptotic(x) - s); - - } - - public static DoubleDouble LogGammaFromSeries (DoubleDouble x) { - - Debug.Assert(x >= 1.5); - - DoubleDouble s = DoubleDouble.Zero; - while (x > 2.5) { - x -= DoubleDouble.One; - s += DoubleDouble.Log(x); - } - - DoubleDouble y = x - 2.0; - Debug.Assert(DoubleDouble.Abs(y) <= 0.5); - return (LogGammaTwoPlus(y) + s); - } - - public static DoubleDouble ZetaMinusOne (int n) { - // For n < 16, needs more than 255 terms. Look into using - // Euler-Maclauren to accelerate. - DoubleDouble s = DoubleDouble.Zero; - for (int k = 2; k < 255; k++) { - DoubleDouble s_old = s; - s += DoubleDouble.Pow(k, -n); - if (s == s_old) { - return s; - } - } - throw new NonconvergenceException(); - } - - public static DoubleDouble LambdaMinusOne (int n) { - DoubleDouble s = DoubleDouble.Zero; - for (int k = 3; k < 512; k += 2) { - DoubleDouble s_old = s; - s += DoubleDouble.Pow(k, -n); - if (s == s_old) { - return s; - } - } - throw new NonconvergenceException(); - } - - private static DoubleDouble[] InitializeZetaMinusOne () { - DoubleDouble[] zetaMinusOne = new DoubleDouble[64]; - zetaMinusOne[0] = -1.5; - zetaMinusOne[1] = Double.PositiveInfinity; - zetaMinusOne[2] = new DoubleDouble("0.64493406684822643647241516664602519"); - zetaMinusOne[3] = new DoubleDouble("0.20205690315959428539973816151144999"); - zetaMinusOne[4] = new DoubleDouble("0.082323233711138191516003696541167903"); - zetaMinusOne[5] = new DoubleDouble("0.036927755143369926331365486457034168"); - zetaMinusOne[6] = new DoubleDouble("0.017343061984449139714517929790920528"); - zetaMinusOne[7] = new DoubleDouble("8.3492773819228268397975498497967596E-3"); - zetaMinusOne[8] = new DoubleDouble("4.0773561979443393786852385086524653E-3"); - zetaMinusOne[9] = new DoubleDouble("2.0083928260822144178527692324120605E-3"); - zetaMinusOne[10] = new DoubleDouble("9.9457512781808533714595890031901701E-4"); - zetaMinusOne[11] = new DoubleDouble("4.9418860411946455870228252646993647E-4"); - zetaMinusOne[12] = new DoubleDouble("2.4608655330804829863799804773967096E-4"); - zetaMinusOne[13] = new DoubleDouble("1.2271334757848914675183652635739571E-4"); - zetaMinusOne[14] = new DoubleDouble("6.1248135058704829258545105135333747E-5"); - zetaMinusOne[15] = new DoubleDouble("3.0588236307020493551728510645062588E-5"); - return zetaMinusOne; - } - - private static readonly DoubleDouble[] zetaMinusOne = InitializeZetaMinusOne(); - - private static DoubleDouble ZetaSeries (DoubleDouble x) { - DoubleDouble s = 0.0; - DoubleDouble xMinus = -x; - DoubleDouble xPower = xMinus; - for (int k = 2; k < zetaMinusOne.Length; k++) { - DoubleDouble s_old = s; - xPower *= xMinus; - // If not yet computed, compute next \zeta - 1 value. - if (zetaMinusOne[k] == DoubleDouble.Zero) { - // Technically this is not thread-safe, because assignment is not atomic for non-native structs. - // But at worst we are filling in the same value from two different threads, so this would only - // be a problem if intermediate values are neither start nor end values in some circumstances. - zetaMinusOne[k] = ZetaMinusOne(k); - } - s += zetaMinusOne[k] * xPower / k; - if (s == s_old) { - return (s); - } - } - throw new NonconvergenceException(); - } - - public static DoubleDouble LogGammaTwoPlus (DoubleDouble x) { - return (DoubleDouble.One - AdvancedDoubleDoubleMath.EulerGamma) * x + ZetaSeries(x); - } - - private static readonly DoubleDouble halfLogTwoPi = 0.5 * DoubleDouble.Log(2.0 * DoubleDouble.Pi); - - - - } - - + [TestClass] public class FutureTest { @@ -338,105 +78,7 @@ public void ErfMaxError () { } - [TestMethod] - public void ChengBeta () { - - Random rng = new Random(271828); - - ChengBetaGenerator g = new ChengBetaGenerator(3.0, 1.1); - List sample = new List(); - for (int i = 0; i < 1000000; i++) sample.Add(g.GetNext(rng)); - - BetaDistribution b = new BetaDistribution(3.0, 1.1); - TestResult r = sample.KolmogorovSmirnovTest(b); - - } - - [TestMethod] - public void PoissonBenchmark () { - - int n = 1000000; - double lambda = 10.0; - int m = 100; - - PoissonDistribution d = new PoissonDistribution(lambda); - Random rng; - - // Inversion - Histogram h0 = new Histogram(m); - rng = new Random(1); - Stopwatch s0 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - int k = d.GetRandomValue(rng); - h0.Add(k); - } - s0.Stop(); - TestResult r0 = h0.ChiSquaredTest(d); - - - // Multiplication - IDeviateGenerator g1 = new PoissonGeneratorMultiplicative(lambda); - Histogram h1 = new Histogram(m); - rng = new Random(1); - Stopwatch s1 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - int k = g1.GetNext(rng); - h1.Add(k); - } - s1.Stop(); - TestResult r1 = h1.ChiSquaredTest(d); - - // Addition - IDeviateGenerator g2 = new PoissonGeneratorAdditive(lambda); - Histogram h2 = new Histogram(m); - rng = new Random(1); - Stopwatch s2 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - int k = g2.GetNext(rng); - h2.Add(k); - } - s2.Stop(); - TestResult r2 = h2.ChiSquaredTest(d); - - // Tabulation - IDeviateGenerator gt = new PoissonGeneratorTabulated(lambda); - Histogram ht = new Histogram(m); - rng = new Random(1); - Stopwatch st = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - int k = gt.GetNext(rng); - ht.Add(k); - } - st.Stop(); - TestResult rt = ht.ChiSquaredTest(d); - - - // NR - IDeviateGenerator g3 = new PoissonGeneratorNR(lambda); - Histogram h3 = new Histogram(m); - rng = new Random(1); - Stopwatch s3 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - int k = g3.GetNext(rng); - h3.Add(k); - } - s3.Stop(); - TestResult r3 = h3.ChiSquaredTest(d); - - // PTRS - IDeviateGenerator g4 = new PoissonGeneratorPTRS(lambda); - Histogram h4 = new Histogram(m); - rng = new Random(1); - Stopwatch s4 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - int k = g4.GetNext(rng); - h4.Add(k); - } - s4.Stop(); - TestResult r4 = h4.ChiSquaredTest(d); - - } - + [TestMethod] public void EigenExample () { @@ -554,136 +196,6 @@ public void SkewnessVariance () { } - [TestMethod] - public void PoissonLargeMeanDeviates () { - - int n = 100000; - List sample = new List(n); - - PoissonDistribution d = new PoissonDistribution(100); - Random rng = new Random(314159265); - for (int i = 0; i < n; i++) { - sample.Add(d.GetRandomValue(rng)); - } - - Histogram h = new Histogram(200); - foreach (int k in sample) h.Add(k); - TestResult ht = h.ChiSquaredTest(d); - - TestResult ct = sample.ChiSquaredTest(d); - - } - - private int PoissonDeviateMultiplication (double expMinusLambda, Random rng) { - double t0 = expMinusLambda; - int k = 0; - double t = rng.NextDouble(); - while (t > t0) { - k++; - t *= rng.NextDouble(); - } - return (k); - } - - private int PoissonDeviateAddition (double lambda, Random rng) { - int k = 0; - double t = rng.NextDouble(); - while (t < lambda) { - k++; - t += rng.NextDouble(); - } - return (k); - } - - private int PoissonDeviateNR (double lambda, Random rng) { - - double sqlam = Math.Sqrt(lambda); - double loglam = Math.Log(lambda); - - while (true) { - double u = 0.64 * rng.NextDouble(); - double v = -0.68 + 1.28 * rng.NextDouble(); - int k = (int) Math.Floor(sqlam * (v / u) + lambda + 0.5); - if (k < 0) continue; - double u2 = u * u; - double lfac = AdvancedIntegerMath.LogFactorial(k); - double p = sqlam * Math.Exp(-lambda + k * loglam - lfac); - if (u2 < p) return (k); - } - - } - - private int PoissonDeviatePTRS (double lambda, Random rng) { - - double lnmu = Math.Log(lambda); - double b = 0.931 + 2.53 * Math.Sqrt(lambda); - double a = -0.059 + 0.02483 * b; - double vr = 0.9277 - 3.6224 / (b - 2.0); - - while (true) { - double u = rng.NextDouble() - 0.5; - double v = rng.NextDouble(); - - double us = 0.5 - Math.Abs(u); - int k = (int) Math.Floor((2.0 * a / us + b) * u + lambda + 0.43); - if (us > 0.07 && v < vr) return (k); - if (k < 0) continue; - if (us < 0.013 && v > us) continue; - double ai = 1.1239 + 1.1328 / (b - 3.4); - if (Math.Log(v * ai / (a / (us * us) + b)) <= -lambda + k * lnmu - AdvancedIntegerMath.LogFactorial(k)) { - return (k); - } - } - - } - - private int PoissonDeviatePRUA (double lambda, Random rng) { - double c = 1.715527770; // 2.0 * Math.Sqrt(2.0 / Math.E); // 1.7155 - double d = 0.898916162; // 3.0 - 2.0 * Math.Sqrt(3.0 / Math.E); // 0.8989 - double a = lambda + 0.5; - double r = Math.Sqrt(a); - double h = c * r + d; - double logLambda = Math.Log(lambda); - - while (true) { - double u = rng.NextDouble(); - double v = rng.NextDouble(); - double x = a + h * (v - 0.5) / u; - if (x < 0.0) continue; - // test x >= b - int k = (int) Math.Floor(x); - // fast accept - // fast reject - double t = k * logLambda - AdvancedIntegerMath.LogFactorial(k) - lambda; - if (u * u < Math.Exp(t)) return (k); - } - - } - - [TestMethod] - public void ComparePoissonGenerators () { - - int n = 10000000; - double lambda = 15.0; - Random rng = new Random(1); - - PoissonDistribution p = new PoissonDistribution(lambda); - - int s0 = 0; - Stopwatch t0 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - s0 += p.GetRandomValue(rng); - } - t0.Stop(); - - int s1 = 0; - Stopwatch t1 = Stopwatch.StartNew(); - for (int i = 0; i < n; i++) { - s1 += PoissonDeviateNR(lambda, rng); - } - t1.Stop(); - } - [TestMethod] public void NonlinearRegressionLinearRegressionAgreementTest () { @@ -703,33 +215,6 @@ public void NonlinearRegressionLinearRegressionAgreementTest () { } - public void JacobiPeriod () { - - double k = 0.875; - - int n = 0; - double a = 1.0; - double b = Math.Sqrt(1.0 - k * k); - double c2 = a * a - b * b; - double s = c2; - while (a != b) { - - double ap = (a + b) / 2.0; - double bp = Math.Sqrt(a * b); - double c2p = MoreMath.Sqr(c2 / 4.0 / ap); - - n++; - a = ap; - b = bp; - c2 = c2p; - s += MoreMath.Pow(2.0, n) * c2; - - } - - // This works, but s~1 for k~1, so don't use to compute E for k~1. - - } - [TestMethod] public void TwoWayAnova () { @@ -3688,232 +3173,7 @@ private void FourthPowerWithDerivative (double x, out double f, out double fp) { fp = 4.0 * Math.Pow(x, 3); } - /* PRIME FACTORIZATION */ - - private void FermatFactor (int n) { - - int c = 0; - int s = (int)Math.Floor(Math.Sqrt(n)); - int x = 2 * s + 1; int y = 1; int r = s * s - n; - while (r != 0) { - r += x; x += 2; - c++; - do { - r -= y; y += 2; - c++; - } while (r > 0); - - } - int p = (x - y) / 2; - int q = (x + y - 2) / 2; - if (p == 1) { - //Console.WriteLine(q); - } else { - FermatFactor(p); - FermatFactor(q); - } - } - - private void PollardRhoFactor (int n) { - - int x = 5; int y = 2; int k = 1; int l = 1; - - for (int c = 0; c < 10000; c++) { - //while (true) { - int g = (int) AdvancedIntegerMath.GCF(Math.Abs(y - x), n); - if (g == n) { - Console.WriteLine("n={0}", n); - return; - } else if (g == 1) { - k--; - if (k == 0) { - y = x; - l = 2 * l; - k = l; - } - //Console.WriteLine("{0}^2 mod {1}", x, n); - x = AdvancedIntegerMath.PowMod(x, 2, n) + 1; - if (x == n) x = 0; - } else { - //Console.WriteLine("g={0}", g); - n = n / g; - x = x % n; - y = y % n; - } - } - - - - } - -#if FUTURE - [TestMethod] - public void TestFactor () { - /* - Stopwatch s1 = Stopwatch.StartNew(); - FermatFactor(1157625); - s1.Stop(); Console.WriteLine(s1.ElapsedMilliseconds); - - Stopwatch s2 = Stopwatch.StartNew(); - PollardRhoFactor(37); - s2.Stop(); Console.WriteLine(s2.ElapsedMilliseconds); - */ - // for (3*5*7)^3 = 1157625, Pollard's Rho fails to factor 125 = 5^3 - - int c = 0; - Stopwatch s1 = Stopwatch.StartNew(); - - //for (int i = 1; i < 1000000; i+=2) { - int i = 220211; - List factors = AdvancedIntegerMath.Factor(i); - //FermatFactor(i); - - int m = 1; - foreach (Factor factor in factors) { - Console.WriteLine(factor); - if (!AdvancedIntegerMath.IsPrime(factor.Value)) { - c++; - Console.WriteLine("for {0}, factor {1} is not prime", i, factor.Value); - } - //Console.WriteLine("{0} {1}", factor.Value, factor.Multiplicity); - m *= (int)MoreMath.Pow(factor.Value, factor.Multiplicity); - } - if (m != i) { - Console.WriteLine("for {0}, factors do not multiply to number", i); - } - //Assert.IsTrue(m == i); - //Console.WriteLine(m); - - // } - s1.Stop(); Console.WriteLine(s1.ElapsedMilliseconds); - - Console.WriteLine(c); - - } -#endif -#if FUTURE - [TestMethod] - public void BigAdd () { - - BigFloat a = new BigFloat(new byte[] { 2, 2, 5 }, 0); - Console.WriteLine(a); - BigFloat b = new BigFloat(new byte[] { 3, 3, 3 }, 0); - Console.WriteLine(b); - BigFloat c = a + b; - Console.WriteLine(c); - - } -#endif - - [TestMethod] - public void MultiRootTest () { - - // here is a system with a zeros at (3,4) and (1,0) - Func f = delegate (double[] u) { - double x = u[0]; double y = u[1]; - double a = 2.0 * y - x; - double b = (x * x + x * (y * y - 2.0) - 4.0 * y) / (x + 4.0); - double c = Math.Sqrt(x * x + y * y); - return (new double[] { b - a, c - a }); - }; - - FunctionMath.FindZero(f, new double[] { 1.0, 1.0 }); - - } - - /* INVERSE BETA */ - - public static double ApproximateInverseBetaSeries (double a, double b, double P) { - - double bigB = AdvancedMath.Beta(a, b); - - double z = Math.Pow(a * P * bigB, 1.0 / a); - - double z2 = (b - 1.0) / (a + 1.0) * z; - - double z3 = z2 * (a * a + 3.0 * b * a - a + 5.0 * b - 4.0) / (a + 1.0) / (a + 2.0) / 2.0 * z; - - double z4 = z2 * (a * a * a * a + (6.0 * b - 1.0) * a * a * a + (b + 2.0) * (8.0 * b - 5.0) * a * a + - (33.0 * b * b - 30.0 * b + 4.0) * a + b * (31.0 * b - 47.0) + 18.0) / MoreMath.Pow(a + 1.0, 2) / (a + 2.0) / (a + 3.0) / 3.0 * z * z; - - Console.WriteLine("z={0} z2={1} z3={2} z4={3}", z, z2, z3, z4); - - return (z * (1.0 + z2 + z3 + z4)); - - } - - public static double ApproximateInverseBeta (double a, double b, double P) { - - if (P > 0.5) { - return (1.0 - ApproximateInverseBeta(b, a, 1.0 - P)); - } else { - double bigB = AdvancedMath.Beta(a, b); - double z = Math.Pow(a * P * bigB, 1.0 / a); - double z2 = (b - 1.0) / (a + 1.0) * z; - if (z2 < 0.25) { - double z3 = z2 * (a * a + 3.0 * b * a - a + 5.0 * b - 4.0) / (a + 1.0) / (a + 2.0) / 2.0 * z; - return (z * (1.0 + z2 + z3)); - } else { - throw new NotImplementedException(); - } - } - - } - - public static double RefineInverseBeta (double a, double b, double P, double x) { - - for (int i = 0; i < 8; i++) { - double x_old = x; - double y = AdvancedMath.LeftRegularizedBeta(a, b, x) - P; - double yp = Math.Pow(x, a - 1.0) * Math.Pow(1.0 - x, b - 1.0) / AdvancedMath.Beta(a, b); - double dx = -y / yp; - x += dx; - if (x == x_old) return (x); - } - return (x); - } - - [TestMethod] - public void TestBeta () { - - double a = 200.0; double b = 200.0; double P = 1.0E-5; - double x1 = ApproximateInverseBetaSeries(a, b, P); - if ((0.0 < x1) && (x1 < 1.0)) { - Console.WriteLine("x1 {0} {1}", x1, AdvancedMath.LeftRegularizedBeta(a, b, x1)); - } - - double x2 = 1.0 - ApproximateInverseBetaSeries(b, a, 1.0 - P); - if ((0.0 < x2) && (x2 < 1.0)) { - Console.WriteLine("x2 {0} {1}", x2, AdvancedMath.LeftRegularizedBeta(a, b, x2)); - } - - //x1 = RefineInverseBeta(a, b, P, x1); - //Console.WriteLine("{0} {1}", x1, AdvancedMath.LeftRegularizedBeta(a, b, x1)); - - NormalDistribution N = new NormalDistribution(); - double m = a / (a + b); double s = Math.Sqrt(a * b / (a + b + 1.0)) / (a + b); - double x3 = m + s * N.InverseLeftProbability(P); - if ((0.0 < x3) && (x3 < 1.0)) { - Console.WriteLine("x3 {0} {1}", x3, AdvancedMath.LeftRegularizedBeta(a, b, x3)); - } - - //Console.WriteLine(AdvancedMath.Beta(a, b, 0.35) / AdvancedMath.Beta(a, b)); - //Console.WriteLine(AdvancedMath.Beta(a, b, 0.40) / AdvancedMath.Beta(a, b)); - //Console.WriteLine(AdvancedMath.Beta(a, b, 0.45) / AdvancedMath.Beta(a, b)); - - } - - private static void WriteMatrix (AnyRectangularMatrix A) { - for (int r = 0; r < A.RowCount; r++) { - for (int c = 0; c < A.ColumnCount; c++) { - Console.Write("{0} ", A[r, c]); - } - Console.WriteLine(); - } - } - } - #endif - } +} \ No newline at end of file diff --git a/Test/MultivariateSampleTest.cs b/Test/MultivariateSampleTest.cs index e19a211..6814dab 100644 --- a/Test/MultivariateSampleTest.cs +++ b/Test/MultivariateSampleTest.cs @@ -226,14 +226,18 @@ public void MultivariateLinearRegressionSimple () { Assert.IsTrue(newResult.Intercept.ConfidenceInterval(0.99).ClosedContains(a)); // The residuals should be compatible with the model predictions + double ssr = 0.0; for (int i = 0; i < table.Rows.Count; i++) { FrameRow row = table.Rows[i]; double x0 = (double) row["x0"]; double x1 = (double) row["x1"]; double yp = newResult.Predict(x0, x1).Value; double y = (double) row["y"]; - Assert.IsTrue(TestUtilities.IsNearlyEqual(newResult.Residuals[i], y - yp)); + double z = y - yp; + Assert.IsTrue(TestUtilities.IsNearlyEqual(newResult.Residuals[i], z)); + ssr += z * z; } + Assert.IsTrue(TestUtilities.IsNearlyEqual(newResult.SumOfSquaredResiduals, ssr)); } @@ -618,20 +622,29 @@ public void PrincipalComponentAnalysis () { Assert.IsTrue(pca.Count == sample.Count); // check that the PCs behave as expected + Assert.IsTrue(pca.Components.Count == pca.Dimension); for (int i = 0; i < pca.Dimension; i++) { PrincipalComponent pc = pca.Components[i]; Assert.IsTrue(pc.Index == i); Assert.IsTrue(pc.Analysis == pca); Assert.IsTrue(TestUtilities.IsNearlyEqual(pc.Weight * pc.NormalizedVector, pc.ScaledVector())); + Assert.IsTrue(pca.MinimumDimension(pc.CumulativeVarianceFraction) == i + 1); + } + + // Check enumerator, and verify that variance fractions behave as expected. + int count = 0; + double cumulative = 0.0; + double previous = Double.PositiveInfinity; + foreach (PrincipalComponent pc in pca.Components) { + Assert.IsTrue(pc.Index == count); + count++; Assert.IsTrue((0.0 <= pc.VarianceFraction) && (pc.VarianceFraction <= 1.0)); - if (i == 0) { - Assert.IsTrue(pc.VarianceFraction == pc.CumulativeVarianceFraction); - } else { - PrincipalComponent ppc = pca.Components[i - 1]; - Assert.IsTrue(pc.VarianceFraction <= ppc.VarianceFraction); - Assert.IsTrue(TestUtilities.IsNearlyEqual(ppc.CumulativeVarianceFraction + pc.VarianceFraction, pc.CumulativeVarianceFraction)); - } + Assert.IsTrue(pc.VarianceFraction <= previous); + previous = pc.VarianceFraction; + cumulative += pc.VarianceFraction; + Assert.IsTrue(TestUtilities.IsNearlyEqual(cumulative, pc.CumulativeVarianceFraction)); } + Assert.IsTrue(count == pca.Components.Count); // express the sample in terms of principal components MultivariateSample csample = pca.TransformedSample(); diff --git a/Test/SampleTest.cs b/Test/SampleTest.cs index a923347..62889e1 100644 --- a/Test/SampleTest.cs +++ b/Test/SampleTest.cs @@ -241,6 +241,19 @@ public void OneSampleProperties () { Assert.IsTrue(sample.Median() == sample[0]); } + [TestMethod] + public void SampleVarianceRelationships () { + double[] sample = new double[] { 1.0, 2.0, 4.0 }; + + double m = sample.Mean(); + double sumOfSqareDeviations = 0.0; + for (int i = 0; i < sample.Length; i++) sumOfSqareDeviations += MoreMath.Sqr(sample[i] - m); + + Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.Variance(), sumOfSqareDeviations / sample.Length)); + Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.StandardDeviation(), Math.Sqrt(sumOfSqareDeviations / sample.Length))); + Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.CorrectedStandardDeviation(), Math.Sqrt(sumOfSqareDeviations / (sample.Length - 1)))); + } + [TestMethod] public void LognormalFit () { diff --git a/Test/StreamingSampleSummaryTest.cs b/Test/StreamingSampleSummaryTest.cs index 3127f85..1185d45 100644 --- a/Test/StreamingSampleSummaryTest.cs +++ b/Test/StreamingSampleSummaryTest.cs @@ -26,6 +26,7 @@ public void StreamSampleSummaryAgreement () { Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.Variance(), summary.Variance)); Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.PopulationMean(), summary.PopulationMean)); Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.PopulationVariance(), summary.PopulationVariance)); + Assert.IsTrue(TestUtilities.IsNearlyEqual(sample.PopulationStandardDeviation(), summary.PopulationStandardDeviation)); Assert.IsTrue(sample.Minimum() == summary.Minimum); Assert.IsTrue(sample.Maximum() == summary.Maximum); diff --git a/Test/Test.csproj b/Test/Test.csproj index 8f2a9b2..0dab870 100644 --- a/Test/Test.csproj +++ b/Test/Test.csproj @@ -9,7 +9,7 @@ Properties Test Test - v4.5.1 + v4.6.1 512 {3AC096D0-A1C2-E12C-1390-A8335801FDAB};{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC} 15.0 @@ -19,6 +19,7 @@ UnitTest + true @@ -133,6 +134,9 @@ + + PreserveNewest +