C# Gaussian.GetMeanAndVariance方法代码示例

		// logw1 = N(mx;m1,vx+v1) phi((mx1 - m2)/sqrt(vx1+v2))
		// a1 = N(mx1;m2,vx1+v2)/phi
		internal static void ComputeStats(Gaussian max, Gaussian a, Gaussian b, out double logz,
			out double logw1, out double a1, out double vx1, out double mx1,
			out double logw2, out double a2, out double vx2, out double mx2)
		{
			double m1, v1, m2, v2;
			a.GetMeanAndVariance(out m1, out v1);
			b.GetMeanAndVariance(out m2, out v2);
			if (max.IsPointMass) {
				vx1 = 0.0;
				mx1 = max.Point;
				vx2 = 0.0;
				mx2 = max.Point;
				if (b.IsPointMass) {
					if (b.Point > max.Point) throw new AllZeroException();
					else if (b.Point == max.Point) {
						// the factor reduces to the constraint (max.Point > a)
						logw1 = Double.NegativeInfinity;
						logw2 = MMath.NormalCdfLn((max.Point - m1)/Math.Sqrt(v1));
						logz = logw2;
						a1 = 0;
						a2 = Math.Exp(Gaussian.GetLogProb(max.Point, m1, v1) - logw2);
						return;
					} else {
						// b.Point < max.Point
						// the factor reduces to the constraint (a == max.Point)
						throw new NotImplementedException();
					}
				} else if (a.IsPointMass) throw new NotImplementedException();
			} else {
				if (a.IsPointMass) {
					vx1 = 0.0;
					mx1 = a.Point;
				} else {
					vx1 = 1.0 / (max.Precision + a.Precision);
					mx1 = vx1 * (max.MeanTimesPrecision + a.MeanTimesPrecision);
				}
				if (b.IsPointMass) {
					vx2 = 0.0;
					mx2 = b.Point;
				} else {
					vx2 = 1.0 / (max.Precision + b.Precision);
					mx2 = vx2 * (max.MeanTimesPrecision + b.MeanTimesPrecision);
				}
			}
			logw1 = max.GetLogAverageOf(a);
			double logPhi1 = MMath.NormalCdfLn((mx1 - m2) / Math.Sqrt(vx1 + v2));
			logw1 += logPhi1;

			logw2 = max.GetLogAverageOf(b);
			double logPhi2 = MMath.NormalCdfLn((mx2 - m1) / Math.Sqrt(vx2 + v1));
			logw2 += logPhi2;

			logz = MMath.LogSumExp(logw1, logw2);

			double logN1 = Gaussian.GetLogProb(mx1, m2, vx1 + v2);
			a1 = Math.Exp(logN1 - logPhi1);
			double logN2 = Gaussian.GetLogProb(mx2, m1, vx2 + v1);
			a2 = Math.Exp(logN2 - logPhi2);
		}

		public static Gaussian SampleAverageConditional(Gaussian sample, [SkipIfUniform] Gaussian mean, [Proper] Gamma precision, [Fresh] Gamma q)
		{
			if (mean.IsUniform() || sample.IsPointMass) return Gaussian.Uniform();
			if (precision.IsPointMass) return GaussianOp.SampleAverageConditional(mean, precision.Point);
			if (q.IsPointMass) throw new Exception();
			double mm, vm;
			mean.GetMeanAndVariance(out mm, out vm);
			if (sample.IsUniform()) {
				if (precision.Shape <= 1.0) {
					sample = Gaussian.FromNatural(1e-20, 1e-20);  // try to proceed instead of throwing an exception
					//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.");
				} else return Gaussian.FromMeanAndVariance(mm, vm + precision.GetMeanInverse());
			}
			double mx, vx;
			sample.GetMeanAndVariance(out mx, out vx);
			double m = mx-mm;
			double v = vx+vm;
			Gaussian sampleMarginal;
			if (true) {
				Gaussian qx = Qx(Gaussian.FromMeanAndVariance(m, v), precision, Gaussian.FromMeanAndPrecision(0, 1));
				double y = qx.GetMean();
				double[] dlogfx = dlogfxs(y, precision);
				double delta = 0.5*dlogfx[2]/(qx.Precision*qx.Precision);
				// sampleMean can be computed directly as:
				//double dlogz = dlogfx[0] + delta/v;
				//double sampleMean = mx + vx*dlogz;
				double my = y + delta;
				double vy = 1/qx.Precision + 4*delta*delta + 0.5*dlogfx[3]/(qx.Precision*qx.Precision*qx.Precision);
				if (vy < 0) throw new Exception();
				Gaussian yMsg = new Gaussian(my, vy)/(new Gaussian(m, v));
				Gaussian sampleMsg = DoublePlusOp.SumAverageConditional(yMsg, mean);
				sampleMarginal = sampleMsg*sample;
				if (!sampleMarginal.IsProper()) throw new Exception();
				if (sampleMarginal.IsPointMass) throw new Exception();
				//Console.WriteLine("{0} {1}", sampleMean, sampleMarginal);
			} else {
				double a = q.Shape;
				double b = q.Rate;
				double x = a/b;
				if (sample.IsPointMass || mean.IsPointMass) {
					double denom = 1 + x*v;
					double denom2 = denom*denom;
					Gaussian y = sample*mean;
					double my, vy;
					y.GetMeanAndVariance(out my, out vy);
					// 1-1/denom = x*v/denom
					// sampleMean = E[ (x*(mx*vm + mm*vx) + mx)/denom ]
					//            = E[ (mx*vm + mm*vx)/v*(1-1/denom) + mx/denom ]
					//            = E[ (mx/vx + mm/vm)/(1/vx + 1/vm)*(1-1/denom) + mx/denom ]
					// sampleVariance = var((my*(1-1/denom) + mx/denom)) + E[ (r*vx*vm + vx)/denom ]
					//                = (mx-my)^2*var(1/denom) + E[ vx*vm/v*(1-1/denom) + vx/denom ]
					double[] g = new double[] { 1/denom, -v/denom2, 2*v*v/(denom2*denom), -6*v*v*v/(denom2*denom2) };
					double [] dlogf = dlogfs(x, m, v);
					double edenom, vdenom;
					LaplaceMoments(q, g, dlogf, out edenom, out vdenom);
					double sampleMean = mx*edenom + my*(1 - edenom);
					double diff = mx-my;
					double sampleVariance = vx*edenom + vy*(1 - edenom) + diff*diff*vdenom;
					sampleMarginal = Gaussian.FromMeanAndVariance(sampleMean, sampleVariance);
				} else {
					// 1 - samplePrec*mPrec/denom = x*yprec/denom
					// sampleMean = E[ (x*ymprec + sampleMP*mPrec)/denom ]
					//            = E[ (1 - samplePrec*mPrec/denom)*ymprec/yprec + sampleMP*mPrec/denom ]
					// sampleVariance = var((1 - samplePrec*mPrec/denom)*ymprec/yprec + sampleMP*mPrec/denom) + E[ (x+mPrec)/denom ]
					//                = (sampleMP*mPrec - samplePrec*mPrec*ymprec/yprec)^2*var(1/denom) + E[ (1-samplePrec*mPrec/denom)/yprec + mPrec/denom ]
					double yprec = sample.Precision + mean.Precision;
					double ymprec = sample.MeanTimesPrecision + mean.MeanTimesPrecision;
					double denom = sample.Precision*mean.Precision + x*yprec;
					double denom2 = denom*denom;
					double[] g = new double[] { 1/denom, -yprec/denom2, 2*yprec*yprec/(denom2*denom), -6*yprec*yprec*yprec/(denom2*denom2) };
					double [] dlogf = dlogfs(x, m, v);
					double edenom, vdenom;
					LaplaceMoments(q, g, dlogf, out edenom, out vdenom);
					double sampleMean = sample.MeanTimesPrecision*mean.Precision*edenom + ymprec/yprec*(1-sample.Precision*mean.Precision*edenom);
					double diff = sample.MeanTimesPrecision*mean.Precision - sample.Precision*mean.Precision*ymprec/yprec;
					double sampleVariance = mean.Precision*edenom + (1 - sample.Precision*mean.Precision*edenom)/yprec + diff*diff*vdenom;
					sampleMarginal = Gaussian.FromMeanAndVariance(sampleMean, sampleVariance);
				}
			}
			Gaussian result = new Gaussian();
			result.SetToRatio(sampleMarginal, sample, GaussianOp.ForceProper);
			if (modified2) {
				if(!mean.IsPointMass) throw new Exception();
				result.SetMeanAndPrecision(mm, q.GetMean());
			} else if (modified && !sample.IsUniform()) {
				// heuristic method to avoid improper messages:
				// the message's mean must be E[mean] (regardless of context) and the variance is chosen to match the posterior mean when multiplied by context
				double sampleMean = sampleMarginal.GetMean();
				if (sampleMean != mm) {
					result.Precision = (sample.MeanTimesPrecision-sampleMean*sample.Precision)/(sampleMean - mm);
					if (result.Precision < 0) throw new Exception("internal: sampleMean is not between sample.Mean and mean.Mean");
					result.MeanTimesPrecision = result.Precision*mm;
				}
			}
			//if (!result.IsProper()) throw new Exception();
			if (result.Precision < -0.001) throw new Exception();
			if (double.IsNaN(result.MeanTimesPrecision) || double.IsNaN(result.Precision)) throw new Exception("result is nan");
			return result;
		}

 public static double AverageLogFactor(Bernoulli sample, Gaussian logOdds)
 {
     // This is the non-conjugate VMP update using the Saul and Jordan (1999) bound.
     double m, v;
     logOdds.GetMeanAndVariance(out m, out v);
     double a = 0.5;
     // TODO: use a buffer to store the value of 'a', so it doesn't need to be re-optimised each time.
     for (int iter = 0; iter < 10; iter++)
     {
         double aOld = a;
         a = MMath.Logistic(m + (1 - 2 * a) * v * 0.5);
         if (Math.Abs(a - aOld) < 1e-8) break;
     }
     return sample.GetProbTrue() * m - .5 * a * a * v - MMath.Log1PlusExp(m + (1 - 2 * a) * v * 0.5);
 }

		/// <summary>
		/// Find a proposal distribution - we will use this to set the limits of
		/// integration
		/// </summary>
		/// <param name="y">Gaussian distribution for y</param>
		/// <param name="a">Gaussian distribution for a</param>
		/// <param name="b">Beta distribution for b</param>
		/// <returns>A proposal distribution</returns>
		public static Gaussian ProposalDistribution(Gaussian y, Gaussian a, Beta b)
		{
			double ma, my, va, vy;
			y.GetMeanAndVariance(out my, out vy);
			a.GetMeanAndVariance(out ma, out va);
			Gaussian g;
			// If ma is zero, just build the proposal distribution on the beta
			if (ma == 0.0)
				g = Gaussian.Uniform();
			else
				g = Gaussian.FromMeanAndVariance(my/ma, (vy+va) / ma);
			return LogisticProposalDistribution(b, g);
		}

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

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

		public static Gaussian AAverageConditional(Gaussian Product, Gaussian A, Gaussian B, Gaussian to_A)
		{
			Gaussian Apost = A*to_A;
			double mx,vx;
			Product.GetMeanAndVariance(out mx, out vx);
			double mb,vb;
			B.GetMeanAndVariance(out mb, out vb);
			double ma,va;
			bool usePost = false;
			if (usePost) Apost.GetMeanAndVariance(out ma, out va);
			else A.GetMeanAndVariance(out ma, out va);
			double ahat = ma;
			ahat = Apost.GetMean();
			double denom = vx + ahat*ahat*vb;
			double diff = mx - ahat*mb;
			double dd = 2*ahat*vb;
			double ddd = 2*vb;
			double n = diff*diff;
			double dn = -2*diff*mb;
			double ddn = 2*mb*mb;
			double dlogf = -0.5*(dd+dn)/denom + 0.5*n*dd/(denom*denom);
			double ddlogf = -0.5*(ddd+ddn)/denom + 0.5*(dd*dd+2*dd*dn+n*ddd)/(denom*denom) - n*dd*dd/(denom*denom*denom);
			if (usePost) {
				double m0,v0;
				A.GetMeanAndVariance(out m0, out v0);
				dlogf += (m0 - ahat)/v0;
				ddlogf += 1/va - 1/v0;
				// at a fixed point, dlogf2 = 0 so (ahat - m0)/v0 = dlogf, ahat = m0 + v0*dlogf
				// this is the same fixed point condition as Laplace
			}
			// Ef'/Ef =approx dlogf*f(ahat)/f(ahat) = dlogf
			double mnew = ma + va*dlogf;
			double vnew = va*(1 + va*(ddlogf - dlogf*dlogf));
			double rnew = Math.Max(1/vnew, A.Precision);
			return Gaussian.FromNatural(rnew*mnew - A.MeanTimesPrecision, rnew - A.Precision);
		}

		public static Gaussian ProductAverageConditional(Gaussian Product, Gaussian A, Gaussian B, Gaussian to_A, Gaussian to_B)
		{
			if (Product.IsUniform()) {
				double ma,va;
				A.GetMeanAndVariance(out ma, out va);
				double mb,vb;
				B.GetMeanAndVariance(out mb, out vb);
				// this is the limiting value of below - it excludes va*vb
				return Gaussian.FromMeanAndVariance(ma*mb, ma*ma*vb + mb*mb*va);
			}
			Gaussian Apost = A*to_A;
			Gaussian Bpost = B*to_B;
			double ahat = Apost.GetMean();
			double bhat = Bpost.GetMean();
			double gx = -(Product.MeanTimesPrecision - ahat*bhat*Product.Precision);
			double gxx = -Product.Precision;
			double gax = bhat*Product.Precision;
			double gbx = ahat*Product.Precision;
			double gaa = -bhat*bhat*Product.Precision;
			double gbb = -ahat*ahat*Product.Precision;
			double gab = Product.MeanTimesPrecision - 2*ahat*bhat*Product.Precision;
			if (modified) {
				double gaaa = 0;
				double gaaab = 0;
				double gaaaa = 0;
				double gaab = -2*bhat*Product.Precision;
				double gaabb = -2*Product.Precision;
				double gabb = -2*ahat*Product.Precision;
				double gabbb = 0;
				double gbbb = 0;
				double gbbbb = 0;
				double adiff = A.Precision - gaa;
				double bdiff = B.Precision - gbb;
				double h = adiff*bdiff;
				double ha = -gaaa*bdiff - gabb*adiff;
				double hb = -gaab*bdiff - gbbb*adiff;
				double haa = -gaaaa*bdiff - gaabb*adiff + 2*gaaa*gabb;
				double hab = -gaaab*bdiff - gabbb*adiff + gaaa*gbbb + gabb*gaab;
				double hbb = -gaabb*bdiff - gbbbb*adiff + 2*gaab*gbbb;
				if (offDiagonal) {
					h += -gab*gab;
					ha += -2*gab*gaab;
					hb += -2*gab*gabb;
					haa += -2*gaab*gaab - 2*gab*gaaab;
					hab += -2*gabb*gaab - 2*gab*gaabb;
					hbb += -2*gabb*gabb - 2*gab*gabbb;
				}
				double logha = ha/h;
				double loghb = hb/h;
				double loghaa = haa/h - logha*logha;
				double loghab = hab/h - logha*loghb;
				double loghbb = hbb/h - loghb*loghb;
				gaa -= 0.5*loghaa;
				gab -= 0.5*loghab;
				gbb -= 0.5*loghbb;
			}
			double cb = gab/(B.Precision - gbb);
			double dax = (gax + gbx*cb)/(A.Precision - (gaa + gab*cb));
			double dbx = gbx/(B.Precision - gbb) + cb*dax;
			double dlogz = gx;
			double ddlogz = gxx + gax*dax + gbx*dbx;
			return GaussianProductOp_Laplace.GaussianFromAlphaBeta(Product, dlogz, -ddlogz, true);
		}

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

		public static double LogAverageFactor([SkipIfUniform] Gaussian Product, Gaussian A, Gaussian B, Gaussian to_product)
		{
			if (A.IsPointMass) return GaussianProductEvidenceOp.LogAverageFactor(Product, B, A.Point, to_product);
			if (B.IsPointMass) return GaussianProductEvidenceOp.LogAverageFactor(Product, A, B.Point, to_product);
			if (Product.IsUniform()) return 0.0;
			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;
			for (int i = 0; i <= GaussianProductOp.QuadratureNodeCount; i++) {
				double a = (2.0 * i) / GaussianProductOp.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 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;
			}
			double inc = 2.0/GaussianProductOp.QuadratureNodeCount;
			return Math.Log(z*inc);
		}

		public static Vector GradientVector(Gaussian Product, Gaussian A, Gamma B)
		{
			double rateb = B.Rate;
			double scaleb = 1.0/rateb;
			double shapeb = B.Shape;
			double ma, va;
			A.GetMeanAndVariance(out ma, out va);

			Matrix C = new Matrix(2, 2);
			C[0, 0] = scaleb;
			C[0, 1] = (2*shapeb+1)*scaleb*scaleb;
			C[1, 0] = -shapeb*scaleb*scaleb;
			C[1, 1] = -2*shapeb*(shapeb+1)*scaleb*scaleb*scaleb;
			Vector v = new Vector(2);
			v[0] = ma*Product.MeanTimesPrecision;
			v[1] = -0.5*(ma*ma+va)*Product.Precision;
			return C*v;
		}

		/// <summary>
		/// Helper method for computing average log factor
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'.</param>
		/// <param name="mean">Constant value for 'mean'.</param>
		/// <param name="precision_Elogx">Expected log value of the incoming message from 'precision'</param>
		/// <param name="precision_Ex">Expected value of incoming message from 'precision'</param>
		/// <returns>Computed average log factor</returns>
		private static double ComputeAverageLogFactor(Gaussian sample, double mean, double precision_Elogx, double precision_Ex)
		{
			if (double.IsPositiveInfinity(precision_Ex)) throw new ArgumentException("precision is infinite");
			double sampleMean, sampleVariance;
			sample.GetMeanAndVariance(out sampleMean, out sampleVariance);
			double diff = sampleMean - mean;
			return -MMath.LnSqrt2PI + 0.5 * (precision_Elogx
											- precision_Ex * (diff * diff + sampleVariance));
		}

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