Merge pull request #157 from dotnet/minka

Improved numerical accuracy of GaussianOp, DoubleIsBetweenOp, and MMath.NormalCdf2

Infer.sln

@@ -1,6 +1,6 @@
 Microsoft Visual Studio Solution File, Format Version 12.00
-# Visual Studio 15
-VisualStudioVersion = 15.0.27428.2005
+# Visual Studio Version 16
+VisualStudioVersion = 16.0.29020.237
 MinimumVisualStudioVersion = 10.0.40219.1
 Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{A181C943-2E01-454D-9008-2E3C53AA09CC}"
 	ProjectSection(SolutionItems) = preProject

src/Compiler/Infer/Transforms/VariableTransform.cs

@@ -431,7 +431,7 @@ namespace Microsoft.ML.Probabilistic.Compiler.Transforms
                 IExpression arrayExpr = imie.Arguments[0];
                 IVariableDeclaration ivd = Recognizer.GetVariableDeclaration(arrayExpr);
                 // restrict to IVariableReferenceExpression for simplicity
-                if (ivd != null && arrayExpr is IVariableReferenceExpression)
+                if (ivd != null && !context.InputAttributes.Has<PointEstimate>(ivd) && arrayExpr is IVariableReferenceExpression)
                 {
                     variablesExcludingVariableFactor.Add(ivd);
                 }

src/Learners/Recommender/MatchboxRecommenderInternal/NativeDataFormatMatchboxRecommender.cs

@@ -964,10 +964,13 @@ namespace Microsoft.ML.Probabilistic.Learners.MatchboxRecommenderInternal
                         communityTrainingAlgorithmPool.Release(currentBatchAlgorithm, instanceData.Ratings.Count);
                     });
             }
-            catch (AggregateException ae)
+            catch (AggregateException ex)
             {
-                // Propagate up the first in the list of exceptions
-                throw ae.InnerException;
+                // To make the Visual Studio debugger stop at the inner exception, check "Enable Just My Code" in Debug->Options.
+                // throw InnerException while preserving stack trace
+                // https://stackoverflow.com/questions/57383/in-c-how-can-i-rethrow-innerexception-without-losing-stack-trace
+                System.Runtime.ExceptionServices.ExceptionDispatchInfo.Capture(ex.InnerException).Throw();
+                throw;
             }
         }
 

src/Runtime/Core/Maths/ExtendedDouble.cs

@@ -0,0 +1,151 @@
+using System;
+using System.Collections.Generic;
+using System.Reflection;
+using System.Text;
+using Microsoft.ML.Probabilistic.Collections;
+using Microsoft.ML.Probabilistic.Utilities;
+
+namespace Microsoft.ML.Probabilistic.Core.Maths
+{
+    /// <summary>
+    /// Represents a number as Mantissa * exp(Exponent).
+    /// </summary>
+    public class ExtendedDouble
+    {
+        public readonly double Mantissa, Exponent;
+
+        public ExtendedDouble(double mantissa, double exponent)
+        {
+            this.Mantissa = mantissa;
+            this.Exponent = exponent;
+        }
+
+        public override bool Equals(object obj)
+        {
+            if (obj is ExtendedDouble that)
+            {
+                return (this.Mantissa == that.Mantissa) && (this.Exponent == that.Exponent);
+            }
+            else return false;
+        }
+
+        public override int GetHashCode()
+        {
+            return Hash.Combine(Mantissa.GetHashCode(), Exponent);
+        }
+
+        public override string ToString()
+        {
+            return $"{Mantissa:r}*exp({Exponent:r})";
+        }
+
+        public static ExtendedDouble Zero()
+        {
+            return new ExtendedDouble(0, 0);
+        }
+
+        public static ExtendedDouble PositiveInfinity()
+        {
+            return new ExtendedDouble(double.PositiveInfinity, 0);
+        }
+
+        public static ExtendedDouble NaN()
+        {
+            return new ExtendedDouble(double.NaN, 0);
+        }
+
+        public static ExtendedDouble FromDouble(double value)
+        {
+            return new ExtendedDouble(value, 0);
+        }
+
+        public double ToDouble()
+        {
+            return Mantissa * System.Math.Exp(Exponent);
+        }
+
+        public double Log()
+        {
+            return Exponent + System.Math.Log(Mantissa);
+        }
+
+        public ExtendedDouble Max(double minimum)
+        {
+            return new ExtendedDouble(System.Math.Max(minimum, Mantissa), Exponent);
+        }
+
+        public ExtendedDouble MultiplyExp(double logarithm)
+        {
+            return new ExtendedDouble(Mantissa, Exponent + logarithm);
+        }
+
+        public static ExtendedDouble operator*(ExtendedDouble x, ExtendedDouble y)
+        {
+            return new ExtendedDouble(x.Mantissa * y.Mantissa, x.Exponent + y.Exponent);
+        }
+
+        public static ExtendedDouble operator /(ExtendedDouble x, ExtendedDouble y)
+        {
+            return new ExtendedDouble(x.Mantissa / y.Mantissa, x.Exponent - y.Exponent);
+        }
+
+        public static ExtendedDouble operator *(ExtendedDouble x, double y)
+        {
+            return new ExtendedDouble(x.Mantissa * y, x.Exponent);
+        }
+
+        public static ExtendedDouble operator /(ExtendedDouble x, double y)
+        {
+            return new ExtendedDouble(x.Mantissa / y, x.Exponent);
+        }
+
+        public static ExtendedDouble operator+(ExtendedDouble x, ExtendedDouble y)
+        {
+            if (x.Mantissa == 0)
+            {
+                return y;
+            }
+            else if (y.Mantissa == 0)
+            {
+                return x;
+            }
+            else if (y.Exponent > x.Exponent)
+            {
+                if (double.IsInfinity(x.Mantissa)) throw new ArgumentOutOfRangeException(nameof(x), x, "x is infinite");
+                return new ExtendedDouble(x.Mantissa * System.Math.Exp(x.Exponent - y.Exponent) + y.Mantissa, y.Exponent);
+            }
+            else
+            {
+                if (double.IsInfinity(y.Mantissa)) throw new ArgumentOutOfRangeException(nameof(y), y, "y is infinite");
+                return new ExtendedDouble(x.Mantissa + y.Mantissa * System.Math.Exp(y.Exponent - x.Exponent), x.Exponent);
+            }
+        }
+
+        public static ExtendedDouble operator-(ExtendedDouble x)
+        {
+            return new ExtendedDouble(-x.Mantissa, x.Exponent);
+        }
+
+        public static ExtendedDouble operator -(ExtendedDouble x, ExtendedDouble y)
+        {
+            if (x.Mantissa == 0)
+            {
+                return -y;
+            }
+            else if (y.Mantissa == 0)
+            {
+                return x;
+            }
+            else if (y.Exponent > x.Exponent)
+            {
+                if (double.IsInfinity(x.Mantissa)) throw new ArgumentOutOfRangeException(nameof(x), x, "x is infinite");
+                return new ExtendedDouble(x.Mantissa * System.Math.Exp(x.Exponent - y.Exponent) - y.Mantissa, y.Exponent);
+            }
+            else
+            {
+                if (double.IsInfinity(y.Mantissa)) throw new ArgumentOutOfRangeException(nameof(y), y, "y is infinite");
+                return new ExtendedDouble(x.Mantissa - y.Mantissa * System.Math.Exp(y.Exponent - x.Exponent), x.Exponent);
+            }
+        }
+    }
+}

src/Runtime/Core/Maths/MeanVarianceAccumulator.cs

@@ -47,18 +47,17 @@ namespace Microsoft.ML.Probabilistic.Math
             // new_mean = mean + w/(count + w)*(x - mean)
             // new_var = count/(count + w)*(var + w/(count+w)*(x-mean)*(x-mean)')
             Count += weight;
-            if (Count == 0) return;
-            double diff;
-            if(x == Mean) diff = 0;  // avoid subtracting infinities
-            else diff = x - Mean;
             if (weight == Count)
             {
                 // avoid numerical overflow
-                Mean += diff;
+                Mean = x;
                 Variance = 0;
             }
-            else
+            else if(weight != 0)
             {
+                double diff;
+                if (x == Mean) diff = 0;  // avoid subtracting infinities
+                else diff = x - Mean;
                 double s = weight / Count;
                 if (double.IsInfinity(Mean))
                     Mean = s * x + (1 - s) * Mean;

src/Runtime/Core/Maths/Region.cs

@@ -43,7 +43,7 @@ namespace Microsoft.ML.Probabilistic.Math
             double logVolume = 0;
             for (int i = 0; i < Dimension; i++)
             {
-                logVolume += Math.Log(Upper[i] - Lower[i]);
+                logVolume += Math.Log(Math.Max(1e-10, Upper[i] - Lower[i]));
             }
             return logVolume;
         }

src/Runtime/Factors/GaussianOp.cs

@@ -238,6 +238,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             else if (sample.Precision == 0)
             {
                 // for large vx, Z =approx N(mx; mm, vx+vm+E[1/prec])
+                if (mean.Precision == 0) return mean;
                 double mm, mv;
                 mean.GetMeanAndVariance(out mm, out mv);
                 // NOTE: this error may happen because sample didn't receive any message yet under the schedule.
@@ -334,6 +335,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                 double[] nodes = new double[QuadratureNodeCount];
                 double[] logWeights = new double[nodes.Length];
                 Gamma precMarginal = precision * to_precision;
+                if (precMarginal.IsPointMass) return SampleAverageConditional(mean, precMarginal.Point);
                 precMarginal = GaussianOp_Laplace.Q(sample, mean, precision, precMarginal);
                 QuadratureNodesAndWeights(precMarginal, nodes, logWeights);
                 if (!to_precision.IsUniform())
@@ -362,14 +364,14 @@ namespace Microsoft.ML.Probabilistic.Factors
                         Assert.IsTrue(prec > 0);
                         if ((prec > lowerBound && prec < upperBound) != precisionIsBetween)
                             continue;
-                        double ivr = prec / (v*prec + 1);
+                        double ivr = prec / (v * prec + 1);
                         double dlogf = -m * ivr;
                         double ddlogf = dlogf * dlogf - ivr;
                         double lp = 0.5 * Math.Log(ivr) - 0.5 * m * m * ivr;
                         if (double.IsPositiveInfinity(lp))
                             throw new Exception("overflow");
                         if (!double.IsNegativeInfinity(lp) && shift == 0)
-                            shift = logWeights[i]+lp;
+                            shift = logWeights[i] + lp;
                         double f = Math.Exp(logWeights[i] + lp - shift);
                         Z += f;
                         sum1 += dlogf * f;
@@ -384,7 +386,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                         throw new Exception("Quadrature found zero mass");
                     }
                     double alpha = sum1 / Z;
-                    double beta = alpha*alpha -sum2 / Z;
+                    double beta = alpha * alpha - sum2 / Z;
                     return GaussianOp.GaussianFromAlphaBeta(sample, alpha, beta, GaussianOp.ForceProper);
                 }
                 else
@@ -415,7 +417,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                         }
                         double lp = Gaussian.GetLogProb(xm, mm, xv + mv + 1.0 / prec);
                         if (i == 0)
-                            shift = logWeights[i]+lp;
+                            shift = logWeights[i] + lp;
                         double f = Math.Exp(logWeights[i] + lp - shift);
                         est.Add(Gaussian.FromMeanAndVariance(newMean, newVar), f);
                     }
@@ -489,43 +491,6 @@ namespace Microsoft.ML.Probabilistic.Factors
             return SampleAverageConditional(mean, sample, precision, to_precision);
         }
 
-#if false
-    /// <summary>
-    /// EP message to 'precision'
-    /// </summary>
-    /// <param name="sample">Constant value for 'sample'.</param>
-    /// <param name="mean">Incoming message from 'mean'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
-    /// <param name="precision">Incoming message from 'precision'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
-    /// <returns>The outgoing EP message to the 'precision' argument</returns>
-    /// <remarks><para>
-    /// The outgoing message is a distribution matching the moments of 'precision' as the random arguments are varied.
-    /// The formula is <c>proj[p(precision) sum_(mean) p(mean) factor(sample,mean,precision)]/p(precision)</c>.
-    /// </para></remarks>
-    /// <exception cref="ImproperMessageException"><paramref name="mean"/> is not a proper distribution</exception>
-    /// <exception cref="ImproperMessageException"><paramref name="precision"/> is not a proper distribution</exception>
-    public static Gamma PrecisionAverageConditional(double sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision, Gamma to_precision)
-    {
-      return PrecisionAverageConditional(Gaussian.PointMass(sample), mean, precision, to_precision);
-    }
-    /// <summary>
-    /// EP message to 'precision'
-    /// </summary>
-    /// <param name="sample">Incoming message from 'sample'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
-    /// <param name="mean">Constant value for 'mean'.</param>
-    /// <param name="precision">Incoming message from 'precision'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
-    /// <returns>The outgoing EP message to the 'precision' argument</returns>
-    /// <remarks><para>
-    /// The outgoing message is a distribution matching the moments of 'precision' as the random arguments are varied.
-    /// The formula is <c>proj[p(precision) sum_(sample) p(sample) factor(sample,mean,precision)]/p(precision)</c>.
-    /// </para></remarks>
-    /// <exception cref="ImproperMessageException"><paramref name="sample"/> is not a proper distribution</exception>
-    /// <exception cref="ImproperMessageException"><paramref name="precision"/> is not a proper distribution</exception>
-    public static Gamma PrecisionAverageConditional([SkipIfUniform] Gaussian sample, double mean, [SkipIfUniform] Gamma precision, Gamma to_precision)
-    {
-      return PrecisionAverageConditional(sample, Gaussian.PointMass(mean), precision, to_precision);
-    }
-#endif
-
         public static Gamma PrecisionAverageConditional_slow([SkipIfUniform] Gaussian sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision)
         {
             return PrecisionAverageConditional(sample, mean, precision);
@@ -533,24 +498,32 @@ namespace Microsoft.ML.Probabilistic.Factors
 
         public static Gamma PrecisionAverageConditional_Point(double ym, double yv, double precision)
         {
+            // xdlogf = 0.5*(1/(yv*r+1) - r*ym*ym/(yv*r+1)^2) = 0.5*(yv*r+1-r*ym*ym)/(yv*r+1)^2
+            // x2ddlogf = (-yv*r-0.5+r*ym*ym)/(yv*r+1)^2 - r*ym*ym/(yv*r+1)^3
+            // xdlogf + x2ddlogf = (-0.5*yv*r + 0.5*r*ym*ym)/(yv*r+1)^2 - r*ym*ym/(yv*r+1)^3
+            // as r->0: -0.5*r*(yv + ym*ym)
+            if (precision == 0) return Gamma.FromShapeAndRate(1.5, 0.5 * (yv + ym * ym));
             // point mass case
             // f(r) = N(xm;mm, xv+mv+1/r)
             // log f(r) = -0.5*log(xv+mv+1/r) - 0.5*(xm-mm)^2/(xv+mv+1/r)
             // (log f)' = (-0.5/(yv + 1/r) + 0.5*ym^2/(yv+1/r)^2)*(-1/r^2)
             // (log f)'' = (-0.5/(yv + 1/r) + 0.5*ym^2/(yv+1/r)^2)*(2/r^3) + (0.5/(yv+1/r)^2 - ym^2/(yv+1/r)^3)*(1/r^4)
-            double v = 1 / precision;
-            double v2 = v * v;
-            double denom = 1 / (yv + v);
-            double ymdenom = ym * denom;
-            double ym2denom2 = ymdenom * ymdenom;
-            double dlogf = (-0.5 * denom + 0.5 * ym2denom2) * (-v2);
-            // This method is slightly more accurate.
+            // r (log f)' = 0.5/(yv*r + 1) - 0.5*r*ym^2/(yv*r+1)^2
+            // r^2 (log f)'' = -1/(yv*r + 1) + r*ym^2/(yv*r+1)^2) + 0.5/(yv*r+1)^2 - r*ym^2/(yv*r+1)^3
+            double vdenom = 1 / (yv * precision + 1);
+            double ymvdenom = ym * vdenom;
+            double ymvdenom2 = precision * ymvdenom * ymvdenom;
+            //dlogf = (-0.5 * denom + 0.5 * ym2denom2) * (-v2);
             //dlogf = 0.5 * (1 - ym * ymdenom) * denom * v2;
             //dlogf = 0.5 * (v - ym * ym/(yv*precision+1))/(yv*precision + 1);
-            //dlogf = 0.5 * (yv+v - ym * ym) / (yv * precision + 1) / (yv * precision + 1);
-            double ddlogf = dlogf * (-2 * v) + (0.5 * denom - ym2denom2) * denom * v2 * v2;
+            //double dlogf = 0.5 * (v * vdenom - ymvdenom2);
+            //double xdlogf = precision * dlogf;
+            double xdlogf = 0.5 * (vdenom - ymvdenom2);
+            //double ddlogf = dlogf * (-2 * v) + (0.5 * denom - ym2denom2) * denom * v2 * v2;
+            //double ddlogf = v * (-2 * dlogf + (0.5 * v * vdenom - ym*ym*vdenom*vdenom) * vdenom);
+            double x2ddlogf = -2 * xdlogf + (0.5 * vdenom - ymvdenom2) * vdenom;
             bool checkDerivatives = false;
-            if(checkDerivatives)
+            if (checkDerivatives)
             {
                 double delta = precision * 1e-6;
                 double logf = Gaussian.GetLogProb(0, ym, yv + 1 / precision);
@@ -561,12 +534,53 @@ namespace Microsoft.ML.Probabilistic.Factors
                 double ulp = MMath.Ulp(logf);
                 if (logfd - logf > ulp && logf - logfd2 > ulp)
                 {
+                    double v = 1 / precision;
+                    double dlogf = v * xdlogf;
+                    double ddlogf = v * v * x2ddlogf;
                     Console.WriteLine($"dlogf={dlogf} check={dlogf2} ddlogf={ddlogf} check={ddlogf2}");
                     if (Math.Abs(dlogf2 - dlogf) > 1e-4) throw new Exception();
                     if (Math.Abs(ddlogf2 - ddlogf) > 1e-4) throw new Exception();
                 }
             }
-            return Gamma.FromDerivatives(precision, dlogf, ddlogf, ForceProper);
+            return GammaFromDerivatives(precision, xdlogf, x2ddlogf, ForceProper);
+        }
+
+        /// <summary>
+        /// Construct a Gamma distribution whose pdf has the given derivatives at a point.
+        /// </summary>
+        /// <param name="x">Must be positive</param>
+        /// <param name="xdLogP">Desired derivative of log-density at x, times x</param>
+        /// <param name="x2ddLogP">Desired second derivative of log-density at x, times x*x</param>
+        /// <param name="forceProper">If true and both derivatives cannot be matched by a proper distribution, match only the first.</param>
+        /// <returns></returns>
+        internal static Gamma GammaFromDerivatives(double x, double xdLogP, double x2ddLogP, bool forceProper)
+        {
+            if (x <= 0)
+                throw new ArgumentException("x <= 0");
+            double a = -x2ddLogP;
+            if (forceProper)
+            {
+                if (xdLogP < 0)
+                {
+                    if (a < 0)
+                        a = 0;
+                }
+                else
+                {
+                    double amin = xdLogP;
+                    if (a < amin)
+                        a = amin;
+                }
+            }
+            double b = (a - xdLogP) / x;
+            if (forceProper)
+            {
+                // correct roundoff errors that might make b negative
+                b = Math.Max(b, 0);
+            }
+            if (double.IsNaN(a) || double.IsNaN(b))
+                throw new InferRuntimeException($"result is NaN.  x={x}, xdlogf={xdLogP}, x2ddlogf={x2ddLogP}");
+            return Gamma.FromShapeAndRate(a + 1, b);
         }
 
         /// <include file='FactorDocs.xml' path='factor_docs/message_op_class[@name="GaussianOp"]/message_doc[@name="PrecisionAverageConditional(Gaussian, Gaussian, Gamma)"]/*'/>
@@ -586,16 +600,16 @@ namespace Microsoft.ML.Probabilistic.Factors
             {
                 result.SetToUniform();
             }
-            else if (!precision.IsProper())
-            {
-                // improper prior
-                throw new ImproperMessageException(precision);
-            }
             else if (precision.IsPointMass)
             {
                 // must return a sensible value since precision could be initialized to a point mass.
                 return PrecisionAverageConditional_Point(ym, yv, precision.Point);
             }
+            else if (!precision.IsProper())
+            {
+                // improper prior
+                throw new ImproperMessageException(precision);
+            }
             else
             {
                 // use quadrature to integrate over the precision
@@ -643,7 +657,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                     Z = mva.Count;
                     if (double.IsNaN(Z))
                         throw new Exception("Z is nan");
-                    if (Z > 2)
+                    if (Z > 2 && mva.Variance > 0)
                         break;
                     // one quadrature node dominates the answer.  must re-try with a different weight function.
                     double delta = (argmax == 0) ? (nodes[argmax + 1] - nodes[argmax]) : (nodes[argmax] - nodes[argmax - 1]);
@@ -665,7 +679,6 @@ namespace Microsoft.ML.Probabilistic.Factors
                         {
                             double v = 1.0 / nodes[i];
                             double denom = 1 / (yv + v);
-                            double v2 = v * v;
                             double denom2 = denom * denom;
                             double dlogf1 = -0.5 * denom + 0.5 * ym2 * denom2;
                             // dlogfr = r f'/f
@@ -673,19 +686,19 @@ namespace Microsoft.ML.Probabilistic.Factors
                             double dlogf2 = (0.5 - ym2 * denom) * denom2;
                             // ddfrr = r^2 f''/f
                             // f''/f = d(f'/f) + (f'/f)^2
-                            double ddfrr = dlogfr * dlogfr + dlogf2 * v2 + (2 * v) * dlogf1;
+                            double ddfrr = dlogfr * dlogfr + dlogf2 * v * v + (2 * v) * dlogf1;
                             sum1 += dlogfr * f;
                             sum2 += ddfrr * f;
                         }
                         else  // nodes[i] is small
                         {
                             double r = nodes[i];
-                            double vdenom = 1 / (r*yv + 1);
+                            double vdenom = 1 / (r * yv + 1);
                             double v2denom2 = vdenom * vdenom;
-                            double vdlogf1 = -0.5 * vdenom + 0.5 * ym2 * v2denom2*r;
+                            double vdlogf1 = -0.5 * vdenom + 0.5 * ym2 * v2denom2 * r;
                             // dlogfr = r f'/f
                             double dlogfr = -vdlogf1;
-                            double v2dlogf2 = (0.5 - ym2 * vdenom*r) * v2denom2;
+                            double v2dlogf2 = (0.5 - ym2 * vdenom * r) * v2denom2;
                             // ddfrr = r^2 f''/f
                             // f''/f = d(f'/f) + (f'/f)^2
                             double ddfrr = dlogfr * dlogfr + v2dlogf2 + 2 * vdlogf1;
@@ -715,17 +728,6 @@ namespace Microsoft.ML.Probabilistic.Factors
                     }
                 }
             }
-#if KeepLastMessage
-      if (LastPrecisionMessage != null) {
-        if (Stepsize != 1 && Stepsize != 0) {
-          LastPrecisionMessage.SetToPower(LastPrecisionMessage, 1 - Stepsize);
-          result.SetToPower(result, Stepsize);
-          result.SetToProduct(result, LastPrecisionMessage);
-        }
-      }
-      // FIXME: this is not entirely safe since the caller could overwrite result.
-      LastPrecisionMessage = result;
-#endif
             return result;
         }
 
@@ -796,17 +798,6 @@ namespace Microsoft.ML.Probabilistic.Factors
         /// <param name="logWeights">Place to put the log-weights</param>
         public static void QuadratureNodesAndWeights(Gamma precision, double[] nodes, double[] logWeights)
         {
-#if KeepLastMessage
-      if (LastPrecisionMessage != null) {
-        Gamma PrecisionPosterior = precision * LastPrecisionMessage;
-        Quadrature.GammaNodesAndWeights(PrecisionPosterior.Precision, PrecisionPosterior.PrecisionOverMean, nodes, weights);
-        // modify the weights to include q(prec)/Ga(prec;a,b)
-        for (int i = 0; i < weights.Length; i++) {
-          weights[i] *= Math.Exp(precision.EvaluateLn(nodes[i]) - Gamma.EvaluateLn(nodes[i], PrecisionPosterior.Precision, PrecisionPosterior.PrecisionOverMean));
-        }
-        return;
-      }
-#endif
             Quadrature.GammaNodesAndWeights(precision.Shape - 1, precision.Rate, nodes, logWeights);
         }
 
@@ -1217,7 +1208,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                     for (int i = 0; i < n; i++)
                     {
                         double logr = logrmin + i * inc;
-                        double ivr = 1/(v + Math.Exp(-logr));
+                        double ivr = 1 / (v + Math.Exp(-logr));
                         double dlogf = -m * ivr;
                         double ddlogf = dlogf * dlogf - ivr;
                         double diff = LogLikelihoodRatio(logr, logr0, m, v, a, b);
@@ -1235,7 +1226,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                         sum2 += ddlogf * f;
                     }
                     double alpha = sum1 / Z;
-                    double beta = alpha*alpha -sum2 / Z;
+                    double beta = alpha * alpha - sum2 / Z;
                     return GaussianOp.GaussianFromAlphaBeta(sample, alpha, beta, GaussianOp.ForceProper);
                 }
                 else
@@ -1460,16 +1451,16 @@ namespace Microsoft.ML.Probabilistic.Factors
         public static void GetIntegrationBoundsForPrecision(double m, double v, double a, double b, out double logrmin, out double logrmax, out double logrmode)
         {
             // compute stationary points and inflection points
-            double v2 = v*v;
+            double v2 = v * v;
             double m2 = m * m;
-            double[] coeffs = { -b*v2, (a*v2 - 2*b*v), (0.5*v-0.5*m2+2*a*v-b), (0.5+a) };
+            double[] coeffs = { -b * v2, (a * v2 - 2 * b * v), (0.5 * v - 0.5 * m2 + 2 * a * v - b), (0.5 + a) };
             List<double> stationaryPoints;
-            Predicate<double> isPositive = (r => r>0.0);
+            Predicate<double> isPositive = (r => r > 0.0);
             GetRealRoots(coeffs, out stationaryPoints, isPositive);
             //Console.WriteLine("stationary points = {0}", StringUtil.CollectionToString(stationaryPoints, " "));
             stationaryPoints = stationaryPoints.ConvertAll(x => Math.Log(x));
-            double v3 = v*v2;
-            double[] coeffs2 = { -b*v3, -b*3*v2, (-b*3*v-0.5*v2+0.5*m2*v), (-0.5*v-b-0.5*m2) };
+            double v3 = v * v2;
+            double[] coeffs2 = { -b * v3, -b * 3 * v2, (-b * 3 * v - 0.5 * v2 + 0.5 * m2 * v), (-0.5 * v - b - 0.5 * m2) };
             List<double> inflectionPoints;
             GetRealRoots(coeffs2, out inflectionPoints, isPositive);
             //Console.WriteLine("inflection points = {0}", StringUtil.CollectionToString(inflectionPoints, " "));
@@ -1512,7 +1503,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             double vx = v + Math.Exp(-logx);
             double vx0 = v + Math.Exp(-logx0);
             double diff = MMath.DifferenceOfExp(logx, logx0);
-            return a * (logx - logx0) - b * diff - 0.5 * Math.Log(vx/vx0) - 0.5 * m*m * (1 / vx - 1 / vx0);
+            return a * (logx - logx0) - b * diff - 0.5 * Math.Log(vx / vx0) - 0.5 * m * m * (1 / vx - 1 / vx0);
         }
 
         /// <summary>
@@ -1655,7 +1646,7 @@ namespace Microsoft.ML.Probabilistic.Factors
         {
             if (double.IsInfinity(lowerBound) || double.IsInfinity(upperBound))
                 throw new ArgumentException("infinite bound");
-            while(true)
+            while (true)
             {
                 double x = (lowerBound + upperBound) / 2;
                 double f = func(x);
@@ -1739,7 +1730,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                     break;
                 }
             }
-            int n = coeffs.Count-1 - firstNonZero;
+            int n = coeffs.Count - 1 - firstNonZero;
             if (n <= 0)
             {
                 rootsReal = new double[0];
@@ -1747,7 +1738,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                 return;
             }
             Matrix m = new Matrix(n, n);
-            for (int i = 0; i < n-1; i++)
+            for (int i = 0; i < n - 1; i++)
             {
                 m[i + 1, i] = 1;
             }
@@ -2073,7 +2064,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                             dlogfss = dlogfs(x, m, v);
                             double c1 = 0.5 * x * x2 * (dlogfss[1] - a / x2);
                             double c2 = (dlogfss[0] + a / x - b) + c1 / x2;
-                            if (c1 < 0 && c2 < 0)
+                            if (c1 < 0 && c2 < 0 && c1 > double.MinValue && c2 > double.MinValue)
                             {
                                 x = Math.Sqrt(c1 / c2);
                                 found = true;
@@ -2134,7 +2125,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                 }
                 if (double.IsNaN(x))
                     throw new Exception("x is nan");
-                //Console.WriteLine("{0}: {1}", iter, x);
+                //System.Diagnostics.Trace.WriteLine($"{iter}: {x}");
                 if (MMath.AbsDiff(oldx, x, 1e-10) < 1e-10)
                 {
                     x = QReinitialize(sample, mean, precision, x);
@@ -2154,7 +2145,18 @@ namespace Microsoft.ML.Probabilistic.Factors
             double xdlogf = xdlogfss[0];
             double xxddlogf = xdlogfss[1];
             a = precision.Shape - xxddlogf;
-            b = precision.Rate - (xdlogf + xxddlogf)/x;
+            if (x == 0)
+            {
+                // xdlogf = 0.5/(v*x+1) - 0.5*x*m*m/(v*x+1)^2
+                // xxddlogf = -1/(v*x+1) + x*m*m/(v*x+1)^2 + 0.5/(v*x+1)^2 - x*m*m/(v*x+1)^3
+                // xdlogf + xxddlogf = -0.5*x*v/(v*x+1)^2 + 0.5*x*m*m/(v*x+1)^2 - x*m*m/(v*x+1)^3
+                // limit x->0: x*(-0.5*v - 0.5*m*m)
+                b = precision.Rate + 0.5 * (v + m * m);
+            }
+            else
+            {
+                b = precision.Rate - (xdlogf + xxddlogf) / x;
+            }
             if (a <= 0 || b <= 0)
                 throw new InferRuntimeException();
             if (double.IsNaN(a) || double.IsNaN(b))
@@ -2207,22 +2209,25 @@ namespace Microsoft.ML.Probabilistic.Factors
         {
             if (double.IsPositiveInfinity(v))
                 return new double[4];
+            if (x == 0)
+                return new double[] { double.PositiveInfinity, double.NegativeInfinity, double.PositiveInfinity, double.NegativeInfinity };
             // log f(x) = -0.5*log(v+1/x) -0.5*m^2/(v+1/x)
-            double m2 = m * m;
             double x2 = x * x;
             double x3 = x * x2;
             double x4 = x * x3;
-            double p = 1 / (v + 1 / x);
-            double p2 = p * p;
-            double p3 = p * p2;
-            double dlogf1 = -0.5 * p + 0.5 * m2 * p2;
-            double dlogf = dlogf1 * (-1 / x2);
-            double ddlogf1 = 0.5 * p2 - m2 * p3;
-            double ddlogf = dlogf1 * 2 / x3 + ddlogf1 / x4;
-            double dddlogf1 = -p3 + 3 * m2 * p * p3;
-            double dddlogf = dlogf1 * (-6) / x4 + ddlogf1 * (-6) / (x * x4) + dddlogf1 * (-1) / (x2 * x4);
-            double d4logf1 = 3 * p * p3 - 12 * m2 * p2 * p3;
-            double d4logf = dlogf1 * (24) / (x2 * x3) + ddlogf1 * 36 / (x3 * x3) + dddlogf1 * (12) / (x4 * x3) + d4logf1 / (x4 * x4);
+            double pix = 1 / (v * x + 1);
+            double pix2 = pix / x;
+            double p2ix2 = pix * pix;
+            double p3ix3 = p2ix2 * pix;
+            double m2p2ix2 = p2ix2 * m * m;
+            double dlogf1ix2 = -0.5 * pix2 + 0.5 * m2p2ix2;
+            double dlogf = dlogf1ix2 * (-1);
+            double ddlogf1ix3 = 0.5 * pix * pix2 - m2p2ix2 * pix;
+            double ddlogf = (dlogf1ix2 * 2 + ddlogf1ix3) / x;
+            double dddlogf1ix4 = -p2ix2 * pix2 + 3 * m2p2ix2 * p2ix2;
+            double dddlogf = (dlogf1ix2 * (-6) + ddlogf1ix3 * (-6) + dddlogf1ix4 * (-1)) / x2;
+            double d4logf1ix5 = 3 * pix2 * p3ix3 - 12 * m2p2ix2 * p3ix3;
+            double d4logf = (dlogf1ix2 * 24 + ddlogf1ix3 * 36 + dddlogf1ix4 * 12 + d4logf1ix5) / x3;
             return new double[] { dlogf, ddlogf, dddlogf, d4logf };
         }
 
@@ -2238,9 +2243,9 @@ namespace Microsoft.ML.Probabilistic.Factors
             if (double.IsPositiveInfinity(v))
                 return new double[4];
             // log f(x) = -0.5*log(v+1/x) -0.5*m^2/(v+1/x)
-            double m2 = m * m;
-            if (x*v > 1)
+            if (x * v > 1 && false)
             {
+                double m2 = m * m;
                 double x2 = x * x;
                 double x3 = x * x2;
                 double x4 = x * x3;
@@ -2259,19 +2264,22 @@ namespace Microsoft.ML.Probabilistic.Factors
             }
             else  // x is small
             {
+                double m2x = x * m * m;
                 double x2 = x * x;
                 double x3 = x * x2;
                 double x4 = x * x3;
-                double p = 1 / (v*x + 1);
+                double p = 1 / (v * x + 1);
                 double p2 = p * p;
                 double p3 = p * p2;
-                double ixdlogf1 = -0.5 * p + 0.5 * m2 * p2*x;
+                // xdlogf = 0.5/(v*x+1) - 0.5*x*m*m/(v*x+1)^2
+                double ixdlogf1 = -0.5 * p + 0.5 * m2x * p2;
                 double xdlogf = -ixdlogf1;
-                double ix2ddlogf1 = 0.5 * p2 - m2 * p3*x;
+                double ix2ddlogf1 = 0.5 * p2 - m2x * p3;
+                // xxddlogf = -1/(v*x+1) + x*m*m/(v*x+1)^2 + 0.5/(v*x+1)^2 - x*m*m/(v*x+1)^3
                 double xxddlogf = ixdlogf1 * 2 + ix2ddlogf1;
-                double ix3dddlogf1 = -p3 + 3 * m2 * p * p3*x;
+                double ix3dddlogf1 = -p3 + 3 * m2x * p * p3;
                 double xxxdddlogf = ixdlogf1 * (-6) + ix2ddlogf1 * (-6) + ix3dddlogf1 * (-1);
-                double ix4d4logf1 = 3 * p * p3 - 12 * m2 * p2 * p3*x;
+                double ix4d4logf1 = 3 * p * p3 - 12 * m2x * p2 * p3;
                 double x4d4logf = ixdlogf1 * (24) + ix2ddlogf1 * 36 + ix3dddlogf1 * (12) + ix4d4logf1;
                 return new double[] { xdlogf, xxddlogf, xxxdddlogf, x4d4logf };
             }
@@ -2549,7 +2557,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             Gaussian diffPost = diffPrior * diffLike;
             double md, vd;
             diffPost.GetMeanAndVariance(out md, out vd);
-            result.Rate = 0.5 * (vd + md*md);
+            result.Rate = 0.5 * (vd + md * md);
             return result;
         }
 

src/Runtime/Factors/IsBetween.cs

@@ -5,6 +5,7 @@
 namespace Microsoft.ML.Probabilistic.Factors
 {
     using System;
+    using System.Diagnostics;
     using Microsoft.ML.Probabilistic.Distributions;
     using Microsoft.ML.Probabilistic.Math;
     using Microsoft.ML.Probabilistic.Factors.Attributes;
@@ -19,6 +20,8 @@ namespace Microsoft.ML.Probabilistic.Factors
         /// Static flag to force a proper distribution
         /// </summary>
         public static bool ForceProper = true;
+        public static double LowPrecisionThreshold = 0; // 1e-8;
+        public static double LargeMeanThreshold = 2.5e3;
 
         //-- TruncatedGaussian bounds ------------------------------------------------------------------------------
 
@@ -651,8 +654,11 @@ namespace Microsoft.ML.Probabilistic.Factors
         //-- Random bounds --------------------------------------------------------------------------------
 
         // Compute the mean and variance of (X-L) and (U-X)
-        internal static void GetDiffMeanAndVariance(Gaussian X, Gaussian L, Gaussian U, out double yl, out double yu, out double r, out double invSqrtVxl,
-                                                    out double invSqrtVxu)
+        // sqrtomr2 is sqrt(1-r*r) with high accuracy.
+        internal static void GetDiffMeanAndVariance(Gaussian X, Gaussian L, Gaussian U, out double yl, out double yu, out double r, 
+            out double sqrtomr2,
+            out double invSqrtVxl,
+            out double invSqrtVxu)
         {
             double mx, vx, ml, vl, mu, vu;
             X.GetMeanAndVariance(out mx, out vx);
@@ -691,14 +697,19 @@ namespace Microsoft.ML.Probabilistic.Factors
             if (X.IsPointMass)
             {
                 r = 0.0;
+                sqrtomr2 = 1;
             }
             else
             {
                 //r = -vx * invSqrtVxl * invSqrtVxu;
-                // This is a more accurate way to compute the above.
-                r = -1 / Math.Sqrt(1 + vl / vx) / Math.Sqrt(1 + vu / vx);
+                // This formula ensures r is between -1 and 1.
+                //r = -1 / Math.Sqrt(1 + vl / vx) / Math.Sqrt(1 + vu / vx);
+                //r = -vx / Math.Sqrt((vx + vl) * (vx + vu));
+                r = Math.Max(-1, Math.Min(1, -vx / Math.Sqrt(vx + vl) / Math.Sqrt(vx + vu)));
                 if (r < -1 || r > 1)
                     throw new Exception("Internal: r is outside [-1,1]");
+                double omr2 = ((vl + vu) * vx + vl * vu) / (vx + vl) / (vx + vu);
+                sqrtomr2 = Math.Sqrt(omr2);
             }
         }
 
@@ -733,15 +744,21 @@ namespace Microsoft.ML.Probabilistic.Factors
                     return Double.NegativeInfinity;
                 }
             }
+            else if (!X.IsProper())
+            {
+                return double.NegativeInfinity;
+            }
             else
             {
                 // at this point, X is not uniform
-                double yl, yu, r, invSqrtVxl, invSqrtVxu;
-                GetDiffMeanAndVariance(X, L, U, out yl, out yu, out r, out invSqrtVxl, out invSqrtVxu);
-                double logp = MMath.NormalCdfLn(yl, yu, r);
-                if (logp > 0)
+                double yl, yu, r, sqrtomr2, invSqrtVxl, invSqrtVxu;
+                GetDiffMeanAndVariance(X, L, U, out yl, out yu, out r, out sqrtomr2, out invSqrtVxl, out invSqrtVxu);
+                //Trace.WriteLine($"yl={yl:r} yu={yu:r} r={r:r} sqrtomr2={sqrtomr2:r}");
+                var prob = MMath.NormalCdf(yl, yu, r, sqrtomr2);
+                double logProb = prob.Log();
+                if (logProb > 0)
                     throw new Exception("LogProbBetween is positive");
-                return logp;
+                return logProb;
             }
         }
 
@@ -813,7 +830,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             Gaussian result = new Gaussian();
             if (isBetween.IsUniform())
                 return result;
-            if (X.IsUniform())
+            if (X.Precision == 0)
             {
                 if (upperBound.IsUniform() || lowerBound.IsUniform())
                 {
@@ -849,7 +866,6 @@ namespace Microsoft.ML.Probabilistic.Factors
             }
             else
             {
-                //double logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
                 double d_p = 2 * isBetween.GetProbTrue() - 1;
                 if (lowerBound.IsPointMass)
                 {
@@ -900,29 +916,80 @@ namespace Microsoft.ML.Probabilistic.Factors
                         }
                     }
                 }
-                double yl, yu, r, invSqrtVxl, invSqrtVxu;
-                GetDiffMeanAndVariance(X, lowerBound, upperBound, out yl, out yu, out r, out invSqrtVxl, out invSqrtVxu);
+                if (X.IsPointMass)
+                {
+                    if (upperBound.IsPointMass && upperBound.Point < X.Point)
+                    {
+                        // The constraint reduces to (lowerBound < upperBound), which is (lowerBound - upperBound < 0)
+                        Gaussian shifted_F = DoublePlusOp.AAverageConditional(lowerBound, upperBound.Point);
+                        Gaussian shifted_B = IsPositiveOp.XAverageConditional(false, shifted_F);
+                        return DoublePlusOp.SumAverageConditional(shifted_B, upperBound.Point);
+                    }
+                    else
+                    {
+                        // The constraint reduces to (lowerBound < X), which is (lowerBound - X < 0)
+                        Gaussian shifted_F = DoublePlusOp.AAverageConditional(lowerBound, X.Point);
+                        Gaussian shifted_B = IsPositiveOp.XAverageConditional(false, shifted_F);
+                        return DoublePlusOp.SumAverageConditional(shifted_B, X.Point);
+                    }
+                }
+                bool precisionWasZero = AdjustXPrecision(isBetween, ref X, lowerBound, upperBound, ref logZ, 1e-0);
+                double yl, yu, r, sqrtomr2, invSqrtVxl, invSqrtVxu;
+                GetDiffMeanAndVariance(X, lowerBound, upperBound, out yl, out yu, out r, out sqrtomr2, out invSqrtVxl, out invSqrtVxu);
+                bool useLogZRatio = (r > smallR) && (logZ < smallLogZ);
+                double logZRatio = useLogZRatio ? MMath.NormalCdfRatioLn(yl, yu, r, sqrtomr2) : 0;
                 // if we get here, we know that -1 < r <= 0 and invSqrtVxl is finite
                 // since lowerBound is not uniform and X is not uniform, invSqrtVxl > 0
-                // yl is always finite.  yu may be +/-infinity, in which case r = 0.
-                double logPhiL = Gaussian.GetLogProb(yl, 0, 1) + MMath.NormalCdfLn((yu - r * yl) / Math.Sqrt(1 - r * r));
-                double alphaL;
-                if (r == 0)
+                // yl is always finite.  yu may be +/-infinity.
+                double alphaL, betaL;
+                GetAlpha(X, lowerBound, upperBound, logZ, logZRatio, d_p, yl, yu, r, sqrtomr2, invSqrtVxl, invSqrtVxu, true, out alphaL, false, out double alphaU, out double alphaX, out double ylInvSqrtVxlPlusAlphaX, out double yuInvSqrtVxuMinusAlphaX);
+                if (d_p == 1 && useLogZRatio && yl < 0 && yu < 0)
                 {
-                    alphaL = -d_p * invSqrtVxl / MMath.NormalCdfRatio(yl);
+                    double invZRatio = Math.Exp(-logZRatio);
+                    double yuryl = (yu - r * yl) / sqrtomr2;
+                    double Ryuryl = MMath.NormalCdfRatio(yuryl);
+                    double ylryu = (yl - r * yu) / sqrtomr2;
+                    double Rylryu = MMath.NormalCdfRatio(ylryu);
+                    // alphaL = -invSqrtVxl*R(yuryl)/ZRatio
+                    // alphaL*(alphaL - yl*invSqrtVxl) = invSqtVxl^2*R(yuryl)/ZRatio*(R(yuryl)/ZRatio + yl)
+                    // beta = q * invSqrtVxl^2 / ZRatio
+                    //double q = Ryuryl * (Ryuryl * invZRatio + yl) + r / Math.Sqrt(omr2);
+                    double R1yuryl = MMath.NormalCdfMomentRatio(1, yuryl);
+                    double R1ylryu = MMath.NormalCdfMomentRatio(1, ylryu);
+                    // This is an asymptotic approximation of phi1(x,y,r)/phi_r(x,y,r)/sqrtomr2 where phi1 is the integral of phi wrt x
+                    double intZRatio = (r * Rylryu + Ryuryl - r * sqrtomr2 * R1ylryu * yl) / (yl * yl + 1);
+                    // Substitute Math.Exp(logZRatio) = (intZRatio - r * Rylryu - Ryuryl) / yl
+                    //double u2 = Ryuryl * (intZRatio - r * Rylryu) * invZRatio + r / Math.Sqrt(omr2);
+                    // r*yuryl + ylryu = yl*sqrtomr2
+                    // yuryl = (R2 - 1)/Ryuryl
+                    // Substitute r/sqrtomr2 = r*yl/(r*yuryl + ylryu) = r*yl/(r*(R2yuryl - 1)/Ryuryl + (R2ylryu - 1)/Rylryu)
+                    // = r*yl*Ryuryl*Rylryu/(r*(R2yuryl - 1)*Rylryu + (R2ylryu - 1)*Ryuryl)
+                    double w = r * R1yuryl * Rylryu + R1ylryu * Ryuryl;
+                    //double u3 = Ryuryl * (intZRatio - r * Rylryu) * invZRatio + r * yl * Ryuryl * Rylryu / (w - (r * Rylryu + Ryuryl));
+                    //double u4 = Ryuryl * (intZRatio - r * Rylryu + r * yl * Rylryu * Math.Exp(logZRatio) / (w - (r * Rylryu + Ryuryl))) * invZRatio;
+                    //u4 = Ryuryl * (intZRatio - r * Rylryu + r * Rylryu * (intZRatio - r * Rylryu - Ryuryl) / (w - (r * Rylryu + Ryuryl))) * invZRatio;
+                    //u4 = Ryuryl * (intZRatio + r * Rylryu * (intZRatio - w) / (w - (r * Rylryu + Ryuryl))) * invZRatio;
+                    double q = Ryuryl * (intZRatio * (w - Ryuryl) - r * Rylryu * w) / (w - (r * Rylryu + Ryuryl)) * invZRatio;
+                    betaL = q * invSqrtVxl * invSqrtVxl * invZRatio;
+                    if (double.IsNaN(betaL)) throw new Exception("betaL is NaN");
                 }
                 else
                 {
-                    alphaL = -d_p * invSqrtVxl * Math.Exp(logPhiL - logZ);
-                }
-                // (mx - ml) / (vl + vx) = yl*invSqrtVxl
-                double betaL = alphaL * (alphaL - yl * invSqrtVxl);
-                if (r > -1 && r != 0)
-                {
-                    double logPhiR = -2 * MMath.LnSqrt2PI - 0.5 * Math.Log(1 - r * r) - 0.5 * (yl * yl + yu * (yu - 2 * r * yl)) / (1 - r * r);
-                    double c = d_p * r * Math.Exp(logPhiR - logZ);
-                    betaL += c * invSqrtVxl * invSqrtVxl;
+                    // (mx - ml) / (vl + vx) = yl*invSqrtVxl
+                    betaL = alphaL * (alphaL - yl * invSqrtVxl);
+                    // This formula is slightly more accurate when r == -1 && yl > 0 && yu < 0:
+                    // betaU = alphaU * invSqrtVxu * (MMath.NormalCdfMomentRatio(1, yu) - yu * ncrl / delta) / ZoverPhiU;
+                    // because yu * ZoverPhiU + 1 = ncru*yu+1 - yu*ncrl / delta
+                    if (r > -1 && r != 0 && !precisionWasZero)
+                    {
+                        double omr2 = sqrtomr2 * sqrtomr2;
+                        double logPhiR = GetLogPhiR(X, lowerBound, upperBound, yl, yu, r, omr2, logZ, logZRatio);
+                        double c = d_p * r * Math.Exp(logPhiR);
+                        betaL += c * invSqrtVxl * invSqrtVxl;
+                        if (double.IsNaN(betaL)) throw new Exception("betaL is NaN");
+                    }
                 }
+                //Trace.WriteLine($"alpha = {alphaL} beta = {betaL} yl = {yl} yu = {yu} r = {r}");
                 return GaussianOp.GaussianFromAlphaBeta(lowerBound, alphaL, betaL, ForceProper);
             }
             if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision))
@@ -943,122 +1010,12 @@ namespace Microsoft.ML.Probabilistic.Factors
             [SkipIfUniform] Bernoulli isBetween, [RequiredArgument] Gaussian X, [RequiredArgument] Gaussian lowerBound, [RequiredArgument] Gaussian upperBound,
             double logZ)
         {
-            Gaussian result = new Gaussian();
-            if (isBetween.IsUniform())
-                return result;
-            if (X.IsUniform())
-            {
-                if (lowerBound.IsUniform() || upperBound.IsUniform())
-                {
-                    result.SetToUniform();
-                }
-                else if (isBetween.IsPointMass && isBetween.Point)
-                {
-                    double ml, vl, mu, vu;
-                    lowerBound.GetMeanAndVariance(out ml, out vl);
-                    upperBound.GetMeanAndVariance(out mu, out vu);
-                    double vlu = vl + vu;
-                    double alpha = Math.Exp(Gaussian.GetLogProb(ml, mu, vlu) - MMath.NormalCdfLn((mu - ml) / Math.Sqrt(vlu)));
-                    double alphaU = 1.0 / (mu - ml + vlu * alpha);
-                    double betaU = alphaU * (alphaU - alpha);
-                    result.SetMeanAndVariance(mu + vu * alphaU, vu - vu * vu * betaU);
-                    result.SetToRatio(result, upperBound);
-                }
-                else
-                    throw new NotImplementedException();
-            }
-            else if (upperBound.IsUniform())
-            {
-                if (isBetween.IsPointMass && !isBetween.Point)
-                {
-                    // lowerBound <= upperBound <= X
-                    // upperBound is not a point mass so upperBound==X is impossible
-                    return XAverageConditional(Bernoulli.PointMass(true), upperBound, lowerBound, X, logZ);
-                }
-                else
-                {
-                    result.SetToUniform();
-                }
-            }
-            else
-            {
-                //double logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
-                double d_p = 2 * isBetween.GetProbTrue() - 1;
-                if (upperBound.IsPointMass)
-                {
-                    if (double.IsPositiveInfinity(upperBound.Point)) return Gaussian.Uniform();
-                    if (X.IsPointMass)
-                    {
-                        if (upperBound.Point > X.Point) return Gaussian.Uniform();
-                        else return X;
-                    }
-                    if (lowerBound.IsPointMass)
-                    {
-                        // r = -1 case
-                        // X is not uniform or point mass
-                        // f(U) = d_p (NormalCdf((U-mx)*sqrtPrec) - NormalCdf((L-mx)*sqrtPrec)) + const.
-                        // dlogf/dU = d_p N(U;mx,vx)/f
-                        // ddlogf = -dlogf^2 + dlogf*(mx-U)/vx
-                        double U = upperBound.Point;
-                        double mx = X.GetMean();
-                        if (mx > U && d_p == 1)
-                        {
-                            // Z = MMath.NormalCdf(sqrtPrec*(U-mx)) - MMath.NormalCdf(sqrtPrec*(L-mx))
-                            // Z/X.GetProb(U)*sqrtPrec = MMath.NormalCdfRatio(sqrtPrec*(U-mx)) - MMath.NormalCdfRatio(sqrtPrec*(L-mx))*X.GetProb(L)/X.GetProb(U)
-                            // X.GetProb(L)/X.GetProb(U) = Math.Exp(-X.MeanTimesPrecision*(U-L) + 0.5*(U*U - L*L)*X.Precision) =approx 0
-                            double sqrtPrec = Math.Sqrt(X.Precision);
-                            double L = lowerBound.Point;
-                            double Umx = sqrtPrec * (U - mx);
-                            double UCdfRatio = MMath.NormalCdfRatio(Umx);
-                            double LCdfRatio = MMath.NormalCdfRatio(sqrtPrec * (L - mx));
-                            double LPdfRatio = Math.Exp((U - L) * (-X.MeanTimesPrecision + 0.5 * (U + L) * X.Precision));
-                            double CdfRatioDiff = UCdfRatio - LCdfRatio * LPdfRatio;
-                            double dlogf = d_p / CdfRatioDiff * sqrtPrec;
-                            // ddlogf = -dlogf * (d_p * sqrtPrec / CdfRatioDiff + sqrtPrec*Umx)
-                            //        = -dlogf * sqrtPrec * (d_p / CdfRatioDiff + Umx)
-                            //        = -dlogf * sqrtPrec * (d_p + Umx * CdfRatioDiff) / CdfRatioDiff
-                            //        = -dlogf * sqrtPrec * (d_p - 1 - Umx * LCdfRatio * LPdfRatio + NormalCdfMomentRatio(1,Umx)) / CdfRatioDiff
-                            double ddlogf = -dlogf * dlogf * (Umx * (CdfRatioDiff - UCdfRatio) + MMath.NormalCdfMomentRatio(1, Umx));
-                            return Gaussian.FromDerivatives(U, dlogf, ddlogf, ForceProper);
-                            // this is equivalent
-                            // U - dlogf/ddlogf = U - 1/(-dlogf + (X.MeanTimesPrecision - U * X.Precision)) =approx mx
-                            //return Gaussian.FromMeanAndPrecision(U - dlogf / ddlogf, -ddlogf);
-                        }
-                        else
-                        {
-                            double dlogf = d_p * Math.Exp(X.GetLogProb(U) - logZ);
-                            double ddlogf = dlogf * (-dlogf + (X.MeanTimesPrecision - U * X.Precision));
-                            return Gaussian.FromDerivatives(U, dlogf, ddlogf, ForceProper);
-                        }
-                    }
-                }
-                double yl, yu, r, invSqrtVxl, invSqrtVxu;
-                GetDiffMeanAndVariance(X, lowerBound, upperBound, out yl, out yu, out r, out invSqrtVxl, out invSqrtVxu);
-                // if we get here, -1 < r <= 0 and invSqrtVxu is finite
-                // since upperBound is not uniform and X is not uniform, invSqrtVxu > 0
-                // yu is always finite.  yl may be infinity, in which case r = 0.
-                double logPhiU = Gaussian.GetLogProb(yu, 0, 1) + MMath.NormalCdfLn((yl - r * yu) / Math.Sqrt(1 - r * r));
-                double alphaU;
-                if (r == 0)
-                {
-                    alphaU = d_p * invSqrtVxu / MMath.NormalCdfRatio(yu);
-                }
-                else
-                {
-                    alphaU = d_p * invSqrtVxu * Math.Exp(logPhiU - logZ);
-                }
-                // (mu - mx) / (vx + vu) = yu*invSqrtVxu
-                double betaU = alphaU * (alphaU + yu * invSqrtVxu);
-                if (r > -1 && r != 0)
-                {
-                    double logPhiR = -2 * MMath.LnSqrt2PI - 0.5 * Math.Log(1 - r * r) - 0.5 * (yu * yu + yl * (yl - 2 * r * yu)) / (1 - r * r);
-                    double c = d_p * r * Math.Exp(logPhiR - logZ);
-                    betaU += c * invSqrtVxu * invSqrtVxu;
-                }
-                return GaussianOp.GaussianFromAlphaBeta(upperBound, alphaU, betaU, ForceProper);
-            }
-            if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision))
-                throw new InferRuntimeException($"{nameof(result)} is NaN.  {nameof(isBetween)}={isBetween}, {nameof(X)}={X}, {nameof(lowerBound)}={lowerBound}, {nameof(upperBound)}={upperBound}, {nameof(logZ)}={logZ}");
+            Gaussian result = LowerBoundAverageConditional(isBetween,
+                Gaussian.FromNatural(-X.MeanTimesPrecision, X.Precision),
+                Gaussian.FromNatural(-upperBound.MeanTimesPrecision, upperBound.Precision),
+                Gaussian.FromNatural(-lowerBound.MeanTimesPrecision, lowerBound.Precision),
+                logZ);
+            result.MeanTimesPrecision *= -1;
             return result;
         }
 
@@ -1089,7 +1046,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             {
                 return upperBound;
             }
-            else if (X.IsUniform())
+            else if (X.Precision == 0)
             {
                 if (lowerBound.IsUniform() || upperBound.IsUniform() ||
                     (lowerBound.IsPointMass && Double.IsInfinity(lowerBound.Point)) ||
@@ -1119,64 +1076,369 @@ namespace Microsoft.ML.Probabilistic.Factors
             }
             else
             {
-                // X is not a point mass or uniform
-                //double logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
+                // X is not uniform
+                bool precisionWasZero = AdjustXPrecision(isBetween, ref X, lowerBound, upperBound, ref logZ);
                 if (Double.IsNegativeInfinity(logZ))
                     throw new AllZeroException();
                 double d_p = 2 * isBetween.GetProbTrue() - 1;
-                double yl, yu, r, invSqrtVxl, invSqrtVxu;
-                GetDiffMeanAndVariance(X, lowerBound, upperBound, out yl, out yu, out r, out invSqrtVxl, out invSqrtVxu);
-                // r == -1 iff lowerBound and upperBound are point masses
-                // since X is not a point mass, invSqrtVxl is finite, invSqrtVxu is finite
-                double alphaL = 0.0;
-                if (X.IsPointMass && !Double.IsInfinity(yl))
-                {
-                    alphaL = -d_p * invSqrtVxl / MMath.NormalCdfRatio(yl);
-                }
-                else if (!lowerBound.IsUniform() && !Double.IsInfinity(yl))
-                {
-                    // since X and lowerBound are not both uniform, invSqrtVxl > 0
-                    double logPhiL = Gaussian.GetLogProb(yl, 0, 1);
-                    if (r > -1)
-                        logPhiL += MMath.NormalCdfLn((yu - r * yl) / Math.Sqrt(1 - r * r));
-                    alphaL = -d_p * invSqrtVxl * Math.Exp(logPhiL - logZ);
-                    // TODO: make this case smoothly blend into the X.IsPointMass case
-                }
-                double alphaU = 0.0;
-                if (X.IsPointMass && !Double.IsInfinity(yu))
-                {
-                    alphaU = d_p * invSqrtVxu / MMath.NormalCdfRatio(yu);
-                }
-                else if (!upperBound.IsUniform() && !Double.IsInfinity(yu))
-                {
-                    // since X and upperBound are not both uniform, invSqrtVxu > 0
-                    double logPhiU = Gaussian.GetLogProb(yu, 0, 1);
-                    if (r > -1)
-                        logPhiU += MMath.NormalCdfLn((yl - r * yu) / Math.Sqrt(1 - r * r));
-                    alphaU = d_p * invSqrtVxu * Math.Exp(logPhiU - logZ);
-                }
-                double alphaX = -alphaL - alphaU;
+                double yl, yu, r, sqrtomr2, invSqrtVxl, invSqrtVxu;
+                GetDiffMeanAndVariance(X, lowerBound, upperBound, out yl, out yu, out r, out sqrtomr2, out invSqrtVxl, out invSqrtVxu);
+                bool useLogZRatio = (r > smallR) && (logZ < smallLogZ);
+                double logZRatio = useLogZRatio ? MMath.NormalCdfRatioLn(yl, yu, r, sqrtomr2) : 0;
+                double alphaL, alphaU, alphaX;
+                GetAlpha(X, lowerBound, upperBound, logZ, logZRatio, d_p, yl, yu, r, sqrtomr2, invSqrtVxl, invSqrtVxu, true, out alphaL, true, out alphaU, out alphaX, out double ylInvSqrtVxlPlusAlphaX, out double yuInvSqrtVxuMinusAlphaX);
+                //alphaX = -alphaL - alphaU;
                 // (mx - ml) / (vl + vx) = yl*invSqrtVxl
-                double betaX = alphaX * alphaX;
+                // To improve accuracy here, could rewrite betaX to be relative to X.Precision, or solve for posterior on X directly as in uniform case.
+                // betaX = alphaX * alphaX - alphaL * (yl * invSqrtVxl) + alphaU * (yu * invSqrtVxu)
+                //       = (-alphaL - alphaU)*alphaX - alphaL * (yl * invSqrtVxl) + alphaU * (yu * invSqrtVxu)
+                //       = - alphaL * (yl * invSqrtVxl + alphaX) + alphaU * (yu * invSqrtVxu - alphaX)
+                double betaX = 0;
                 if (!Double.IsInfinity(yl))
                 {
-                    betaX -= alphaL * (yl * invSqrtVxl);
+                    // if yl is infinity then alphaL == 0
+                    betaX -= alphaL * ylInvSqrtVxlPlusAlphaX;
                 }
                 if (!Double.IsInfinity(yu))
                 {
-                    betaX += alphaU * (yu * invSqrtVxu);
+                    // if yu is infinity then alphaU == 0
+                    betaX += alphaU * yuInvSqrtVxuMinusAlphaX;
                 }
-                if (r > -1 && r != 0 && !Double.IsInfinity(yl) && !Double.IsInfinity(yu))
+                if (r > -1 && r != 0 && !Double.IsInfinity(yl) && !Double.IsInfinity(yu) && !precisionWasZero)
                 {
-                    double logPhiR = -2 * MMath.LnSqrt2PI - 0.5 * Math.Log(1 - r * r) - 0.5 * (yl * yl + yu * yu - 2 * r * yl * yu) / (1 - r * r);
-                    double c = d_p * r * Math.Exp(logPhiR - logZ);
+                    double omr2 = sqrtomr2 * sqrtomr2;
+                    double logPhiR = GetLogPhiR(X, lowerBound, upperBound, yl, yu, r, omr2, logZ, logZRatio);
+                    double c = d_p * r * Math.Exp(logPhiR);
                     betaX += c * (-2 * X.Precision + invSqrtVxl * invSqrtVxl + invSqrtVxu * invSqrtVxu);
                 }
+                if (TraceAlpha)
+                    Trace.WriteLine($"yu = {yu} yl = {yl} r = {r} alphaX={alphaX}, alphaL={alphaL}, alphaU={alphaU}, betaX={betaX}, ylInvSqrtVxlPlusAlphaX = {ylInvSqrtVxlPlusAlphaX}, yuInvSqrtVxuMinusAlphaX = {yuInvSqrtVxuMinusAlphaX}");
                 return GaussianOp.GaussianFromAlphaBeta(X, alphaX, betaX, ForceProper);
             }
             return result;
         }
 
+        public static bool TraceAlpha;
+        private const double smallR = -1;
+        private const double smallLogZ = -1e4;
+
+        // Invariant to swapping yu and yl
+        private static double GetLogPhiR(Gaussian X, Gaussian lowerBound, Gaussian upperBound, double yl, double yu, double r, double omr2, double logZ, double logZRatio)
+        {
+            bool useLogZRatio = (r > smallR) && (logZ < smallLogZ);
+            double logPhiR = -0.5 * Math.Log(omr2);
+            if (useLogZRatio) return logPhiR - logZRatio;
+            logPhiR += -2 * MMath.LnSqrt2PI;
+            if (r > smallR)
+            {
+                if (Math.Abs(yu) > Math.Abs(yl))
+                    logPhiR -= 0.5 * (yl * yl + yu * (yu - 2 * r * yl)) / omr2;
+                else
+                    logPhiR -= 0.5 * (yu * yu + yl * (yl - 2 * r * yu)) / omr2;
+            }
+            else
+            {
+                // yu = (mu-mx)*invSqrtVxu
+                // r = -vx*invSqrtVxu*invSqrtVxl
+                // (yu * yu + yl * (yl - 2 * r * yu)) =
+                // ((mu-mx)*(mu-mx)/(vx+vu) + (mx-ml)*invSqrtVxl * ((mx-ml)*invSqrtVxl - 2 * r * (mu-mx)*invSqrtVxu))
+                // ((mu-mx)*(mu-mx)/(vx+vu) + (mx-ml) * ((mx-ml) + 2 *vx * invSqrtVxu* (mu-mx)*invSqrtVxu)/(vx+vl))
+                // ((mu-mx)*(mu-mx)/(vx+vu) + (mx-ml) * ((mx-ml) + 2 *vx * (mu-mx)/(vx+vu))/(vx+vl))
+                // (mu-mx)*(mu-mx)/(vx+vu) + (mx-ml) * ((mx-ml)*(vx+vu) + (2*vx*mu-2*vx*mx))/(vx+vu)/(vx+vl)
+                // ((mu-mx)*(mu-mx) + (mx-ml) * ((mx-ml)*(vx+vu) + (2*vx*mu-2*vx*mx))/(vx+vl))/(vx+vu)
+                // ((mu-mx)*(mu-mx) + (mx-ml) * ((mx-ml)*vu + (2*mu-ml-mx)*vx)/(vx+vl))/(vx+vu)
+                // ((mu-mx)*(mu-mx)*(vx+vl) + (mx-ml) * ((mx-ml)*vu + (2*mu-ml-mx)*vx))/(vx+vl)/(vx+vu)
+                // (((mu-mx)*(mu-mx) + (mx-ml)*(2*mu-ml-mx))*vx + (mu-mx)*(mu-mx)*vl + (mx-ml)*(mx-ml)*vu)/(vx+vl)/(vx+vu)
+                // ((mu-ml)*(mu-ml)*vx + (mu-mx)*(mu-mx)*vl + (mx-ml)*(mx-ml)*vu)/(vx+vl)/(vx+vu)
+                // (1-r*r) = 1 - vx*vx/(vx+vl)/(vx+vu) = (vx*vl + vx*vu + vl*vu)/(vx+vl)/(vx+vu)
+                double mx, vx, ml, vl, mu, vu;
+                X.GetMeanAndVariance(out mx, out vx);
+                lowerBound.GetMeanAndVariance(out ml, out vl);
+                upperBound.GetMeanAndVariance(out mu, out vu);
+                double dmul = mu - ml;
+                double dmux = mu - mx; // yu/invSqrtVxu
+                double dmxl = mx - ml;
+                logPhiR -= 0.5 * (dmul * dmul * vx + dmux * dmux * vl + dmxl * dmxl * vu) / (vx * vl + vx * vu + vl * vu);
+            }
+            if (double.IsNaN(logPhiR)) throw new Exception($"logPhiR is NaN for X={X}, lowerBound={lowerBound}, upperBound={upperBound}, yl={yl}, yu={yu}, r={r}");
+            return logPhiR - logZ;
+        }
+
+        private static bool AdjustXPrecision(Bernoulli isBetween, ref Gaussian X, Gaussian lowerBound, Gaussian upperBound, ref double logZ, double precisionScaling = 1)
+        {
+            bool precisionWasZero = false;
+            double minPrecision = precisionScaling * LowPrecisionThreshold * Math.Min(lowerBound.Precision, upperBound.Precision);
+            // minPrecision cannot be infinity here
+            if (double.IsPositiveInfinity(minPrecision)) throw new Exception("minPrecision is infinity");
+            //Trace.WriteLine($"X.Precision={X.Precision}, minPrecision={minPrecision}");
+            // When X.Precision is small, r is close to -1.
+            if (X.Precision < minPrecision)
+            {
+                double mx, vx, ml, vl, mu, vu;
+                X.GetMeanAndVariance(out mx, out vx);
+                lowerBound.GetMeanAndVariance(out ml, out vl);
+                upperBound.GetMeanAndVariance(out mu, out vu);
+                double maxMean = Math.Max(Math.Abs(ml), Math.Abs(mu)) * LargeMeanThreshold;
+                //Trace.WriteLine($"mx={Math.Abs(mx)}, maxMean={maxMean}");
+                if (Math.Abs(mx) > maxMean)
+                {
+                    // in this case, logZ is inaccurate, making all further computations inaccurate.
+                    // To prevent this, we make X.Precision artificially larger.
+                    //Trace.WriteLine("precisionWasZero");
+                    precisionWasZero = true;
+                    //X.SetMeanAndPrecision(mx, minPrecision);
+                    X.Precision = minPrecision;
+                    logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
+                }
+            }
+
+            return precisionWasZero;
+        }
+
+        private static void GetAlpha(Gaussian X, Gaussian lowerBound, Gaussian upperBound, double logZ, double logZRatio,
+            double d_p, double yl, double yu, double r, double sqrtomr2, double invSqrtVxl, double invSqrtVxu,
+            bool getAlphaL, out double alphaL, bool getAlphaU, out double alphaU, out double alphaX, out double ylInvSqrtVxlPlusAlphaX, out double yuInvSqrtVxuMinusAlphaX)
+        {
+            double mx, vx, ml, vl, mu, vu;
+            X.GetMeanAndVariance(out mx, out vx);
+            lowerBound.GetMeanAndVariance(out ml, out vl);
+            upperBound.GetMeanAndVariance(out mu, out vu);
+            bool useLogZRatio = (r > smallR) && (logZ < smallLogZ);
+            double rPlus1 = MMath.GetRPlus1(r, sqrtomr2);
+            double omr2 = sqrtomr2 * sqrtomr2;
+            alphaL = 0.0;
+            if (r == 0 && !Double.IsInfinity(yl))
+            {
+                alphaL = -d_p * invSqrtVxl / MMath.NormalCdfRatio(yl);
+            }
+            else if (!lowerBound.IsUniform() && !Double.IsInfinity(yl))
+            {
+                // since X and lowerBound are not both uniform, invSqrtVxl > 0
+                double logPhiL = 0;
+                if (!useLogZRatio)
+                    logPhiL += Gaussian.GetLogProb(yl, 0, 1);
+                if (r > -1)
+                {
+                    double yuryl;
+                    if (r > smallR)
+                    {
+                        yuryl = MMath.GetXMinusRY(yu, yl, r, omr2) / sqrtomr2;
+                    }
+                    else
+                    {
+                        // yu - r * yl = 
+                        // (mu-mx)*invSqrtVxu + vx*invSqrtVxu*invSqrtVxl*(mx-ml)*invSqrtVxl =
+                        // (mu-mx + (mx-ml)*vx/(vx+vl))*invSqrtVxu =
+                        // (mx*vl/(vx+vl) + mu - ml*vx/(vx+vu))*invSqrtVxl
+                        yuryl = (mx * vl / (vx + vl) + mu - ml * vx / (vx + vl)) / Math.Sqrt(vx * vl + vx * vu + vl * vu) / invSqrtVxl;
+                    }
+                    if (useLogZRatio)
+                    {
+                        logPhiL = MMath.NormalCdfRatioLn(yuryl);
+                    }
+                    else
+                    {
+                        logPhiL += MMath.NormalCdfLn(yuryl);
+                    }
+                    if (logPhiL > double.MaxValue) throw new Exception();
+                    //Trace.WriteLine($"yuryl = {yuryl}, invSqrtVxl = {invSqrtVxl} useLogZRatio = {useLogZRatio}");
+                }
+                alphaL = -d_p * invSqrtVxl * Math.Exp(logPhiL - (useLogZRatio ? logZRatio : logZ));
+            }
+            alphaU = 0.0;
+            if (r == 0 && !Double.IsInfinity(yu))
+            {
+                alphaU = d_p * invSqrtVxu / MMath.NormalCdfRatio(yu);
+            }
+            else if (!upperBound.IsUniform() && !Double.IsInfinity(yu))
+            {
+                // since X and upperBound are not both uniform, invSqrtVxu > 0
+                double logPhiU = 0;
+                if (!useLogZRatio)
+                    logPhiU = Gaussian.GetLogProb(yu, 0, 1);
+                if (r > -1)
+                {
+                    double ylryu;
+                    if (r > smallR)
+                    {
+                        ylryu = MMath.GetXMinusRY(yl, yu, r, omr2) / sqrtomr2;
+                    }
+                    else
+                    {
+                        // yl - r * yu = 
+                        // (mx-ml)*invSqrtVxl + vx*invSqrtVxu*invSqrtVxl*(mu-mx)*invSqrtVxu =
+                        // (mx-ml + (mu-mx)*vx/(vx+vu))*invSqrtVxl =
+                        // (mx*vu/(vx+vu) - ml + mu*vx/(vx+vu))*invSqrtVxl
+                        // Sqrt(1 - r * r) = Sqrt(vx*vl + vx*vu + vl*vu)*invSqrtVxl*invSqrtVxu
+                        ylryu = (mx * vu / (vx + vu) - ml + mu * vx / (vx + vu)) / Math.Sqrt(vx * vl + vx * vu + vl * vu) / invSqrtVxu;
+                    }
+                    if (useLogZRatio)
+                        logPhiU = MMath.NormalCdfRatioLn(ylryu);
+                    else
+                        logPhiU += MMath.NormalCdfLn(ylryu);
+                    //Trace.WriteLine($"ylryu = {ylryu}, invSqrtVxu = {invSqrtVxu}");
+                }
+                alphaU = d_p * invSqrtVxu * Math.Exp(logPhiU - (useLogZRatio ? logZRatio : logZ));
+            }
+            alphaX = -alphaL - alphaU;
+            if (useLogZRatio && Math.Abs(alphaX) < Math.Abs(alphaL) * 1e-6)
+            {
+                // In this regime, we can assume NormalCdfRatio(x) == -1/x (since x < -1e8)
+                // alphaX = -alphaL - alphaU
+                // = d_p * Math.Exp(-logZRatio) * (invSqrtVxl * NormalCdfRatio(yuryl) - invSqrtVxu * NormalCdfRatio(ylryu))
+                // = d_p * Math.Exp(-logZRatio) * (-invSqrtVxl / yuryl + invSqrtVxu / ylryu)
+                //double q = -1 / ((mx + mu) * vl + (mu - ml) * vx) + 1 / ((mx - ml) * vu + (mu - ml) * vx);
+                double q2 = ((ml - mx) * vu + (mx + mu) * vl) / ((mx + mu) * vl + (mu - ml) * vx) / ((mx - ml) * vu + (mu - ml) * vx);
+                alphaX = d_p * Math.Exp(-logZRatio) * Math.Sqrt(vx * vl + vx * vu + vl * vu) * q2;
+            }
+            ylInvSqrtVxlPlusAlphaX = yl * invSqrtVxl + alphaX;
+            yuInvSqrtVxuMinusAlphaX = yu * invSqrtVxu - alphaX;
+            double invSqrtVxlMinusInvSqrtVxu = invSqrtVxl - invSqrtVxu;
+            if (Math.Abs(invSqrtVxlMinusInvSqrtVxu) < invSqrtVxl * 1e-10)
+            {
+                // Use Taylor expansion:
+                // f(y) = 1/sqrt(vx+y)
+                // f'(y) = -0.5/(vx+y)^(3/2)
+                // f(vl) - f(vu) =approx f(vu) + (vl-vu)*f'(vu)
+                invSqrtVxlMinusInvSqrtVxu = (vl - vu) * (-0.5) * (invSqrtVxu * invSqrtVxu * invSqrtVxu);
+            }
+            if (d_p == 1)
+            {
+                double ylInvSqrtVxlPlusAlphaX2 = ylInvSqrtVxlPlusAlphaX;
+                //   yl * invSqrtVxl + alphaX 
+                // = yl * invSqrtVxl - alphaL - alphaU
+                // = invSqrtVxl * (yl + Zx / Z) - invSqrtVxu * Zy / Z
+                // = invSqrtVxl * (yl + Zx / Z - Zy / Z) + (invSqrtVxl - invSqrtVxu) * Zy /Z
+                //   yl + Zx / Z - Zy / Z 
+                // = (yl*Z + Zx - Zy)/Z
+                // = (intZ - Zx - r*Zy + Zx - Zy)/Z
+                // = intZ/Z - (1+r)*Zy/Z
+                double intZOverZ = MMath.NormalCdfIntegralRatio(yl, yu, r, sqrtomr2);
+                if (r == 0)
+                    ylInvSqrtVxlPlusAlphaX = invSqrtVxl * intZOverZ - alphaU;
+                else
+                    ylInvSqrtVxlPlusAlphaX = invSqrtVxl * intZOverZ + (invSqrtVxlMinusInvSqrtVxu - rPlus1 * invSqrtVxl) / invSqrtVxu * alphaU;
+                //ylInvSqrtVxlPlusAlphaX = -invSqrtVxl / yl;
+                if (TraceAlpha)
+                    Trace.WriteLine($"ylInvSqrtVxlPlusAlphaX = {ylInvSqrtVxlPlusAlphaX} replaces {ylInvSqrtVxlPlusAlphaX2} (intZOverZ = {intZOverZ}, alphaU = {alphaU}, r = {r}, yl = {yl})");
+                if (double.IsNaN(ylInvSqrtVxlPlusAlphaX)) throw new Exception("ylInvSqrtVxlPlusAlphaX is NaN");
+                // alphaL = -invSqrtVxl * Zx / Z
+                // yuInvSqrtVxuMinusAlphaX
+                // = yu * invSqrtVxu - alphaX
+                // = yu * invSqrtVxu + alphaU + alphaL
+                // = invSqrtVxu * (yu + Zy / Z) - invSqrtVxl * Zx / Z
+                // = invSqrtVxu * (yu + Zy / Z - Zx / Z) - (invSqrtVxl - invSqrtVxu) * Zx / Z
+                //   yu + Zy / Z - Zx / Z 
+                // = (yu*Z + Zy - Zx)/Z
+                // = (intZ2 - Zy - r*Zx + Zy - Zx)/Z
+                // = intZ2/Z - (1+r)*Zx/Z
+                double yuInvSqrtVxuMinusAlphaX2 = yuInvSqrtVxuMinusAlphaX;
+                double intZ2OverZ = MMath.NormalCdfIntegralRatio(yu, yl, r, sqrtomr2);
+                if (r == 0)
+                    yuInvSqrtVxuMinusAlphaX = invSqrtVxu * intZ2OverZ + alphaL;
+                else
+                    yuInvSqrtVxuMinusAlphaX = invSqrtVxu * intZ2OverZ + (invSqrtVxlMinusInvSqrtVxu + rPlus1 * invSqrtVxu) / invSqrtVxl * alphaL;
+                if (TraceAlpha)
+                    Trace.WriteLine($"yuInvSqrtVxuMinusAlphaX = {yuInvSqrtVxuMinusAlphaX} replaces {yuInvSqrtVxuMinusAlphaX2} (intZ2OverZ = {intZ2OverZ}, alphaL = {alphaL})");
+                //if (Math.Abs(yuInvSqrtVxuMinusAlphaX - yuInvSqrtVxuMinusAlphaX2) > 1e-2) throw new Exception();
+                if (double.IsNaN(yuInvSqrtVxuMinusAlphaX)) throw new Exception("yuInvSqrtVxuMinusAlphaX is NaN");
+            }
+            // TODO: make this case smoothly blend into the X.IsPointMass case
+            if (string.Empty.Length > 0)//(r == -1 || (!Double.IsInfinity(yl) && !Double.IsInfinity(yu) && logZ == MMath.NormalCdfLn(yl, yu, -1)))
+            {
+                if ((yl > 0 && yu < 0) || (yl < 0 && yu > 0))
+                {
+                    Trace.WriteLine("special case for r == -1");
+                    // yu + yl = (mu-mx)/sqrt(vx+vu) + (mx-ml)/sqrt(vx+vl) = mu/sqrt(vx+vu) - ml/sqrt(vx+vl) + mx*(1/sqrt(vx+vl) - 1/sqrt(vx+vu))
+                    // This is more accurate than (yu+yl)
+                    double yuPlusyl = mu * invSqrtVxu - ml * invSqrtVxl + mx * invSqrtVxlMinusInvSqrtVxu;
+                    // double delta = Math.Exp(-0.5 * (yu * yu - yl * yl));
+                    // This is more accurate than above.
+                    double deltaMinus1 = MMath.ExpMinus1(-0.5 * yuPlusyl * (yu - yl));
+                    double delta = 1 + deltaMinus1;
+                    bool deltaMinus1IsSmall = Math.Abs(deltaMinus1) < 1e-8;
+                    double ZoverPhiL, ZoverPhiU;
+                    if (yl > 0 && yu < 0)
+                    {
+                        // This code inlines the following formula for Z:
+                        // Z = NormalCdf(yu) - NormalCdf(-yl)
+                        double ncru = MMath.NormalCdfRatio(yu);
+                        double ncrl = MMath.NormalCdfRatio(-yl);
+                        if (deltaMinus1IsSmall)
+                        {
+                            double ncruMinusNcrl = MMath.NormalCdfRatioDiff(-yl, yuPlusyl);
+                            ZoverPhiL = ncruMinusNcrl + ncru * deltaMinus1;
+                            ZoverPhiU = ZoverPhiL / delta;
+                        }
+                        else
+                        {
+                            ZoverPhiL = ncru * delta - ncrl;
+                            ZoverPhiU = ncru - ncrl / delta;
+                        }
+                    }
+                    else // (yl < 0 && yu > 0)
+                    {
+                        // This code inlines the following formula for Z:
+                        // Z = NormalCdf(yl) - NormalCdf(-yu)
+                        double ncru = MMath.NormalCdfRatio(-yu);
+                        double ncrl = MMath.NormalCdfRatio(yl);
+                        if (deltaMinus1IsSmall)
+                        {
+                            double ncrlMinusNcru = MMath.NormalCdfRatioDiff(-yu, yuPlusyl);
+                            ZoverPhiL = ncrlMinusNcru - ncru * deltaMinus1;
+                            ZoverPhiU = ZoverPhiL / delta;
+                        }
+                        else
+                        {
+                            ZoverPhiL = ncrl - ncru * delta;
+                            ZoverPhiU = ncrl / delta - ncru;
+                        }
+                    }
+                    alphaL = -d_p * invSqrtVxl / ZoverPhiL;
+                    //alphaU = d_p * invSqrtVxu * delta / ZoverPhiL;
+                    alphaU = d_p * invSqrtVxu / ZoverPhiU;
+                    if (deltaMinus1IsSmall)
+                    {
+                        // alphaX = d_p/ZoverPhiL * (invSqrtVxl - invSqrtVxu * delta)
+                        alphaX = d_p / ZoverPhiL * (invSqrtVxlMinusInvSqrtVxu - invSqrtVxu * deltaMinus1);
+                    }
+                    else
+                    {
+                        alphaX = -alphaL - alphaU;
+                    }
+                    if (d_p == 1)
+                    {
+                        //   yl * invSqrtVxl + alphaX 
+                        // = yl * invSqrtVxl - alphaL - alphaU
+                        // = invSqrtVxl * (yl + 1/ZoverPhiL) - invSqrtVxu / ZoverPhiU
+                        // = invSqrtVxl * (yl + 1/ZoverPhiL - 1/ZoverPhiU) + (invSqrtVxl - invSqrtVxu) / ZoverPhiU
+                        // yl + 1/ZoverPhiL - 1/ZoverPhiU = (yl*ZRatio + Rylryu - Ryuryl)/ZRatio
+                        // = (intZRatio - Rylryu - r*Ryuryl + Rylryu - Ryuryl)/ZRatio
+                        // = intZRatio/ZRatio = intZ/Z
+                        double ylInvSqrtVxlPlusAlphaX2 = ylInvSqrtVxlPlusAlphaX;
+                        double intZOverZ = MMath.NormalCdfIntegralRatio(yl, yu, -1);
+                        ylInvSqrtVxlPlusAlphaX = invSqrtVxl * intZOverZ + invSqrtVxlMinusInvSqrtVxu / ZoverPhiU;
+                        Trace.WriteLine($"ylInvSqrtVxlPlusAlphaX = {ylInvSqrtVxlPlusAlphaX} replaces {ylInvSqrtVxlPlusAlphaX2} (intZOverZ = {intZOverZ}, alphaU = {alphaU})");
+                        if (double.IsNaN(ylInvSqrtVxlPlusAlphaX)) throw new Exception("ylInvSqrtVxlPlusAlphaX is NaN");
+                        //   yu * invSqrtVxu - alphaX
+                        // = yu * invSqrtVxu + alphaL + alphaU
+                        // = invSqrtVxu * (yu + 1/ZoverPhiU) - invSqrtVxl / ZoverPhiL
+                        // = invSqrtVxu * (yu + 1/ZoverPhiU - 1/ZoverPhiL) - (invSqrtVxl - invSqrtVxu) / ZoverPhiL
+                        // yu + 1/ZoverPhiU - 1/ZoverPhiL = (yu*ZRatio + Ryuryl - Rylryu)/ZRatio
+                        // = (intZ2Ratio - r*Rylryu - Ryuryl + Ryuryl - Rylryu)/ZRatio
+                        // = intZ2/Z
+                        double yuInvSqrtVxuMinusAlphaX2 = yuInvSqrtVxuMinusAlphaX;
+                        double intZ2OverZ = MMath.NormalCdfIntegralRatio(yu, yl, -1);
+                        yuInvSqrtVxuMinusAlphaX = invSqrtVxu * intZ2OverZ - invSqrtVxlMinusInvSqrtVxu / ZoverPhiL;
+                        Trace.WriteLine($"yuInvSqrtVxuMinusAlphaX = {yuInvSqrtVxuMinusAlphaX} replaces {yuInvSqrtVxuMinusAlphaX2} (intZ2OverZ = {intZ2OverZ}, alphaL = {alphaL})");
+                        if (double.IsNaN(yuInvSqrtVxuMinusAlphaX)) throw new Exception("yuInvSqrtVxuMinusAlphaX is NaN");
+                    }
+                    if (double.IsNaN(alphaL)) throw new Exception("alphaL is NaN");
+                    if (double.IsNaN(alphaU)) throw new Exception("alphaU is NaN");
+                }
+            }
+        }
+
 #if false
     /// <summary>
     /// EP message to 'isBetween'

src/Runtime/Factors/IsPositive.cs

@@ -137,9 +137,7 @@ namespace Microsoft.ML.Probabilistic.Factors
                         return (new Gaussian(mp, vp)) / x;
                     }
                 }
-                if (x.IsUniform())
-                    return Gaussian.Uniform();
-                throw new ImproperMessageException(x);
+                return Gaussian.Uniform();
             }
             else if (prec < 0)
             {

src/Runtime/Factors/JaggedSubarray.cs

@@ -363,7 +363,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             for (int j = 0; j < indices_i_Count; j++)
             {
                 DistributionType value = result[j];
-                value.SetToRatio(marginal[indices_i[j]], items[j]);
+                value.SetToRatio(marginal[indices_i[j]], items[j], GaussianOp.ForceProper);
                 result[j] = value;
             }
             return result;
@@ -904,7 +904,7 @@ namespace Microsoft.ML.Probabilistic.Factors
         /// <typeparam name="ItemType">The type of a sub-array.</typeparam>
         public static ItemType ItemsAverageConditional<ArrayType, DistributionType, ItemType>(
             [Indexed, Cancels] ItemType items, // items dependency must be ignored for Sequential schedule
-            [IgnoreDependency] ArrayType array,
+            ArrayType array,
             [SkipIfAllUniform] ArrayType to_marginal,
             IList<IList<int>> indices,
             int resultIndex,
@@ -922,7 +922,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             for (int j = 0; j < indices_i_Count; j++)
             {
                 DistributionType value = result[j];
-                value.SetToRatio(to_marginal[indices_i[j]], items[j]);
+                value.SetToRatio(to_marginal[indices_i[j]], items[j], GaussianOp.ForceProper);
                 result[j] = value;
             }
             return result;
@@ -936,7 +936,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             ArrayType result)
             where ArrayType : SettableToRatio<ArrayType>
         {
-            result.SetToRatio(to_marginal, array);
+            result.SetToRatio(to_marginal, array, GaussianOp.ForceProper);
             return result;
         }
 

src/Runtime/Factors/PlusDouble.cs

@@ -7,6 +7,7 @@ namespace Microsoft.ML.Probabilistic.Factors
     using System;
     using Microsoft.ML.Probabilistic.Distributions;
     using Microsoft.ML.Probabilistic.Factors.Attributes;
+    using Microsoft.ML.Probabilistic.Math;
 
     /// <include file='FactorDocs.xml' path='factor_docs/message_op_class[@name="DoublePlusOp"]/doc/*'/>
     [FactorMethod(typeof(Factor), "Plus", typeof(double), typeof(double), Default = true)]
@@ -186,8 +187,13 @@ namespace Microsoft.ML.Probabilistic.Factors
             if (b.IsPointMass)
                 return SumAverageConditional(a, b.Point);
             if (a.IsUniform() || b.IsUniform()) return Gaussian.Uniform();
+            if (a.Precision == b.Precision)
+            {
+                // This formula works even when Precision == 0
+                return Gaussian.FromNatural(MMath.Average(a.MeanTimesPrecision, b.MeanTimesPrecision), a.Precision / 2);
+            }
             double prec = a.Precision + b.Precision;
-            if (prec <= 0)
+            if (prec == 0)
                 throw new ImproperDistributionException(a.IsProper() ? b : a);
             return Gaussian.FromNatural((a.MeanTimesPrecision * b.Precision + b.MeanTimesPrecision * a.Precision) / prec, a.Precision * b.Precision / prec);
         }
@@ -222,8 +228,13 @@ namespace Microsoft.ML.Probabilistic.Factors
             if (b.IsPointMass)
                 return AAverageConditional(Sum, b.Point);
             if (Sum.IsUniform() || b.IsUniform()) return Gaussian.Uniform();
+            if (Sum.Precision == b.Precision)
+            {
+                // This formula works even when Precision == 0
+                return Gaussian.FromNatural(MMath.Average(Sum.MeanTimesPrecision, -b.MeanTimesPrecision), Sum.Precision / 2);
+            }
             double prec = Sum.Precision + b.Precision;
-            if (prec <= 0)
+            if (prec == 0)
                 throw new ImproperDistributionException(Sum.IsProper() ? b : Sum);
             return Gaussian.FromNatural((Sum.MeanTimesPrecision * b.Precision - b.MeanTimesPrecision * Sum.Precision) / prec, Sum.Precision * b.Precision / prec);
         }

test/TestApp/TestApp.csproj

@@ -14,6 +14,7 @@
     <RootNamespace>TestApp</RootNamespace>
     <Configurations>Debug;Release;DebugFull;DebugCore;ReleaseFull;ReleaseCore</Configurations>
     <StartupObject>TestApp.Program</StartupObject>
+    <DelaySign>false</DelaySign>
   </PropertyGroup>
   <Choose>
     <When Condition="'$(Configuration)'=='DebugFull' OR '$(Configuration)'=='ReleaseFull'">
@@ -53,6 +54,10 @@
     <Optimize>true</Optimize>
   </PropertyGroup>
 
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='DebugFull|AnyCPU'">
+    <PlatformTarget>AnyCPU</PlatformTarget>
+  </PropertyGroup>
+
   <ItemGroup>
     <ProjectReference Include="..\..\src\Compiler\Compiler.csproj" />
     <ProjectReference Include="..\..\src\Csoft\Csoft.csproj" />

test/Tests/Core/RandomTest.cs

@@ -421,14 +421,13 @@ namespace Microsoft.ML.Probabilistic.Tests
         }
 
         [Fact]
-                public void RandInt_ThrowsIfEqualArguments()
+        public void RandInt_ThrowsIfEqualArguments()
         {
-Assert.Throws<ArgumentException>(() =>
-{
-            Rand.Int(9, 9);
-        
-});
-}
+            Assert.Throws<ArgumentException>(() =>
+            {
+                Rand.Int(9, 9);        
+            });
+        }
 
         [Fact]
         public void MeanVarianceAccumulator_Add_Infinity()
@@ -439,5 +438,16 @@ Assert.Throws<ArgumentException>(() =>
             Assert.True(double.IsPositiveInfinity(mva.Mean));
             Assert.True(double.IsPositiveInfinity(mva.Variance));
         }
+
+        [Fact]
+        public void MeanVarianceAccumulator_Add_ZeroWeight()
+        {
+            MeanVarianceAccumulator mva = new MeanVarianceAccumulator();
+            mva.Add(4.5, 0.0);
+            Assert.Equal(4.5, mva.Mean);
+            mva.Add(4.5);
+            mva.Add(double.PositiveInfinity, 0.0);
+            Assert.Equal(4.5, mva.Mean);
+        }
     }
 }

test/Tests/Core/SpecialFunctionsTest.cs

@@ -5,12 +5,14 @@
 using System;
 using System.Collections.Generic;
 using System.Diagnostics;
+using System.Linq;
 using Xunit;
 using Assert = Microsoft.ML.Probabilistic.Tests.AssertHelper;
 using Microsoft.ML.Probabilistic.Utilities;
 using Microsoft.ML.Probabilistic.Collections;
 using Microsoft.ML.Probabilistic.Distributions;
 using Microsoft.ML.Probabilistic.Math;
+using Microsoft.ML.Probabilistic.Core.Maths;
 
 namespace Microsoft.ML.Probabilistic.Tests
 {
@@ -19,11 +21,10 @@ namespace Microsoft.ML.Probabilistic.Tests
     {
         public const double TOLERANCE = 1e-15;
 
-        public delegate double MathFcn(double x);
-
-        public delegate double MathFcn2(double x1, double x2);
-
-        public delegate double MathFcn3(double x1, double x2, double x3);
+        public delegate double MathFcn(double arg);
+        public delegate double MathFcn2(double arg1, double arg2);
+        public delegate double MathFcn3(double arg1, double arg2, double arg3);
+        public delegate double MathFcn4(double arg1, double arg2, double arg3, double arg4);
 
         public struct Test
         {
@@ -372,7 +373,6 @@ x = mpf("1e-4")
                     {1e-2, 0.0016708416805754216545690286 },
                     {1e-3, 0.000166708341668055753993 },
                     {1e-4, 0.0000166670833416668055575397},
-                    /*
                     {-1e-4, -0.0000166662500083331944464},
                     {-1e-3, -0.000166625008331944642832344 },
                     {-1e-2, -0.0016625083194642609402281996344222792 },
@@ -386,7 +386,6 @@ x = mpf("1e-4")
                     {-6e-1, -0.08663434418325990936539 },
                     {-1, -0.132120558828557678404476229838539},
                     {Double.NegativeInfinity, -0.5},
-                    */
                     {Double.PositiveInfinity, Double.PositiveInfinity},
                     {Double.NaN, Double.NaN}
                 };
@@ -1036,70 +1035,10 @@ gammainc(mpf('1'),mpf('1'),mpf('inf'),regularized=True)
             CheckFunctionValues("GammaUpper", MMath.GammaUpper, gammaUpper_pairs);
         }
 
-        private static void CheckFunctionValues(string name, MathFcn fcn, double[,] pairs, double assertTolerance = 1e-11)
-        {
-            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
-        }
-
-        private static void CheckFunctionValues(string name, MathFcn2 fcn, double[,] pairs, double assertTolerance = 1e-11)
-        {
-            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
-        }
-
-        private static void CheckFunctionValues(string name, MathFcn3 fcn, double[,] pairs, double assertTolerance = 1e-11)
-        {
-            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
-        }
-
-        private static void CheckFunctionValues(string name, Delegate fcn, double[,] pairs, double assertTolerance)
-        {
-            Vector x = Vector.Zero(pairs.GetLength(1) - 1);
-            object[] args = new object[x.Count];
-            for (int i = 0; i < pairs.GetLength(0); i++)
-            {
-                for (int k = 0; k < x.Count; k++)
-                {
-                    x[k] = pairs[i, k];
-                    args[k] = x[k];
-                }
-                double fx = pairs[i, x.Count];
-                double result = (double)Util.DynamicInvoke(fcn, args);
-                if (!double.IsNaN(result) && System.Math.Sign(result) != System.Math.Sign(fx) && fx != 0)
-                {
-                    string strMsg = string.Format("{0}({1})\t has wrong sign (result = {2})", name, x.ToString("r"), result);
-                    Console.WriteLine(strMsg);
-                    Assert.True(false, strMsg);
-                }
-                double err = System.Math.Abs(result - fx);
-                if (fx != 0)
-                    err /= System.Math.Abs(fx);
-                if (Double.IsInfinity(fx))
-                {
-                    err = (result == fx) ? 0 : Double.PositiveInfinity;
-                }
-                if (Double.IsNaN(fx))
-                {
-                    err = Double.IsNaN(result) ? 0 : Double.NaN;
-                }
-                if (err < TOLERANCE)
-                {
-                    Console.WriteLine("{0}({1})\t ok", name, x.ToString("r"));
-                }
-                else
-                {
-                    string strMsg = String.Format("{0}({1})\t wrong by {2}", name, x.ToString("r"), err.ToString("g2"));
-                    Console.WriteLine(strMsg);
-                    //Console.WriteLine("{0}({1}) = {2} should be {3}", name, x.ToString("r"), result.ToString("r"), fx.ToString("r"));
-                    if (err > assertTolerance || double.IsNaN(err))
-                        Assert.True(false, strMsg);
-                }
-            }
-        }
-
         [Fact]
         public void NormalCdfRatioTest()
         {
-            double r0 = MMath.Sqrt2PI/2;
+            double r0 = MMath.Sqrt2PI / 2;
             Assert.Equal(r0, MMath.NormalCdfRatio(0));
             Assert.Equal(r0, MMath.NormalCdfRatio(6.6243372842224754E-170));
             Assert.Equal(r0, MMath.NormalCdfRatio(-6.6243372842224754E-170));
@@ -1829,21 +1768,23 @@ exp(x*x/4)*pcfu(0.5+n,-x)
         }
 
         /// <summary>
-        /// Used to debug MMath.NormalCdfLn
+        /// Used to tune MMath.NormalCdfLn.  The best tuning minimizes the number of messages printed.
         /// </summary>
         internal void NormalCdf2Test2()
         {
+            // Call both routines now to speed up later calls.
             MMath.NormalCdf(-2, -2, -0.5);
             NormalCdf_Quadrature(-2, -2, -0.5);
             Stopwatch watch = new Stopwatch();
-            double xmin = -1;
-            double xmax = 0;
-            double n = 10;
+            double xmin = 0.1;
+            double xmax = 0.1;
+            double n = 20;
             double xinc = (xmax - xmin) / (n - 1);
             for (int xi = 0; xi < n; xi++)
             {
+                if (xinc == 0 && xi > 0) break;
                 double x = xmin + xi * xinc;
-                double ymin = x;
+                double ymin = -System.Math.Abs(x)*10;
                 double ymax = -ymin;
                 double yinc = (ymax - ymin) / (n - 1);
                 for (int yi = 0; yi < n; yi++)
@@ -1851,8 +1792,14 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                     double y = ymin + yi * yinc;
                     double rmin = -0.999999;
                     double rmax = -0.000001;
-                    rmin = -0.6;
-                    rmax = -0.5;
+                    rmin = MMath.NextDouble(-1);
+                    rmax = -1 + 1e-6;
+                    //rmax = -0.5;
+                    //rmax = 0.1;
+                    //rmax = -0.58;
+                    //rmax = -0.9;
+                    //rmin = -0.5;
+                    rmax = 1;
                     double rinc = (rmax - rmin) / (n - 1);
                     for (int ri = 0; ri < n; ri++)
                     {
@@ -1862,7 +1809,7 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                         double result1 = double.NaN;
                         try
                         {
-                            result1 = MMath.NormalCdfLn(x, y, r);
+                            result1 = MMath.NormalCdfIntegral(x, y, r);
                         }
                         catch
                         {
@@ -1881,8 +1828,9 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                         {
                         }
                         long ticks2 = watch.ElapsedTicks;
-                        if (double.IsNaN(result1) || ticks > 10 * ticks2)
-                            Trace.WriteLine($"({x},{y},{r}): {good} {ticks} {ticks2} {result1} {result2}");
+                        bool overtime = ticks > 10 * ticks2;
+                        if (double.IsNaN(result1) /*|| overtime*/)
+                            Trace.WriteLine($"({x:r},{y:r},{r:r},{x-r*y}): {good} {ticks} {ticks2} {result1} {result2}");
                     }
                 }
             }
@@ -1905,7 +1853,7 @@ exp(x*x/4)*pcfu(0.5+n,-x)
         }
 
         // Used to debug MMath.NormalCdf
-        internal static void NormalCdf2Test3()
+        internal void NormalCdf2Test3()
         {
             double x, y, r;
             bool first = true;
@@ -1935,17 +1883,27 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                 x = -1.5;
                 y = 1.5;
                 r = -0.49;
+
+                x = -1.6450031341281908;
+                y = 1.2645625117080999;
+                r = -0.054054238344620031;
+
+                x = -0.5;
+                y = -0.5;
+                r = 0.001;
+
                 Console.WriteLine(1 - r * r);
 
                 Console.WriteLine("NormalCdfBrute: {0}", NormalCdfBrute(0, x, y, r));
                 Console.WriteLine("NormalCdf_Quadrature: {0}", NormalCdf_Quadrature(x, y, r));
                 //Console.WriteLine("{0}", NormalCdfAlt2(x, y, r));
                 //Console.WriteLine("NormalCdfAlt: {0}", NormalCdfAlt(x, y, r));
-                //Console.WriteLine("NormalCdfTaylor: {0}", NormalCdfTaylor(x, y, r));
-                Console.WriteLine("NormalCdfConFrac3: {0}", NormalCdfConFrac3(x, y, r));
+                //Console.WriteLine("NormalCdfTaylor: {0}", MMath.NormalCdfRatioTaylor(x, y, r));
+                //Console.WriteLine("NormalCdfConFrac3: {0}", NormalCdfConFrac3(x, y, r));
                 //Console.WriteLine("NormalCdfConFrac4: {0}", NormalCdfConFrac4(x, y, r));
                 //Console.WriteLine("NormalCdfConFrac5: {0}", NormalCdfConFrac5(x, y, r));
                 Console.WriteLine("MMath.NormalCdf: {0}", MMath.NormalCdf(x, y, r));
+                Console.WriteLine("MMath.NormalCdfLn: {0}", MMath.NormalCdfLn(x, y, r));
                 for (int i = 1; i < 50; i++)
                 {
                     //Console.WriteLine("{0}: {1}", i, NormalCdfBrute(i, x, y, r));
@@ -1977,10 +1935,21 @@ exp(x*x/4)*pcfu(0.5+n,-x)
             double r2 = -0.95, y2 = 100;
             double xy1 = System.Math.Max(0.0, MMath.NormalCdf(x) + MMath.NormalCdf(y) - 1);
 
+            Assert.True(0 < MMath.NormalCdf(6.8419544775976187E-08, -5.2647906596206016E-08, -1, 3.1873689658872377E-10).Mantissa);
+
             // In sage:
             // integral(1/(2*pi*sqrt(1-t*t))*exp(-(x*x+y*y-2*t*x*y)/(2*(1-t*t))),t,-1,r).n(digits=200);
             double[,] normalcdf2_pairs = new double[,]
                 {
+                    { -0.025562923054653991, 0.015700461401932524, -0.99999824265582826, 9.5694701738831278E-12 },
+                    { -1.9315113095095533, 1.8359512109009408, 0.52222045656167426, 0.026703981978683056 },
+                    { -16, -16.8, 0.5, 0 },
+                    { -1E-19, 1.1000000000000002E-19, -1, 0.398942280401433e-20 },
+                    { -20, 21, -1, 2.75362411532694E-89 },
+                    { -9945842.6290678252, 9945822.06800303, -0.98919581102480514, 0 },
+                    { -312498.36862450332, 312498.298221121, -0.999989333908269, 0 },
+                    { -1.143307502481119, 0.49983795334601688, -0.65, 0.0306849096621242 },
+                    { -1.143307502481119, 0.64555140288768986, -0.53867143302788545, 0.049792920736062646 },
                     {-1.5, 1.5, -0.49, 0.04889082718786 },
                     {-1.5, 1.5, -0.50, 0.048484159276754 },
                     {-1.5, 1.5, -0.51, 0.0480706563619095 },
@@ -1995,14 +1964,18 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                     {-10,-0.5,-0.5, 5.94859187434255E-34},
                     {10,1,0.5, System.Math.Exp(-0.17275377902345) },
                     {-17.16232660642066, 13.435028842544401, -0.67419986246324193, 2.5244995486068114E-66 },
-                    {-2,1.3,-0.5, 0.0125543619945072},
+                    { -1.5235871609407412, 0.069296909625515962, -0.027053117898400668, 0.0323181622672724 },
+                    { -1.64500313412819, 1.2645625117081, -0.05405423834462, 0.0437792307578264 },
+                    { -2,1.3,-0.5, 0.0125543619945072},
                     {-2,1.5,-0.5, 0.0145466302871907},
+                    {-3,3,-0.404040595959596, 0.0013030754478730753 },
                     {x, y, 0, MMath.NormalCdf(x)*MMath.NormalCdf(y)},
                     {x, y, 1, MMath.NormalCdf(System.Math.Min(x, y))},
                     {x, y, -1, xy1},
                     {0.6, -1, -1, 0.0},
                     {-1, 0.6, -1, 0.0},
                     {-1, -1, -1, 0.0},
+                    {-0.5,0.5,-0.4, 0.160231262969512 },
                     {x, y, 0.95,0.0817971046701358 + xy1 },
                     {x,y,-0.95, 2.08560436479985e-17 + xy1},
                     {x, y, 0.9, 0.0726860333050785 + xy1},
@@ -2032,6 +2005,13 @@ exp(x*x/4)*pcfu(0.5+n,-x)
             // log(integrate(1/(2*pi*sqrt(1-t*t))*exp(-(2*1.5*1.5 + 2*t*1.5*1.5)/(2*(1-t*t))),t,-1,-0.5))
             double[,] normalcdfln2_pairs = new double[,]
                 {
+                    { -20, -21, 0.5, -287.59599647373 },
+                    { -14, -14.700000000000001, 0.5, -143.75347547431542 },
+                    { -1E-19, 1.1000000000000002E-19, -1, System.Math.Log(0.398942280401433e-20) },
+                    { -20, 21, -1, System.Math.Log(2.75362411532694E-89) },
+                    { -9945842.6290678252, 9945822.06800303, -0.98919581102480514, -49459892801108.422 },
+                    { -312498.36862450332, 312498.298221121, -0.999989333908269, -48827615210.059265 },
+                    { -1.5235871609407412, 0.069296909625515962, -0.027053117898400668, -3.43212590748818 },
                     { -63, 63, -0.4637494637494638, -1989.56232505372569173398 },
                     {-0.1,-0.1, -0.999999, -10018.30269574383335648},
                     {-0.1,-0.1, -0.99999, -1014.8501635587853745},
@@ -2049,9 +2029,16 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                     {-1, 0.6, -1, double.NegativeInfinity},
                     {-1, -1, -1, double.NegativeInfinity},
                     {1, y2, r2, System.Math.Log(0.841344746068543)},
+                    {-0.01,0.00175879396984925,-0.590954773869347, -1.9123400906813246 },
+                    {-0.5,0.5,-0.4, -1.83113711376467 },
+                    {-0.5,0.464824120603015,-0.523618090452261, -1.990034060610423},
+                    {-0.6,0.6,-0.4, -1.8931285572857055 },
+                    {-0.6,0.6,-0.488442211055276, -1.977994514649918 },
+                    {-0.7,0.7,-0.456281407035176, -2.0179010527218026 },
                     {-1.999999999999, 0, -9.9999999999950042E-07, -4.47633340779343 },
                     {-2.9,-2.9,-0.49, System.Math.Log(3.55239065152386E-10)},
                     {-2, -2, -0.9, -47.0355622406161041990406706314832},
+                    {-2,2,-0.484848515151515, -3.9692750279327953 },
                     {-3, -3.5, -0.72, -44.175239867478986},
                     {-4.5, -4.5, -0.6, -57.1114368213654},
                     //{ 8,y2,r2, 0 },
@@ -2064,6 +2051,314 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                     {x, double.PositiveInfinity, 0.5, MMath.NormalCdfLn(x)}
                 };
             CheckFunctionValues("NormalCdfLn2", new MathFcn3(MMath.NormalCdfLn), normalcdfln2_pairs);
+
+            double[,] normalcdfRatioLn2_pairs = new double[,]
+                {
+                    { 379473319.22020543, -379473319.22020543, -0.99999999999999978, 4.0991122317442558E-10, -21.61306377550903 },
+                    { -9945842.6290678252, 9945822.06800303, -0.98919581102480514, 0.14660029826360518, 268535429229.67789 },
+                    { -312498.36862450332, 312498.298221121, -0.999989333908269, 0.0046186653587789346, 249505.2149799816 },
+                };
+            CheckFunctionValues("NormalCdfRatioLn2", new MathFcn4(MMath.NormalCdfRatioLn), normalcdfRatioLn2_pairs);
+
+            double[,] normalcdfIntegral_pairs = new double[,]
+                {
+                    { -10, -100000000000, -0.1, 0 },
+                    { 1, -100, -0.010000000000000002, 0 },
+                    { 4.3993699214502984E-08, 2.1162549143007975E-09, -1, 4.2410115039363254E-16 },
+                    { 54.144028629649441, 16.352706352205796, -1, 54.144028629649441 },
+                    { double.PositiveInfinity, -57.9771953893382, -1, double.PositiveInfinity },
+                    { -x, double.PositiveInfinity, -1, MMath.NormalCdfMomentRatio(1,-x)*System.Math.Exp(Gaussian.GetLogProb(-x,0,1)) },
+                    { 1E-08, 1E-08, -0.999999999999999, 0 },
+                    { -1, -8.9473684210526319, -0.999999999999999, 0 },
+                    { 19.073484197181429, -19.073488459566978, -0.99999999999996048, 4.37349630514655E-147 },
+                    { -0.4999, 0.5, -0.9999, 1.7769677765658425E-05 },
+                    { 0.021034851174404436, -0.37961242087533614, -0.999999997317639, 0 },
+                    { -0.013170888687821042, 0.013170891631143039, -1, 1.7278974097756908E-18 },
+                    { -2, -2, 1, 0.008490702616829637 },
+                    { 2, 2, 1, 2.0084907026168297 },
+                    { 2, 1.9, 1, 2.0081826951426729 },
+                    { 0.73080682429779031, -0.57951293845369067, -0.99999152863368379, 0.0037455902014599447 },
+                    { -0.94102098773740084, -1.2461486442846208, 0.5240076921033775, 0.035369263379357981 },
+                    { 0.89626183425208061, -0.94178051490858161, -0.77953292732402091, 0.03003800734061967 },
+                    { 1.144658872358864, -1.2215859551105948, -0.81617290357740435, 0.017564280718115784 },
+                    { double.PositiveInfinity, -0.40225579098340325, -0.4697418841876283, double.PositiveInfinity },
+                    { -1.08E+31, 790.80368892437889, -0.94587440643473975, 0 },
+                    { -1.081776354231671E+31, 790.80368892437889, -0.94587440643473975, 0 },
+                    { 790.80368892437889, -108177636171.16806, -0.94587440643473975, 0 },
+                    { double.PositiveInfinity, 1, -0.5, double.PositiveInfinity },
+                    { double.NegativeInfinity, 0.87287156094396956, -0.87287156094396956, 0 },
+                    { -0.1, 0.5, -0.1, 0.22610911461618027 }, // 0.233962568815606
+                    { -0.33333333333333331, -1.5, 0.16666666666666666, 0.02599549974965194 }, // 0.022270228761468
+                    { -0.33333333333333331, 1.5, -0.16666666666666666, 0.22824061454942 }, // 0.2319658855376
+                    { 1003, -1001, 0, 0 },
+                    { -1001, double.PositiveInfinity, 0, MMath.NormalCdfMomentRatio(1, -1001)*System.Math.Exp(Gaussian.GetLogProb(-1001,0,1)) },
+                    { -0.5, double.PositiveInfinity, 0.1, MMath.NormalCdfMomentRatio(1, -0.5)*System.Math.Exp(Gaussian.GetLogProb(-0.5,0,1)) },
+                    { 0.5, 0.1, 0.9999, 0.666826770860778 },
+                    { -2499147.006377392, 2499147.273918618, -1, 0 },
+                    { -8, 8.4, -1, 6.3985926111712448E-17 }, // confrac converges rapidly
+                    { 50, 0.5, -1, 34.221057736936238 }, // 461
+                    { 5, 0.5, -1, 3.1052470330674216 }, // 77
+                    { 0.5, 0.5, -1, 0.19146246127401312 }, // 29
+                    { -0.4999, 0.5, -1, 1.7603559714981849E-09 },
+                    { -0.499, 0.5, -1, 1.7606199115326087E-07 },
+                    { -0.49, 0.5, -1, 1.7632494692391017E-05 },
+                    { 0.1, 0.5, -0.9, 0.10185469947014977 },
+                    { 0.1, 0.5, -0.9999, 0.068051457100492874 },
+                    { 0.1, 0.5, -1, 0.06801625056781653 },
+                    { 0.1, -0.5, -1, 0 },
+                    { 0.1, -0.5, -0.9999, 2.5567149109314317E-183 },
+                    { 0.1, -0.5, 0.9999, 0.38288387410422181 },
+                    { 0.5, 0.1, -1, 0.0707579285628088 },
+                    { 0.5, 0.1, -0.9999, 0.0707976238175565 },
+                    { -0.5, 0.1, -1, 0 },
+                    { -0.5, 0.1, -0.9999, 2.5563333094766041E-183 },
+                    { -0.1, 0.5, -1, 0.0297237583130139 },
+                    { -0.1, 0.5, -0.9999, 0.029758964845705 },
+                    { -0.1, -0.5, 0.9999, 0.32117636635902441 },
+                    { -0.1, -1, 0.9999, 0.22608100205354575 },
+                };
+            CheckFunctionValues("NormalCdfIntegral", new MathFcn3(MMath.NormalCdfIntegral), normalcdfIntegral_pairs);
+
+            double[,] normalcdfIntegralRatio_pairs = new double[,]
+                {
+                    { -39062.492380206008, 39062.501110681893, -0.99999983334056686, 2.5600004960154716E-05 },
+                    { -2499147.006377392, 2499147.273918618, -1, 4.001365255616626E-07 }, // 4.0013422697176734E-07 4.0009884107980741E-07
+                    { -0.12449767319913681, double.PositiveInfinity, -0.99757523672528092, 0.75429721821181472 },
+                    { -213393529.2046707, 213393529.2046707, -1, 0 },
+                    { -824.43680216388009, -23300.713731480908, -0.99915764591723821, 6.9859450855259114E-08 },
+                    { 790.80368892437889, -1081777102.2326407, -0.94587440643473975, 1.0293108592794333E-10 },
+                    { 790.80368892437889, -1081776354979.6719, -0.94587440643473975, 1.0293107755790882E-13 },
+                    { -40, 42, -1, 0.0249688472072654 },
+                    { 0.1, 0.5, -0.9, 0.426708959089942 },
+                    { 0.1, 0.5, -0.9999, 0.29422529836665412 },
+                    { double.NegativeInfinity, 0, -0.1, 0 },
+                    { 0, double.NegativeInfinity, -0.1, 0 },
+                };
+            CheckFunctionValues("NormalCdfIntegralRatio", new MathFcn3(MMath.NormalCdfIntegralRatio), normalcdfIntegralRatio_pairs);
+        }
+
+        [Fact]
+        public void NormalCdfIntegralTest()
+        {
+            Assert.True(0 <= MMath.NormalCdfIntegral(1.9018718309533485E+77, -1.9018718309533485E+77, -1, 8.17880416082724E-79).Mantissa);
+            Assert.True(0 <= MMath.NormalCdfIntegral(213393529.2046707, -213393529.2046707, -1, 7.2893668811495072E-10).Mantissa);
+            Assert.True(0 < MMath.NormalCdfIntegral(-0.42146853220760722, 0.42146843802130329, -0.99999999999999989, 6.2292398855983019E-09).Mantissa);
+ 
+            foreach (var x in OperatorTests.Doubles())
+            {
+                foreach (var y in OperatorTests.Doubles())
+                {
+                    foreach (var r in OperatorTests.Doubles().Where(d => d >= -1 && d <= 1))
+                    {
+                        MMath.NormalCdfIntegral(x, y, r);
+                    }
+                }
+            }
+        }
+
+        internal void NormalCdfIntegralTest2()
+        {
+            double x = 0.0093132267868981222;
+            double y = -0.0093132247056551785;
+            double r = -1;
+            y = -2499147.006377392;
+            x = 2499147.273918618;
+            //MMath.TraceConFrac = true;
+            //MMath.TraceConFrac2 = true;
+            for (int i = 0; i < 100; i++)
+            {
+                //x = 2.1 * (i + 1);
+                //y = -2 * (i + 1);
+                //x = -2 * (i + 1);
+                //y = 2.1 * (i + 1);
+                //x = -System.Math.Pow(10, -i);
+                //y = -x * 1.1;
+                x = -0.33333333333333331;
+                y = -1.5;
+                r = 0.16666666666666666;
+                x = -0.4999;
+                y = 0.5;
+                x = -0.1;
+                y = 0.5;
+                r = -0.1;
+
+                x = -824.43680216388009;
+                y = -23300.713731480908;
+                r = -0.99915764591723821;
+                x = -0.94102098773740084;
+                x = 1 + i * 0.01;
+                y = 2;
+                r = 1;
+
+                x = 0.021034851174404436;
+                y = -0.37961242087533614;
+                //x = -0.02;
+                //y += -1;
+                //x -= -1;
+                r = -1 + System.Math.Pow(10, -i);
+
+                //x = i * 0.01;
+                //y = -1;
+                //r = -1 + 1e-8;
+
+                // 1.81377005549484E-40 with exponent
+                // flipped is 1.70330340479022E-40
+                //x = -1;
+                //y = -8.9473684210526319;
+                //x = System.Math.Pow(10, -i);
+                //y = x;
+                //r = -0.999999999999999;
+
+                //x = -0.94102098773740084;
+                //y = -1.2461486442846208;
+                //r = 0.5240076921033775;
+
+                x = 790.80368892437889;
+                y = -1081776354979.6719;
+                y = -System.Math.Pow(10, i);
+                r = -0.94587440643473975;
+
+                x = -39062.492380206008;
+                y = 39062.501110681893;
+                r = -0.99999983334056686;
+
+                //x = -2;
+                //y = 1.5789473684210522;
+                //r = -0.78947368421052622;
+
+                //x = -1.1;
+                //y = -1.1;
+                //r = 0.052631578947368474;
+
+                //x = 0.001;
+                //y = -0.0016842105263157896;
+                //r = -0.4;
+
+                //x = 0.1;
+                //x = 2000;
+                //y = -2000;
+                //r = -0.99999999999999989;
+
+                x = double.MinValue;
+                y = double.MinValue;
+                r = 0.1;
+
+
+                Trace.WriteLine($"(x,y,r) = {x:r}, {y:r}, {r:r}");
+
+                double intZOverZ;
+                try
+                {
+                    intZOverZ = MMath.NormalCdfIntegralRatio(x, y, r);
+                }
+                catch
+                {
+                    intZOverZ = double.NaN;
+                }
+                Trace.WriteLine($"intZOverZ = {intZOverZ:r}");
+
+                double intZ0 = NormalCdfIntegralBasic(x, y, r);
+                double intZ1 = 0; // NormalCdfIntegralFlip(x, y, r);
+                double intZr = 0;// NormalCdfIntegralBasic2(x, y, r);
+                ExtendedDouble intZ;
+                double sqrtomr2 = System.Math.Sqrt((1 - r) * (1 + r));
+                try
+                {
+                    intZ = MMath.NormalCdfIntegral(x, y, r, sqrtomr2);
+                }
+                catch
+                {
+                    intZ = ExtendedDouble.NaN();
+                }
+                //double intZ = intZ0;
+                Trace.WriteLine($"intZ = {intZ:r} {intZ.ToDouble():r} {intZ0:r} {intZ1:r} {intZr:r}");
+                if (intZ.Mantissa < 0) throw new Exception();
+                //double intZ2 = NormalCdfIntegralBasic(y, x, r);
+                //Trace.WriteLine($"intZ2 = {intZ2} {r*intZ}");
+                double Z = MMath.NormalCdf(x, y, r);
+                if (Z < 0) throw new Exception();
+            }
+        }
+
+        private double NormalCdfIntegralFlip(double x, double y, double r)
+        {
+            double logProbX = Gaussian.GetLogProb(x, 0, 1);
+            return -MMath.NormalCdfIntegral(x, -y, -r) + x * MMath.NormalCdf(x) + System.Math.Exp(logProbX);
+        }
+
+        private double NormalCdfIntegralTaylor(double x, double y, double r)
+        {
+            double omr2 = 1 - r * r;
+            double sqrtomr2 = System.Math.Sqrt(omr2);
+            double ymrx = y / sqrtomr2;
+            double dx0 = MMath.NormalCdf(0, y, r);
+            double ddx0 = System.Math.Exp(Gaussian.GetLogProb(0, 0, 1) + MMath.NormalCdfLn(ymrx));
+            // \phi_{xx} &= -x \phi_x - r \phi_r
+            double dddx0 = -r * System.Math.Exp(Gaussian.GetLogProb(0, 0, 1) + Gaussian.GetLogProb(ymrx, 0, 1));
+            Trace.WriteLine($"dx0 = {dx0} {ddx0} {dddx0}");
+            return MMath.NormalCdfIntegral(0, y, r) + x * dx0 + 0.5 * x * x * ddx0 + 1.0 / 6 * x * x * x * dddx0;
+        }
+
+        private double NormalCdfIntegralBasic2(double x, double y, double r)
+        {
+            double omr2 = 1 - r * r;
+            double sqrtomr2 = System.Math.Sqrt(omr2);
+            double ymrx = (y - r * x) / sqrtomr2;
+            double xmry = (x - r * y) / sqrtomr2;
+            double func(double t)
+            {
+                return (y - r * x + r * t) * System.Math.Exp(Gaussian.GetLogProb(t, x, 1) + MMath.NormalCdfLn(ymrx + r * t / sqrtomr2));
+            }
+            func(0);
+            double func2(double t)
+            {
+                double ymrxt = ymrx + r * t / sqrtomr2;
+                return sqrtomr2 * System.Math.Exp(Gaussian.GetLogProb(t, x, 1) + Gaussian.GetLogProb(ymrxt, 0, 1)) * (MMath.NormalCdfMomentRatio(1, ymrxt) - 1);
+            }
+            func2(0);
+            //return -MMath.NormalCdf(x, y, r) * (y / r - x) + Integrate(func2) / r;
+            double func3(double t)
+            {
+                double xmryt = xmry + r * t / sqrtomr2;
+                return sqrtomr2 * System.Math.Exp(Gaussian.GetLogProb(t, y, 1) + Gaussian.GetLogProb(xmryt, 0, 1)) * MMath.NormalCdfMomentRatio(1, xmryt);
+            }
+            //double Z = MMath.NormalCdf(x, y, r, out double exponent);
+            double Z3 = Integrate(func3);
+            //return System.Math.Exp(exponent)*(-Z * (y / r - x) - omr2 / r * MMath.NormalCdfRatio(xmry)) + Z3/r;
+            return Z3;
+        }
+
+        private static double Integrate(Func<double, double> func)
+        {
+            double sum = 0;
+            var ts = EpTests.linspace(0, 1, 100000);
+            double inc = ts[1] - ts[0];
+            for (int i = 0; i < ts.Length; i++)
+            {
+                double t = ts[i];
+                double term = func(t);
+                if (i == 0 || i == ts.Length - 1) term /= 2;
+                sum += term * inc;
+            }
+            return sum;
+        }
+
+        private double NormalCdfIntegralBasic(double x, double y, double r)
+        {
+            double omr2 = 1 - r * r;
+            double sqrtomr2 = System.Math.Sqrt(omr2);
+            double ymrx = (y - r * x) / sqrtomr2;
+            double xmry = (x - r * y) / sqrtomr2;
+            // should use this whenever x > 0 and Rymrx >= Rxmry (y-r*x >= x-r*y implies y*(1+r) >= x*(1+r) therefore y >= x)
+            // we need a special routine to compute 2nd half without cancellation and without dividing by phir
+            // what about x > y > 0?
+            //double t = MMath.NormalCdfIntegral(-x, y, -r) + x * MMath.NormalCdf(y) + r * System.Math.Exp(Gaussian.GetLogProb(y, 0, 1));
+            //Console.WriteLine(t);
+            double phix = System.Math.Exp(Gaussian.GetLogProb(x, 0, 1) + MMath.NormalCdfLn(ymrx));
+            double phiy = System.Math.Exp(Gaussian.GetLogProb(y, 0, 1) + MMath.NormalCdfLn(xmry));
+            //Trace.WriteLine($"phix = {phix} phiy = {phiy}");
+            return x * MMath.NormalCdf(x, y, r) + phix + r * phiy;
+            //return y * MMath.NormalCdf(x, y, r) + r * System.Math.Exp(Gaussian.GetLogProb(x, 0, 1) + MMath.NormalCdfLn(ymrx)) + System.Math.Exp(Gaussian.GetLogProb(y, 0, 1) + MMath.NormalCdfLn(xmry));
         }
 
         [Fact]
@@ -2496,10 +2791,10 @@ exp(x*x/4)*pcfu(0.5+n,-x)
                 {
                     //Console.WriteLine("denom/dfact = {0}", denom / dfact);
                     double result = -numer / denom;
-                    //Console.WriteLine("iter {0}: {1}", i, result);
+                    Console.WriteLine("iter {0}: {1}", i, result);
                     if (double.IsInfinity(result) || double.IsNaN(result))
                         throw new Exception();
-                    if (result == rprev)
+                    if (MMath.AreEqual(result, rprev))
                         return result;
                     rprev = result;
                 }
@@ -23607,5 +23902,69 @@ exp(x*x/4)*pcfu(0.5+n,-x)
             Assert.Equal(double.MaxValue, MMath.Median(new[] { double.MaxValue, double.MaxValue }));
             Assert.Equal(3, MMath.Median(new[] { double.MaxValue, 3, double.MinValue }));
         }
+
+        private static void CheckFunctionValues(string name, MathFcn fcn, double[,] pairs, double assertTolerance = 1e-11)
+        {
+            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
+        }
+
+        private static void CheckFunctionValues(string name, MathFcn2 fcn, double[,] pairs, double assertTolerance = 1e-11)
+        {
+            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
+        }
+
+        private static void CheckFunctionValues(string name, MathFcn3 fcn, double[,] pairs, double assertTolerance = 1e-11)
+        {
+            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
+        }
+
+        private static void CheckFunctionValues(string name, MathFcn4 fcn, double[,] pairs, double assertTolerance = 1e-11)
+        {
+            CheckFunctionValues(name, (Delegate)fcn, pairs, assertTolerance);
+        }
+
+        private static void CheckFunctionValues(string name, Delegate fcn, double[,] pairs, double assertTolerance)
+        {
+            Vector x = Vector.Zero(pairs.GetLength(1) - 1);
+            object[] args = new object[x.Count];
+            for (int i = 0; i < pairs.GetLength(0); i++)
+            {
+                for (int k = 0; k < x.Count; k++)
+                {
+                    x[k] = pairs[i, k];
+                    args[k] = x[k];
+                }
+                double fx = pairs[i, x.Count];
+                double result = (double)Util.DynamicInvoke(fcn, args);
+                if (!double.IsNaN(result) && System.Math.Sign(result) != System.Math.Sign(fx) && fx != 0)
+                {
+                    string strMsg = $"{name}({x:r})\t has wrong sign (result = {result:r})";
+                    Trace.WriteLine(strMsg);
+                    Assert.True(false, strMsg);
+                }
+                double err = System.Math.Abs(result - fx);
+                if (fx != 0)
+                    err /= System.Math.Abs(fx);
+                if (Double.IsInfinity(fx))
+                {
+                    err = (result == fx) ? 0 : Double.PositiveInfinity;
+                }
+                if (Double.IsNaN(fx))
+                {
+                    err = Double.IsNaN(result) ? 0 : Double.NaN;
+                }
+                if (err < TOLERANCE)
+                {
+                    Trace.WriteLine($"{name}({x:r})\t ok");
+                }
+                else
+                {
+                    string strMsg = $"{name}({x:r})\t wrong by {err.ToString("g2")} (result = {result:r})";
+                    Trace.WriteLine(strMsg);
+                    if (err > assertTolerance || double.IsNaN(err))
+                        Assert.True(false, strMsg);
+                }
+            }
+        }
     }
 }

test/Tests/Operators/GateOpTests.cs

@@ -21,10 +21,10 @@ namespace Microsoft.ML.Probabilistic.Tests
         [Fact]
         public void BernoulliEnterTest()
         {
-            this.DoBernoulliEnterTest(0.3, 0.8, 0.4);
-            this.DoBernoulliEnterTest(0.7, 0.8, 0.2);
-            this.DoBernoulliEnterTest(0.7, 0.8, 0.0);
-            this.DoBernoulliEnterTest(0.7, 0.8, 1.0);
+            DoBernoulliEnterTest(0.3, 0.8, 0.4);
+            DoBernoulliEnterTest(0.7, 0.8, 0.2);
+            DoBernoulliEnterTest(0.7, 0.8, 0.0);
+            DoBernoulliEnterTest(0.7, 0.8, 1.0);
         }
 
         /// <summary>

test/Tests/Operators/GaussianOpTests.cs

@@ -23,7 +23,7 @@ namespace Microsoft.ML.Probabilistic.Tests
         {
             Gaussian uniform = new Gaussian();
             Gaussian X0 = Gaussian.FromMeanAndVariance(3, 0.5);
-            Gaussian Mean0 = Gaussian.FromMeanAndVariance(7, 1.0/3);
+            Gaussian Mean0 = Gaussian.FromMeanAndVariance(7, 1.0 / 3);
             Gamma Precision0 = Gamma.FromShapeAndScale(3, 3);
 
             // Fixed precision
@@ -219,6 +219,17 @@ namespace Microsoft.ML.Probabilistic.Tests
             Console.WriteLine("Q2 = {0}", watch.ElapsedTicks);
         }
 
+        [Fact]
+        public void GaussianOp_Laplace_Q_Test()
+        {
+            Gaussian sample = Gaussian.FromNatural(5.3861033232682936E-79, 2.901010900892175E-157);
+            Gaussian mean = Gaussian.FromNatural(-2.7232954231713977, 0.074308384738968308);
+            Gamma precision = Gamma.FromShapeAndRate(656.04827139518625, 1.4379651227587877E+159);
+            Gamma q = Gamma.FromShapeAndRate(1.1680712992725464, 2.0252330545334344E+155);
+            // Fails in 32-bit
+            GaussianOp.MeanAverageConditional(sample, mean, precision, q);
+        }
+
         [Fact]
         public void GaussianOpX()
         {
@@ -241,7 +252,13 @@ namespace Microsoft.ML.Probabilistic.Tests
                 xActual = GaussianOp.SampleAverageConditional(X, Mean, Precision, to_precision);
             }
 
-            X = Gaussian.FromNatural(0.1559599323109816, 8.5162535450918462);
+            X = Gaussian.FromNatural(-2.7793306963303595, 0.050822473645365768);
+            Mean = Gaussian.FromNatural(-5.9447032851878134E-09, 3.2975231004586637E-204);
+            Precision = Gamma.FromShapeAndRate(318.50907574398883, 9.6226982361933746E+205);
+            to_precision = Gamma.PointMass(0);
+            xActual = GaussianOp.SampleAverageConditional(X, Mean, Precision, to_precision);
+			
+			X = Gaussian.FromNatural(0.1559599323109816, 8.5162535450918462);
             Mean = Gaussian.PointMass(0.57957597647840942);
             Precision = Gamma.FromShapeAndRate(7.8308812008325587E+30, 8.2854255911709925E+30);
             to_precision = Gamma.FromShapeAndRate(1.4709139487775529, 0.14968339171493822);
@@ -352,7 +369,7 @@ namespace Microsoft.ML.Probabilistic.Tests
         public void GaussianOpMean()
         {
             Gaussian result = new Gaussian();
-            Gaussian X0 = Gaussian.FromMeanAndVariance(7, 1.0/3);
+            Gaussian X0 = Gaussian.FromMeanAndVariance(7, 1.0 / 3);
             Gaussian Mean0 = Gaussian.FromMeanAndVariance(3, 0.5);
             Gamma Precision0 = Gamma.FromShapeAndScale(3, 3);
 
@@ -369,32 +386,40 @@ namespace Microsoft.ML.Probabilistic.Tests
             Mean = Gaussian.PointMass(0);
             Precision = Gamma.FromShapeAndRate(3, 1);
             Gaussian xPostExpected = Gaussian.FromMeanAndVariance(0.178378819440295, 0.365796599498963);
-            Console.WriteLine(GaussianOp.SampleAverageConditional_slow(X, Mean, Precision)*X);
-            Assert.True(GaussianOp.SampleAverageConditional_slow(X, Mean, Precision).MaxDiff(xPostExpected/X) < 5e-7);
-            Console.WriteLine(GaussianOp_Slow.SampleAverageConditional(X, Mean, Precision)*X);
-            Assert.True(GaussianOp_Slow.SampleAverageConditional(X, Mean, Precision).MaxDiff(xPostExpected/X) < 5e-7);
+            Console.WriteLine(GaussianOp.SampleAverageConditional_slow(X, Mean, Precision) * X);
+            Assert.True(GaussianOp.SampleAverageConditional_slow(X, Mean, Precision).MaxDiff(xPostExpected / X) < 5e-7);
+            Console.WriteLine(GaussianOp_Slow.SampleAverageConditional(X, Mean, Precision) * X);
+            Assert.True(GaussianOp_Slow.SampleAverageConditional(X, Mean, Precision).MaxDiff(xPostExpected / X) < 5e-7);
         }
 
         /// <summary>
         /// Test that the operator behaves correctly when sample has large variance.
         /// Here we see that the message.Rate is non-monotonic in the sample variance, which doesn't seem right.
         /// </summary>
+        [Fact]
         [Trait("Category", "ModifiesGlobals")]
-        internal void GaussianOpPrecision3()
+        public void GaussianOpPrecision_IsMonotonicInSampleVariance()
         {
             using (TestUtils.TemporarilyAllowGaussianImproperMessages)
             {
                 Gaussian mean = Gaussian.PointMass(0);
-                Gamma precision = Gamma.FromShapeAndRate(2, 10);
-                for (int i = -10; i < 10; i++)
+                for (int logRate = 0; logRate < 310; logRate++)
                 {
-                    Gaussian sample = Gaussian.FromMeanAndPrecision(0, System.Math.Pow(10, -i));
-                    Gamma precMsg = GaussianOp.PrecisionAverageConditional(sample, mean, precision);
-                    //precMsg = GaussianOp_Laplace.PrecisionAverageConditional_slow(sample, mean, precision);
-                    //Gamma precMsg2 = GaussianOp_Slow.PrecisionAverageConditional(sample, mean, precision);
-                    //Console.WriteLine("{0}: {1} should be {2}", sample, precMsg, precMsg2);
-                    Gamma post = precMsg * precision;
-                    Console.WriteLine("{0}: {1} post = {2}", sample, precMsg.Rate, post.Rate);
+                    Gamma precision = Gamma.FromShapeAndRate(300, System.Math.Pow(10, logRate));
+                    double previousRate = double.PositiveInfinity;
+                    for (int i = 0; i < 310; i++)
+                    {
+                        Gaussian sample = Gaussian.FromMeanAndPrecision(0, System.Math.Pow(10, -i));
+                        Gamma precMsg = GaussianOp.PrecisionAverageConditional(sample, mean, precision);
+                        //precMsg = GaussianOp_Laplace.PrecisionAverageConditional_slow(sample, mean, precision);
+                        //Gamma precMsg2 = GaussianOp_Slow.PrecisionAverageConditional(sample, mean, precision);
+                        //Console.WriteLine("{0}: {1} should be {2}", sample, precMsg, precMsg2);
+                        Gamma post = precMsg * precision;
+                        //Trace.WriteLine($"{sample}: {precMsg.Rate} post = {post.Rate}");
+                        if (i >= logRate)
+                            Assert.True(precMsg.Rate <= previousRate);
+                        previousRate = precMsg.Rate;
+                    }
                 }
             }
         }
@@ -406,6 +431,24 @@ namespace Microsoft.ML.Probabilistic.Tests
             Gaussian X, Mean;
             Gamma Precision;
 
+            X = Gaussian.FromNatural(-1.5098177152950143E-09, 1.061649960537027E-168);
+            Mean = Gaussian.FromNatural(-3.6177299471249587, 0.11664740799025652);
+            Precision = Gamma.FromShapeAndRate(306.39423695125572, 1.8326832031565403E+170);
+            precMsg = Gamma.PointMass(0);
+            precMsg2 = GaussianOp_Slow.PrecisionAverageConditional(X, Mean, Precision);
+            Assert.True(precMsg.MaxDiff(precMsg2) < 1e-4);
+            precMsg2 = GaussianOp.PrecisionAverageConditional(X, Mean, Precision);
+            Assert.True(precMsg.MaxDiff(precMsg2) < 1e-4);
+
+            X = Gaussian.FromNatural(-0.55657497231637854, 6.6259783218464713E-141);
+            Mean = Gaussian.FromNatural(-2.9330116542965374, 0.07513822741674292);
+            Precision = Gamma.FromShapeAndRate(308.8184220331475, 4.6489382805051884E+142);
+            precMsg = Gamma.FromShapeAndRate(1.5000000000000628, 3.5279086383286634E+279);
+            precMsg2 = GaussianOp_Slow.PrecisionAverageConditional(X, Mean, Precision);
+            Assert.True(precMsg.MaxDiff(precMsg2) < 1e-4);
+            precMsg2 = GaussianOp.PrecisionAverageConditional(X, Mean, Precision);
+            Assert.True(precMsg.MaxDiff(precMsg2) < 1e-4);
+
             X = Gaussian.FromNatural(0, 1.0705890985886898E-153);
             Mean = Gaussian.PointMass(0);
             Precision = Gamma.FromShapeAndRate(1.6461630749684018, 1.0021354807958952E+153);
@@ -428,7 +471,7 @@ namespace Microsoft.ML.Probabilistic.Tests
             Assert.True(!double.IsNaN(precMsg.Rate));
 
             Gaussian X0 = Gaussian.FromMeanAndVariance(3, 0.5);
-            Gaussian Mean0 = Gaussian.FromMeanAndVariance(7, 1.0/3);
+            Gaussian Mean0 = Gaussian.FromMeanAndVariance(7, 1.0 / 3);
             Gamma Precision0 = Gamma.FromShapeAndScale(3, 3);
 
             precMsg = GaussianOp_Slow.PrecisionAverageConditional(Gaussian.FromNatural(0.010158033515400506, 0.0041117304509528533),
@@ -532,7 +575,7 @@ namespace Microsoft.ML.Probabilistic.Tests
                     }
                 }
                 count++;
-                if(count % 100 == 0)
+                if (count % 100 == 0)
                     Console.WriteLine("{0}", count);
             }
         }

test/Tests/Operators/OperatorTests.cs

@@ -25,14 +25,14 @@ namespace Microsoft.ML.Probabilistic.Tests
         [Fact]
         public void AverageTest()
         {
-            foreach(var a in Doubles())
+            foreach (var a in Doubles())
             {
-                foreach(var b in Doubles())
+                foreach (var b in Doubles())
                 {
-                    if (double.IsNaN(a+b)) continue;
+                    if (double.IsNaN(a + b)) continue;
                     double midpoint = MMath.Average(a, b);
-                    Assert.True(midpoint >= System.Math.Min(a,b));
-                    Assert.True(midpoint <= System.Math.Max(a,b));
+                    Assert.True(midpoint >= System.Math.Min(a, b));
+                    Assert.True(midpoint <= System.Math.Max(a, b));
                 }
             }
         }
@@ -1180,25 +1180,6 @@ namespace Microsoft.ML.Probabilistic.Tests
             Assert.True(IsPositiveOp.XAverageConditional(false, new Gaussian(-127, 11)).Equals(uniform));
             Assert.True(IsPositiveOp.XAverageConditional(true, new Gaussian(-1e5, 10)).IsProper());
             Assert.True(IsPositiveOp.XAverageConditional(false, new Gaussian(1e5, 10)).IsProper());
-
-            try
-            {
-                IsPositiveOp.XAverageConditional(true, Gaussian.FromNatural(1, 0));
-                Assert.True(false, "Did not throw exception");
-            }
-            catch (ImproperMessageException)
-            {
-                Console.WriteLine("Correctly threw ImproperMessageException");
-            }
-            try
-            {
-                IsPositiveOp.XAverageConditional(false, Gaussian.FromNatural(-1, 0));
-                Assert.True(false, "Did not throw exception");
-            }
-            catch (ImproperMessageException)
-            {
-                Console.WriteLine("Correctly threw ImproperMessageException");
-            }
         }
 
         [Fact]
@@ -1993,18 +1974,18 @@ zL = (L - mx)*sqrt(prec)
                 if (trial == 0)
                     upperBound = Gaussian.PointMass(1);
                 else
-                    upperBound = new Gaussian(1, 1);
-                Console.WriteLine("upperBound = {0}", upperBound);
+                    upperBound = new Gaussian(1, 1e-8);
+                Console.WriteLine($"upperBound = {upperBound}");
                 Gaussian lowerBound = Gaussian.PointMass(-1);
                 Gaussian result2 = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(isBetween, x, lowerBound, upperBound);
-                Console.WriteLine("{0}: {1}", lowerBound, result2);
-                for (int i = 8; i < 30; i++)
+                Console.WriteLine($"{lowerBound}: {result2}");
+                for (int i = 6; i < 30; i++)
                 {
                     double v = System.Math.Pow(0.1, i);
                     lowerBound = Gaussian.FromMeanAndVariance(-1, v);
                     Gaussian result = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(isBetween, x, lowerBound, upperBound);
                     double error = result.MaxDiff(result2);
-                    Console.WriteLine("{0}: {1} {2}", lowerBound, result, error);
+                    Console.WriteLine($"{lowerBound}: {result} {error}");
                     Assert.True(error < 1e-6);
                 }
             }
@@ -2019,20 +2000,21 @@ zL = (L - mx)*sqrt(prec)
             Gaussian lowerBound;
             for (int trial = 0; trial < 2; trial++)
             {
-                if (trial == 1)
+                if (trial == 0)
                     lowerBound = Gaussian.PointMass(-1);
                 else
-                    lowerBound = new Gaussian(-1, 1);
+                    lowerBound = new Gaussian(-1, 1e-8);
+                Console.WriteLine($"lowerBound = {lowerBound}");
                 Gaussian upperBound = Gaussian.PointMass(1);
                 Gaussian result2 = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(isBetween, x, lowerBound, upperBound);
-                Console.WriteLine("{0}: {1}", upperBound, result2);
-                for (int i = 8; i < 30; i++)
+                Console.WriteLine($"{upperBound}: {result2}");
+                for (int i = 6; i < 300; i++)
                 {
                     double v = System.Math.Pow(0.1, i);
                     upperBound = Gaussian.FromMeanAndVariance(1, v);
                     Gaussian result = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(isBetween, x, lowerBound, upperBound);
                     double error = result.MaxDiff(result2);
-                    Console.WriteLine("{0}: {1} {2}", lowerBound, result, error);
+                    Console.WriteLine($"{upperBound}: {result} {error}");
                     Assert.True(error < 1e-6);
                 }
             }
@@ -2653,6 +2635,27 @@ zL = (L - mx)*sqrt(prec)
                 (a.Precision + b.Precision - GetSumError(a.Precision, b.Precision));
         }
 
+        [Fact]
+        public void GaussianIsBetween_CRRR_UncertainXTest()
+        {
+            Gaussian lowerBound = Gaussian.FromNatural(17.111433288915187, 0.66938434508155731);
+            Gaussian upperBound = Gaussian.FromNatural(7.7694959349146462, 0.49485730932861044);
+            lowerBound = Gaussian.FromNatural(-54.625321469620474, 24.217302563266891);
+            upperBound = Gaussian.FromNatural(54.406085493265941, 24.077474982519977);
+            double previousLogProb = double.PositiveInfinity;
+            for (int i = 5; i < 200; i++)
+            {
+                Gaussian X = Gaussian.FromNatural(2.2204460492503131E-16, System.Math.Pow(0.1, i));
+                X = Gaussian.FromMeanAndPrecision(-60, System.Math.Pow(0.1, i));
+                double logProb = DoubleIsBetweenOp.LogProbBetween(X, lowerBound, upperBound);
+                Trace.WriteLine($"X={X}: logProb={logProb}");
+                Assert.True(!double.IsNaN(logProb));
+                Assert.True(!double.IsInfinity(logProb));
+                Assert.True(logProb < previousLogProb);
+                previousLogProb = logProb;
+            }
+        }
+
         // Test that the operator behaves correctly for arguments with small variance
         [Fact]
         public void GaussianIsBetween_PointX_Test()
@@ -2775,11 +2778,245 @@ weight * (tau + alphaX) + alphaX
             Assert.True(MaxUlpDiff(expected, result2) <= 5);
         }
 
+        [Fact]
+        [Trait("Category", "OpenBug")]
+        public void GaussianIsBetweenCRRR_IsMonotonicInX()
+        {
+            Gaussian previousToX = new Gaussian();
+            Gaussian previousToLowerBound = new Gaussian();
+            Gaussian previousToUpperBound = new Gaussian();
+            Gaussian lowerBound = Gaussian.FromNatural(-200, 100);
+            Gaussian upperBound = Gaussian.FromNatural(255, 147);
+            Bernoulli isBetween = Bernoulli.PointMass(true);
+            double xMeanMaxUlpError = 0;
+            double xPrecisionMaxUlpError = 0;
+            double uMeanMaxUlpError = 0;
+            double uPrecisionMaxUlpError = 0;
+            double lMeanMaxUlpError = 0;
+            for (int i = 10; i < 3000; i++)
+            {
+                Gaussian X = Gaussian.FromNatural(1.2, System.Math.Pow(10, -i * 0.1));
+                var toX = DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, X, lowerBound, upperBound);
+                var toLowerBound = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(isBetween, X, lowerBound, upperBound);
+                var toUpperBound = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(isBetween, X, lowerBound, upperBound);
+                Trace.WriteLine($"{X}: {toX} {toLowerBound} {toUpperBound}");
+                if (!previousToX.IsUniform())
+                {
+                    if (previousToX.GetMean() > toX.GetMean())
+                    {
+                        xMeanMaxUlpError = System.Math.Max(xMeanMaxUlpError, UlpDiff(previousToX.GetMean(), toX.GetMean()));
+                        //Assert.True(xMeanMaxUlpError < 1e13);
+                    }
+                    if (previousToX.Precision > toX.Precision)
+                    {
+                        xPrecisionMaxUlpError = System.Math.Max(xPrecisionMaxUlpError, UlpDiff(previousToX.Precision, toX.Precision));
+                        //Assert.True(xPrecisionMaxUlpError <= 0);
+                    }
+                    if (previousToLowerBound.GetMean() < toLowerBound.GetMean())
+                    {
+                        lMeanMaxUlpError = System.Math.Max(lMeanMaxUlpError, UlpDiff(previousToLowerBound.GetMean(), toLowerBound.GetMean()));
+                    }
+                    //Assert.True(previousToLowerBound.GetVariance() <= toLowerBound.GetVariance());
+                    if (previousToUpperBound.GetMean() > toUpperBound.GetMean())
+                    {
+                        uMeanMaxUlpError = System.Math.Max(uMeanMaxUlpError, UlpDiff(previousToUpperBound.GetMean(), toUpperBound.GetMean()));
+                    }
+                    if (previousToUpperBound.Precision < toUpperBound.Precision)
+                    {
+                        uPrecisionMaxUlpError = System.Math.Max(uPrecisionMaxUlpError, UlpDiff(previousToUpperBound.Precision, toUpperBound.Precision));
+                    }
+                }
+                previousToX = toX;
+                previousToLowerBound = toLowerBound;
+                previousToUpperBound = toUpperBound;
+            }
+            Trace.WriteLine($"xMeanMaxUlpError = {xMeanMaxUlpError}, xPrecisionMaxUlpError = {xPrecisionMaxUlpError}, uMeanMaxUlpError = {uMeanMaxUlpError}, uPrecisionMaxUlpError = {uPrecisionMaxUlpError}, lMeanMaxUlpError = {lMeanMaxUlpError}");
+            // TODO: tighten these thresholds
+            Assert.True(xMeanMaxUlpError < 1e15);
+            Assert.True(xPrecisionMaxUlpError < 1e15);
+            Assert.True(uPrecisionMaxUlpError < 1);
+            Assert.True(uMeanMaxUlpError < 1);
+            Assert.True(lMeanMaxUlpError < 1);
+        }
+
+        [Fact]
+        [Trait("Category", "OpenBug")]
+        public void GaussianIsBetweenCRRR_IsMonotonicInX2()
+        {
+            Gaussian previousToX = new Gaussian();
+            Gaussian previousToLowerBound = new Gaussian();
+            Gaussian previousToUpperBound = new Gaussian();
+            Gaussian upperBound = Gaussian.FromNatural(200, 100);
+            Gaussian lowerBound = Gaussian.FromNatural(-255, 147);
+            Bernoulli isBetween = Bernoulli.PointMass(true);
+            DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, Gaussian.FromNatural(-1.2, 6.3095734448019427E-17), lowerBound, upperBound);
+            double xMeanMaxUlpError = 0;
+            double xPrecisionMaxUlpError = 0;
+            double lMeanMaxUlpError = 0;
+            double lPrecisionMaxUlpError = 0;
+            double uMeanMaxUlpError = 0;
+            for (int i = 10; i < 3000; i++)
+            {
+                Gaussian X = Gaussian.FromNatural(-1.2, System.Math.Pow(10, -i * 0.1));
+                //Gaussian toXExpected = new Gaussian(-0.3047, 0.5397);
+                //SolveAlphaBeta(X, toXExpected, out double alpha, out double beta);
+                //Trace.WriteLine($"expected alpha = {alpha}, beta = {beta}");
+                var toX = DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, X, lowerBound, upperBound);
+                var toLowerBound = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(isBetween, X, lowerBound, upperBound);
+                var toUpperBound = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(isBetween, X, lowerBound, upperBound);
+                Trace.WriteLine($"{X}: {toX} {toLowerBound} {toUpperBound}");
+                if (!previousToX.IsUniform())
+                {
+                    if (previousToX.GetMean() < toX.GetMean())
+                    {
+                        xMeanMaxUlpError = System.Math.Max(xMeanMaxUlpError, UlpDiff(previousToX.GetMean(), toX.GetMean()));
+                        //Assert.True(xMeanMaxUlpError < 1e11);
+                    }
+                    if (previousToX.Precision > toX.Precision)
+                    {
+                        xPrecisionMaxUlpError = System.Math.Max(xPrecisionMaxUlpError, UlpDiff(previousToX.Precision, toX.Precision));
+                        //Assert.True(xPrecisionMaxUlpError <= 0);
+                    }
+                    if (previousToLowerBound.GetMean() < toLowerBound.GetMean())
+                    {
+                        lMeanMaxUlpError = System.Math.Max(lMeanMaxUlpError, UlpDiff(previousToLowerBound.GetMean(), toLowerBound.GetMean()));
+                    }
+                    if (previousToLowerBound.Precision < toLowerBound.Precision)
+                    {
+                        lPrecisionMaxUlpError = System.Math.Max(lPrecisionMaxUlpError, UlpDiff(previousToLowerBound.Precision, toLowerBound.Precision));
+                    }
+                    if (previousToUpperBound.GetMean() > toUpperBound.GetMean())
+                    {
+                        uMeanMaxUlpError = System.Math.Max(uMeanMaxUlpError, UlpDiff(previousToUpperBound.GetMean(), toUpperBound.GetMean()));
+                    }
+                    //Assert.True(previousToUpperBound.GetVariance() <= toUpperBound.GetVariance());
+                }
+                previousToX = toX;
+                previousToLowerBound = toLowerBound;
+                previousToUpperBound = toUpperBound;
+            }
+            Trace.WriteLine($"xMeanMaxUlpError = {xMeanMaxUlpError}, xPrecisionMaxUlpError = {xPrecisionMaxUlpError}, lMeanMaxUlpError = {lMeanMaxUlpError}, lPrecisionMaxUlpError = {lPrecisionMaxUlpError}, uMeanMaxUlpError={uMeanMaxUlpError}");
+            // TODO: tighten these thresholds
+            Assert.True(xMeanMaxUlpError < 1e15);
+            Assert.True(xPrecisionMaxUlpError < 1e15);
+            Assert.True(lMeanMaxUlpError < 1);
+            Assert.True(lPrecisionMaxUlpError < 1);
+            Assert.True(uMeanMaxUlpError < 1);
+        }
+
+        private static void SolveAlphaBeta(Gaussian prior, Gaussian msg, out double alpha, out double beta)
+        {
+            if (prior.IsPointMass)
+            {
+                beta = msg.Precision;
+                alpha = msg.MeanTimesPrecision - prior.Point * beta;
+                return;
+            }
+            if (msg.IsPointMass)
+            {
+                beta = prior.Precision;
+                alpha = (msg.Point - prior.GetMean()) * beta;
+                return;
+            }
+            beta = 1 / (1 / msg.Precision + 1 / prior.Precision);
+            double weight = beta / (prior.Precision - beta);
+            // weight*tau + (weight + 1)*alpha = msg.MeanTimesPrecision
+            alpha = (msg.MeanTimesPrecision - weight * prior.MeanTimesPrecision) / (1 + weight);
+        }
+
+        [Fact]
+        [Trait("Category", "OpenBug")]
+        public void GaussianIsBetweenCRRR_NegativeUpperBoundTest()
+        {
+            Gaussian X = Gaussian.FromNatural(813.982758311301, 1.0594806725507477);
+            Gaussian previousToX = new Gaussian();
+            Gaussian previousXpost = new Gaussian();
+            double tolerance = 1e-10;
+            for (int i = 8; i < 100; i++)
+            {
+                // seems like answer should always be Gaussian(m/v=-814, 1/v=0)
+                Gaussian upperBound = Gaussian.FromNatural(-System.Math.Pow(10, i), 9);
+                Gaussian toX = DoubleIsBetweenOp.XAverageConditional_Slow(Bernoulli.PointMass(true), X, Gaussian.PointMass(0), upperBound);
+                Gaussian Xpost = X * toX;
+                Trace.WriteLine($"{upperBound}: {toX} {toX.MeanTimesPrecision} {Xpost}");
+                // lowerBound is point mass, so cannot be violated.
+                Assert.True(Xpost.GetMean() >= 0 - tolerance);
+                if (!previousToX.IsUniform())
+                {
+                    // upperBound is decreasing, so posterior mean should be decreasing.
+                    Assert.True(Xpost.GetMean() <= previousXpost.GetMean() + tolerance);
+                }
+                previousToX = toX;
+                previousXpost = Xpost;
+            }
+        }
+
+        [Fact]
+        public void GaussianIsBetweenCRRR_LowerBoundTest()
+        {
+            Gaussian X = Gaussian.FromNatural(898.71395259259464, 1.4308788553248037);
+            Gaussian lowerBound = Gaussian.FromNatural(17028358.45574614, 9);
+            Gaussian upperBound = Gaussian.FromNatural(412820.08287991461, 423722.55474045349);
+            for (int i = -10; i <= 0; i++)
+            {
+                lowerBound = Gaussian.FromNatural(17028358.45574614*System.Math.Pow(2,i), 9);
+                Gaussian toLowerBound = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), X, lowerBound, upperBound);
+                Trace.WriteLine($"{lowerBound}: {toLowerBound.MeanTimesPrecision} {toLowerBound.Precision}");
+                Assert.False(toLowerBound.IsPointMass);
+            }
+        }
+
+        [Fact]
+        [Trait("Category", "OpenBug")]
+        public void GaussianIsBetweenCRRR_SmallXPrecisionTest()
+        {
+            Gaussian lowerBound = Gaussian.FromNatural(-102.3311202057678, 91.572320438929935);
+            Gaussian upperBound = Gaussian.FromNatural(102.27224205502382, 91.541070478258376);
+            foreach (var mean in new[] { 1e7, -1e7 })
+            {
+                Gaussian toLowerBoundPrev = Gaussian.FromNatural(double.MaxValue, double.MaxValue);
+                Gaussian toXPrev = Gaussian.FromNatural(double.MaxValue, double.MaxValue);
+                double xMeanTimesPrecisionMaxUlpError = 0;
+                double lowerBoundMeanTimesPrecisionMaxUlpError = 0;
+                for (int i = 0; i < 200; i++)
+                {
+                    Gaussian X = Gaussian.FromMeanAndPrecision(mean, System.Math.Pow(2, -i*1-20));
+                    Gaussian toX = DoubleIsBetweenOp.XAverageConditional_Slow(Bernoulli.PointMass(true), X, lowerBound, upperBound);
+                    Gaussian toLowerBound = toLowerBoundPrev;// DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), X, lowerBound, upperBound);
+                    Trace.WriteLine($"{i} {X}: {toX.MeanTimesPrecision:r} {toX.Precision:r} {toLowerBound.MeanTimesPrecision:r} {toLowerBound.Precision:r}");
+                    Assert.False(toLowerBound.IsPointMass);
+                    if ((mean > 0 && toLowerBound.MeanTimesPrecision > toLowerBoundPrev.MeanTimesPrecision) ||
+                        (mean < 0 && toLowerBound.MeanTimesPrecision < toLowerBoundPrev.MeanTimesPrecision))
+                    {
+                        lowerBoundMeanTimesPrecisionMaxUlpError = System.Math.Max(lowerBoundMeanTimesPrecisionMaxUlpError, UlpDiff(toLowerBound.MeanTimesPrecision, toLowerBoundPrev.MeanTimesPrecision));
+                        //Assert.True(meanTimesPrecisionMaxUlpError < 5);
+                    }
+                    if ((mean > 0 && toX.MeanTimesPrecision > toXPrev.MeanTimesPrecision) ||
+                        (mean < 0 && toX.MeanTimesPrecision < toXPrev.MeanTimesPrecision))
+                    {
+                        xMeanTimesPrecisionMaxUlpError = System.Math.Max(xMeanTimesPrecisionMaxUlpError, UlpDiff(toX.MeanTimesPrecision, toXPrev.MeanTimesPrecision));
+                        //Assert.True(xMeanTimesPrecisionMaxUlpError < 1e12);
+                    }
+                    toLowerBoundPrev = toLowerBound;
+                    toXPrev = toX;
+                }
+                Trace.WriteLine($"xMeanTimesPrecisionMaxUlpError = {xMeanTimesPrecisionMaxUlpError} lowerBoundMeanTimesPrecisionMaxUlpError = {lowerBoundMeanTimesPrecisionMaxUlpError}");
+                Assert.True(xMeanTimesPrecisionMaxUlpError < 1e12);
+                Assert.True(lowerBoundMeanTimesPrecisionMaxUlpError < 1e4);
+            }
+        }
+
         [Fact]
         public void GaussianIsBetweenTest2()
         {
-            Assert.True(!double.IsNaN(DoubleIsBetweenOp.XAverageConditional(Bernoulli.PointMass(true), new Gaussian(1, 2), double.PositiveInfinity, double.PositiveInfinity).MeanTimesPrecision));
             Bernoulli isBetween = new Bernoulli(1);
+            Assert.False(double.IsNaN(DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, Gaussian.FromNatural(0.9106071714590378, 5.9521837280027985E-11), Gaussian.FromNatural(-49.9894026120194, 107.30343404076896), Gaussian.FromNatural(49.051818445888259, 107.26846525506932)).MeanTimesPrecision));
+            Assert.False(double.IsNaN(DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, Gaussian.FromNatural(2.2204460492503131E-16, 6.9388939039072284E-18), Gaussian.FromNatural(17.111433288915187, 0.66938434508155731), Gaussian.FromNatural(7.7694959349146462, 0.49485730932861044)).MeanTimesPrecision));
+            DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, Gaussian.FromNatural(0, 0.0038937777431664196), Gaussian.PointMass(double.NegativeInfinity), Gaussian.FromNatural(-1.6, 0.8));
+            DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(isBetween, Gaussian.FromNatural(1.156293233217532E-25, 6.162975822039154E-33), Gaussian.FromNatural(-102.3311202057678, 91.572320438929935), Gaussian.FromNatural(102.27224205502382, 91.541070478258376));
+            Assert.False(double.IsNaN(DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, Gaussian.FromNatural(980.18122429721575, 1.409544490082087), Gaussian.FromNatural(17028174.685026139, 837.26675043005957), Gaussian.FromNatural(412820.4122154137, 423722.54499249317)).MeanTimesPrecision));
+            Assert.False(double.IsNaN(DoubleIsBetweenOp.XAverageConditional_Slow(isBetween, Gaussian.FromNatural(-5.3548456213550253E-41, 4.61370960061741E-81), Gaussian.FromNatural(-15848812.635800883, 13451.362337266379), Gaussian.FromNatural(-22204349.280881952, 389690.00236138358)).MeanTimesPrecision));
+            Assert.False(double.IsNaN(DoubleIsBetweenOp.XAverageConditional(isBetween, new Gaussian(1, 2), double.PositiveInfinity, double.PositiveInfinity).MeanTimesPrecision));
             Gaussian x = new Gaussian(0, 1);
             Gaussian lowerBound = new Gaussian(1, 8);
             Gaussian upperBound = new Gaussian(3, 3);
@@ -2894,7 +3131,7 @@ weight * (tau + alphaX) + alphaX
                 double m, v;
                 Gaussian x2 = new Gaussian(-lowerBounds[i], 1);
                 Xpost.SetToProduct(x2, IsPositiveOp.XAverageConditional(true, x2));
-                Console.WriteLine(Xpost);
+                //Console.WriteLine(Xpost);
                 Xpost.GetMeanAndVariance(out m, out v);
                 Xpost.SetToProduct(x, DoubleIsBetweenOp.XAverageConditional(true, x, lowerBounds[i], Double.PositiveInfinity));
                 Assert.True(Xpost.MaxDiff(new Gaussian(m + lowerBounds[i], v)) < 1e-5 / v);
@@ -2902,15 +3139,27 @@ weight * (tau + alphaX) + alphaX
                 //Assert.True(Xpost.MaxDiff(new Gaussian(m+lowerBounds[i], v)) < 1e-8);
             }
 
-            //Gaussian upperBound = new Gaussian(3,4);
-            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0.0), lowerBound, Gaussian.PointMass(0.0));
-            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0.0), lowerBound, Gaussian.PointMass(1.0));
-            Assert.True(Lpost.MaxDiff(IsPositiveOp.XAverageConditional(false, lowerBound)) < 1e-10);
+
+            Gaussian Lexpected = IsPositiveOp.XAverageConditional(false, lowerBound);
+            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.FromMeanAndVariance(0, 1e-10), lowerBound, Gaussian.FromMeanAndVariance(-1, 1));
+            Assert.True(Lpost.MaxDiff(Lexpected) < 1e-10);
+            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0), lowerBound, Gaussian.FromMeanAndVariance(-1, 1));
+            Assert.True(Lpost.MaxDiff(Lexpected) < 1e-10);
+            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(1), lowerBound, Gaussian.PointMass(0));
+            Assert.True(Lpost.MaxDiff(Lexpected) < 1e-10);
+            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0), lowerBound, Gaussian.PointMass(0));
+            Assert.True(Lpost.MaxDiff(Lexpected) < 1e-10);
+            Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0), lowerBound, Gaussian.PointMass(1));
+            Assert.True(Lpost.MaxDiff(Lexpected) < 1e-10);
             //Lpost = DoubleIsBetweenOp.LowerBoundAverageConditional(true,Gaussian.Uniform(),lowerBound,0);
             //Assert.True(Lpost.MaxDiff(IsPositiveOp.XAverageConditional(false,lowerBound)) < 1e-3);
-            Upost = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(-1), Gaussian.PointMass(0.0), upperBound);
-            Upost = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0.0), Gaussian.PointMass(0.0), upperBound);
-            Assert.True(Upost.MaxDiff(IsPositiveOp.XAverageConditional(true, upperBound)) < 1e-10);
+            Gaussian Uexpected = IsPositiveOp.XAverageConditional(true, upperBound);
+            Upost = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(-1), Gaussian.PointMass(0), upperBound);
+            Assert.True(Upost.MaxDiff(Uexpected) < 1e-10);
+            Upost = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0), Gaussian.PointMass(0), upperBound);
+            Assert.True(Upost.MaxDiff(Uexpected) < 1e-10);
+            Upost = DoubleIsBetweenOp.UpperBoundAverageConditional_Slow(Bernoulli.PointMass(true), Gaussian.PointMass(0), Gaussian.PointMass(-1), upperBound);
+            Assert.True(Upost.MaxDiff(Uexpected) < 1e-10);
             //Upost = DoubleIsBetweenOp.UpperBoundAverageConditional(true,Gaussian.Uniform(),0,upperBound);
             //Assert.True(Upost.MaxDiff(IsPositiveOp.XAverageConditional(true,upperBound)) < 1e-3);
         }

test/Tests/Properties/AssemblyInfo.cs

@@ -43,5 +43,5 @@ using System.Runtime.InteropServices;
 [assembly: Microsoft.ML.Probabilistic.Factors.Attributes.HasMessageFunctions]
 [assembly:
     InternalsVisibleTo(
-        "ML.Probabilistic.TestApp,PublicKey=0024000004800000940000000602000000240000525341310004000001000100551f07a755a3e3f2901fa321ab631d13d6192b4e6ac9c87279500f49d6635cde6902587752eff20402f46f6ea9c3d80e827580a799840aaab9a49b1d2597e4c1798ee93c5cb66851e9d22f4d6e8110571f4a2e59f1d760f7be04fb10e7dc43ee7ed2831907731427b9815c5fe7f4888f9933ee7a1ad5d1f293fd8ab834fac1be"
+        "Microsoft.ML.Probabilistic.TestApp,PublicKey=0024000004800000940000000602000000240000525341310004000001000100551f07a755a3e3f2901fa321ab631d13d6192b4e6ac9c87279500f49d6635cde6902587752eff20402f46f6ea9c3d80e827580a799840aaab9a49b1d2597e4c1798ee93c5cb66851e9d22f4d6e8110571f4a2e59f1d760f7be04fb10e7dc43ee7ed2831907731427b9815c5fe7f4888f9933ee7a1ad5d1f293fd8ab834fac1be"
         )]