C# Beta.IsUniform方法代码示例

		/// <summary>
		/// Find the Laplace approximation for Beta(Logistic(x)) * Gaussian(x))
		/// </summary>
		/// <param name="beta">Beta distribution</param>
		/// <param name="gauss">Gaussian distribution</param>
		/// <returns>A proposal distribution</returns>
		public static Gaussian LogisticProposalDistribution(Beta beta, Gaussian gauss)
		{
			if (beta.IsUniform())
				return new Gaussian(gauss);

			// if gauss is uniform, m,p = 0 below, and the following code will just ignore the Gaussian
			// and do a Laplace approximation for Beta(Logistic(x))

			double c = beta.TrueCount-1;
			double d = beta.FalseCount-1;
			double m = gauss.GetMean();
			double p = gauss.Precision;
			// We want to find the mode of
			// ln(g(x)) = c.ln(f(x)) + d.ln(1 - f(x)) - 0.5p((x - m)^2) + constant
			// First deriv:
			// h(x) = (ln(g(x))' = c.(1 - f(x)) - d.f(x) - p(x-m)
			// Second deriv:
			// h'(x) = (ln(g(x))' = -(c+d).f'(x) - p
			// Use Newton-Raphson to find unique root of h(x).
			// g(x) is log-concave so Newton-Raphson should converge quickly.
			// Set the initial point by projecting beta
			// to a Gaussian and taking the mean of the product:
			double bMean, bVar;
			beta.GetMeanAndVariance(out bMean, out bVar);
			Gaussian prod = new Gaussian();
			double invLogisticMean = Math.Log(bMean) - Math.Log(1.0-bMean);
			prod.SetToProduct(Gaussian.FromMeanAndVariance(invLogisticMean, bVar), gauss);
			double xnew = prod.GetMean();
			double x=0, fx, dfx, hx, dhx=0;
			int maxIters = 100; // Should only need a handful of iters
			int cnt = 0;
			do {
				x = xnew;
				fx = MMath.Logistic(x);
				dfx = fx * (1.0-fx);
				// Find the root of h(x)
				hx = c * (1.0 - fx) - d * fx - p*(x-m);
				dhx = -(c+d)*dfx - p;
				xnew = x - (hx / dhx); // The Newton step
				if (Math.Abs(x - xnew) < 0.00001)
					break;
			} while (++cnt < maxIters);
			if (cnt >= maxIters)
				throw new ApplicationException("Unable to find proposal distribution mode");
			return Gaussian.FromMeanAndPrecision(x, -dhx);
		}

		/// <summary>
		/// EP message to 'logistic'
		/// </summary>
		/// <param name="logistic">Incoming message from 'logistic'.</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>The outgoing EP message to the 'logistic' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'logistic' as the random arguments are varied.
		/// The formula is <c>proj[p(logistic) sum_(x) p(x) factor(logistic,x)]/p(logistic)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="x"/> is not a proper distribution</exception>
		public static Beta LogisticAverageConditional(Beta logistic, [Proper] Gaussian x, Gaussian falseMsg)
		{
			if (x.IsPointMass)
				return Beta.PointMass(MMath.Logistic(x.Point));

			if (logistic.IsPointMass || x.IsUniform())
				return Beta.Uniform();

			double m,v;
			x.GetMeanAndVariance(out m, out v);
			if ((logistic.TrueCount == 2 && logistic.FalseCount == 1) ||
				  (logistic.TrueCount == 1 && logistic.FalseCount == 2) ||
				   logistic.IsUniform()) {
				// shortcut for the common case
				// result is a Beta distribution satisfying:
				// int_p to_p(p) p dp = int_x sigma(x) qnoti(x) dx
				// int_p to_p(p) p^2 dp = int_x sigma(x)^2 qnoti(x) dx
				// the second constraint can be rewritten as:
				// int_p to_p(p) p (1-p) dp = int_x sigma(x) (1 - sigma(x)) qnoti(x) dx
				// the constraints are the same if we replace p with (1-p)
				double mean = MMath.LogisticGaussian(m, v);
				// meanTF = E[p] - E[p^2]
				double meanTF = MMath.LogisticGaussianDerivative(m, v);
				double meanSquare = mean - meanTF;
				return Beta.FromMeanAndVariance(mean, meanSquare - mean*mean);
			} else {
				// stabilized EP message
				// choose a normalized distribution to_p such that:
				// int_p to_p(p) qnoti(p) dp = int_x qnoti(sigma(x)) qnoti(x) dx
				// int_p to_p(p) p qnoti(p) dp = int_x qnoti(sigma(x)) sigma(x) qnoti(x) dx
				double logZ = LogAverageFactor(logistic, x, falseMsg) + logistic.GetLogNormalizer(); // log int_x logistic(sigma(x)) N(x;m,v) dx
				Gaussian post = XAverageConditional(logistic, falseMsg) * x;
				double mp,vp;
				post.GetMeanAndVariance(out mp, out vp);
				double tc1 = logistic.TrueCount-1;
				double fc1 = logistic.FalseCount-1;
				double Ep;
				if (tc1+fc1 == 0) {
					Beta logistic1 = new Beta(logistic.TrueCount+1, logistic.FalseCount);
					double logZp = LogAverageFactor(logistic1, x, falseMsg) + logistic1.GetLogNormalizer();
					Ep = Math.Exp(logZp - logZ);
				} else {
					// Ep = int_p to_p(p) p qnoti(p) dp / int_p to_p(p) qnoti(p) dp
					// mp = m + v (a - (a+b) Ep)
					Ep = (tc1 - (mp - m)/v)/(tc1+fc1);
				}
				return BetaFromMeanAndIntegral(Ep, logZ, tc1, fc1);
			}
		}