C# PositiveDefiniteMatrix.QuadraticForm方法代码示例

		/// <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">Constant value for 'A'.</param>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="MeanOfB">Buffer 'MeanOfB'.</param>
		/// <param name="CovarianceOfB">Buffer 'CovarianceOfB'.</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_(B) p(B) factor(Sum,A,B)]</c>.
		/// </para><para>
		/// Uses John Winn's rule for deterministic factors.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian SumAverageLogarithm(bool[] A, [SkipIfUniform] VectorGaussian B, Vector MeanOfB, PositiveDefiniteMatrix CovarianceOfB)
		{
			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.FromArray(A.Select(x => x?1.0:0.0).ToArray());
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			result.SetMeanAndVariance(ma.Inner(MeanOfB), CovarianceOfB.QuadraticForm(ma));
			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;
		}

        private static Gaussian ShapeAverageConditional(
            Vector point, Bernoulli label, Gaussian shapeX, Gaussian shapeY, PositiveDefiniteMatrix shapeOrientation, bool resultForXCoord)
        {
            if (shapeX.IsPointMass && shapeY.IsPointMass)
            {
                double labelProbTrue = label.GetProbTrue();
                double labelProbFalse = 1.0 - labelProbTrue;
                double probDiff = labelProbTrue - labelProbFalse;

                Vector shapeLocation = Vector.FromArray(shapeX.Point, shapeY.Point);
                Vector diff = point - shapeLocation;
                Vector orientationTimesDiff = shapeOrientation * diff;
                Matrix orientationTimesDiffOuter = orientationTimesDiff.Outer(orientationTimesDiff);

                double factorValue = Math.Exp(-0.5 * shapeOrientation.QuadraticForm(diff));
                double funcValue = factorValue * probDiff + labelProbFalse;

                Vector dFunc = probDiff * factorValue * orientationTimesDiff;
                Vector dLogFunc = 1.0 / funcValue * dFunc;
                Matrix ddLogFunc =
                    ((orientationTimesDiffOuter + shapeOrientation) * factorValue * funcValue - orientationTimesDiffOuter * probDiff * factorValue * factorValue)
                    * (probDiff / (funcValue * funcValue));

                double x = resultForXCoord ? shapeX.Point : shapeY.Point;
                double d = resultForXCoord ? dLogFunc[0] : dLogFunc[1];
                double dd = resultForXCoord ? ddLogFunc[0, 0] : ddLogFunc[1, 1];
                return Gaussian.FromDerivatives(x, d, dd, forceProper: true);
            }
            else if (!shapeX.IsPointMass && !shapeY.IsPointMass)
            {
                VectorGaussian shapeLocationTimesFactor = ShapeLocationTimesFactor(point, shapeX, shapeY, shapeOrientation);
                double labelProbFalse = label.GetProbFalse();
                double shapeLocationWeight = labelProbFalse;
                double shapeLocationTimesFactorWeight =
                    Math.Exp(shapeLocationTimesFactor.GetLogNormalizer() - shapeX.GetLogNormalizer() - shapeY.GetLogNormalizer() - 0.5 * shapeOrientation.QuadraticForm(point)) *
                    (1 - 2 * labelProbFalse);

                var projectionOfSum = new Gaussian();
                projectionOfSum.SetToSum(
                    shapeLocationWeight,
                    resultForXCoord ? shapeX : shapeY,
                    shapeLocationTimesFactorWeight,
                    shapeLocationTimesFactor.GetMarginal(resultForXCoord ? 0 : 1));
                Gaussian result = new Gaussian();
                result.SetToRatio(projectionOfSum, resultForXCoord ? shapeX : shapeY);

                return result;
            }
            else
            {
                throw new NotSupportedException();
            }
        }

 public static double LogAverageFactor(
     Bernoulli label, Vector point, Gaussian shapeX, Gaussian shapeY, PositiveDefiniteMatrix shapeOrientation)
 {
     VectorGaussian shapeLocationTimesFactor = ShapeLocationTimesFactor(point, shapeX, shapeY, shapeOrientation);
     double labelProbFalse = label.GetProbFalse();
     double normalizerProduct = Math.Exp(
         shapeLocationTimesFactor.GetLogNormalizer() - 0.5 * shapeOrientation.QuadraticForm(point)
         - shapeX.GetLogNormalizer() - shapeY.GetLogNormalizer());
     double averageFactor = labelProbFalse + (1 - 2 * labelProbFalse) * normalizerProduct;
     Debug.Assert(averageFactor > 0);
     return Math.Log(averageFactor);
 }

 public static Bernoulli LabelAverageConditional(
     Vector point, double shapeX, double shapeY, PositiveDefiniteMatrix shapeOrientation)
 {
     Vector shapeLocation = Vector.FromArray(shapeX, shapeY);
     return new Bernoulli(Math.Exp(-0.5 * shapeOrientation.QuadraticForm(point - shapeLocation)));
 }