C# Gaussian.GetMean方法代码示例

		//-- Constant bounds --------------------------------------------------------------------------------

		/// <summary>
		/// The logarithm of the probability that L &lt;= X &lt; U.
		/// </summary>
		/// <param name="X"></param>
		/// <param name="L">Can be negative infinity.</param>
		/// <param name="U">Can be positive infinity.</param>
		/// <returns></returns>
		public static double LogProbBetween(Gaussian X, double L, double U)
		{
			if (L > U) throw new AllZeroException("low > high (" + L + " > " + U + ")");
			if (X.IsPointMass)
			{
				return Factor.IsBetween(X.Point, L, U) ? 0.0 : Double.NegativeInfinity;
			}
			else if (X.IsUniform())
			{
				if (Double.IsNegativeInfinity(L))
				{
					if (Double.IsPositiveInfinity(U)) return 0.0; // always between
					else return -MMath.Ln2;  // between half the time
				}
				else if (Double.IsPositiveInfinity(U)) return -MMath.Ln2;  // between half the time
				else return Double.NegativeInfinity;  // never between two finite numbers
			}
			else
			{
				double sqrtPrec = Math.Sqrt(X.Precision);
				double mx = X.GetMean();
				double pl = MMath.NormalCdfLn(sqrtPrec * (L - mx));  // log(p(x <= L))
				double pu = MMath.NormalCdfLn(sqrtPrec * (U - mx));  // log(p(x <= U))
				if (pl == pu) return Double.NegativeInfinity;
				if (Double.IsNegativeInfinity(pl)) return pu;
				// log(NormalCdf(yu) - NormalCdf(yl)) = NormalCdfLn(yu) + log(1 - NormalCdf(yl)/NormalCdf(yu))
				return pu + MMath.Log1MinusExp(pl - pu);
			}
		}

		/// <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>
		/// VMP message to 'b'
		/// </summary>
		/// <param name="ProductExp">Incoming message from 'productExp'. 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>
		/// <param name="to_B">Previous outgoing message to 'B'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></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, 'productExp' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(a) p(a) log(sum_productExp p(productExp) factor(productExp,a,b)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="ProductExp"/> 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 NonconjugateGaussian BAverageLogarithm([SkipIfUniform] Gaussian ProductExp, [Proper, SkipIfUniform] Gaussian A, [Proper, SkipIfUniform] Gaussian B, NonconjugateGaussian to_B, NonconjugateGaussian result)
        {
            if (B.IsPointMass) return NonconjugateGaussian.Uniform();
            if (ProductExp.IsPointMass) return BAverageLogarithm(ProductExp.Point, A);
            if (!B.IsProper()) throw new ImproperMessageException(B);
            // catch uniform case to avoid 0*Inf
            if (ProductExp.IsUniform()) return NonconjugateGaussian.Uniform();
            double mx, vx, m, v, mz, vz;
            ProductExp.GetMeanAndVariance(out mz, out vz);
            A.GetMeanAndVariance(out mx, out vx);
            B.GetMeanAndVariance(out m, out v);

            //if (mx * mz < 0)
            //{
            //    Console.WriteLine("Warning: mx*mz < 0, setting to uniform");
            //    result.SetToUniform();
            //    return result;
            //}

            double Ex2 = mx * mx + vx;
            double grad2_S_m2 = -2 * Ex2 * Math.Exp(2 * m + 2 * v) / vz + mx * mz * Math.Exp(m + .5 * v) / vz;
            double grad_m = -Ex2 * Math.Exp(2 * m + 2 * v) / vz + mx * mz * Math.Exp(m + .5 * v) / vz;
            double threshold = 10;
            double mf, vf, afm1, bf;
            if (grad2_S_m2 >= -threshold && mx * mz > 0)
            {
                mf = Math.Log(mx * mz / Ex2) - 1.5 * v;
                vf = (mf - m) / grad_m;
            }
            else
            {
                vf = -1 / grad2_S_m2;
                mf = m - grad_m / grad2_S_m2;
            }
            
            Gaussian priorG;
            if (result.IsUniform())
                priorG = B;
            else
            {
                var prior = new NonconjugateGaussian();
                prior.SetToRatio((new NonconjugateGaussian(B)), to_B);
                priorG = prior.GetGaussian();
            }

            result.MeanTimesPrecision = mf / vf ;
            result.Precision = 1 / vf;

            var updatedM = new Gaussian(mf,vf) * priorG;
            m = updatedM.GetMean(); 
            
            double grad_S2_v2 = -2 * Ex2 * Math.Exp(2 * m + 2 * v) / vz + .25 * mx * mz * Math.Exp(m + .5 * v) / vz;
            double grad_S_v = -Ex2 * Math.Exp(2 * m + 2 * v) / vz + .5 * mx * mz * Math.Exp(m + .5 * v) / vz;

            afm1 = -1;
            bf = -1;
            if (grad2_S_m2 >= -threshold)
            {
                afm1 = -v * v * grad_S2_v2;
                bf = -grad_S_v + afm1 / v;
            }

            if ((afm1 < 0 || bf < 0) && mx * mz > 0)
            {
                double v_opt = 2 / 3 * (Math.Log(mx * mz / Ex2 / 2) - m);
                if (v_opt != v)
                {
                    bf = v * grad_S_v / (v_opt - v);
                    afm1 = v_opt * bf;
                }
            }

            if (afm1 < 0 || bf < 0)
            {
                afm1 = -v * v * grad_S2_v2;
                bf = -grad_S_v + afm1 / v;
            }

            if (afm1 < 0 || bf < 0)
            {
                result.Shape = 1;
                result.Rate = 0; 
            }
//.........这里部分代码省略.........

		public static double XInit(Gaussian d)
		{
			return d.GetMean();
		}

 private static VectorGaussian ShapeLocationTimesFactor(
     Vector point, Gaussian shapeX, Gaussian shapeY, PositiveDefiniteMatrix shapeOrientation)
 {
     VectorGaussian shapeLocationDistr = VectorGaussian.FromMeanAndVariance(
         Vector.FromArray(shapeX.GetMean(), shapeY.GetMean()),
         new PositiveDefiniteMatrix(new double[,] { { shapeX.GetVariance(), 0.0 }, { 0.0, shapeY.GetVariance() } }));
     VectorGaussian factorDistribution = VectorGaussian.FromMeanAndPrecision(point, shapeOrientation);
     VectorGaussian result = new VectorGaussian(2);
     result.SetToProduct(shapeLocationDistr, factorDistribution);
     return result;
 }

		public static void LaplaceMoments(Gaussian q, double[] dlogfx, out double m, out double v)
		{
			double vx = 1/q.Precision;
			double delta = 0.5*dlogfx[2]*vx*vx;
			m = q.GetMean() + delta;
			v = vx + 4*delta*delta + 0.5*dlogfx[3]*vx*vx*vx;
			if (v < 0) throw new Exception();
		}

		private static double QxReinitialize(Gaussian y, Gamma precision, double x)
		{
			double init0 = 0;
			double init1 = y.GetMean();
			double a = precision.Shape;
			double b = precision.Rate;
			double a2 = -(a+0.5);
			double logz0 = a2*Math.Log(b + init0*init0/2) + y.GetLogProb(init0);
			double logz1 = a2*Math.Log(b + init1*init1/2) + y.GetLogProb(init1);
			double logz = a2*Math.Log(b + x*x/2) + y.GetLogProb(x);
			if (logz0 > Math.Max(logz1, logz)) return init0;
			else if (logz1 > logz) return init1;
			else return x;
		}

		public static Gaussian Qx(Gaussian y, Gamma precision, Gaussian qx)
		{
			if (y.IsPointMass) return y;
			double x = QxReinitialize(y, precision, qx.GetMean());
			double r = 0;
			for (int iter = 0; iter < 1000; iter++) {
				double oldx = x;
				double[] dlogfs = dlogfxs(x, precision);
				double ddlogf = dlogfs[1];
				r = Math.Max(0, -ddlogf);
				double t = r*x + dlogfs[0];
				x = (t + y.MeanTimesPrecision)/(r + y.Precision);
				if (Math.Abs(oldx - x) < 1e-10) break;
				//Console.WriteLine("{0}: {1}", iter, x);		
				if (iter == 1000-1) throw new Exception("not converging");
				if (iter % 100 == 99) x = QxReinitialize(y, precision, x);
			}
			return Gaussian.FromMeanAndPrecision(x, r + y.Precision);
		}