C# Gaussian.SetToUniform方法代码示例

		public static Gaussian XAverageConditional_Helper([SkipIfUniform] Bernoulli isPositive, [SkipIfUniform, Proper] Gaussian x, bool forceProper)
		{
			Gaussian result = new Gaussian();
			if (x.IsPointMass) {
				result.SetToUniform();
				return result;
			}
			double prec = x.Precision;
			if (prec == 0.0) {
				result.SetToUniform();
				return result;
			} else if (prec < 0) {
				throw new ImproperMessageException(x);
			}
			double sqrtPrec = Math.Sqrt(prec);
			double tau = x.MeanTimesPrecision;
			// m/sqrt(v) = (m/v)/sqrt(1/v)
			double z = tau / sqrtPrec;
			// epsilon = p(b=F)
			// eq (51) in EP quickref
			double alpha;
			if (isPositive.IsPointMass) {
				if (isPositive.Point) {
					if (z < -100) {
						double Y = MMath.NormalCdfRatio(z);
						double d2Y = 2*MMath.NormalCdfMomentRatio(2, z);
						double m2TimesPrec = d2Y/Y;
						Assert.IsTrue(tau != 0);
						double mp = (m2TimesPrec - 1)/tau;
						double vp = m2TimesPrec/prec - mp*mp;
						return (new Gaussian(mp, vp))/x;
					}
					alpha = sqrtPrec / MMath.NormalCdfRatio(z);
				} else {
					if (z > 100) {
						double Y = MMath.NormalCdfRatio(-z);
						// dY = -(d2Y/prec - Y)/(-z)*sqrtPrec
						// dY/Y/prec = -(d2Y/Y/prec/prec - 1/prec)/(-z)*sqrtPrec
						// dY/Y/prec = -(d2Y/Y/prec - 1)/(-tau)
						//double dY = -MMath.NormalCdfMomentRatio(1,-z)*sqrtPrec;
						double d2Y = 2*MMath.NormalCdfMomentRatio(2, -z);
						double m2TimesPrec = d2Y/Y;
						Assert.IsTrue(tau != 0);
						double mp = (m2TimesPrec - 1)/tau;
						double vp = m2TimesPrec/prec - mp*mp;
						return (new Gaussian(mp, vp))/x;
					}
					alpha = -sqrtPrec / MMath.NormalCdfRatio(-z);
				}
			} else {
				//double v = MMath.LogSumExp(isPositive.LogProbTrue + MMath.NormalCdfLn(z), isPositive.LogProbFalse + MMath.NormalCdfLn(-z));
				double v = LogAverageFactor(isPositive, x);
				alpha = sqrtPrec * Math.Exp(-z * z * 0.5 - MMath.LnSqrt2PI - v) * (2 * isPositive.GetProbTrue() - 1);
			}
			// eq (52) in EP quickref (where tau = mnoti/Vnoti)
			double beta = alpha * (alpha + tau);
			double weight = beta / (prec - beta);
			if (forceProper && weight < 0) weight = 0;
			// eq (31) in EP quickref; same as inv(inv(beta)-inv(prec))
			result.Precision = prec * weight;
			// eq (30) in EP quickref times above and simplified
			result.MeanTimesPrecision = weight * (tau + alpha) + alpha;
			if (double.IsNaN(result.Precision) || double.IsNaN(result.MeanTimesPrecision)) throw new ApplicationException("result is nan");
			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;
		}

		/// <summary>
		/// Update the buffer 'falseMsg'
		/// </summary>
		/// <param name="logistic">Incoming message from 'logistic'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="x">Incoming message from 'x'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="falseMsg">Buffer 'falseMsg'.</param>
		/// <returns>New value of buffer 'falseMsg'</returns>
		/// <remarks><para>
		/// 
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="logistic"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="x"/> is not a proper distribution</exception>
		public static Gaussian FalseMsg([SkipIfUniform] Beta logistic, [Proper] Gaussian x, Gaussian falseMsg)
		{
			// falseMsg approximates sigma(-x)
			// logistic(sigma(x)) N(x;m,v)
			// = sigma(x)^(a-1) sigma(-x)^(b-1) N(x;m,v) 
			// = e^((a-1)x) sigma(-x)^(a+b-2) N(x;m,v)
			// = sigma(-x)^(a+b-2) N(x;m+(a-1)v,v) exp((a-1)m + (a-1)^2 v/2)
			// = sigma(-x) (prior)
			// where prior = sigma(-x)^(a+b-3) N(x;m+(a-1)v,v)
			double tc1 = logistic.TrueCount-1;
			double fc1 = logistic.FalseCount-1;
			double m,v;
			x.GetMeanAndVariance(out m, out v);
			if (tc1+fc1 == 0) {
				falseMsg.SetToUniform();
				return falseMsg;
			} else if (tc1+fc1 < 0) {
				// power EP update, using 1/sigma(-x) as the factor
				Gaussian prior = new Gaussian(m + tc1*v, v) * (falseMsg^(tc1+fc1+1));
				double mprior,vprior;
				prior.GetMeanAndVariance(out mprior, out vprior);
				// posterior moments can be computed exactly
				double w = MMath.Logistic(mprior+0.5*vprior);
				Gaussian post = new Gaussian(mprior + w*vprior, vprior*(1 + w*(1-w)*vprior));
				return prior/post;
			} else {
				// power EP update
				Gaussian prior = new Gaussian(m + tc1*v, v) * (falseMsg^(tc1+fc1-1));
				Gaussian newMsg = BernoulliFromLogOddsOp.LogOddsAverageConditional(false, prior);
				//Console.WriteLine("prior = {0}, falseMsg = {1}, newMsg = {2}", prior, falseMsg, newMsg);
				if (true) {
					// adaptive damping scheme
					Gaussian ratio = newMsg/falseMsg;
					if ((ratio.MeanTimesPrecision < 0 && prior.MeanTimesPrecision > 0) ||
						(ratio.MeanTimesPrecision > 0 && prior.MeanTimesPrecision < 0)) {
						// if the update would change the sign of the mean, take a fractional step so that the new prior has exactly zero mean
						// newMsg = falseMsg * (ratio^step)
						// newPrior = prior * (ratio^step)^(tc1+fc1-1)
						// 0 = prior.mp + ratio.mp*step*(tc1+fc1-1)
						double step = -prior.MeanTimesPrecision/(ratio.MeanTimesPrecision*(tc1+fc1-1));
						if (step > 0 && step < 1) {
							newMsg = falseMsg * (ratio^step);
							// check that newPrior has zero mean
							//Gaussian newPrior = prior * ((ratio^step)^(tc1+fc1-1));
							//Console.WriteLine(newPrior);
						}
					}
				}
				return newMsg;
			}
		}

		/// <summary>
		/// EP message to 'x'
		/// </summary>
		/// <param name="isBetween">Incoming message from 'isBetween'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="X">Incoming message from 'x'.</param>
		/// <param name="lowerBound">Constant value for 'lowerBound'.</param>
		/// <param name="upperBound">Constant value for 'upperBound'.</param>
		/// <returns>The outgoing EP message to the 'x' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'x' as the random arguments are varied.
		/// The formula is <c>proj[p(x) sum_(isBetween) p(isBetween) factor(isBetween,x,lowerBound,upperBound)]/p(x)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="isBetween"/> is not a proper distribution</exception>
		public static Gaussian XAverageConditional([SkipIfUniform] Bernoulli isBetween, Gaussian X, double lowerBound, double upperBound)
		{
			Gaussian result = new Gaussian();
			//return XAverageConditional(isBetween, X, Gaussian.PointMass(lowerBound), Gaussian.PointMass(upperBound), result);
			if (X.IsPointMass)
			{
				result.SetToUniform();
			}
			else if (X.IsUniform())
			{
				if (Double.IsInfinity(lowerBound) || Double.IsInfinity(upperBound) ||
					!Double.IsPositiveInfinity(isBetween.LogOdds))
				{
					result.SetToUniform();
				}
				else
				{
					double diff = upperBound - lowerBound;
					result.SetMeanAndVariance((lowerBound + upperBound) / 2, diff * diff / 12);
				}
			}
			else
			{
				// X is not a point mass or uniform
				double mx, vx;
				X.GetMeanAndVariance(out mx, out vx);
				if (double.IsPositiveInfinity(upperBound) && !double.IsInfinity(lowerBound)) {
					Gaussian Xshifted = new Gaussian(mx - lowerBound, vx);
					Gaussian XshiftedPost = Xshifted*IsPositiveOp.XAverageConditional(isBetween, Xshifted);
					double mt,vt;
					XshiftedPost.GetMeanAndVariance(out mt, out vt);
					Gaussian XPost = new Gaussian(mt + lowerBound, vt);
					return XPost/X;
				}
				double logZ = LogAverageFactor(isBetween, X, lowerBound, upperBound);
				double d_p = 2 * isBetween.GetProbTrue() - 1;
				double logPhiL = Gaussian.GetLogProb(lowerBound, mx, vx);
				double alphaL = d_p * Math.Exp(logPhiL - logZ);
				double logPhiU = Gaussian.GetLogProb(upperBound, mx, vx);
				double alphaU = d_p * Math.Exp(logPhiU - logZ);
				double alphaX = alphaL - alphaU;
#if false
				// minka: testing numerical accuracy
				double diff = upperBound - lowerBound;
				double center = (lowerBound+upperBound)/2;
				double logNdiff = diff*(-X.MeanTimesPrecision + center*X.Precision);
				//logNdiff = logPhiL - logPhiU;
				if(logNdiff >= 0) {
					alphaX = d_p*Math.Exp(MMath.LogExpMinus1(logNdiff) + logPhiU-logZ);
				} else {
					alphaX = -d_p*Math.Exp(MMath.LogExpMinus1(-logNdiff) + logPhiL-logZ);
				}
				double m = mx + vx*alphaX;
#endif
				double betaX = alphaX * alphaX;
				if (alphaU != 0.0) betaX += (upperBound - mx) / vx * alphaU;
				if (alphaL != 0.0) betaX -= (lowerBound - mx) / vx * alphaL;
				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;
		}

		/// <summary>
		/// VMP message to 'b'
		/// </summary>
		/// <param name="Difference">Incoming message from 'difference'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="a">Incoming message from 'a'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>The outgoing VMP message to the 'b' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'b'.
		/// Because the factor is deterministic, 'difference' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(a) p(a) log(sum_difference p(difference) factor(difference,a,b)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="Difference"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="a"/> is not a proper distribution</exception>
		public static Gaussian BAverageLogarithm([SkipIfUniform, Proper] Gaussian Difference, [Proper] Gaussian a)
		{
			Gaussian result = new Gaussian();
			if (Difference.IsUniform()) {
				result.SetToUniform();
			} else if (Difference.Precision < 0) {
				throw new ImproperMessageException(Difference);
			} else {
				// p(b|diff,a) = N(E[a] - E[diff], var(diff) )
				double ms, vs;
				double ma = a.GetMean();
				Difference.GetMeanAndVariance(out ms, out vs);
				result.SetMeanAndVariance(ma - ms, vs);
			}
			return result;
		}

		/// <summary>
		/// VMP message to 'A'
		/// </summary>
		/// <param name="Sum">Incoming message from 'Sum'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="b">Incoming message from 'B'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>The outgoing VMP message to the 'A' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'A'.
		/// Because the factor is deterministic, 'Sum' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(B) p(B) log(sum_Sum p(Sum) factor(Sum,A,B)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="Sum"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="b"/> is not a proper distribution</exception>
		public static Gaussian AAverageLogarithm([SkipIfUniform, Proper] Gaussian Sum, [Proper] Gaussian b)
		{
			Gaussian result = new Gaussian();
			if (Sum.IsUniform()) {
				result.SetToUniform();
			} else if (Sum.Precision < 0) {
				throw new ImproperMessageException(Sum);
			} else {
				// p(a|sum,b) = N(E[sum] - E[b], var(sum) )
				double ms, vs;
				double mb = b.GetMean();
				Sum.GetMeanAndVariance(out ms, out vs);
				result.SetMeanAndVariance(ms - mb, vs);
			}
			return result;
		}

		/// <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
//.........这里部分代码省略.........