C# Gamma.GetLogNormalizer方法代码示例

		public static Gamma SampleAverageConditional(Gamma sample, double shape, [SkipIfUniform] Gamma rate, [Fresh] Gamma q)
		{
			if (sample.IsPointMass) return Gamma.Uniform();
			if (rate.IsPointMass) return GammaFromShapeAndRateOp.SampleAverageConditional(shape, rate.Point);
			if (q.IsPointMass) throw new Exception();
			double shape2 = GetShape2(sample, shape); // sample.Shape+shape-1
			double sampleMean, sampleVariance;
			if (sample.Rate == 0) sample = Gamma.FromShapeAndRate(sample.Shape, 1e-10);
			if (sample.Rate == 0) {
				// result.Shape > 1 iff sampleMean^2 > sampleVariance iff shape*(rate.Shape-2) > shape+rate.Shape-1
				// e.g. shape=4, rate.Shape=4
				sampleMean = shape2*rate.GetMeanInverse();
				sampleVariance = shape2*(1+shape2)*rate.GetMeanPower(-2) - sampleMean*sampleMean;
			} else if(true) {
				// Laplace for sample
				Gamma temp = Gamma.FromShapeAndRate(rate.Shape+2, rate.Rate);
				Gamma qy = Q_slow(temp, shape-1, sample);
				double x = qy.GetMean();
				double[] g = new double[] { x, 1, 0, 0 };
				GaussianOp_Laplace.LaplaceMoments(qy, g, dlogfs(x, shape-1, temp), out sampleMean, out sampleVariance);
			} else {
				// Laplace for rate
				// tends to be less accurate than above
				double x = q.GetMean();
				double x2 = x*x;
				double p = 1/(x + sample.Rate);
				double p2 = p*p;
				if (sample.Rate < x) {
					// another approach might be to write 1/(sample.Rate + r) = 1/r - sample.Rate/r/(sample.Rate + r)
					//double a1 = q.Shape - x2*p2;
					//double b1 = q.Rate + sample.Rate*p2;
					double logz = LogAverageFactor_slow(sample, shape, rate) + sample.GetLogNormalizer();
					Gamma sample1 = Gamma.FromShapeAndRate(sample.Shape+1, sample.Rate);
					double logz1 = LogAverageFactor_slow(sample1, shape, rate) + sample1.GetLogNormalizer();
					double pMean = Math.Exp(logz1 - logz);
					sampleMean = shape2*pMean;
					Gamma sample2 = Gamma.FromShapeAndRate(sample.Shape+2, sample.Rate);
					double logz2 = LogAverageFactor_slow(sample2, shape, rate) + sample2.GetLogNormalizer();
					double pMean2 = Math.Exp(logz2 - logz);
					double pVariance = pMean2 - pMean*pMean;
					sampleVariance = shape2*shape2*pVariance + shape2*(pVariance + pMean*pMean);
				} else {
					// sampleMean = E[ shape2/(sample.Rate + r) ]
					// sampleVariance = var(shape2/(sample.Rate + r)) + E[ shape2/(sample.Rate+r)^2 ]
					//                = shape2^2*var(1/(sample.Rate + r)) + shape2*(var(1/(sample.Rate+r)) + (sampleMean/shape2)^2)
					// Note: this is not a good approximation if sample.Rate is small
					double[] g = new double[] { p, -p2, 2*p2*p, -6*p2*p2 };
					double pMean, pVariance;
					GaussianOp_Laplace.LaplaceMoments(q, g, dlogfs(x, shape, sample), out pMean, out pVariance);
					sampleMean = shape2*pMean;
					sampleVariance = shape2*shape2*pVariance + shape2*(pVariance + pMean*pMean);
				}
			}
			Gamma sampleMarginal = Gamma.FromMeanAndVariance(sampleMean, sampleVariance);
			Gamma result = new Gamma();
			result.SetToRatio(sampleMarginal, sample, true);
			if (double.IsNaN(result.Shape) || double.IsNaN(result.Rate)) throw new Exception("result is nan");
			return result;
		}

		public static Gamma ExpAverageConditional(Gamma exp, Gaussian d, double x)
		{
			if (d.IsPointMass) return Gamma.PointMass(Math.Exp(d.Point));
			if (exp.IsPointMass) return Gamma.Uniform();
			double a = exp.Shape;
			double b = exp.Rate;
			double expx = Math.Exp(x);
			double ddlogf = -b*expx;
			double v = 1/(d.Precision - ddlogf);
			double logz = LogAverageFactor(exp, d, x) + exp.GetLogNormalizer();
			Gamma exp1 = Gamma.FromShapeAndRate(a+1,b);
			double x1 = x;
			for (int iter = 0; iter < 100; iter++) {
				double oldx1 = x1;
				x1 = X(exp1, d, x1);
				if (Math.Abs(x1 - oldx1) < 1e-10) break;
			}
			double logz1 = LogAverageFactor(exp1, d, x1) + exp1.GetLogNormalizer();
			Gamma exp2 = Gamma.FromShapeAndRate(a+2, b);
			double x2 = x;
			for (int iter = 0; iter < 100; iter++) {
				double oldx2 = x2;
				x2 = X(exp2, d, x2);
				if (Math.Abs(x2 - oldx2) < 1e-10) break;
			}
			double logz2 = LogAverageFactor(exp2, d, x2) + exp2.GetLogNormalizer();
			double mpost = Math.Exp(logz1 - logz);
			double m2post = Math.Exp(logz2 - logz);
			double vpost = m2post - mpost*mpost;
			double dddlogf = ddlogf;
			if (false) {
				// this is better for large b, worse for small b
				double mpost2 = expx + 0.5*v*expx + 0.5*v*v*expx*dddlogf;
				Console.WriteLine(mpost2);
			}
			if (false) {
				// this is identical to above
				double dlogf = (a-1) + ddlogf;
				double dlogz = dlogf + 0.5*v*v*dddlogf*d.Precision;
				double mpost2 = ((a-1)-dlogz)/b; // expx+0.5*v*v*prec*expx
				Console.WriteLine(mpost2);
			}
			Gamma result = Gamma.FromMeanAndVariance(mpost, vpost);
			result.SetToRatio(result, exp, true);
			return result;
		}

		public static double LogAverageFactor(Gamma sample, double shape, Gamma rate, [Fresh] Gamma q)
		{
			double x = q.GetMean();
			double shape2 = GetShape2(sample, shape);
			double logf = shape*Math.Log(x) - shape2*Math.Log(x + sample.Rate) + 
				MMath.GammaLn(shape2) - MMath.GammaLn(shape) - sample.GetLogNormalizer();
			double logz = logf + rate.GetLogProb(x) - q.GetLogProb(x);
			return logz;
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="exp">Incoming message from 'exp'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <param name="to_d">Previous outgoing message to 'd'.</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_(exp,d) p(exp,d) factor(exp,d))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Gamma exp, Gaussian d, Gaussian to_d)
		{
			if (d.IsPointMass) return LogAverageFactor(exp, d.Point);
			if (d.IsUniform()) return exp.GetLogAverageOf(new Gamma(0, 0));
			if (exp.IsPointMass) return LogAverageFactor(exp.Point, d);
			if (exp.IsUniform()) return 0.0;
			double[] nodes = new double[QuadratureNodeCount];
			double[] weights = new double[QuadratureNodeCount];
			double mD, vD;
			Gaussian dMarginal = d * to_d;
			dMarginal.GetMeanAndVariance(out mD, out vD);
			Quadrature.GaussianNodesAndWeights(mD, vD, nodes, weights);
			if (!to_d.IsUniform()) {
				// modify the weights to include q(y_k)/N(y_k;m,v)
				for (int i = 0; i < weights.Length; i++) {
					weights[i] *= Math.Exp(d.GetLogProb(nodes[i]) - Gaussian.GetLogProb(nodes[i], mD, vD));
				}
			}
			double Z = 0;
			for (int i = 0; i < weights.Length; i++) {
				double y = nodes[i];
				double f = weights[i] * Math.Exp((exp.Shape - 1) * y - exp.Rate * Math.Exp(y));
				Z += f;
			}
			return Math.Log(Z) - exp.GetLogNormalizer();
		}