C# Gaussian.SetMeanAndVariance方法代码示例

		/// <summary>
		/// VMP message to 'product'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="B">Incoming message from 'b'. 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 'product' as the random arguments are varied.
		/// The formula is <c>proj[sum_(b) p(b) factor(product,a,b)]</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian ProductAverageLogarithm(double A, [SkipIfUniform] Beta B)
		{
			double mb, vb;
			B.GetMeanAndVariance(out mb, out vb);	
			Gaussian result = new Gaussian();
			result.SetMeanAndVariance(A * mb, A * A * vb);
			return result;
		}

		/// <summary>
		/// VMP message to 'innerProduct'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</param>
		/// <returns>The outgoing VMP message to the 'innerProduct' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'innerProduct' conditioned on the given values.
		/// </para></remarks>
		public static Gaussian InnerProductAverageLogarithm(Vector A, Vector BMean, PositiveDefiniteMatrix BVariance)
		{
			Gaussian result = new Gaussian();
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			// p(x) = N(a' E[b], a' var(b) a)
			result.SetMeanAndVariance(A.Inner(BMean), BVariance.QuadraticForm(A));
			return result;
		}

		/// <summary>
		/// VMP message to 'Sum'
		/// </summary>
		/// <param name="A">Incoming message from 'A'.</param>
		/// <param name="B">Constant value for 'B'.</param>
		/// <returns>The outgoing VMP message to the 'Sum' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'Sum' as the random arguments are varied.
		/// The formula is <c>proj[sum_(A) p(A) factor(Sum,A,B)]</c>.
		/// </para></remarks>
		public static Gaussian SumAverageLogarithm(DistributionStructArray<Bernoulli, bool> A, [SkipIfUniform] Vector B)
		{
			Gaussian result = new Gaussian();
			// p(x|a,b) = N(E[a]'*E[b], E[b]'*var(a)*E[b] + E[a]'*var(b)*E[a] + trace(var(a)*var(b)))
			Vector ma = Vector.Zero(A.Count);
			Vector va = Vector.Zero(A.Count);
			for (int i = 0; i < A.Count; i++) {
				ma[i] = A[i].GetMean();
				va[i] = A[i].GetVariance();
			}
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			var MeanOfBSquared = Vector.Zero(B.Count);
			MeanOfBSquared.SetToFunction(B, x => x * x);
			result.SetMeanAndVariance(ma.Inner(B), va.Inner(MeanOfBSquared));
			return result;
		}

		/// <summary>
		/// VMP message to 'X'
		/// </summary>
		/// <param name="A">Incoming message from 'A'. Must be a proper distribution.  If all elements are uniform, the result will be uniform.</param>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If all elements are uniform, the result will be uniform.</param>
		/// <param name="MeanOfB">Buffer 'MeanOfB'.</param>
		/// <param name="CovarianceOfB">Buffer 'CovarianceOfB'.</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 'X' as the random arguments are varied.
		/// The formula is <c>proj[sum_(A,B) p(A,B) factor(X,A,B)]</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="A"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian XAverageLogarithm([SkipIfAllUniform] GaussianArray A, [SkipIfAllUniform] VectorGaussian B, Vector MeanOfB, PositiveDefiniteMatrix CovarianceOfB)
		{
			int K = MeanOfB.Count;
			// p(x|a,b) = N(E[a]'*E[b], E[b]'*var(a)*E[b] + E[a]'*var(b)*E[a] + trace(var(a)*var(b)))
			var ma = Vector.Zero(K);
			var va = Vector.Zero(K);
			for (int k = 0; k < K; k++) {
				double m, v;
				A[k].GetMeanAndVariance(out m, out v);
				ma[k] = m;
				va[k] = v;
			}
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			var mbj2 = Vector.Zero(K);
			mbj2.SetToFunction(MeanOfB, x => x * x);
			// slooow
			Gaussian result = new Gaussian();
			result.SetMeanAndVariance(ma.Inner(MeanOfB), va.Inner(mbj2) + CovarianceOfB.QuadraticForm(ma) + va.Inner(CovarianceOfB.Diagonal()));
			if (result.Precision < 0)
				throw new ApplicationException("improper message");

			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>
		/// 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;
		}

//.........这里部分代码省略.........
			} 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
						f = weights[i] * Math.Exp(Gaussian.GetLogProb(xm, mm, xv + mv + 1.0 / nodes[i]));
					double fm = f * newMean;
					double fmm = f * (newVar + newMean * newMean);
					Z += f;
					rmean += fm;
					rvariance += fmm;
				}
				if (Z == 0.0) {
					//throw new Exception("Quadrature failed");
					//Console.WriteLine("Warning: Quadrature found zero mass.  Results may be inaccurate.");
					result.SetToUniform();
					return result;
				}
				double s = 1.0 / Z;
				rmean *= s;
				rvariance = rvariance * s - rmean * rmean;
				if (Double.IsInfinity(rmean)) {
					result.SetToUniform();
				} else {
					result.SetMeanAndVariance(rmean, rvariance);
					if (ForceProper) result.SetToRatioProper(result, sample);
					else result.SetToRatio(result, sample);
				}
			}
			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;
		}

		public static Gaussian AAverageConditional([SkipIfUniform] Gaussian Product, Gaussian A, [SkipIfUniform] Gaussian B)
		{
			if (B.IsPointMass) return AAverageConditional(Product, B.Point);
			if (A.IsPointMass || Product.IsUniform()) return Gaussian.Uniform();
			Gaussian result = new Gaussian();
			// algorithm: quadrature on A from -1 to 1, plus quadrature on 1/A from -1 to 1.
			double mProduct, vProduct;
			Product.GetMeanAndVariance(out mProduct, out vProduct);
			double mA, vA;
			A.GetMeanAndVariance(out mA, out vA);
			double mB, vB;
			B.GetMeanAndVariance(out mB, out vB);
			double z = 0, sumA = 0, sumA2 = 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;
				sumA += a * fA;
				sumA2 += a * a * 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;
				sumA += a * fInvA;
				sumA2 += a * a * fInvA;
			}
			double mean = sumA / z;
			double var = sumA2 / z - mean * mean;
			result.SetMeanAndVariance(mean, var);
			if (ForceProper) result.SetToRatioProper(result, A);
			else result.SetToRatio(result, A);
			return result;
		}

		/// <summary>
		/// VMP message to 'product'
		/// </summary>
		/// <param name="A">Incoming message from 'a'. 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 'product' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'product' as the random arguments are varied.
		/// The formula is <c>proj[sum_(a,b) p(a,b) factor(product,a,b)]</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="A"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian ProductAverageLogarithm([SkipIfUniform] Gaussian A, [SkipIfUniform] Gaussian B)
		{
			Gaussian result = new Gaussian();
			// p(x|a,b) = N(E[a]*E[b], E[b]^2*var(a) + E[a]^2*var(b) + var(a)*var(b))
			double ma, va, mb, vb;
			A.GetMeanAndVariance(out ma, out va);
			B.GetMeanAndVariance(out mb, out vb);
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			result.SetMeanAndVariance(ma * mb, mb * mb * va + ma * ma * vb + va * vb);
			return result;
		}

		/// <summary>
		/// VMP 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="to_d">Previous outgoing message to 'd'.</param>
		/// <returns>The outgoing VMP message to the 'd' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'd' with 'exp' integrated out.
		/// The formula is <c>sum_exp p(exp) factor(exp,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>
		public static Gaussian DAverageLogarithm([Proper] Gamma exp, [Proper, Stochastic] Gaussian d, Gaussian to_d)
		{
			if (exp.IsPointMass) return ExpOp.DAverageLogarithm(exp.Point);

			double m, v;
			d.GetMeanAndVariance(out m, out v);
			Gaussian msg = new Gaussian();
			double mu, s2;
			var prior = d / to_d;
			prior.GetMeanAndVariance(out mu, out s2);
			var z = Vector.Zero(2);
			z[0] = m;
			z[1] = Math.Log(v);
			double startingValue = GradientAndValueAtPoint(mu, s2, z, exp.Shape, exp.Rate, null);
			var s = new BFGS();
			int evalCounter = 0;
			s.MaximumStep = 1e3;
			s.MaximumIterations = 100;
			s.Epsilon = 1e-5;
			s.convergenceCriteria = BFGS.ConvergenceCriteria.Objective;
			z = s.Run(z, 1.0, delegate(Vector y, ref Vector grad) { evalCounter++; return GradientAndValueAtPoint(mu, s2, y, exp.Shape, exp.Rate, grad); });
			m = z[0];
			v = Math.Exp(z[1]);
			to_d.SetMeanAndVariance(m, v);
			to_d.SetToRatio(to_d, prior);
			double endValue = GradientAndValueAtPoint(mu, s2, z, exp.Shape, exp.Rate, null);
			//Console.WriteLine("Went from {0} to {1} in {2} steps, {3} evals", startingValue, endValue, s.IterationsPerformed, evalCounter);
			if (startingValue < endValue)
				Console.WriteLine("Warning: BFGS resulted in an increased objective function");
			return to_d;

			/* ---------------- NEWTON ITERATION VERSION 1 ------------------- 
			double meanTimesPrec, prec;
			d.GetNatural(out meanTimesPrec, out prec);
			Matrix K = new Matrix(2, 2);
			K[0, 0]=1/prec; // d2K by dmu^2
			K[1, 0]=K[0, 1]=-meanTimesPrec/(prec*prec);
			K[1, 1]=meanTimesPrec*meanTimesPrec/Math.Pow(prec, 3)+1/(2*prec*prec);
			double[,,] Kprime = new double[2, 2, 2];
			Kprime[0, 0, 0]=0;
			Kprime[0, 0, 1]=Kprime[0, 1, 0]=Kprime[1, 0, 0]=-1/(prec*prec);
			Kprime[0, 1, 1]=Kprime[1, 1, 0]=Kprime[1, 0, 1]=2*meanTimesPrec/Math.Pow(prec, 3);
			Kprime[1, 1, 1]=-3*meanTimesPrec*meanTimesPrec/Math.Pow(prec, 4)-1/Math.Pow(prec, 3);
			Vector gradS = new Vector(2);
			gradS[0]=(exp.Shape-1)/prec-exp.Rate/prec*Math.Exp((meanTimesPrec+.5)/prec);
			gradS[1]=-(exp.Shape-1)*meanTimesPrec/(prec*prec)+exp.Rate*(meanTimesPrec+.5)/(prec*prec)*Math.Exp((meanTimesPrec+.5)/prec);
			Matrix grad2S = new Matrix(2, 2);
			grad2S[0, 0]=-exp.Rate/(prec*prec)*Math.Exp((meanTimesPrec+.5)/prec);
			grad2S[0, 1]=grad2S[1, 0]=-(exp.Shape-1)/(prec*prec)+exp.Rate*(1/(prec*prec)+(meanTimesPrec+.5)/Math.Pow(prec, 3))*Math.Exp((meanTimesPrec+.5)/prec);
			grad2S[1, 1]=2*(exp.Shape-1)*meanTimesPrec/Math.Pow(prec, 3)-exp.Rate*(meanTimesPrec+.5)/(prec*prec)*(2/prec+(meanTimesPrec+.5)/(prec*prec))*Math.Exp((meanTimesPrec+.5)/prec);
			Vector phi = new Vector(new double[] { result.MeanTimesPrecision, result.Precision });
			Vector gradKL = K*phi-gradS;
			Matrix hessianKL = K - grad2S;
			for (int i=0; i<2; i++)
							for (int j=0; j<2; j++)
											for (int k=0; k<2; k++)
															hessianKL[i, j]+=Kprime[i, j, k]*phi[k];
			double step = 1;
			Vector direction = GammaFromShapeAndRate.twoByTwoInverse(hessianKL)*gradKL;
			Vector newPhi = phi - step * direction;
			result.SetNatural(newPhi[0], newPhi[1]);
			return result;

			---------------- NEWTON ITERATION VERSION 2 ------------------- 
			double mean, variance;
			d.GetMeanAndVariance(out mean, out variance); 
			Gaussian prior = d / result; 
			double mean1, variance1;
			prior.GetMeanAndVariance(out mean1, out variance1); 
			Vector gradKL = new Vector(2);
			gradKL[0]=-(exp.Shape-1)+exp.Rate*Math.Exp(mean+variance/2)+mean/variance1-mean1/variance1;
			gradKL[1]=-1/(2*variance)+exp.Rate*Math.Exp(mean+variance/2)+1/(2*variance1);
			Matrix hessianKL = new Matrix(2, 2);
			hessianKL[0, 0]=exp.Rate*Math.Exp(mean+variance/2)+1/variance1;
			hessianKL[0, 1]=hessianKL[1, 0]=.5*exp.Rate*Math.Exp(mean+variance/2);
			hessianKL[1, 1]=1/(2*variance*variance)+exp.Rate*Math.Exp(mean+variance/2)/4;
			result.GetMeanAndVariance(out mean, out variance);
			if (double.IsInfinity(variance))
							variance=1000;
			Vector theta = new Vector(new double[] { mean, variance });
			theta -= GammaFromShapeAndRate.twoByTwoInverse(hessianKL)*gradKL;
			result.SetMeanAndVariance(theta[0], theta[1]);
			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 'log'
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</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 'log' as the random arguments are varied.
		/// The formula is <c>proj[sum_(d) p(d) factor(log,d)]</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="d"/> is not a proper distribution</exception>
		public static Gaussian LogAverageLogarithm(Gaussian log, [SkipIfUniform] Gamma d, Gaussian result)
		{
			if (d.IsPointMass)
				return LogAverageLogarithm(log, d.Point, result);
			result.SetMeanAndVariance(d.GetMeanLog(), MMath.Trigamma(d.Shape));
			return result;
		}

		// Get click observations for each chunk and label class
		static private Gaussian[][][] getClickObservations(int numLabs, int chunkSize, int[] labels, int[] clicks, int[] exams)
		{

			int nData = labels.Length;
			int numChunks = (nData + chunkSize - 1) / chunkSize;
			Gaussian[][][] chunks = new Gaussian[numChunks][][];
			int[] obsX = new int[numLabs];

			int startChunk = 0;
			int endChunk = 0;
			for (int c = 0; c < numChunks; c++) {
				startChunk = endChunk;
				endChunk = startChunk + chunkSize;
				if (endChunk > nData)
					endChunk = nData;

				int[] labCnts = getLabelCounts(numLabs, labels, startChunk, endChunk);
				chunks[c] = new Gaussian[numLabs][];
				Gaussian[][] currChunk = chunks[c];
				for (int l = 0; l < numLabs; l++) {
					currChunk[l] = new Gaussian[labCnts[l]];
					obsX[l] = 0;
				}

				for (int d = startChunk; d < endChunk; d++) {
					int lab = labels[d];
					int nC = clicks[d];
					int nE = exams[d];
					int nNC = nE - nC;
					double b0 = 1.0 + nC;  // Observations of clicks
					double b1 = 1.0 + nNC;   // Observations of no clicks
					Beta b = new Beta(b0, b1);
					double m, v;
					b.GetMeanAndVariance(out m, out v);
					Gaussian g = new Gaussian();
					g.SetMeanAndVariance(m, v);
					currChunk[lab][obsX[lab]++] = g;
				}
			}
			return chunks;
		}