C# Gamma.SetToRatio方法代码示例

		public static Gamma SampleAverageConditional(Gamma sample, double shape, [SkipIfUniform] Gamma rate)
		{
			if (sample.IsPointMass || rate.Rate == 0) return Gamma.Uniform();
			if (rate.IsPointMass) return SampleAverageConditional(shape, rate.Point);
			double shape1 = AddShapesMinus1(shape, rate.Shape);
			double shape2 = AddShapesMinus1(shape, sample.Shape);
			double sampleMean, sampleVariance;
			if (sample.Rate == 0) sample = Gamma.FromShapeAndRate(sample.Shape, 1e-20);
			if (sample.Rate == 0) {
				sampleMean = shape2*rate.GetMeanInverse();
				sampleVariance = shape2*(1+shape2)*rate.GetMeanPower(-2) - sampleMean*sampleMean;
			} else if(true) {
				// quadrature over sample
				double y, ymin, ymax;
				GetIntegrationBoundsForSample(sample, shape, rate, out y, out ymin, out ymax);
				int n = QuadratureNodeCount;
				double inc = (ymax-ymin)/(n-1);
				shape1 = sample.Shape+shape-2;
				shape2 = shape+rate.Shape;
				double shift = shape1*Math.Log(y) - shape2*Math.Log(y+rate.Rate) - y*sample.Rate;
				MeanVarianceAccumulator mva = new MeanVarianceAccumulator();
				for (int i = 0; i < n; i++) {
					double x = ymin + i*inc;
					double logp = shape1*Math.Log(x) - shape2*Math.Log(x+rate.Rate) - x*sample.Rate;
					//if (i == 0 || i == n-1) Console.WriteLine(logp-shift);
					if ((i == 0 || i == n-1) && (logp-shift > -50)) throw new Exception("invalid integration bounds");
					double p = Math.Exp(logp - shift);
					mva.Add(x, p);
				}
				sampleMean = mva.Mean;
				sampleVariance = mva.Variance;
			} else {
				// quadrature over rate
				// 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)
				double r, rmin, rmax;
				GetIntegrationBoundsForRate(sample, shape, rate, out r, out rmin, out rmax);
				int n = QuadratureNodeCount;
				double inc = (rmax-rmin)/(n-1);
				double shift = shape1*Math.Log(r) - shape2*Math.Log(r+sample.Rate) - r*rate.Rate;
				MeanVarianceAccumulator mva = new MeanVarianceAccumulator();
				for (int i = 0; i < n; i++) {
					double x = rmin + i*inc;
					double logp = shape1*Math.Log(x) - shape2*Math.Log(x+sample.Rate) - x*rate.Rate;
					//if (i == 0 || i == n-1) Console.WriteLine(logp-shift);
					if ((i == 0 || i == n-1) && (logp-shift > -50)) throw new Exception("invalid integration bounds");
					double p = Math.Exp(logp - shift);
					double f = 1/(x + sample.Rate);
					mva.Add(f, p);
				}
				sampleMean = shape2*mva.Mean;
				sampleVariance = shape2*(1+shape2)*mva.Variance + shape2*mva.Mean*mva.Mean;
			}
			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 RateAverageConditional([SkipIfUniform] Gamma sample, double shape, Gamma rate)
		{
			if (rate.IsPointMass) return Gamma.Uniform();
			if (sample.IsPointMass) return RateAverageConditional(sample.Point, shape);
			if (sample.Rate == 0) return Gamma.Uniform();
			double shape1 = AddShapesMinus1(shape, rate.Shape);
			double shape2 = AddShapesMinus1(shape, sample.Shape);
			double rateMean, rateVariance;
			double r, rmin, rmax;
			GetIntegrationBoundsForRate(sample, shape, rate, out r, out rmin, out rmax);
			int n = QuadratureNodeCount;
			double inc = (rmax-rmin)/(n-1);
			double shift = shape1*Math.Log(r) - shape2*Math.Log(r+sample.Rate) - r*rate.Rate;
			MeanVarianceAccumulator mva = new MeanVarianceAccumulator();
			for (int i = 0; i < n; i++) {
				double x = rmin + i*inc;
				double logp = shape1*Math.Log(x) - shape2*Math.Log(x+sample.Rate) - x*rate.Rate;
				if ((i == 0 || i == n-1) && (logp-shift > -50)) throw new Exception("invalid integration bounds");
				double p = Math.Exp(logp - shift);
				mva.Add(x, p);
			}
			rateMean = mva.Mean;
			rateVariance = mva.Variance;
			Gamma rateMarginal = Gamma.FromMeanAndVariance(rateMean, rateVariance);
			Gamma result = new Gamma();
			result.SetToRatio(rateMarginal, rate, true);
			if (double.IsNaN(result.Shape) || double.IsNaN(result.Rate)) throw new Exception("result is nan");
			return result;
		}

		/// <summary>
		/// EP message to 'precision'
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="mean">Incoming message from 'mean'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="precision">Incoming message from 'precision'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>The outgoing EP message to the 'precision' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'precision' as the random arguments are varied.
		/// The formula is <c>proj[p(precision) sum_(sample,mean) p(sample,mean) factor(sample,mean,precision)]/p(precision)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="sample"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="mean"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="precision"/> is not a proper distribution</exception>
		public static Gamma PrecisionAverageConditional([SkipIfUniform] Gaussian sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision)
		{
			if (sample.IsPointMass && mean.IsPointMass)
				return PrecisionAverageConditional(sample.Point, mean.Point);

			Gamma result = new Gamma();
			if (precision.IsPointMass) {
				// The correct answer here is not uniform, but rather a limit.  
				// However it doesn't really matter what we return since multiplication by a point mass 
				// always yields a point mass.
				result.SetToUniform();
			} else if (sample.IsUniform() || mean.IsUniform()) {
				result.SetToUniform();
			} else if (!precision.IsProper()) {
				// improper prior
				throw new ImproperMessageException(precision);
			} else {
				// use quadrature to integrate over the precision
				// see LogAverageFactor
				double xm, xv, mm, mv;
				sample.GetMeanAndVarianceImproper(out xm, out xv);
				mean.GetMeanAndVarianceImproper(out mm, out mv);
				double upperBound = Double.PositiveInfinity;
				if (xv + mv < 0) upperBound = -1.0 / (xv + mv);
				double[] nodes = new double[QuadratureNodeCount];
				double[] weights = new double[nodes.Length];
				QuadratureNodesAndWeights(precision, nodes, weights);
				double Z = 0, rmean = 0, rvariance = 0;
				double shift = 0;
				for (int i = 0; i < nodes.Length; i++) {
					double v = 1.0 / nodes[i] + xv + mv;
					if (v < 0) continue;
					double lp = Gaussian.GetLogProb(xm, mm, v);
					if (shift == 0) shift = lp;
					double f = weights[i] * Math.Exp(lp - shift);
					double fm = f * nodes[i];
					double fmm = fm * nodes[i];
					Z += f;
					rmean += fm;
					rvariance += fmm;
				}

                // Adaptive Clenshaw-Curtis quadrature: gives same results on easy integrals but 
                // still fails ExpFactorTest2
                //Converter<double,double> func = delegate(double y) {
                //    double x = Math.Exp(y); 
                //    double v = 1.0 / x + xv + mv;
                //    if (v < 0) return 0.0;
                //    return Math.Exp(Gaussian.GetLogProb(xm, mm, v) + Gamma.GetLogProb(x, precision.Shape+1, precision.Rate));
                //};
                //double Z2 = BernoulliFromLogOddsOp.AdaptiveClenshawCurtisQuadrature(func, 1, 16, 1e-10);
                //Converter<double, double> func2 = delegate(double y)
                //{
                //    return Math.Exp(y) * func(y); 
                //};
                //double rmean2 = BernoulliFromLogOddsOp.AdaptiveClenshawCurtisQuadrature(func2, 1, 16, 1e-10);
                //Converter<double, double> func3 = delegate(double y)
                //{
                //    return Math.Exp(2*y) * func(y);
                //};
                //double rvariance2 = BernoulliFromLogOddsOp.AdaptiveClenshawCurtisQuadrature(func3, 1, 16, 1e-10);
                //rmean2 = rmean2/ Z2;
                //rvariance2 = rvariance2 / Z2 - rmean2 * rmean2;

				if (Z == 0.0) {
					//throw new Exception("Quadrature failed");
					//Console.WriteLine("Warning: Quadrature found zero mass.  Results may be inaccurate.");
					result.SetToUniform();
					return result;
				}
				double s = 1.0 / Z;
				rmean *= s;
				rvariance = rvariance * s - rmean * rmean;
				if (Double.IsInfinity(rmean)) {
					result.SetToUniform();
				} else {
					result.SetMeanAndVariance(rmean, rvariance);
					if (ForceProper) result.SetToRatioProper(result, precision);
					else result.SetToRatio(result, precision);
				}
			}
#if KeepLastMessage
			if (LastPrecisionMessage != null) {
				if (Stepsize != 1 && Stepsize != 0) {
					LastPrecisionMessage.SetToPower(LastPrecisionMessage, 1 - Stepsize);
					result.SetToPower(result, Stepsize);
//.........这里部分代码省略.........

		public static Gamma PrecisionAverageConditional([SkipIfUniform] Gaussian sample, [SkipIfUniform] Gaussian mean, [Proper] Gamma precision, [Fresh] Gamma q)
		{
			if (sample.IsUniform() || mean.IsUniform()) return Gamma.Uniform();
			double mx, vx;
			sample.GetMeanAndVariance(out mx, out vx);
			double mm, vm;
			mean.GetMeanAndVariance(out mm, out vm);
			double m = mx-mm;
			double v = vx+vm;
			double x = q.GetMean();
			double[] g = new double[] { x, 1, 0, 0 };
			double precMean, precVariance;
			LaplaceMoments(q, g, dlogfs(x, m, v), out precMean, out precVariance);
			if (precMean < 0) throw new Exception("internal: precMean < 0");
			Gamma precMarginal = Gamma.FromMeanAndVariance(precMean, precVariance);
			Gamma result = new Gamma();
			result.SetToRatio(precMarginal, precision, true);
			return result;
		}