C# Gamma.GetMeanLog方法代码示例

		public static Vector CalculateDerivativesNaive(Gamma q)
		{
			Vector gradElogGamma = Vector.Zero(2);
			// Get shape and scale of the distribution
			double a = q.Shape;
			double b = q.Rate;
			// Mean in the transformed domain
			double proposalMean = q.GetMeanLog();
			// Laplace approximation of variance in transformed domain 
			double proposalVariance = 1 / a;

			// Quadrature coefficient
			int n = 11;
			Vector nodes = Vector.Zero(n);
			Vector weights = Vector.Zero(n);
			Quadrature.GaussianNodesAndWeights(proposalMean, proposalVariance, nodes, weights);

			double EXDigamma = 0;
			double ELogGam = 0;
			double ELogXLogGamma = 0;

			// Calculate expectations in x=log(s) space using Gauss-Hermite quadrature
			double logZ = MMath.GammaLn(a) - a * Math.Log(b);
			for (int i = 0; i < n; i++) {
				double x = nodes[i];
				double expx = Math.Exp(x);
				double p = a * x - b * expx - logZ - Gaussian.GetLogProb(x, proposalMean, proposalVariance);
				p = Math.Exp(p) * weights[i];
				EXDigamma += p * (expx * MMath.Digamma(expx));
				ELogGam += p * (MMath.GammaLn(expx));
				ELogXLogGamma += p * (x * MMath.GammaLn(expx));
			}

			// Calculate derivatives
			gradElogGamma[0] = ELogXLogGamma - proposalMean * ELogGam;
			gradElogGamma[1] = -1.0 / b * EXDigamma;
			//gradElogGamma[1] = (ELogGamma(q) - ELogGamma(new Gamma(a + 1, b))) * a / b; 
			return gradElogGamma;
		}

		/// <summary>
		/// Evidence message for VMP
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <returns>Zero</returns>
		/// <remarks><para>
		/// In Variational Message Passing, the evidence contribution of a deterministic factor is zero.
		/// Adding up these values across all factors and variables gives the log-evidence estimate for VMP.
		/// </para></remarks>
		public static double AverageLogFactor(Gaussian log, Gamma d)
		{
			double m, v;
			log.GetMeanAndVariance(out m, out v);
			double Elogd=d.GetMeanLog();
			double Elogd2;
			if (!d.IsPointMass)
				Elogd2 = MMath.Trigamma(d.Shape) + Elogd * Elogd;
			else
				Elogd2 = Math.Log(d.Point) * Math.Log(d.Point);
			return -Elogd2/(2*v)+m*Elogd/v-m*m/(2*v)-MMath.LnSqrt2PI-.5*Math.Log(v);
		}

		/// <summary>
		/// Calculates \int G(x;a,b) LogGamma(x) dx
		/// </summary>
		/// <param name="q">G(x;a,b)</param>
		/// <returns>\int G(x;a,b) LogGamma(x) dx</returns>
		public static double ELogGamma(Gamma q)
		{
			if (q.IsPointMass)
				return MMath.GammaLn(q.Point);
			double a = q.Shape;
			double b = q.Rate;
			// Mean in the transformed domain
			double proposalMean = q.GetMeanLog();
			// Laplace approximation of variance in transformed domain 
			double proposalVariance = 1 / a;

			// Quadrature coefficient
			int n = 11;
			Vector nodes = Vector.Zero(n);
			Vector weights = Vector.Zero(n);
			Quadrature.GaussianNodesAndWeights(proposalMean, proposalVariance, nodes, weights);

			double ELogGamma = 0;
			double logZ = -a * Math.Log(b) + MMath.GammaLn(a);
			// Calculate expectations in x=log(s) space using Gauss-Hermite quadrature
			for (int i = 0; i < n; i++) {
				double x = nodes[i];
				double expx = Math.Exp(x);
				double p = a * x - b * expx - Gaussian.GetLogProb(x, proposalMean, proposalVariance) - logZ;
				p = Math.Exp(p + Math.Log(weights[i]));
				ELogGamma += p * (MMath.GammaLn(expx) + x);
			}

			// Add removed components
			return ELogGamma - proposalMean;
		}

		/// <summary>
		/// Perform the quadrature required for the VMP evidence message
		/// </summary>
		/// <param name="meanQ">Incoming message from m='mean'.</param>
		/// <param name="totalCountQ">Incoming message from s='totalCount'.</param>
		/// <returns>Vector of E[ LogGamma(s*m_k)].</returns>
		/// <remarks><para>
		/// The quadrature over 'totalCount' (which is Gamma-distributed) is 
		/// peformed by a change of variable x=log(s) followed by Gauss-Hermite 
		/// quadrature. The quadrature over m is performed using Gauss-Legendre. 
		/// </para></remarks>
		public static Vector EvidenceMessageExpectations(
				Dirichlet meanQ,
				Gamma totalCountQ)
		{
			// Get shape and scale of the distribution
			double at, bt;
			totalCountQ.GetShapeAndScale(out at, out bt);
			bt = 1 / bt; // want rate not scale

			// Mean in the transformed domain
			double proposalMean = totalCountQ.GetMeanLog();
			// Laplace approximation of variance in transformed domain 
			double proposalVariance = 1 / at;

			// Quadrature coefficient
			int nt = 32;
			Vector nodes = Vector.Zero(nt);
			Vector weights = Vector.Zero(nt);
			Vector expx = Vector.Zero(nt);
			if (!totalCountQ.IsPointMass) {
				Quadrature.GaussianNodesAndWeights(proposalMean, proposalVariance, nodes, weights);
				// Precompute weights for each m slice
				for (int i = 0; i < nt; i++) {
					double x = nodes[i];
					expx[i] = Math.Exp(x);
					double p = at * x - bt * expx[i] - Gaussian.GetLogProb(x, proposalMean, proposalVariance);
					weights[i] *= Math.Exp(p);
				}
			}

			int nm = 20;
			Vector mnodes = Vector.Zero(nm);
			Vector mweight = Vector.Zero(nm);
			Quadrature.UniformNodesAndWeights(0, 1, mnodes, mweight);
			int K = meanQ.Dimension;
			Vector[] mweights = new Vector[K];
			Beta[] mkDist = new Beta[K];
			double[] EELogGamma = new double[K];
			for (int i = 0; i < K; i++) {
				mweights[i] = Vector.Copy(mweight);
				mkDist[i] = new Beta(meanQ.PseudoCount[i], meanQ.TotalCount - meanQ.PseudoCount[i]);
				EELogGamma[i] = 0;
			}

			for (int j = 0; j < nm; j++) {
				double m = mnodes[j];
				double ELogGamma = 0;
				if (totalCountQ.IsPointMass)
					ELogGamma = MMath.GammaLn(m * totalCountQ.Point);
				else {
					// Calculate expectations in x=log(s) space using Gauss-Hermite quadrature
					for (int i = 0; i < nt; i++) {
						double x = nodes[i];
						ELogGamma += weights[i] * (MMath.GammaLn(m * expx[i]) + x);
					}
					// Normalise and add removed components
					double normalisation = Math.Pow(bt, at) / MMath.Gamma(at);
					ELogGamma = normalisation * ELogGamma - proposalMean;
				}
				for (int i = 0; i < K; i++) {
					mweights[i][j] *= Math.Exp(mkDist[i].GetLogProb(m));
					EELogGamma[i] += mweights[i][j] * ELogGamma;
				}
			}
			return Vector.FromArray(EELogGamma);
		}

		/// <summary>
		/// Perform the quadrature required for the Nonconjugate VMP message to 'totalCount'
		/// </summary>
		/// <param name="meanQPseudoCount">Incoming message from 'mean'.</param>
		/// <param name="totalCountQ">Incoming message from 'totalCount'.</param>
		/// <param name="EELogGamma">Array to be filled with E[LogGamma(s*m_k)].</param>
		/// <param name="EELogSLogGamma">Array to be filled with E[Log(s)*LogGamma(s*m_k)].</param>
		/// <param name="EEMSDigamma">Array to be filled with E[s*m_k*Digamma(s*m_k)].</param>
		/// <remarks><para>
		/// All three arrays are calculated simultaneously for efficiency. The quadrature over 
		/// 'totalCount' (which is Gamma-distributed) is peformed by a change of variable x=log(s)
		/// followed by Gauss-Hermite quadrature. The quadrature over m is performed using 
		/// Gauss-Legendre. 
		/// </para></remarks>
		public static void TotalCountMessageExpectations(
				Vector meanQPseudoCount,
				Gamma totalCountQ,
				out double[] EELogGamma,
				out double[] EELogSLogGamma,
				out double[] EEMSDigamma)
		{
			// Get shape and rate of the distribution
			double at = totalCountQ.Shape, bt = totalCountQ.Rate;

			// Mean in the transformed domain
			double proposalMean = totalCountQ.GetMeanLog();
			// Laplace approximation of variance in transformed domain 
			double proposalVariance = 1 / at;

			// Quadrature coefficient
			int nt = 32;
			Vector nodes = Vector.Zero(nt);
			Vector weights = Vector.Zero(nt);
			Vector expx = Vector.Zero(nt);
			if (!totalCountQ.IsPointMass) {
				Quadrature.GaussianNodesAndWeights(proposalMean, proposalVariance, nodes, weights);
				// Precompute weights for each m slice
				for (int i = 0; i < nt; i++) {
					double x = nodes[i];
					expx[i] = Math.Exp(x);
					double p = at * x - bt * expx[i] - Gaussian.GetLogProb(x, proposalMean, proposalVariance);
					weights[i] *= Math.Exp(p);
				}
			}

			int nm = 20;
			Vector mnodes = Vector.Zero(nm);
			Vector mweight = Vector.Zero(nm);
			Quadrature.UniformNodesAndWeights(0, 1, mnodes, mweight);
			int K = meanQPseudoCount.Count;
			Vector[] mweights = new Vector[K];
			Beta[] mkDist = new Beta[K];
			EELogGamma = new double[K];
			EELogSLogGamma = new double[K];
			EEMSDigamma = new double[K];
			double meanQTotalCount = meanQPseudoCount.Sum();
			for (int i = 0; i < K; i++) {
				mweights[i] = Vector.Copy(mweight);
				mkDist[i] = new Beta(meanQPseudoCount[i], meanQTotalCount - meanQPseudoCount[i]);
				EELogGamma[i] = 0;
				EELogSLogGamma[i] = 0;
				EEMSDigamma[i] = 0;
			}

			for (int j = 0; j < nm; j++) {
				double m = mnodes[j];
				double ESDigamma = 0;
				double ELogGamma = 0;
				double ELogSLogGamma = 0;
				if (totalCountQ.IsPointMass) {
					ESDigamma = totalCountQ.Point * MMath.Digamma(m * totalCountQ.Point);
					ELogGamma = MMath.GammaLn(m * totalCountQ.Point);
					ELogSLogGamma = Math.Log(totalCountQ.Point) * ELogGamma;
				} else {
					// Calculate expectations in x=log(s) space using Gauss-Hermite quadrature
					for (int i = 0; i < nt; i++) {
						double x = nodes[i];
						ELogGamma += weights[i] * (MMath.GammaLn(m * expx[i]) + x);
						ESDigamma += weights[i] * (expx[i] * MMath.Digamma(m * expx[i]) + 1);
						ELogSLogGamma += weights[i] * (x * MMath.GammaLn(m * expx[i]) + x * x + x * Math.Log(m));
					}
					// Normalise and add removed components
					double normalisation = Math.Pow(bt, at) / MMath.Gamma(at);
					ELogGamma = normalisation * ELogGamma - proposalMean;
					ELogSLogGamma = normalisation * ELogSLogGamma
                    - (MMath.Trigamma(at) + proposalMean * proposalMean + Math.Log(m) * proposalMean);
					ESDigamma = normalisation * ESDigamma - 1;
				}
				for (int i = 0; i < K; i++) {
					mweights[i][j] *= Math.Exp(mkDist[i].GetLogProb(m));
					EELogGamma[i] += mweights[i][j] * ELogGamma;
					EELogSLogGamma[i] += mweights[i][j] * ELogSLogGamma;
					EEMSDigamma[i] += mweights[i][j] * m * ESDigamma;
				}
			}
		}

		/// <summary>
		/// Perform the quadrature required for the Nonconjugate VMP message to 'mean'
		/// </summary>
		/// <param name="meanQPseudoCount">Incoming message from 'mean'.</param>
		/// <param name="totalCountQ">Incoming message from 'totalCount'.</param>
		/// <param name="EELogGamma">Array to be filled with E[LogGamma(s*m_k)].</param>
		/// <param name="EELogMLogGamma">Array to be filled with E[Log(m_k)*LogGamma(s*m_k)].</param>
		/// <param name="EELogOneMinusMLogGamma">Array to be filled with E[Log(1-m_k)*LogGamma(s*m_k)].</param>
		/// <remarks><para>
		/// All three arrays are calculated simultaneously for efficiency. The quadrature over 
		/// 'totalCount' (which is Gamma-distributed) is peformed by a change of variable x=log(s)
		/// followed by Gauss-Hermite quadrature. The quadrature over m is performed using 
		/// Gauss-Legendre. 
		/// </para></remarks>
		public static void MeanMessageExpectations(
				Vector meanQPseudoCount,
				Gamma totalCountQ,
				out double[] EELogGamma,
				out double[] EELogMLogGamma,
				out double[] EELogOneMinusMLogGamma)
		{
			// Get shape and scale of the distribution
			double at, bt;
			at = totalCountQ.Shape;
			bt = totalCountQ.Rate;

			// Mean in the transformed domain
			double ELogS = totalCountQ.GetMeanLog();
			// Laplace approximation of variance in transformed domain 
			double proposalVariance = 1 / at;

			// Quadrature coefficient
			int nt = 32;
			Vector nodes = Vector.Zero(nt);
			Vector weights = Vector.Zero(nt);
			Vector expx = Vector.Zero(nt);
			if (!totalCountQ.IsPointMass) {
				Quadrature.GaussianNodesAndWeights(ELogS, proposalVariance, nodes, weights);
				// Precompute weights for each m slice
				for (int i = 0; i < nt; i++) {
					double x = nodes[i];
					expx[i] = Math.Exp(x);
					double p = at * x - bt * expx[i] - Gaussian.GetLogProb(x, ELogS, proposalVariance);
					weights[i] *= Math.Exp(p);
				}
			}

			int nm = 20;
			Vector mnodes = Vector.Zero(nm);
			Vector mweight = Vector.Zero(nm);
			Quadrature.UniformNodesAndWeights(0, 1, mnodes, mweight);
			int K = meanQPseudoCount.Count;
			Vector[] mweights = new Vector[K];
			Beta[] mkDist = new Beta[K];
			EELogGamma = new double[K];
			EELogMLogGamma = new double[K];
			EELogOneMinusMLogGamma = new double[K];
			double meanQTotalCount = meanQPseudoCount.Sum();
			for (int i = 0; i < K; i++) {
				mweights[i] = Vector.Copy(mweight);
				mkDist[i] = new Beta(meanQPseudoCount[i], meanQTotalCount - meanQPseudoCount[i]);
				EELogGamma[i] = 0;
				EELogMLogGamma[i] = 0;
				EELogOneMinusMLogGamma[i] = 0;
			}

			double ES = totalCountQ.GetMean();
			double ESLogS = ELogS * ES + 1 / bt;

			for (int j = 0; j < nm; j++) {
				double m = mnodes[j];
				double ELogGamma = 0;
				if (totalCountQ.IsPointMass)
					ELogGamma = MMath.GammaLn(m * totalCountQ.Point);
				else {
					// Calculate expectations in x=log(s) space using Gauss-Hermite quadrature
					for (int i = 0; i < nt; i++)
						ELogGamma += weights[i] * (MMath.GammaLn(m * expx[i]) + nodes[i]);
					// Normalise and add removed components
					double normalisation = Math.Pow(bt, at) / MMath.Gamma(at);
					ELogGamma = normalisation * ELogGamma - ELogS;
				}

				double EELogMLogGammaTemp = Math.Log(m) * (ELogGamma + ELogS + Math.Log(m));
				double EELogOneMinusMLogGammaTemp = Math.Log(1 - m) *
                    (ELogGamma - (.5 * Math.Log(2 * Math.PI) - .5 * ELogS
                    - .5 * Math.Log(m) + m * ESLogS + ES * m * Math.Log(m) - ES * m));
				for (int i = 0; i < K; i++) {
					mweights[i][j] *= Math.Exp(mkDist[i].GetLogProb(m));
					EELogGamma[i] += mweights[i][j] * ELogGamma;
					EELogMLogGamma[i] += mweights[i][j] * EELogMLogGammaTemp;
					EELogOneMinusMLogGamma[i] += mweights[i][j] * EELogOneMinusMLogGammaTemp;
				}
			}
			for (int i = 0; i < K; i++)
				AddAnalyticComponent(
						mkDist[i],
						ELogS,
						ES,
						ESLogS,
//.........这里部分代码省略.........

		/// <summary>
		/// VMP message to 'totalCount'. This functionality is separated out to allow use by BetaOp. 
		/// </summary>
		/// <param name="meanPseudoCount">Pseudocount of incoming message from 'mean'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="totalCount">Incoming message from 'totalCount'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="meanLogProb">E[log(prob)] from incoming message from 'prob'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <remarks><para>
		/// The outgoing message here would not be Dirichlet distributed, so we use Nonconjugate VMP, which
		/// sends the approximate factor ensuring the gradient of the KL wrt to the variational parameters match. 
		/// </para></remarks>
		internal static Gamma TotalCountAverageLogarithmHelper(Vector meanPseudoCount, Gamma totalCount, Vector meanLogProb)
		{
			double[] EELogGamma;
			double[] EELogSLogGamma;
			double[] EEMSDigamma;
			// 2D quadrature
			TotalCountMessageExpectations(
					meanPseudoCount,
					totalCount,
					out EELogGamma,
					out EELogSLogGamma,
					out EEMSDigamma);
			double at = totalCount.Shape;
			double bt = totalCount.Rate;
			// Find required expectations using quadrature
			Vector gradElogGamma = GammaFromShapeAndRateOp.CalculateDerivatives(totalCount);
			Vector gradS = gradElogGamma;
			Vector EM = Vector.Zero(meanPseudoCount.Count);
			EM.SetToProduct(meanPseudoCount, 1.0 / meanPseudoCount.Sum());
			double c = 0;
			for (int k = 0; k < meanPseudoCount.Count; k++) {
				gradS[0] -= EELogSLogGamma[k] - totalCount.GetMeanLog() * EELogGamma[k];
				gradS[1] -= -EEMSDigamma[k] / bt;
				c += EM[k] * meanLogProb[k];
			}
			// Analytic 
			gradS[0] += c / bt;
			gradS[1] -= c * at / (bt * bt);
			Matrix mat = new Matrix(2, 2);
			mat[0, 0] = MMath.Trigamma(at);
			mat[1, 0] = mat[0, 1] = -1 / bt;
			mat[1, 1] = at / (bt * bt);
			Vector v = GammaFromShapeAndRateOp.twoByTwoInverse(mat) * gradS;
			return Gamma.FromShapeAndRate(v[0] + 1, v[1]);
		}