C# Gamma.SetToProduct方法代码示例

//.........这里部分代码省略.........
		/// 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);
					result.SetToProduct(result, LastPrecisionMessage);
				}
			}
			// FIXME: this is not entirely safe since the caller could overwrite result.
			LastPrecisionMessage = result;
#endif
			return result;
		}