C# Gamma.GetLogProb方法代码示例

        private static void ARSTest(int sampleCount = 100000)
        {
            List<double> acceptedSamples = new List<double>();

            double shape = 2;
            double scale = 2;
            Gamma gaussian = new Gamma(shape, scale);
            DiscreteEnvelope envelope = new DiscreteEnvelope(0, 1000, gaussian, new double[] { 1, 5 });

            int rejectedCount = 0;
            while  (acceptedSamples.Count < sampleCount)
            {
                var sampleRegion = envelope.SampleDiscrete();
                double sample = envelope.SampleContinuous(sampleRegion);

                double ratio = Math.Exp(gaussian.GetLogProb(sample) - envelope.GetLogProb(sampleRegion, sample));
                double u = Utils.SRandom.GetUniform();

                if (u < ratio)
                {
                    Console.WriteLine("Sample accepted {0}/{1} : {2}", acceptedSamples.Count + 1, sampleCount, sample);
                    acceptedSamples.Add(sample);
                }
                else
                {
                    Console.WriteLine("Rejected, adding cut at {0}", sample);
                    rejectedCount++;
                    envelope.AddCutPoint(sample);
                }
            }
            double mean = acceptedSamples.Sum() / acceptedSamples.Count;
            double variance = acceptedSamples.Select(s => (s - mean) * (s - mean)).Sum() / (sampleCount - 1);

            Console.WriteLine("Total Rejected = {0}", rejectedCount);
            Console.WriteLine("Sample Mean = {0}, Sample Variance = {1}", mean, variance);
            Console.WriteLine("True Mean = {0},     True Variance = {1}", shape * scale, shape * scale * scale);
        }

		//-- EP -------------------------------------------------------------------------------------------

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="log">Constant value for 'log'.</param>
		/// <param name="d">Incoming message from '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_(d) p(d) factor(log,d))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(double log, Gamma d)
		{
			return d.GetLogProb(Math.Exp(log));
		}

		public static Vector CalculateDerivativesTrapezoid(Gamma q)
		{
			Vector gradElogGamma = Vector.Zero(2);
			// Get shape and scale of the distribution
			double a = q.Shape;
			double b = q.Rate;
			double mean, variance;
			q.GetMeanAndVariance(out mean, out variance);
			double upperBound = 10;

			int n = 10000;
			double ELogGamma = 0, ELogXLogGamma = 0, ExDigamma = 0;
			double inc = upperBound/n;
			for (int i = 0; i < n; i++) {
				double x = inc * (i+1);
				double logp = q.GetLogProb(x);
				double p = Math.Exp(logp);
				double f = p * MMath.GammaLn(x);
				ELogGamma += f;
				ELogXLogGamma += Math.Log(x)*f;
				ExDigamma += x*MMath.Digamma(x)*p;
			}
			ELogGamma *= inc;
			ELogXLogGamma *= inc;
			ExDigamma *= inc;
			gradElogGamma[0] = ELogXLogGamma + (Math.Log(b) - MMath.Digamma(a))*ELogGamma;
			gradElogGamma[1] = -ExDigamma/b;
			return gradElogGamma;
		}

		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">Constant value for '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) p(exp) factor(exp,d))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Gamma exp, double d)
		{
			return exp.GetLogProb(Math.Exp(d));
		}

		public static double LogAverageFactor([SkipIfUniform] Gaussian sample, [SkipIfUniform] Gaussian mean, Gamma precision, [Fresh] Gamma q)
		{
			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 m2 = m*m;
			double x = q.GetMean();
			double logf = -MMath.LnSqrt2PI -0.5*Math.Log(v + 1/x) - 0.5*m2/(v + 1/x);
			double logz = logf + precision.GetLogProb(x) - q.GetLogProb(x);
			return logz;
		}