C# Bernoulli.GetProbTrue方法代码示例

 public static Bernoulli AAverageConditional(Bernoulli b, double penalty)
 {
     double penaltyFactor = Math.Exp(-penalty);
     double bProbTrue = b.GetProbTrue();
     double probTrue = (penaltyFactor * (1 - bProbTrue) + bProbTrue) / (penaltyFactor + 1);
     return new Bernoulli(probTrue);
 }

		/// <summary>
		/// Evidence message for VMP.
		/// </summary>
		/// <param name="sample">Incoming message from sample</param>
		/// <param name="logOdds">Incoming message from logOdds</param>
		/// <returns><c>sum_x marginal(x)*log(factor(x))</c></returns>
		/// <remarks><para>
		/// The formula for the result is <c>int log(f(x)) q(x) dx</c>
		/// where <c>x = (sample,logOdds)</c>.
		/// </para></remarks>
		public static double AverageLogFactor(Bernoulli sample, [Proper, SkipIfUniform] Gaussian logOdds)
		{
			if (logOdds.IsUniform()) return 0.0;
			double m, v;
			logOdds.GetMeanAndVariance(out m, out v);
			double t = Math.Sqrt(m * m + v);
			double s = 2 * sample.GetProbTrue() - 1;  // probTrue - probFalse
			return MMath.LogisticLn(t) + (s * m - t) / 2;
		}

		/// <summary>
		/// VMP message to LogOdds
		/// </summary>
		/// <param name="sample">Incoming message from sample</param>
		/// <param name="logOdds">Incoming message from logOdds</param>
		/// <returns><c>sum_x marginal(x)*log(factor(x))</c></returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the integral of the log-factor times incoming messages, over all arguments except 'logOdds'.
		/// The formula is <c>int log(f(logOdds,x)) q(x) dx</c> where <c>x = (sample)</c>.
		/// </para></remarks>
		public static Gaussian LogOddsAverageLogarithm(Bernoulli sample, [Proper, SkipIfUniform] Gaussian logOdds)
		{
			if (logOdds.IsUniform()) return logOdds;
			double m, v;
			logOdds.GetMeanAndVariance(out m, out v);
			double t = Math.Sqrt(m * m + v);
			double lambda = (t == 0) ? 0.25 : Math.Tanh(t / 2) / (2 * t);
			return Gaussian.FromNatural(sample.GetProbTrue() - 0.5, lambda);
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="isBetween">Incoming message from 'isBetween'.</param>
		/// <param name="X">Incoming message from 'x'.</param>
		/// <param name="lowerBound">Incoming message from 'lowerBound'.</param>
		/// <param name="upperBound">Incoming message from 'upperBound'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(isBetween,x,lowerBound,upperBound) p(isBetween,x,lowerBound,upperBound) factor(isBetween,x,lowerBound,upperBound))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Bernoulli isBetween, Gaussian X, Gaussian lowerBound, Gaussian upperBound)
		{
			if (isBetween.LogOdds == 0.0) return -MMath.Ln2;
			else
			{
#if true
				double logitProbBetween = MMath.LogitFromLog(LogProbBetween(X, lowerBound, upperBound));
				return Bernoulli.LogProbEqual(isBetween.LogOdds, logitProbBetween);
#else
			double d_p = isBetween.GetProbTrue() - isBetween.GetProbFalse();
			return Math.Log(d_p * Math.Exp(LogProbBetween()) + isBetween.GetProbFalse());
#endif
			}
		}

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

        /// <summary>
        /// VMP message to 'sample'.
        /// </summary>
        /// <param name="choice">Incoming message from 'choice'.</param>
        /// <param name="probTrue0">Constant value for 'probTrue0'.</param>
        /// <param name="probTrue1">Constant value for 'probTrue1'.</param>
        /// <returns>The outgoing VMP message to the 'sample' argument.</returns>
        /// <remarks><para>
        /// The outgoing message is the exponential of the integral of the log-factor times incoming messages, over all arguments except 'sample'.
        /// The formula is <c>int log(f(sample,x)) q(x) dx</c> where <c>x = (choice,probTrue0,probTrue1)</c>.
        /// </para></remarks>
        public static Bernoulli SampleAverageLogarithm(Bernoulli choice, double probTrue0, double probTrue1)
		{
            Bernoulli result = new Bernoulli();
			if(choice.IsPointMass) return SampleConditional(choice.Point,probTrue0,probTrue1);
			// log(p(X=true)/p(X=false)) = sum_k p(Y=k) log(ProbTrue[k]/(1-ProbTrue[k]))
			result.LogOdds = choice.GetProbFalse() * MMath.Logit(probTrue0) + choice.GetProbTrue() * MMath.Logit(probTrue1);
			return result;
		}

        /// <summary>
        /// VMP message to 'choice'.
        /// </summary>
        /// <param name="sample">Incoming message from 'sample'.</param>
        /// <param name="probTrue0">Constant value for 'probTrue0'.</param>
        /// <param name="probTrue1">Constant value for 'probTrue1'.</param>
        /// <returns>The outgoing VMP message to the 'choice' argument.</returns>
        /// <remarks><para>
        /// The outgoing message is the exponential of the integral of the log-factor times incoming messages, over all arguments except 'choice'.
        /// The formula is <c>int log(f(choice,x)) q(x) dx</c> where <c>x = (sample,probTrue0,probTrue1)</c>.
        /// </para></remarks>
		public static Bernoulli ChoiceAverageLogarithm(Bernoulli sample, double probTrue0, double probTrue1)
		{
            Bernoulli result = new Bernoulli();
			if(sample.IsPointMass) return ChoiceConditional(sample.Point,probTrue0,probTrue1);
			// p(Y=k) =propto ProbTrue[k]^p(X=true) (1-ProbTrue[k])^p(X=false)
			// log(p(Y=true)/p(Y=false)) = p(X=true)*log(ProbTrue[1]/ProbTrue[0]) + p(X=false)*log((1-ProbTrue[1])/(1-ProbTrue[0]))
			//                           = p(X=false)*(log(ProbTrue[0]/(1-ProbTrue[0]) - log(ProbTrue[1]/(1-ProbTrue[1]))) + log(ProbTrue[1]/ProbTrue[0])
			if (probTrue0 == 0 || probTrue1 == 0) throw new ArgumentException("probTrue is zero");
			result.LogOdds = sample.GetProbTrue() * Math.Log(probTrue1 / probTrue0) + sample.GetProbFalse() * Math.Log((1 - probTrue1) / (1 - probTrue0));
			return result;
		}

        /// <summary>
        /// EP message to 'sample'.
        /// </summary>
        /// <param name="choice">Incoming message from 'choice'.</param>
        /// <param name="probTrue0">Constant value for 'probTrue0'.</param>
        /// <param name="probTrue1">Constant value for 'probTrue1'.</param>
        /// <returns>The outgoing EP message to the 'sample' argument.</returns>
        /// <remarks><para>
        /// The outgoing message is the integral of the factor times incoming messages, over all arguments except 'sample'.
        /// The formula is <c>int f(sample,x) q(x) dx</c> where <c>x = (choice,probTrue0,probTrue1)</c>.
        /// </para></remarks>
        public static Bernoulli SampleAverageConditional(Bernoulli choice, double probTrue0, double probTrue1)
		{
            Bernoulli result = new Bernoulli();
			if(choice.IsPointMass) return SampleConditional(choice.Point,probTrue0,probTrue1);
#if FAST
			result.SetProbTrue(choice.GetProbFalse() * probTrue0 + choice.GetProbTrue() * probTrue1);
#else
			// This method is more numerically stable but slower.
			// let oX = log(p(X)/(1-p(X))
			// let oY = log(p(Y)/(1-p(Y))
			// oX = log( (TT*sigma(oY) + TF*sigma(-oY))/(FT*sigma(oY) + FF*sigma(-oY)) )
			//    = log( (TT*exp(oY) + TF)/(FT*exp(oY) + FF) )
			//    = log( (exp(oY) + TF/TT)/(exp(oY) + FF/FT) ) + log(TT/FT)
			// ay = log(TF/TT)
			// by = log(FF/FT)
			// offset = log(TT/FT)
			if (probTrue0 == 0 || probTrue1 == 0) throw new ArgumentException("probTrue is zero");
			double ay = Math.Log(probTrue0 / probTrue1);
			double by = Math.Log((1 - probTrue0) / (1 - probTrue1));
			double offset = MMath.Logit(probTrue1);
			result.LogOdds = MMath.DiffLogSumExp(choice.LogOdds, ay, by) + offset;
#endif
			return result;
		}

		/// <summary>
		/// 
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'.</param>
		/// <param name="index">Incoming message from 'index'.</param>
		/// <param name="ProbTrue">Constant value for 'probTrue'.</param>
		/// <returns></returns>
		/// <remarks><para>
		/// 
		/// </para></remarks>
		public static double AverageValueLn(Bernoulli sample, Discrete index, double[] ProbTrue)
		{
			double p = 0;
			for (int i = 0; i < ProbTrue.Length; i++)
			{
				p += ProbTrue[i] * index[i];
			}
			double b = sample.GetProbTrue();
			return Math.Log(b * p + (1 - b) * (1 - p));
		}

		/// <summary>
		/// VMP message to 'index'.
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'.</param>
		/// <param name="ProbTrue">Constant value for 'probTrue'.</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 integral of the log-factor times incoming messages, over all arguments except 'index'.
		/// The formula is <c>int log(f(index,x)) q(x) dx</c> where <c>x = (sample,probTrue)</c>.
		/// </para></remarks>
		public static Discrete IndexAverageLogarithm(Bernoulli sample, double[] ProbTrue, Discrete result)
		{
			if (result == default(Discrete)) result = Discrete.Uniform(ProbTrue.Length);
			// E[sum_k I(Y=k) (X*log(ProbTrue[k]) + (1-X)*log(1-ProbTrue[k]))]
			// = sum_k I(Y=k) (p(X=true)*log(ProbTrue[k]) + p(X=false)*log(1-ProbTrue[k]))
			// p(Y=k) =propto ProbTrue[k]^p(X=true) (1-ProbTrue[k])^p(X=false)
			Vector probs = result.GetWorkspace();
			double p = sample.GetProbTrue();
			probs.SetTo(ProbTrue);
			probs.SetToFunction(probs, x => Math.Pow(x, p) * Math.Pow(1.0 - x, 1.0 - p));
			result.SetProbs(probs);
			return result;
		}

		//-- VMP -------------------------------------------------------------------------------------------

		/// <summary>
		/// Evidence message for VMP
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="probTrue">Incoming message from 'probTrue'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>Average of the factor's log-value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>sum_(sample,probTrue) p(sample,probTrue) log(factor(sample,probTrue))</c>.
		/// Adding up these values across all factors and variables gives the log-evidence estimate for VMP.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="sample"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="probTrue"/> is not a proper distribution</exception>
		public static double AverageLogFactor(Bernoulli sample, [Proper] Beta probTrue)
		{
			if (sample.IsPointMass) return AverageLogFactor(sample.Point, probTrue);
			double eLogP, eLog1MinusP;
			probTrue.GetMeanLogs(out eLogP, out eLog1MinusP);
			double p = sample.GetProbTrue();
			return p * eLogP + (1 - p) * eLog1MinusP;
		}

		/// <summary>
		/// VMP message to 'a'.
		/// </summary>
		/// <param name="and">Incoming message from 'and'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="B">Incoming message from 'b'.</param>
		/// <returns>The outgoing VMP message to the 'a' argument.</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the integral of the log-factor times incoming messages, over all arguments except 'a'.
		/// The formula is <c>int log(f(a,x)) q(x) dx</c> where <c>x = (and,b)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="and"/> is not a proper distribution</exception>
		public static Bernoulli AAverageLogarithm([SkipIfUniform] Bernoulli and, Bernoulli B)
		{
			// when 'and' is marginalized, the factor is proportional to exp(A*B*and.LogOdds)
			return Bernoulli.FromLogOdds(and.LogOdds * B.GetProbTrue());
		}

		/// <summary>
		/// VMP message to 'probTrue'
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>The outgoing VMP message to the 'probTrue' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'probTrue'.
		/// The formula is <c>exp(sum_(sample) p(sample) log(factor(sample,probTrue)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="sample"/> is not a proper distribution</exception>
		public static Beta ProbTrueAverageLogarithm(Bernoulli sample)
		{
			// E[x*log(p) + (1-x)*log(1-p)] = E[x]*log(p) + (1-E[x])*log(1-p)
			double ex = sample.GetProbTrue();
			return new Beta(1 + ex, 2 - ex);
		}

        public static Wishart ShapeOrientationAverageConditional(
            Vector point, Bernoulli label, Gaussian shapeX, Gaussian shapeY, Wishart shapeOrientation, Wishart result)
        {
            if (shapeOrientation.IsPointMass && 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 = shapeLocation - point;
                Matrix diffOuter = diff.Outer(diff);
                Matrix orientationTimesDiffOuter = shapeOrientation.Point * diffOuter;
                double trace = orientationTimesDiffOuter.Trace();

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

                PositiveDefiniteMatrix dLogFunc = new PositiveDefiniteMatrix(diffOuter * (-0.5 * probDiff * factorValue / funcValue));
                double xxddLogFunc =
                    -0.5 * probDiff * (-0.5 * labelProbFalse * factorValue * trace * trace / (funcValue * funcValue) + factorValue * trace / funcValue);

                LowerTriangularMatrix cholesky = new LowerTriangularMatrix(2, 2);
                cholesky.SetToCholesky(shapeOrientation.Point);
                PositiveDefiniteMatrix inverse = shapeOrientation.Point.Inverse();
                result.SetDerivatives(cholesky, inverse, dLogFunc, xxddLogFunc, forceProper: true);
                return result;
            }
            else
            {
                throw new NotSupportedException();
            }
        }