/// <summary>
		/// Gibbs message to 'precision'
		/// </summary>
		/// <param name="sample">Constant value for 'sample'.</param>
		/// <param name="mean">Constant value for 'mean'.</param>
		/// <returns>The outgoing Gibbs message to the 'precision' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'precision' conditioned on the given values.
		/// </para></remarks>
		public static Gamma PrecisionConditional(double sample, double mean)
		{
			Gamma result = new Gamma();
			double diff = sample - mean;
			result.Rate = 0.5 * diff * diff;
			result.Shape = 1.5;
			return result;
		}

        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);
        }

		/// <summary>
		/// VMP message to 'totalCount'
		/// </summary>
		/// <param name="mean">Incoming message from 'mean'. Must be a proper distribution.  If 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. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="prob">Incoming message from 'prob'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="to_totalCount">Previous outgoing message to 'TotalCount'.</param>
		/// <returns>The outgoing VMP message to the 'totalCount' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'totalCount'.
		/// The formula is <c>exp(sum_(mean,prob) p(mean,prob) log(factor(prob,mean,totalCount)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="mean"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="totalCount"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="prob"/> is not a proper distribution</exception>
		public static Gamma TotalCountAverageLogarithm([Proper] Beta mean, [Proper] Gamma totalCount, [SkipIfUniform] Beta prob, Gamma to_totalCount)
		{
			double ELogP, ELogOneMinusP;
			prob.GetMeanLogs(out ELogP, out ELogOneMinusP);
			Gamma approximateFactor = DirichletOp.TotalCountAverageLogarithmHelper(
				Vector.FromArray(new double[] { mean.TrueCount, mean.FalseCount }),
				totalCount,
				Vector.FromArray(new double[] { ELogP, ELogOneMinusP }));
			if (damping == 0.0)
				return approximateFactor;
			else
				return (approximateFactor^(1-damping)) * (to_totalCount ^ damping);
		}

        private void Gamma_Click(object sender, EventArgs e)
        {
            Gamma dlg = new Gamma();
            dlg.rouge = dlg.vert = dlg.bleu = 1;

            if (DialogResult.OK == dlg.ShowDialog())
            {
                if (filtres.Gamma(m_Bitmap, dlg.rouge, dlg.vert, dlg.bleu))
                    this.Invalidate();
                    this.Refresh();
            }
        }

		//-- 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));
		}

		/// <summary>
		/// VMP message to 'd'
		/// </summary>
		/// <param name="log">Incoming message from 'log'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="d">Incoming message from 'd'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="to_D">Previous outgoing message to 'D'.</param>
		/// <returns>The outgoing VMP message to the 'd' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'd' with 'log' integrated out.
		/// The formula is <c>sum_log p(log) factor(log,d)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="d"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="log"/> is not a proper distribution</exception>
		public static Gamma DAverageLogarithm([SkipIfUniform] Gaussian log, [SkipIfUniform] Gamma d, Gamma to_D)
		{
			if (log.IsPointMass)
				return DAverageLogarithm(log.Point, d, to_D);
			Vector grad = Vector.Zero(2);
			double meanLog = d.GetMeanLog();
			double m,v;
			log.GetMeanAndVariance(out m, out v);
			grad[0] = -MMath.Tetragamma(d.Shape) / (2 * v) - MMath.Trigamma(d.Shape) / v * (meanLog-m);
			grad[1] = (meanLog - m) / (v * d.Rate);
			Gamma approximateFactor = GammaFromShapeAndRateOp.NonconjugateProjection(d, grad);
			if (damping == 0.0)
				return approximateFactor;
			else
				return (approximateFactor ^ (1 - damping)) * (to_D ^ damping);
		}

		/// <summary>
		/// EP message to 'log'
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'log' as the random arguments are varied.
		/// The formula is <c>proj[p(log) sum_(d) p(d) factor(log,d)]/p(log)</c>.
		/// </para></remarks>
		public static Gaussian LogAverageConditional(Gaussian log, Gamma d, Gaussian result)
		{
			var g = Gamma.FromShapeAndRate(d.Shape + 1, d.Rate);
			return ExpOp.DAverageConditional(g, log, result);
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="d">Constant value for 'd'.</param>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <returns>Logarithm of the factor's contribution the EP model evidence</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(log) p(log) factor(log,d) / sum_log p(log) messageTo(log))</c>.
		/// Adding up these values across all factors and variables gives the log-evidence estimate for EP.
		/// </para></remarks>
		public static double LogEvidenceRatio(double log, Gamma d)
		{
			return LogAverageFactor(log, d);
		}

 public void ValidateMean([Values(0.0, 1.0, 1.0, 10.0, 10.0, 10.0)] double shape, [Values(0.0, 0.1, 1.0, 10.0, 1.0, Double.PositiveInfinity)] double invScale, [Values(Double.NaN, 10.0, 1.0, 1.0, 10.0, 10.0)] double mean)
 {
     var n = new Gamma(shape, invScale);
     Assert.AreEqual(mean, n.Mean);
 }

 public void ValidateMaximum()
 {
     var n = new Gamma(1.0, 1.0);
     Assert.AreEqual(Double.PositiveInfinity, n.Maximum);
 }

 public void ValidateEntropy([Values(0.0, 1.0, 1.0, 10.0, 10.0, 10.0)] double shape, [Values(0.0, 0.1, 1.0, 10.0, 1.0, Double.PositiveInfinity)] double invScale, [Values(Double.NaN, 3.3025850929940456285068402234265387271634735938763824, 1.0, 0.23346908548693395836262094490967812177376750477943892, 2.5360541784809796423806123995940423293748689934081866, 0.0)] double entropy)
 {
     var n = new Gamma(shape, invScale);
     AssertHelpers.AlmostEqual(entropy, n.Entropy, 13);
 }

 public void ValidateDensityLn(
     [Values(0, 0, 0, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10)] int shape, 
     [Values(0.0, 0.0, 0.0, 0.1, 0.1, 0.1, 1.0, 1.0, 1.0, 10.0, 10.0, 10.0, 1.0, 1.0, 1.0, Double.PositiveInfinity, Double.PositiveInfinity, Double.PositiveInfinity)] double invScale, 
     [Values(0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0)] double x, 
     [Values(Double.NegativeInfinity, Double.NegativeInfinity, Double.NegativeInfinity, -2.3025850929940456285068402234265387271634735938763824, -2.402585092994045634057955346552321429281631934330484, -3.3025850929940456285068402234265387271634735938763824, 0.0, -1.0, -10.0, Double.NegativeInfinity, 0.22402344985898722897219667227693591172986563062456522, -69.052710713194601614865880235563786219860220971716511, Double.NegativeInfinity, -13.801827480081469611207717874566706164281149255663166, -2.0785616431350584550457947824074282958712358580042068, Double.NegativeInfinity, Double.NegativeInfinity, Double.PositiveInfinity)] double pdfln)
 {
     var n = new Gamma(shape, invScale);
     AssertHelpers.AlmostEqual(pdfln, n.DensityLn(x), 14);
 }

 public void ValidateDensity(
     [Values(0, 0, 0, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10)] int shape, 
     [Values(0.0, 0.0, 0.0, 0.1, 0.1, 0.1, 1.0, 1.0, 1.0, 10.0, 10.0, 10.0, 1.0, 1.0, 1.0, Double.PositiveInfinity, Double.PositiveInfinity, Double.PositiveInfinity)] double invScale, 
     [Values(0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0)] double x, 
     [Values(0.0, 0.0, 0.0, 0.10000000000000000555111512312578270211815834045410156, 0.090483741803595961836995913651194571475319347018875963, 0.036787944117144234201693506390001264039984687455876246, 1.0, 0.36787944117144232159552377016146086744581113103176804, 0.000045399929762484851535591515560550610237918088866564953, 0.0, 1.2511003572113329898476497894772544708420990097708588, 1.0251532120868705806216092933926141802686541811003037e-30, 0.0, 0.0000010137771196302974029859010421116095333052555418644397, 0.12511003572113329898476497894772544708420990097708601, 0.0, 0.0, Double.PositiveInfinity)] double pdf)
 {
     var n = new Gamma(shape, invScale);
     AssertHelpers.AlmostEqual(pdf, n.Density(x), 14);
 }

 public void ValidateCumulativeDistribution(
     [Values(0, 0, 0, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10)] int shape, 
     [Values(0.0, 0.0, 0.0, 0.1, 0.1, 0.1, 1.0, 1.0, 1.0, 10.0, 10.0, 10.0, 1.0, 1.0, 1.0, Double.PositiveInfinity, Double.PositiveInfinity, Double.PositiveInfinity)] double invScale, 
     [Values(0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0, 0.0, 1.0, 10.0)] double x, 
     [Values(0.0, 0.0, 0.0, 0.0, 0.095162581964040431858607615783064404690935346242622848, 0.63212055882855767840447622983853913255418886896823196, 0.0, 0.63212055882855767840447622983853913255418886896823196, 0.99995460007023751514846440848443944938976208191113396, 0.0, 0.54207028552814779168583514294066541824736464003242184, 0.99999999999999999999999999999988746526039157266114706, 0.0, 0.00000011142547833872067735305068724025236288094949815466035, 0.54207028552814779168583514294066541824736464003242184, 0.0, 0.0, 1.0)] double cdf)
 {
     var n = new Gamma(shape, invScale);
     AssertHelpers.AlmostEqual(cdf, n.CumulativeDistribution(x), 14);
 }

 public void SetShapeFailsWithNegativeShape()
 {
     var n = new Gamma(1.0, 1.0);
     Assert.Throws<ArgumentOutOfRangeException>(() => n.Shape = -1.0);
 }

 public void CanCreateGamma([Values(0.0, 1.0, 1.0, 10.0, 10.0, 10.0)] double shape, [Values(0.0, 0.1, 1.0, 10.0, 1.0, Double.PositiveInfinity)] double invScale)
 {
     var n = new Gamma(shape, invScale);
     Assert.AreEqual(shape, n.Shape);
     Assert.AreEqual(invScale, n.InvScale);
 }

 public void ValidateMedian()
 {
     var n = new Gamma(0.0, 0.0);
     Assert.Throws<NotSupportedException>(() => { var median = n.Median; });
 }

		/// <summary>
		/// EP message to 'd'
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <param name="to_log">Previous outgoing message to 'log'.</param>
		/// <returns>The outgoing EP message to the 'd' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'd' as the random arguments are varied.
		/// The formula is <c>proj[p(d) sum_(log) p(log) factor(log,d)]/p(d)</c>.
		/// </para></remarks>
		public static Gamma DAverageConditional(Gaussian log, Gamma d, Gaussian to_log)
		{
			var g = Gamma.FromShapeAndRate(d.Shape + 1, d.Rate);
			return ExpOp.ExpAverageConditional(g, log, to_log);
		}

 public void ValidateMinimum()
 {
     var n = new Gamma(1.0, 1.0);
     Assert.AreEqual(0.0, n.Minimum);
 }

 public void ValidateMode([Values(0.0, 1.0, 1.0, 10.0, 10.0, 10.0)] double shape, [Values(0.0, 0.1, 1.0, 10.0, 1.0, Double.PositiveInfinity)] double invScale, [Values(Double.NaN, 0.0, 0.0, 0.9, 9.0, 10.0)] double mode)
 {
     var n = new Gamma(shape, invScale);
     Assert.AreEqual(mode, n.Mode);
 }