C# Gaussian.IsUniform方法代码示例

		//-- Constant bounds --------------------------------------------------------------------------------

		/// <summary>
		/// The logarithm of the probability that L &lt;= X &lt; U.
		/// </summary>
		/// <param name="X"></param>
		/// <param name="L">Can be negative infinity.</param>
		/// <param name="U">Can be positive infinity.</param>
		/// <returns></returns>
		public static double LogProbBetween(Gaussian X, double L, double U)
		{
			if (L > U) throw new AllZeroException("low > high (" + L + " > " + U + ")");
			if (X.IsPointMass)
			{
				return Factor.IsBetween(X.Point, L, U) ? 0.0 : Double.NegativeInfinity;
			}
			else if (X.IsUniform())
			{
				if (Double.IsNegativeInfinity(L))
				{
					if (Double.IsPositiveInfinity(U)) return 0.0; // always between
					else return -MMath.Ln2;  // between half the time
				}
				else if (Double.IsPositiveInfinity(U)) return -MMath.Ln2;  // between half the time
				else return Double.NegativeInfinity;  // never between two finite numbers
			}
			else
			{
				double sqrtPrec = Math.Sqrt(X.Precision);
				double mx = X.GetMean();
				double pl = MMath.NormalCdfLn(sqrtPrec * (L - mx));  // log(p(x <= L))
				double pu = MMath.NormalCdfLn(sqrtPrec * (U - mx));  // log(p(x <= U))
				if (pl == pu) return Double.NegativeInfinity;
				if (Double.IsNegativeInfinity(pl)) return pu;
				// log(NormalCdf(yu) - NormalCdf(yl)) = NormalCdfLn(yu) + log(1 - NormalCdf(yl)/NormalCdf(yu))
				return pu + MMath.Log1MinusExp(pl - pu);
			}
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="isPositive">Incoming message from 'isPositive'.</param>
		/// <param name="x">Incoming message from 'x'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(isPositive,x) p(isPositive,x) factor(isPositive,x))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Bernoulli isPositive, Gaussian x)
		{
			Bernoulli to_isPositive = IsPositiveAverageConditional(x);
			return isPositive.GetLogAverageOf(to_isPositive);
#if false
			// Z = p(b=T) p(x > 0) + p(b=F) p(x <= 0)
			//   = p(b=F) + (p(b=T) - p(b=F)) p(x > 0)
			if (x.IsPointMass) {
				return Factor.IsPositive(x.Point) ? isPositive.GetLogProbTrue() : isPositive.GetLogProbFalse();
			} else if(x.IsUniform()) {
				return Bernoulli.LogProbEqual(isPositive.LogOdds,0.0);
			} else {
				// m/sqrt(v) = (m/v)/sqrt(1/v)
				double z = x.MeanTimesPrecision / Math.Sqrt(x.Precision);
				if (isPositive.IsPointMass) {
					return isPositive.Point ? MMath.NormalCdfLn(z) : MMath.NormalCdfLn(-z);
				} else {
					return MMath.LogSumExp(isPositive.GetLogProbTrue() + MMath.NormalCdfLn(z), isPositive.GetLogProbFalse() + MMath.NormalCdfLn(-z));
				}
			}
#endif
		}

		public static Gaussian ProductAverageConditional2(Gaussian Product, Gaussian A, Gaussian B, Gaussian to_A)
		{
			if (Product.IsUniform()) return GaussianProductVmpOp.ProductAverageLogarithm(A, B);
			double mx,vx;
			Product.GetMeanAndVariance(out mx, out vx);
			double ma,va;
			A.GetMeanAndVariance(out ma, out va);
			Gaussian Apost = A*to_A;
			double ahat = Apost.GetMean();
			double mb,vb;
			B.GetMeanAndVariance(out mb, out vb);
			double denom = 1/(vx + ahat*ahat*vb);
			double diff = mx-ahat*mb;
			double dlogz = -diff*denom;
			double dd = ahat*vb;
			double q = denom*(mb + 2*dd*diff*denom);
			double ddd = vb;
			double n = diff*diff;
			double dn = -diff*mb;
			double ddn = mb*mb;
			double dda = denom*(-(ddd+ddn) + denom*(2*dd*(dd+2*dn)+n*ddd - denom*4*n*dd*dd));
			double da = va*q/(1 - va*dda);
			double ddlogz = -denom + q*da;
			return GaussianProductOp_Laplace.GaussianFromAlphaBeta(Product, dlogz, -ddlogz, true);
		}

		//-- Random bounds --------------------------------------------------------------------------------

		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)
		{
			double mx, vx, ml, vl, mu, vu;
			X.GetMeanAndVariance(out mx, out vx);
			L.GetMeanAndVariance(out ml, out vl);
			U.GetMeanAndVariance(out mu, out vu);
			if (X.IsPointMass && L.IsPointMass)
			{
				invSqrtVxl = Double.PositiveInfinity;
				yl = (X.Point >= L.Point) ? Double.PositiveInfinity : Double.NegativeInfinity;
			}
			else if (L.IsUniform())
			{
				invSqrtVxl = 0.0;
				yl = 0;
			}
			else
			{
				invSqrtVxl = 1.0 / Math.Sqrt(vx + vl);
				yl = (mx - ml) * invSqrtVxl;
			}
			if (X.IsPointMass && U.IsPointMass)
			{
				invSqrtVxu = Double.PositiveInfinity;
				yu = (X.Point < U.Point) ? Double.PositiveInfinity : Double.NegativeInfinity;
			}
			else if (U.IsUniform())
			{
				invSqrtVxu = 0.0;
				yu = 0;
			}
			else
			{
				invSqrtVxu = 1.0 / Math.Sqrt(vx + vu);
				yu = (mu - mx) * invSqrtVxu;
			}
			if (X.IsPointMass)
			{
				r = 0.0;
			}
			else
			{
				r = -vx * invSqrtVxl * invSqrtVxu;
			}
		}

		public static Gaussian XAverageConditional([SkipIfUniform] Bernoulli isBetween, Gaussian X, Gaussian lowerBound, Gaussian upperBound)
		{
			if (lowerBound.IsPointMass && upperBound.IsPointMass) return XAverageConditional(isBetween, X, lowerBound.Point, upperBound.Point);
			Gaussian result = new Gaussian();
			if(isBetween.IsUniform()) return result;
			if (X.IsPointMass)
			{
				result.SetToUniform(); // TODO: return the limiting distribution
			}
			else if (X.IsUniform())
			{
				if (lowerBound.IsUniform() || upperBound.IsUniform() ||
					(lowerBound.IsPointMass && Double.IsInfinity(lowerBound.Point)) ||
					(upperBound.IsPointMass && Double.IsInfinity(upperBound.Point)) ||
					!Double.IsPositiveInfinity(isBetween.LogOdds))
				{
					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; // vlu > 0
					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);
					double munew = mu + vu * alphaU;
					double mlnew = ml - vl * alphaU;
					double vunew = vu - vu * vu * betaU;
					double vlnew = vl - vl * vl * betaU;
					double diff = munew - mlnew;
					result.SetMeanAndVariance((munew + mlnew) / 2, diff * diff / 12 + (vunew + vlnew + vu * vl * betaU) / 3);
				}
				else throw new NotImplementedException();
			}
			else
			{
				// X is not a point mass or uniform
				double logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
				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 (!lowerBound.IsUniform() && !Double.IsInfinity(yl))
				{
					// since X and lowerBound are not both uniform, invSqrtVxl > 0
					double logPhiL = Math.Log(invSqrtVxl) + Gaussian.GetLogProb(yl, 0, 1);
					if (r > -1) logPhiL += MMath.NormalCdfLn((yu - r * yl) / Math.Sqrt(1 - r * r));
					alphaL = -d_p * Math.Exp(logPhiL - logZ);
				}
				double alphaU = 0.0;
				if (!upperBound.IsUniform() && !Double.IsInfinity(yu))
				{
					// since X and upperBound are not both uniform, invSqrtVxu > 0
					double logPhiU = Math.Log(invSqrtVxu) + Gaussian.GetLogProb(yu, 0, 1);
					if (r > -1) logPhiU += MMath.NormalCdfLn((yl - r * yu) / Math.Sqrt(1 - r * r));
					alphaU = d_p * Math.Exp(logPhiU - logZ);
				}
				double alphaX = -alphaL - alphaU;
				// (mx - ml) / (vl + vx) = yl*invSqrtVxl
				double betaX = alphaX * alphaX;
				if (!Double.IsInfinity(yl))
				{
					betaX -= alphaL * (yl * invSqrtVxl);
				}
				if (!Double.IsInfinity(yu))
				{
					betaX += alphaU * (yu * invSqrtVxu);
				}
				if (r > -1 && r != 0 && !Double.IsInfinity(yl) && !Double.IsInfinity(yu))
				{
					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);
					betaX += c * (-2 * X.Precision + invSqrtVxl * invSqrtVxl + invSqrtVxu * invSqrtVxu);
				}
				double weight = betaX / (X.Precision - betaX);
				if (ForceProper && weight < 0) weight = 0;
				result.Precision = weight * X.Precision;
				result.MeanTimesPrecision = weight * (X.MeanTimesPrecision + alphaX) + alphaX;
				if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision)) throw new ApplicationException("result is NaN");
			}
			return result;
		}

		public static Gaussian ProductAverageConditional(Gaussian Product, [SkipIfUniform] Gaussian A, [SkipIfUniform] Gaussian B)
		{
			if (A.IsPointMass) return ProductAverageConditional(A.Point, B);
			if (B.IsPointMass) return ProductAverageConditional(A, B.Point);
			if (Product.IsPointMass) return Gaussian.Uniform();
			if (Product.IsUniform()) return GaussianProductVmpOp.ProductAverageLogarithm(A, B);
			double mA, vA;
			A.GetMeanAndVariance(out mA, out vA);
			double mB, vB;
			B.GetMeanAndVariance(out mB, out vB);
			double mProduct, vProduct;
			Product.GetMeanAndVariance(out mProduct, out vProduct);
			// algorithm: quadrature on A from -1 to 1, plus quadrature on 1/A from -1 to 1.
			double z = 0, sumX = 0, sumX2 = 0;
			for (int i = 0; i <= QuadratureNodeCount; i++) {
				double a = (2.0 * i) / QuadratureNodeCount - 1;
				double logfA = Gaussian.GetLogProb(mProduct, a * mB, vProduct + a * a * vB) + Gaussian.GetLogProb(a, mA, vA);
				double fA = Math.Exp(logfA);

				z += fA;
				double b = (mB * vProduct + a * mProduct * vB) / (vProduct + a * a * vB);
				double b2 = b * b + (vProduct * vB) / (vProduct + a * a * vB);
				double x = a * b;
				double x2 = a * a * b2;
				sumX += x * fA;
				sumX2 += x2 * fA;

				double invA = a;
				a = 1.0 / invA;
				double logfInvA = Gaussian.GetLogProb(mProduct * invA, mB, vProduct * invA * invA + vB) + Gaussian.GetLogProb(a, mA, vA) - Math.Log(Math.Abs(invA + Double.Epsilon));
				double fInvA = Math.Exp(logfInvA);
				z += fInvA;
				b = (mB * vProduct + a * mProduct * vB) / (vProduct + a * a * vB);
				b2 = b * b + (vProduct * vB) / (vProduct + a * a * vB);
				x = a * b;
				x2 = a * a * b2;
				sumX += x * fInvA;
				sumX2 += x2 * fInvA;
			}
			double mean = sumX / z;
			double var = sumX2 / z - mean * mean;
			Gaussian result = Gaussian.FromMeanAndVariance(mean, var);
			if (ForceProper) result.SetToRatioProper(result, Product);
			else result.SetToRatio(result, Product);
			return result;
		}

		public static Gaussian SampleAverageConditional(Gaussian sample, [SkipIfUniform] Gaussian mean, [Proper] Gamma variance)
		{
			if (variance.IsPointMass) return SampleAverageConditional(mean, variance.Point);
			if (mean.IsUniform()) return Gaussian.Uniform();
			double mx,vx,mm,vm;
			double a,b;
			sample.GetMeanAndVariance(out mx, out vx);
			mean.GetMeanAndVariance(out mm, out vm);
			a = variance.Shape;
			b = variance.Rate;
			if (sample.IsUniform()) return new Gaussian(mm, vm + a/b);
			double c = Math.Sqrt(2*b);
			double m = c*(mx-mm);
			double v = c*c*(vx+vm);
			double mu,vu;
			if(a < 1) throw new ArgumentException("variance.shape < 1");
			double m2p = m*m/v;
			if (false && Math.Abs(m*m/v) > 30) {
				double invV = 1/v;
				double invSqrtV = Math.Sqrt(invV);
				double mPlus = (m-v)*invSqrtV;
				double mMinus = (-m-v)*invSqrtV;
				double Zplus = MMath.NormalCdfRatio(mPlus);
				double Zminus = MMath.NormalCdfRatio(mMinus);
				double Z0 = Zminus + Zplus;
				double alpha = (Zminus - Zplus)/Z0;
				double beta = -(1 - 2/Z0*invSqrtV - alpha*alpha);
				alpha *= c;
				beta *= c*c;
				double weight = beta / (sample.Precision - beta);
				Gaussian result = new Gaussian();
				result.Precision = sample.Precision * weight;
				result.MeanTimesPrecision = weight * (sample.MeanTimesPrecision + alpha) + alpha;
				if (!result.IsProper()) throw new ImproperMessageException(result);
				if (double.IsNaN(result.Precision) || double.IsNaN(result.MeanTimesPrecision)) throw new ApplicationException("result is nan");
				return result;
			} else {
				//double logZ;
				//if(a == 1) LaplacianTimesGaussianMoments(m, v, out logZ, out mu, out vu);
				VarianceGammaTimesGaussianMoments5(a, m, v, out mu, out vu);
				//Console.WriteLine("mu = {0}, vu = {1}", mu, vu);
				//VarianceGammaTimesGaussianMoments(a, m, v, out mu, out vu);
				//Console.WriteLine("mu = {0}, vu = {1}", mu, vu);
				double r = vx/(vx+vm);
				double mp = r*(mu/c + mm) + (1-r)*mx;
				double vp = r*r*vu/(c*c) + r*vm;
				if (Double.IsNaN(vp)) throw new Exception("vp is NaN");
				Gaussian result = Gaussian.FromMeanAndVariance(mp, vp)/sample;
				if (!result.IsProper()) throw new ImproperMessageException(result);
				return result;
			}
		}

		public static Gaussian UpperBoundAverageConditional([SkipIfUniform] Bernoulli isBetween, Gaussian X, Gaussian lowerBound, Gaussian upperBound)
		{
			Gaussian result = new Gaussian();
			if(isBetween.IsUniform()) return result;
			if (upperBound.IsPointMass)
			{
				result.SetToUniform(); // TODO: return the limiting distribution
			}
			else 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(true, upperBound, lowerBound, X);
				}
				else
				{
					result.SetToUniform();
				}
			}
			else
			{
				double logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
				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);
				// since upperBound is not a point mass, -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 = Math.Log(invSqrtVxu) + Gaussian.GetLogProb(yu, 0, 1) + MMath.NormalCdfLn((yl - r * yu) / Math.Sqrt(1 - r * r));
				double alphaU = d_p * Math.Exp(logPhiU - logZ);
				// (mu - mx) / (vx + vu) = yu*invSqrtVxu
				double betaU = alphaU * (alphaU + yu * invSqrtVxu);
				if (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;
				}
				double weight = betaU / (upperBound.Precision - betaU);
				if (ForceProper && weight < 0) weight = 0;
				result.Precision = weight * upperBound.Precision;
				result.MeanTimesPrecision = weight * (upperBound.MeanTimesPrecision + alphaU) + alphaU;
			}
			if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision)) throw new ApplicationException("result is NaN");
			return result;
		}

		// Perform one update of Q
		private static Gamma QUpdate(Gaussian sample, Gaussian mean, Gamma precision, Gamma q)
		{
			if (sample.IsUniform() || mean.IsUniform() || precision.IsPointMass) return precision;
			double mx, vx;
			sample.GetMeanAndVariance(out mx, out vx);
			double mm, vm;
			mean.GetMeanAndVariance(out mm, out vm);
			double m = mx-mm;
			double v = vx+vm;
			double m2 = m*m;
			double a = q.Shape;
			double b = q.Rate;
			if (b==0 || q.IsPointMass) {
				a = precision.Shape;
				// this guess comes from solving dlogf=0 for x
				double guess = m2-v;
				b = Math.Max(precision.Rate, a*guess);
			}
			double x = a/b;
			double x2 = x*x;
			double x3 = x*x2;
			double x4 = x*x3;
			double p = 1/(v+1/x);
			double dlogf1 = -0.5*p + 0.5*m2*p*p;
			double dlogf = dlogf1*(-1/x2);
			double ddlogf1 = 0.5*p*p - m2*p*p*p;
			double ddlogf = dlogf1*2/x3 + ddlogf1/x4;
			b = precision.Rate - (dlogf + x*ddlogf);
			if (b < 0) {
				if (q.Rate == precision.Rate || true) return QUpdate(sample, mean, precision, Gamma.FromShapeAndRate(precision.Shape, precision.Shape*(m2-v)));
			}
			a = precision.Shape - x2*ddlogf;
			if (a <= 0) a = b*precision.Shape/(precision.Rate - dlogf);
			if (a <= 0 || b <= 0) throw new Exception();
			if (double.IsNaN(a) || double.IsNaN(b)) throw new Exception("result is nan");
			return Gamma.FromShapeAndRate(a, b);
		}

		/// <summary>
		/// EP message to 'exp'
		/// </summary>
		/// <param name="exp">Incoming message from 'exp'.</param>
		/// <param name="d">Incoming message from 'd'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="to_d">Previous outgoing message to 'd'.</param>
		/// <returns>The outgoing EP message to the 'exp' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'exp' as the random arguments are varied.
		/// The formula is <c>proj[p(exp) sum_(d) p(d) factor(exp,d)]/p(exp)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="d"/> is not a proper distribution</exception>
		public static Gamma ExpAverageConditional(Gamma exp, [Proper] Gaussian d, Gaussian to_d)
		{
			if (d.IsPointMass) return Gamma.PointMass(Math.Exp(d.Point));
			if (d.IsUniform()) return Gamma.FromShapeAndRate(0, 0);
			if (exp.IsPointMass) {
				// Z = int_y delta(x - exp(y)) N(y; my, vy) dy 
				//   = int_u delta(x - u) N(log(u); my, vy)/u du
				//   = N(log(x); my, vy)/x
				// logZ = -log(x) -0.5/vy*(log(x)-my)^2
				// dlogZ/dx = -1/x -1/vy*(log(x)-my)/x
				// d2logZ/dx2 = -dlogZ/dx/x -1/vy/x^2
				// log Ga(x;a,b) = (a-1)*log(x) - bx
				// dlogGa/dx = (a-1)/x - b
				// d2logGa/dx2 = -(a-1)/x^2
				// match derivatives and solve for (a,b)
				double shape = (1 + d.GetMean() - Math.Log(exp.Point)) * d.Precision;
				double rate = d.Precision / exp.Point;
				return Gamma.FromShapeAndRate(shape, rate);
			}
			if (exp.IsUniform()) return ExpAverageLogarithm(d);

			if (to_d.IsUniform() && exp.Shape > 1) {
				to_d = new Gaussian(MMath.Digamma(exp.Shape - 1) - Math.Log(exp.Rate), MMath.Trigamma(exp.Shape - 1));
			}

			double mD, vD;
			Gaussian dMarginal = d * to_d;
			dMarginal.GetMeanAndVariance(out mD, out vD);
			double Z = 0;
			double sumy = 0;
			double sumexpy = 0;

			if (vD < 1e-6) {
				double m, v;
				d.GetMeanAndVariance(out m, out v);
				return Gamma.FromLogMeanAndMeanLog(m + v / 2.0, m);
			}

			//if (vD < 10)
			if (true) {
				// Use Gauss-Hermite quadrature
				double[] nodes = new double[QuadratureNodeCount];
				double[] weights = new double[QuadratureNodeCount];

				Quadrature.GaussianNodesAndWeights(mD, vD, nodes, weights);
				for (int i = 0; i < weights.Length; i++) {
					weights[i] = Math.Log(weights[i]);
				}
				if (!to_d.IsUniform()) {
					// modify the weights to include q(y_k)/N(y_k;m,v)
					for (int i = 0; i < weights.Length; i++) {
						weights[i] += d.GetLogProb(nodes[i]) - dMarginal.GetLogProb(nodes[i]);
					}
				}

				double maxLogF = Double.NegativeInfinity;
				// f(x,y) = Ga(exp(y); shape, rate) = exp(y*(shape-1) -rate*exp(y))
				// Z E[x] = int_y int_x x Ga(x;a,b) delta(x - exp(y)) N(y;my,vy) dx dy
				//        = int_y exp(y) Ga(exp(y);a,b) N(y;my,vy) dy
				// Z E[log(x)] = int_y y Ga(exp(y);a,b) N(y;my,vy) dy
				for (int i = 0; i < weights.Length; i++) {
					double y = nodes[i];
					double logf = weights[i] + (exp.Shape - 1) * y - exp.Rate * Math.Exp(y);
					if (logf > maxLogF) maxLogF = logf;
					weights[i] = logf;
				}
				for (int i = 0; i < weights.Length; i++) {
					double y = nodes[i];
					double f = Math.Exp(weights[i] - maxLogF);
					double f_y = f * y;
					double fexpy = f * Math.Exp(y);
					Z += f;
					sumy += f_y;
					sumexpy += fexpy;
				}
			} else {
				Converter<double, double> p = delegate(double y) {
					return d.GetLogProb(y) + (exp.Shape - 1) * y - exp.Rate * Math.Exp(y);
				};
				double sc = Math.Sqrt(vD);
				double offset = p(mD);
				Z = Quadrature.AdaptiveClenshawCurtis(z => Math.Exp(p(sc * z + mD) - offset), 1, 16, 1e-6);
				sumy = Quadrature.AdaptiveClenshawCurtis(z => (sc * z + mD) * Math.Exp(p(sc * z + mD) - offset), 1, 16, 1e-6);
				sumexpy = Quadrature.AdaptiveClenshawCurtis(z => Math.Exp(sc * z + mD + p(sc * z + mD) - offset), 1, 16, 1e-6);
			}
			if (Z == 0) throw new ApplicationException("Z==0");
			double s = 1.0 / Z;
			if (Double.IsPositiveInfinity(s)) throw new ApplicationException("s is -inf");
//.........这里部分代码省略.........

		/// <summary>
		/// EP message to 'd'
		/// </summary>
		/// <param name="exp">Incoming message from 'exp'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="d">Incoming message from 'd'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'd' as the random arguments are varied.
		/// The formula is <c>proj[p(d) sum_(exp) p(exp) factor(exp,d)]/p(d)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="exp"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="d"/> is not a proper distribution</exception>
		//internal static Gaussian DAverageConditional_slow([SkipIfUniform] Gamma exp, [Proper] Gaussian d)
		//{
		//  Gaussian to_d = exp.Shape<=1 || exp.Rate==0 ? 
		//            Gaussian.Uniform() 
		//            : new Gaussian(MMath.Digamma(exp.Shape-1) - Math.Log(exp.Rate), MMath.Trigamma(exp.Shape));
		//  //var to_d = Gaussian.Uniform();
		//  for (int i = 0; i < QuadratureIterations; i++) {
		//    to_d = DAverageConditional(exp, d, to_d);
		//  }
		//  return to_d;
		//}
		// to_d does not need to be Fresh. it is only used for quadrature proposal.
		public static Gaussian DAverageConditional([SkipIfUniform] Gamma exp, [Proper] Gaussian d, Gaussian result)
		{
			if (exp.IsUniform() || d.IsPointMass) return Gaussian.Uniform();
			if (exp.IsPointMass) return DAverageConditional(exp.Point);
			if (exp.Rate < 0) throw new ImproperMessageException(exp);
			if (d.IsUniform()) {
				if (exp.Shape <= 1) throw new ArgumentException("The posterior has infinite variance due to input of Exp distributed as "+d+" and output of Exp distributed as "+exp+" (shape <= 1)");
				// posterior for d is a shifted log-Gamma distribution:
				// exp((a-1)*d - b*exp(d)) =propto exp(a*(d+log(b)) - exp(d+log(b)))
				// we find the Gaussian with same moments.
				// u = d+log(b)
				// E[u] = digamma(a-1)
				// E[d] = E[u]-log(b) = digamma(a-1)-log(b)
				// var(d) = var(u) = trigamma(a-1)
				double lnRate = Math.Log(exp.Rate);
				return new Gaussian(MMath.Digamma(exp.Shape - 1) - lnRate, MMath.Trigamma(exp.Shape - 1));
			}
			// We use moment matching to find the best Gaussian message.
			// The moments are computed via quadrature.
			// Z = int_y f(x,y) q(y) dy =approx sum_k w_k f(x,y_k) q(y_k)/N(y_k;m,v) 
			// f(x,y) = Ga(exp(y); shape, rate) = exp(y*(shape-1) -rate*exp(y))
			double[] nodes = new double[QuadratureNodeCount];
			double[] weights = new double[QuadratureNodeCount];
			double moD, voD;
			d.GetMeanAndVariance(out moD, out voD);
			double mD, vD;
			if (result.IsUniform() && exp.Shape > 1)
				result = new Gaussian(MMath.Digamma(exp.Shape - 1) - Math.Log(exp.Rate), MMath.Trigamma(exp.Shape - 1));
			Gaussian dMarginal = d * result;
			dMarginal.GetMeanAndVariance(out mD, out vD);
			Quadrature.GaussianNodesAndWeights(mD, vD, nodes, weights);
			if (!result.IsUniform()) {
				// modify the weights to include q(y_k)/N(y_k;m,v)
				for (int i = 0; i < weights.Length; i++) {
					weights[i] *= Math.Exp(d.GetLogProb(nodes[i]) - Gaussian.GetLogProb(nodes[i], mD, vD));
				}
			}
			double Z = 0;
			double sumy = 0;
			double sumy2 = 0;
			double maxLogF = Double.NegativeInfinity;
			for (int i = 0; i < weights.Length; i++) {
				double y = nodes[i];
				double logf = Math.Log(weights[i]) + (exp.Shape - 1) * y - exp.Rate * Math.Exp(y);
				if (logf > maxLogF) maxLogF = logf;
				weights[i] = logf;
			}
			for (int i = 0; i < weights.Length; i++) {
				double y = nodes[i];
				double f = Math.Exp(weights[i] - maxLogF);
				double f_y = f * y;
				double fyy = f_y * y;
				Z += f;
				sumy += f_y;
				sumy2 += fyy;
			}
			if (Z == 0) return Gaussian.Uniform();
			double s = 1.0 / Z;
			double mean = sumy * s;
			double var = sumy2 * s - mean * mean;
			// TODO: explain this
			if (var <= 0.0) {
				double quadratureGap = 0.1;
				var = 2 * vD * quadratureGap * quadratureGap;
			}
			result = new Gaussian(mean, var);
			result.SetToRatio(result, d, ForceProper);
			if (result.Precision < -1e10) throw new ApplicationException("result has negative precision");
			if (Double.IsPositiveInfinity(result.Precision)) throw new ApplicationException("result is point mass");
			if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision)) throw new ApplicationException("result is nan");
			return result;
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="exp">Incoming message from 'exp'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <param name="to_d">Previous outgoing message to 'd'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(exp,d) p(exp,d) factor(exp,d))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Gamma exp, Gaussian d, Gaussian to_d)
		{
			if (d.IsPointMass) return LogAverageFactor(exp, d.Point);
			if (d.IsUniform()) return exp.GetLogAverageOf(new Gamma(0, 0));
			if (exp.IsPointMass) return LogAverageFactor(exp.Point, d);
			if (exp.IsUniform()) return 0.0;
			double[] nodes = new double[QuadratureNodeCount];
			double[] weights = new double[QuadratureNodeCount];
			double mD, vD;
			Gaussian dMarginal = d * to_d;
			dMarginal.GetMeanAndVariance(out mD, out vD);
			Quadrature.GaussianNodesAndWeights(mD, vD, nodes, weights);
			if (!to_d.IsUniform()) {
				// modify the weights to include q(y_k)/N(y_k;m,v)
				for (int i = 0; i < weights.Length; i++) {
					weights[i] *= Math.Exp(d.GetLogProb(nodes[i]) - Gaussian.GetLogProb(nodes[i], mD, vD));
				}
			}
			double Z = 0;
			for (int i = 0; i < weights.Length; i++) {
				double y = nodes[i];
				double f = weights[i] * Math.Exp((exp.Shape - 1) * y - exp.Rate * Math.Exp(y));
				Z += f;
			}
			return Math.Log(Z) - exp.GetLogNormalizer();
		}

		/// <summary>
		/// EP message to 'sample'
		/// </summary>
		/// <param name="sample">Incoming message from '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 'sample' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'sample' as the random arguments are varied.
		/// The formula is <c>proj[p(sample) sum_(mean,precision) p(mean,precision) factor(sample,mean,precision)]/p(sample)</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 Gaussian SampleAverageConditional(Gaussian sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision)
		{
			Gaussian result = new Gaussian();
			if (precision.IsPointMass) {
				return SampleAverageConditional(mean, precision.Point);
			} else if (sample.IsUniform()) {
				// for large vx, Z =approx N(mx; mm, vx+vm+E[1/prec])
				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.
				// Need to make the scheduler smarter to avoid this.
				if(precision.Shape <= 1.0) throw new ArgumentException("The posterior has infinite variance due to precision distributed as "+precision+" (shape <= 1).  Try using a different prior for the precision, with shape > 1.");
				return Gaussian.FromMeanAndVariance(mm, mv + precision.GetMeanInverse());
			} else if (mean.IsUniform() || precision.IsUniform()) {
				result.SetToUniform();
			} else if (sample.IsPointMass) {
				// The correct answer here is not uniform, but rather a limit.  
				// However it doesn't really matter what we return since multiplication by a point mass 
				// always yields a point mass.
				result.SetToUniform();
			} else if (!precision.IsProper()) {
				throw new ImproperMessageException(precision);
			} else {
				// The formula is int_prec int_mean N(x;mean,1/prec) p(x) p(mean) p(prec) =
				// int_prec N(x; mm, mv + 1/prec) p(x) p(prec) =
				// int_prec N(x; new xm, new xv) N(xm; mm, mv + xv + 1/prec) p(prec)
				// Let R = Prec/(Prec + mean.Prec)
				// new xv = inv(R*mean.Prec + sample.Prec)
				// new xm = xv*(R*mean.PM + sample.PM)

				// In the case where sample and mean are improper distributions, 
				// we must only consider values of prec for which (new xv > 0).
				// This happens when R*mean.Prec > -sample.Prec
				// As a function of Prec, R*mean.Prec has a singularity at Prec=-mean.Prec
				// This function is greater than a threshold when Prec is sufficiently small or sufficiently large.
				// Therefore we construct an interval of Precs to exclude from the integration.
				double xm, xv, mm, mv;
				sample.GetMeanAndVarianceImproper(out xm, out xv);
				mean.GetMeanAndVarianceImproper(out mm, out mv);
				double lowerBound = 0;
				double upperBound = Double.PositiveInfinity;
				bool precisionIsBetween = true;
				if (mean.Precision >= 0) {
					if (sample.Precision < -mean.Precision) throw new ImproperMessageException(sample);
					//lowerBound = -mean.Precision * sample.Precision / (mean.Precision + sample.Precision);
					lowerBound = -1.0 / (xv + mv);
				} else {  // mean.Precision < 0
					if (sample.Precision < 0) {
						precisionIsBetween = true;
						lowerBound = -1.0 / (xv + mv);
						upperBound = -mean.Precision;
					} else if (sample.Precision < -mean.Precision) {
						precisionIsBetween = true;
						lowerBound = 0;
						upperBound = -mean.Precision;
					} else {
						// in this case, the precision should NOT be in this interval.
						precisionIsBetween = false;
						lowerBound = -mean.Precision;
						lowerBound = -1.0 / (xv + mv);
					}
				}
				double[] nodes = new double[QuadratureNodeCount];
				double[] weights = new double[nodes.Length];
				QuadratureNodesAndWeights(precision, nodes, weights);
				double Z = 0, rmean = 0, rvariance = 0;
				for (int i = 0; i < nodes.Length; i++) {
					double newVar, newMean;
					Assert.IsTrue(nodes[i] > 0);
					if ((nodes[i] > lowerBound && nodes[i] < upperBound) != precisionIsBetween) continue;
					// the following works even if sample is uniform. (sample.Precision == 0)
					if (mean.IsPointMass) {
						// take limit mean.Precision -> Inf
						newVar = 1.0 / (nodes[i] + sample.Precision);
						newMean = newVar * (nodes[i] * mean.Point + sample.MeanTimesPrecision);
					} else {
						// mean.Precision < Inf
						double R = nodes[i] / (nodes[i] + mean.Precision);
						newVar = 1.0 / (R * mean.Precision + sample.Precision);
						newMean = newVar * (R * mean.MeanTimesPrecision + sample.MeanTimesPrecision);
					}

					double f;
					// If p(x) is uniform, xv=Inf and the term N(xm; mm, mv + xv + 1/prec) goes away
					if (sample.IsUniform())
						f = weights[i];
					else
//.........这里部分代码省略.........

		/// <summary>
		/// EP message to 'sample'
		/// </summary>
		/// <param name="sample">Incoming message from '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 'sample' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'sample' as the random arguments are varied.
		/// The formula is <c>proj[p(sample) sum_(mean,precision) p(mean,precision) factor(sample,mean,precision)]/p(sample)</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 Gaussian SampleAverageConditional(Gaussian sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision, Gamma to_precision)
		{
			if (sample.IsUniform() && precision.Shape <= 1.0) sample = Gaussian.FromNatural(1e-20, 1e-20);
			if (precision.IsPointMass) {
				return SampleAverageConditional(mean, precision.Point);
			} else if (sample.IsUniform()) {
				// for large vx, Z =approx N(mx; mm, vx+vm+E[1/prec])
				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.
				// Need to make the scheduler smarter to avoid this.
				if (precision.Shape <= 1.0) throw new ArgumentException("The posterior has infinite variance due to precision distributed as "+precision+" (shape <= 1).  Try using a different prior for the precision, with shape > 1.");
				return Gaussian.FromMeanAndVariance(mm, mv + precision.GetMeanInverse());
			} else if (mean.IsUniform() || precision.IsUniform()) {
				return Gaussian.Uniform();
			} else if (sample.IsPointMass) {
				// The correct answer here is not uniform, but rather a limit.  
				// However it doesn't really matter what we return since multiplication by a point mass 
				// always yields a point mass.
				return Gaussian.Uniform();
			} else if (!precision.IsProper()) {
				throw new ImproperMessageException(precision);
			} else {
				// The formula is int_prec int_mean N(x;mean,1/prec) p(x) p(mean) p(prec) =
				// int_prec N(x; mm, mv + 1/prec) p(x) p(prec) =
				// int_prec N(x; new xm, new xv) N(xm; mm, mv + xv + 1/prec) p(prec)
				// Let R = Prec/(Prec + mean.Prec)
				// new xv = inv(R*mean.Prec + sample.Prec)
				// new xm = xv*(R*mean.PM + sample.PM)

				// In the case where sample and mean are improper distributions, 
				// we must only consider values of prec for which (new xv > 0).
				// This happens when R*mean.Prec > -sample.Prec
				// As a function of Prec, R*mean.Prec has a singularity at Prec=-mean.Prec
				// This function is greater than a threshold when Prec is sufficiently small or sufficiently large.
				// Therefore we construct an interval of Precs to exclude from the integration.
				double xm, xv, mm, mv;
				sample.GetMeanAndVarianceImproper(out xm, out xv);
				mean.GetMeanAndVarianceImproper(out mm, out mv);
				double lowerBound = 0;
				double upperBound = Double.PositiveInfinity;
				bool precisionIsBetween = true;
				if (mean.Precision >= 0) {
					if (sample.Precision < -mean.Precision) throw new ImproperMessageException(sample);
					//lowerBound = -mean.Precision * sample.Precision / (mean.Precision + sample.Precision);
					lowerBound = -1.0 / (xv + mv);
				} else {  // mean.Precision < 0
					if (sample.Precision < 0) {
						precisionIsBetween = true;
						lowerBound = -1.0 / (xv + mv);
						upperBound = -mean.Precision;
					} else if (sample.Precision < -mean.Precision) {
						precisionIsBetween = true;
						lowerBound = 0;
						upperBound = -mean.Precision;
					} else {
						// in this case, the precision should NOT be in this interval.
						precisionIsBetween = false;
						lowerBound = -mean.Precision;
						lowerBound = -1.0 / (xv + mv);
					}
				}
				double[] nodes = new double[QuadratureNodeCount];
				double[] logWeights = new double[nodes.Length];
				Gamma precMarginal = precision*to_precision;
				QuadratureNodesAndWeights(precMarginal, nodes, logWeights);
				if (!to_precision.IsUniform()) {
					// modify the weights
					for (int i = 0; i < logWeights.Length; i++) {
						logWeights[i] += precision.GetLogProb(nodes[i]) - precMarginal.GetLogProb(nodes[i]);
					}
				}
				GaussianEstimator est = new GaussianEstimator();
				double shift = 0;
				for (int i = 0; i < nodes.Length; i++) {
					double newVar, newMean;
					Assert.IsTrue(nodes[i] > 0);
					if ((nodes[i] > lowerBound && nodes[i] < upperBound) != precisionIsBetween) continue;
					// the following works even if sample is uniform. (sample.Precision == 0)
					if (mean.IsPointMass) {
						// take limit mean.Precision -> Inf
						newVar = 1.0 / (nodes[i] + sample.Precision);
						newMean = newVar * (nodes[i] * mean.Point + sample.MeanTimesPrecision);
					} else {
						// mean.Precision < Inf
						double R = nodes[i] / (nodes[i] + mean.Precision);
						newVar = 1.0 / (R * mean.Precision + sample.Precision);
//.........这里部分代码省略.........

		public static Gamma Q(Gaussian sample, Gaussian mean, Gamma precision, Gamma q)
		{
			if (precision.IsPointMass || sample.IsUniform() || mean.IsUniform()) return precision;
			double mx, vx;
			sample.GetMeanAndVariance(out mx, out vx);
			double mm, vm;
			mean.GetMeanAndVariance(out mm, out vm);
			double m = mx-mm;
			double v = vx+vm;
			if (double.IsInfinity(v)) return precision;
			double m2 = m*m;
			double[] dlogfss;
			double x = q.GetMean();
			double a = precision.Shape;
			double b = precision.Rate;
			if (double.IsPositiveInfinity(x)) x = (a + Math.Sqrt(a*a + 2*b/v))/(2*b);
			// TODO: check initialization
			if (q.IsUniform()) x = 0;
			//x = QReinitialize(sample, mean, precision, x);
			for (int iter = 0; iter < 1000; iter++) {
				double oldx = x;
				if (true) {
					// want to find x that optimizes -0.5*log(v + 1/x) - 0.5*m2/(v + 1/x) + a*log(x) - b*x
					// we fit a lower bound, then optimize the bound.
					// this is guaranteed to improve x.
					double logSlope = precision.Shape;
					double slope = -precision.Rate;
					double denom = v*x+1;
					// 1/(v+1/x) <= 1/(v+1/x0) + (x-x0)/(v*x0+1)^2
					slope += -0.5*m2/(denom*denom);
					if (v*x < 1) {
						// log(v+1/x)  = log(v*x+1) - log(x)
						// log(v*x+1) <= log(v*x0+1) + (x-x0)*v/(v*x0+1)
						logSlope += 0.5;
						slope += -0.5*v/denom;
						x = -logSlope/slope;
						// at x=0 the update is x = (a+0.5)/(b + 0.5*(m2+v))
						// at x=inf the update is x = (a+0.5)/b
					} else {
						// if v*x > 1:
						// log(v+1/x) <= log(v+1/x0) + (1/x - 1/x0)/(v + 1/x0)
						double invSlope = -0.5*x/denom;
						// solve for the maximum of logslope*log(r)+slope*r+invslope./r
						// logslope/r + slope - invslope/r^2 = 0
						// logslope*r + slope*r^2 - invslope = 0
						//x = (-logSlope - Math.Sqrt(logSlope*logSlope + 4*invSlope*slope))/(2*slope);
						double c = 0.5*logSlope/slope;
						double d = invSlope/slope;
						// note c < 0 always
						x = Math.Sqrt(c*c + d) - c;
						// at x=inf, invSlope=-0.5/v so the update is x = (ax + sqrt(ax*ax + 2*b/v))/(2*b)
					}
					if (x < 0) throw new Exception("x < 0");
				} else {
					dlogfss = dlogfs(x, m, v);
					if (true) {
						x = (precision.Shape - x*x*dlogfss[1])/(precision.Rate - dlogfss[0] - x*dlogfss[1]);
					} else if (true) {
						double delta = dlogfss[0] - precision.Rate;
						x *= Math.Exp(-(delta*x + precision.Shape)/(delta*x + dlogfss[1]*x));
					} else {
						x = precision.Shape/(precision.Rate - dlogfss[0]);
					}
					if (x < 0) throw new Exception("x < 0");
				}
				if (double.IsNaN(x)) throw new Exception("x is nan");
				//Console.WriteLine("{0}: {1}", iter, x);
				if (MMath.AbsDiff(oldx, x, 1e-10) < 1e-10) {
					x = QReinitialize(sample, mean, precision, x);
					if (MMath.AbsDiff(oldx, x, 1e-10) < 1e-10) break;
				}
				if (iter == 1000-1) throw new Exception("not converging");
			}
			//x = r0;
			dlogfss = dlogfs(x, m, v);
			double dlogf = dlogfss[0];
			double ddlogf = dlogfss[1];
			a = precision.Shape - x*x*ddlogf;
			b = precision.Rate - (dlogf + x*ddlogf);
			if (a <= 0 || b <= 0) throw new Exception();
			if (double.IsNaN(a) || double.IsNaN(b)) throw new Exception("result is nan");
			return Gamma.FromShapeAndRate(a, b);
		}