C# Gaussian.SetToRatioProper方法代码示例

		/// <summary>
		/// EP message to 'a'
		/// </summary>
		/// <param name="max">Incoming message from 'max'. 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>
		/// <param name="b">Incoming message from 'b'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>The outgoing EP message to the 'a' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'a' as the random arguments are varied.
		/// The formula is <c>proj[p(a) sum_(max,b) p(max,b) factor(max,a,b)]/p(a)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="max"/> is not a proper distribution</exception>
		/// <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 AAverageConditional([SkipIfUniform] Gaussian max, [Proper] Gaussian a, [Proper] Gaussian b)
		{
			if (a.IsPointMass) return Gaussian.Uniform();
			if (max.IsUniform()) return Gaussian.Uniform();
			if (!a.IsProper()) throw new ImproperMessageException(a);
			if (!b.IsProper()) throw new ImproperMessageException(b);
			double m1, v1, m2, v2;
			a.GetMeanAndVariance(out m1, out v1);
			b.GetMeanAndVariance(out m2, out v2);

#if true
			double logw1, a1, vx1, mx1;
			double logw2, a2, vx2, mx2;
			double logz;
			ComputeStats(max, a, b, out logz, out logw1, out a1, out vx1, out mx1,
				out logw2, out a2, out vx2, out mx2);
			double w1 = Math.Exp(logw1 - logz);
			double w2 = Math.Exp(logw2 - logz);
			// the posterior is a mixture model with weights exp(logw1-logz), exp(logw2-logz) and distributions
			// N(a; mx1, vx1) phi((a - m2)/sqrt(v2)) / phi((mx1 - m2)/sqrt(vx1 + v2))
			// N(a; m1, v1) phi((mx2 - a)/sqrt(vx2)) / phi((mx2 - m1)/sqrt(vx2 + v1))
			// the moments of the posterior are computed via the moments of these two components.
			double mc1 = mx1 + a1 * vx1;
			a2 = -a2;
			double mc2 = m1 + a2 * v1;
			double m = w1 * mc1 + w2 * mc2;
			double r1 = (mx1 - m2) / (vx1 + v2);
			double r2 = (mx2 - m1) / (vx2 + v1);
			double beta1 = a1 * (a1 + r1);
			double beta2 = a2 * (a2 - r2);
			double vc1 = vx1 * (1 - vx1 * beta1);
			double vc2 = v1 * (1 - v1 * beta2);
			double diff = mc1 - mc2;
			double v = w1 * vc1 + w2 * vc2 + w1 * w2 * diff * diff;
			Gaussian result = new Gaussian(m, v);
			if (ForceProper) result.SetToRatioProper(result, a);
			else result.SetToRatio(result, a);
#else
			double logw1, logPhi1, logu1, vx1, mx1;
			double logw2, logPhi2, logu2, vx2, mx2;
			double logz;
			ComputeStats(max, a, b, out logz, out logw1, out logPhi1, out logu1, out vx1, out mx1,
				out logw2, out logPhi2, out logu2, out vx2, out mx2);
			if (a.IsPointMass || max.IsPointMass || b.IsPointMass) throw new NotImplementedException();
			//double ra = (mx-ma)/(vx+va);
			double ra = (max.MeanTimesPrecision * a.Precision - a.MeanTimesPrecision * max.Precision) / (a.Precision + max.Precision);
			//double rb = (mx-mb)/(vx+vb);
			double rb = (max.MeanTimesPrecision * b.Precision - b.MeanTimesPrecision * max.Precision) / (b.Precision + max.Precision);
			double vxinva = a.Precision / (a.Precision + max.Precision);
			double vxinvb = b.Precision / (b.Precision + max.Precision);
			double mx, vx;
			max.GetMeanAndVariance(out mx, out vx);
			double invxa = 1.0 / (vx + v1);
			double invxb = 1.0 / (vx + v2);
			double alpha = Math.Exp(logw1 + logPhi1 - logz) * ra + Math.Exp(logw1 + logu1 - logz) * vxinva - Math.Exp(logw2 + logu2 - logz);
			double dvx = Math.Exp(logw1 + logPhi1 - logz) * (ra * ra - invxa) + Math.Exp(logw1 + logu1 - logz) * (
				2 * ra * vxinva - (mx1 - m2) / (vx1 + v2) * vxinva * vxinva) - Math.Exp(logw2 + logu2 - logz) * (mx2 - m1) / (vx2 + v1);
			double beta = alpha * alpha - dvx;

			double tau = a.MeanTimesPrecision;
			double prec = a.Precision;
			double weight = beta / (prec - beta);
			if (ForceProper && weight < 0) weight = 0;
			Gaussian result = new Gaussian();
			// 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;
#endif
			if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision)) throw new Exception("result is NaN");
			return result;
		}

		/// <summary>
		/// EP message to 'max'
		/// </summary>
		/// <param name="max">Incoming message from 'max'.</param>
		/// <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 EP message to the 'max' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'max' as the random arguments are varied.
		/// The formula is <c>proj[p(max) sum_(a,b) p(a,b) factor(max,a,b)]/p(max)</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 MaxAverageConditional(Gaussian max, [Proper] Gaussian a, [Proper] Gaussian b)
		{
			if (max.IsPointMass) return Gaussian.Uniform();
			// the following code works correctly even if max is uniform or improper.
			if (!a.IsProper()) throw new ImproperMessageException(a);
			if (!b.IsProper()) throw new ImproperMessageException(b);
			double m1, v1, m2, v2;
			a.GetMeanAndVariance(out m1, out v1);
			b.GetMeanAndVariance(out m2, out v2);

			double logw1, a1, vx1, mx1;
			double logw2, a2, vx2, mx2;
			double logz;
			ComputeStats(max, a, b, out logz, out logw1, out a1, out vx1, out mx1,
				out logw2, out a2, out vx2, out mx2);
			double w1 = Math.Exp(logw1 - logz);
			double w2 = Math.Exp(logw2 - logz);
			// the posterior is a mixture model with weights exp(logw1-logz), exp(logw2-logz) and distributions
			// N(x; mx1, vx1) phi((x - m2)/sqrt(v2)) / phi((mx1 - m2)/sqrt(vx1 + v2))
			// the moments of the posterior are computed via the moments of these two components.
			double mc1 = mx1 + a1 * vx1;
			double mc2 = mx2 + a2 * vx2;
			double m = w1 * mc1 + w2 * mc2;
			double r1 = (mx1 - m2) / (vx1 + v2);
			double r2 = (mx2 - m1) / (vx2 + v1);
			double beta1 = a1 * (a1 + r1);
			double beta2 = a2 * (a2 + r2);
			double vc1 = vx1 * (1 - vx1 * beta1);
			double vc2 = vx2 * (1 - vx2 * beta2);
			double diff = mc1 - mc2;
			double v = w1 * vc1 + w2 * vc2 + w1 * w2 * diff * diff;
			Gaussian result = new Gaussian(m, v);
			if (ForceProper) result.SetToRatioProper(result, max);
			else result.SetToRatio(result, max);
			if (Double.IsNaN(result.Precision) || Double.IsNaN(result.MeanTimesPrecision)) throw new Exception("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;
		}

		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>
		/// 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()) {
				// 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;
			if (var <= 0.0) {
				double quadratureGap = 0.1;
				var = 2 * vD * quadratureGap * quadratureGap;
			}
			result = new Gaussian(mean, var);
			if (ForceProper) result.SetToRatioProper(result, d);
			else result.SetToRatio(result, d);
			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;
		}