C# Gamma.IsUniform方法代码示例

		/// <summary>
		/// EP message to 'exp'
		/// </summary>
		/// <param name="exp">Incoming message from 'exp'.</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 EP message to the 'exp' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'exp' as the random arguments are varied.
		/// The formula is <c>proj[p(exp) sum_(d) p(d) factor(exp,d)]/p(exp)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="d"/> is not a proper distribution</exception>
		public static Gamma ExpAverageConditional(Gamma exp, [Proper] Gaussian d, Gaussian to_d)
		{
			if (d.IsPointMass) return Gamma.PointMass(Math.Exp(d.Point));
			if (d.IsUniform()) return Gamma.FromShapeAndRate(0, 0);
			if (exp.IsPointMass) {
				// Z = int_y delta(x - exp(y)) N(y; my, vy) dy 
				//   = int_u delta(x - u) N(log(u); my, vy)/u du
				//   = N(log(x); my, vy)/x
				// logZ = -log(x) -0.5/vy*(log(x)-my)^2
				// dlogZ/dx = -1/x -1/vy*(log(x)-my)/x
				// d2logZ/dx2 = -dlogZ/dx/x -1/vy/x^2
				// log Ga(x;a,b) = (a-1)*log(x) - bx
				// dlogGa/dx = (a-1)/x - b
				// d2logGa/dx2 = -(a-1)/x^2
				// match derivatives and solve for (a,b)
				double shape = (1 + d.GetMean() - Math.Log(exp.Point)) * d.Precision;
				double rate = d.Precision / exp.Point;
				return Gamma.FromShapeAndRate(shape, rate);
			}
			if (exp.IsUniform()) return ExpAverageLogarithm(d);

			if (to_d.IsUniform() && exp.Shape > 1) {
				to_d = new Gaussian(MMath.Digamma(exp.Shape - 1) - Math.Log(exp.Rate), MMath.Trigamma(exp.Shape - 1));
			}

			double mD, vD;
			Gaussian dMarginal = d * to_d;
			dMarginal.GetMeanAndVariance(out mD, out vD);
			double Z = 0;
			double sumy = 0;
			double sumexpy = 0;

			if (vD < 1e-6) {
				double m, v;
				d.GetMeanAndVariance(out m, out v);
				return Gamma.FromLogMeanAndMeanLog(m + v / 2.0, m);
			}

			//if (vD < 10)
			if (true) {
				// Use Gauss-Hermite quadrature
				double[] nodes = new double[QuadratureNodeCount];
				double[] weights = new double[QuadratureNodeCount];

				Quadrature.GaussianNodesAndWeights(mD, vD, nodes, weights);
				for (int i = 0; i < weights.Length; i++) {
					weights[i] = Math.Log(weights[i]);
				}
				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] += d.GetLogProb(nodes[i]) - dMarginal.GetLogProb(nodes[i]);
					}
				}

				double maxLogF = Double.NegativeInfinity;
				// f(x,y) = Ga(exp(y); shape, rate) = exp(y*(shape-1) -rate*exp(y))
				// Z E[x] = int_y int_x x Ga(x;a,b) delta(x - exp(y)) N(y;my,vy) dx dy
				//        = int_y exp(y) Ga(exp(y);a,b) N(y;my,vy) dy
				// Z E[log(x)] = int_y y Ga(exp(y);a,b) N(y;my,vy) dy
				for (int i = 0; i < weights.Length; i++) {
					double y = nodes[i];
					double logf = weights[i] + (exp.Shape - 1) * y - exp.Rate * Math.Exp(y);
					if (logf > maxLogF) maxLogF = logf;
					weights[i] = logf;
				}
				for (int i = 0; i < weights.Length; i++) {
					double y = nodes[i];
					double f = Math.Exp(weights[i] - maxLogF);
					double f_y = f * y;
					double fexpy = f * Math.Exp(y);
					Z += f;
					sumy += f_y;
					sumexpy += fexpy;
				}
			} else {
				Converter<double, double> p = delegate(double y) {
					return d.GetLogProb(y) + (exp.Shape - 1) * y - exp.Rate * Math.Exp(y);
				};
				double sc = Math.Sqrt(vD);
				double offset = p(mD);
				Z = Quadrature.AdaptiveClenshawCurtis(z => Math.Exp(p(sc * z + mD) - offset), 1, 16, 1e-6);
				sumy = Quadrature.AdaptiveClenshawCurtis(z => (sc * z + mD) * Math.Exp(p(sc * z + mD) - offset), 1, 16, 1e-6);
				sumexpy = Quadrature.AdaptiveClenshawCurtis(z => Math.Exp(sc * z + mD + p(sc * z + mD) - offset), 1, 16, 1e-6);
			}
			if (Z == 0) throw new ApplicationException("Z==0");
			double s = 1.0 / Z;
			if (Double.IsPositiveInfinity(s)) throw new ApplicationException("s is -inf");
//.........这里部分代码省略.........

		// Perform one update of Q
		public static Gamma Q(Gamma sample, double shape, Gamma rate, Gamma q)
		{
			if (sample.IsUniform() || rate.IsPointMass) return rate;
			double a = q.Shape;
			double b = q.Rate;
			if (b==0) {
				a = rate.Shape;
				b = rate.Rate;
				if (shape > 0) {
					// this guess comes from solving dlogf=0 for x
					double guess = shape*sample.GetMeanInverse();
					b = Math.Max(rate.Rate, a/guess);
				}
			}
			double x = a/b;
			double x2 = x*x;
			double[] dlogfss = dlogfs(x, shape, sample);
			double dlogf = dlogfss[0];
			double ddlogf = dlogfss[1];
			b = rate.Rate - (dlogf + x*ddlogf);
			a = rate.Shape - x2*ddlogf;
			if (a <= 0) a = b*rate.Shape/(rate.Rate - dlogf);
			if (a <= 0 || b <= 0) throw new Exception();
			if (double.IsNaN(a) || double.IsNaN(b)) throw new Exception("result is nan");
			return Gamma.FromShapeAndRate(a, b);
		}

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

		/// <summary>
		/// Evidence message for EP
		/// </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>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(sample,mean,precision) p(sample,mean,precision) factor(sample,mean,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 double LogAverageFactor([SkipIfUniform] Gaussian sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision, Gamma to_precision)
		{
			if (precision.IsPointMass) return LogAverageFactor(sample, mean, precision.Point);
			if (precision.IsUniform()) return Double.PositiveInfinity;
			if (sample.IsPointMass && mean.IsPointMass) return LogAverageFactor(sample.Point, mean.Point, precision);
			if (sample.IsUniform() || mean.IsUniform()) return 0.0;
			// this code works even if sample and mean are point masses, but not if any variable is uniform.
			double xm, xv, mm, mv;
			sample.GetMeanAndVariance(out xm, out xv);
			mean.GetMeanAndVariance(out mm, out mv);
			// use quadrature to integrate over the precision
			double[] nodes = new double[QuadratureNodeCount];
			double[] logWeights = new double[nodes.Length];
			Gamma precMarginal = precision*to_precision;
			QuadratureNodesAndWeights(precMarginal, nodes, logWeights);
			for (int i = 0; i < nodes.Length; i++) {
				logWeights[i] = logWeights[i] + Gaussian.GetLogProb(xm, mm, xv + mv + 1.0 / nodes[i]);
			}
			if (!to_precision.IsUniform()) {
				// modify the weights
				for (int i = 0; i < logWeights.Length; i++) {
					logWeights[i] += precision.GetLogProb(nodes[i]) - precMarginal.GetLogProb(nodes[i]);
				}
			}
			return MMath.LogSumExp(logWeights);
		}

		/// <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, Gamma to_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 {
				//to_precision.SetToUniform(); // TEMPORARY
				// 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[] logWeights = new double[nodes.Length];
				Gamma precMarginal = precision*to_precision;
				QuadratureNodesAndWeights(precMarginal, nodes, logWeights);
				if (!to_precision.IsUniform()) {
					// modify the weights
					for (int i = 0; i < logWeights.Length; i++) {
						logWeights[i] += precision.GetLogProb(nodes[i]) - precMarginal.GetLogProb(nodes[i]);
					}
				}
				double shift = 0;
				MeanVarianceAccumulator mva = new MeanVarianceAccumulator();
				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 = Math.Exp(logWeights[i] + lp - shift);
					mva.Add(nodes[i], f);
				}

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

				double Z = mva.Count;
				if (double.IsNaN(Z)) throw new Exception("Z is nan");
				if (Z == 0.0) {
					throw new Exception("Quadrature found zero mass");
					//Console.WriteLine("Warning: Quadrature found zero mass.  Results may be inaccurate.");
					result.SetToUniform();
					return result;
				}
				double rmean = mva.Mean;
				double rvariance = mva.Variance;
				if (rvariance == 0) throw new Exception("Quadrature found zero variance");
				if (Double.IsInfinity(rmean)) {
					result.SetToUniform();
				} else {
					result.SetMeanAndVariance(rmean, rvariance);
					result.SetToRatio(result, precision, ForceProper);
				}
			}
//.........这里部分代码省略.........

		/// <summary>
		/// EP message to 'sample'
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'.</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 'sample' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'sample' as the random arguments are varied.
		/// The formula is <c>proj[p(sample) sum_(mean,precision) p(mean,precision) factor(sample,mean,precision)]/p(sample)</c>.
		/// </para></remarks>
		/// <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 Gaussian SampleAverageConditional(Gaussian sample, [SkipIfUniform] Gaussian mean, [SkipIfUniform] Gamma precision, Gamma to_precision)
		{
			if (sample.IsUniform() && precision.Shape <= 1.0) sample = Gaussian.FromNatural(1e-20, 1e-20);
			if (precision.IsPointMass) {
				return SampleAverageConditional(mean, precision.Point);
			} else if (sample.IsUniform()) {
				// for large vx, Z =approx N(mx; mm, vx+vm+E[1/prec])
				double mm,mv;
				mean.GetMeanAndVariance(out mm, out mv);
				// NOTE: this error may happen because sample didn't receive any message yet under the schedule.
				// Need to make the scheduler smarter to avoid this.
				if (precision.Shape <= 1.0) throw new ArgumentException("The posterior has infinite variance due to precision distributed as "+precision+" (shape <= 1).  Try using a different prior for the precision, with shape > 1.");
				return Gaussian.FromMeanAndVariance(mm, mv + precision.GetMeanInverse());
			} else if (mean.IsUniform() || precision.IsUniform()) {
				return Gaussian.Uniform();
			} else if (sample.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.
				return Gaussian.Uniform();
			} else if (!precision.IsProper()) {
				throw new ImproperMessageException(precision);
			} else {
				// The formula is int_prec int_mean N(x;mean,1/prec) p(x) p(mean) p(prec) =
				// int_prec N(x; mm, mv + 1/prec) p(x) p(prec) =
				// int_prec N(x; new xm, new xv) N(xm; mm, mv + xv + 1/prec) p(prec)
				// Let R = Prec/(Prec + mean.Prec)
				// new xv = inv(R*mean.Prec + sample.Prec)
				// new xm = xv*(R*mean.PM + sample.PM)

				// In the case where sample and mean are improper distributions, 
				// we must only consider values of prec for which (new xv > 0).
				// This happens when R*mean.Prec > -sample.Prec
				// As a function of Prec, R*mean.Prec has a singularity at Prec=-mean.Prec
				// This function is greater than a threshold when Prec is sufficiently small or sufficiently large.
				// Therefore we construct an interval of Precs to exclude from the integration.
				double xm, xv, mm, mv;
				sample.GetMeanAndVarianceImproper(out xm, out xv);
				mean.GetMeanAndVarianceImproper(out mm, out mv);
				double lowerBound = 0;
				double upperBound = Double.PositiveInfinity;
				bool precisionIsBetween = true;
				if (mean.Precision >= 0) {
					if (sample.Precision < -mean.Precision) throw new ImproperMessageException(sample);
					//lowerBound = -mean.Precision * sample.Precision / (mean.Precision + sample.Precision);
					lowerBound = -1.0 / (xv + mv);
				} else {  // mean.Precision < 0
					if (sample.Precision < 0) {
						precisionIsBetween = true;
						lowerBound = -1.0 / (xv + mv);
						upperBound = -mean.Precision;
					} else if (sample.Precision < -mean.Precision) {
						precisionIsBetween = true;
						lowerBound = 0;
						upperBound = -mean.Precision;
					} else {
						// in this case, the precision should NOT be in this interval.
						precisionIsBetween = false;
						lowerBound = -mean.Precision;
						lowerBound = -1.0 / (xv + mv);
					}
				}
				double[] nodes = new double[QuadratureNodeCount];
				double[] logWeights = new double[nodes.Length];
				Gamma precMarginal = precision*to_precision;
				QuadratureNodesAndWeights(precMarginal, nodes, logWeights);
				if (!to_precision.IsUniform()) {
					// modify the weights
					for (int i = 0; i < logWeights.Length; i++) {
						logWeights[i] += precision.GetLogProb(nodes[i]) - precMarginal.GetLogProb(nodes[i]);
					}
				}
				GaussianEstimator est = new GaussianEstimator();
				double shift = 0;
				for (int i = 0; i < nodes.Length; i++) {
					double newVar, newMean;
					Assert.IsTrue(nodes[i] > 0);
					if ((nodes[i] > lowerBound && nodes[i] < upperBound) != precisionIsBetween) continue;
					// the following works even if sample is uniform. (sample.Precision == 0)
					if (mean.IsPointMass) {
						// take limit mean.Precision -> Inf
						newVar = 1.0 / (nodes[i] + sample.Precision);
						newMean = newVar * (nodes[i] * mean.Point + sample.MeanTimesPrecision);
					} else {
						// mean.Precision < Inf
						double R = nodes[i] / (nodes[i] + mean.Precision);
						newVar = 1.0 / (R * mean.Precision + sample.Precision);
//.........这里部分代码省略.........

		public static Gamma Q(Gaussian sample, Gaussian mean, Gamma precision, Gamma q)
		{
			if (precision.IsPointMass || sample.IsUniform() || mean.IsUniform()) return precision;
			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;
			if (double.IsInfinity(v)) return precision;
			double m2 = m*m;
			double[] dlogfss;
			double x = q.GetMean();
			double a = precision.Shape;
			double b = precision.Rate;
			if (double.IsPositiveInfinity(x)) x = (a + Math.Sqrt(a*a + 2*b/v))/(2*b);
			// TODO: check initialization
			if (q.IsUniform()) x = 0;
			//x = QReinitialize(sample, mean, precision, x);
			for (int iter = 0; iter < 1000; iter++) {
				double oldx = x;
				if (true) {
					// want to find x that optimizes -0.5*log(v + 1/x) - 0.5*m2/(v + 1/x) + a*log(x) - b*x
					// we fit a lower bound, then optimize the bound.
					// this is guaranteed to improve x.
					double logSlope = precision.Shape;
					double slope = -precision.Rate;
					double denom = v*x+1;
					// 1/(v+1/x) <= 1/(v+1/x0) + (x-x0)/(v*x0+1)^2
					slope += -0.5*m2/(denom*denom);
					if (v*x < 1) {
						// log(v+1/x)  = log(v*x+1) - log(x)
						// log(v*x+1) <= log(v*x0+1) + (x-x0)*v/(v*x0+1)
						logSlope += 0.5;
						slope += -0.5*v/denom;
						x = -logSlope/slope;
						// at x=0 the update is x = (a+0.5)/(b + 0.5*(m2+v))
						// at x=inf the update is x = (a+0.5)/b
					} else {
						// if v*x > 1:
						// log(v+1/x) <= log(v+1/x0) + (1/x - 1/x0)/(v + 1/x0)
						double invSlope = -0.5*x/denom;
						// solve for the maximum of logslope*log(r)+slope*r+invslope./r
						// logslope/r + slope - invslope/r^2 = 0
						// logslope*r + slope*r^2 - invslope = 0
						//x = (-logSlope - Math.Sqrt(logSlope*logSlope + 4*invSlope*slope))/(2*slope);
						double c = 0.5*logSlope/slope;
						double d = invSlope/slope;
						// note c < 0 always
						x = Math.Sqrt(c*c + d) - c;
						// at x=inf, invSlope=-0.5/v so the update is x = (ax + sqrt(ax*ax + 2*b/v))/(2*b)
					}
					if (x < 0) throw new Exception("x < 0");
				} else {
					dlogfss = dlogfs(x, m, v);
					if (true) {
						x = (precision.Shape - x*x*dlogfss[1])/(precision.Rate - dlogfss[0] - x*dlogfss[1]);
					} else if (true) {
						double delta = dlogfss[0] - precision.Rate;
						x *= Math.Exp(-(delta*x + precision.Shape)/(delta*x + dlogfss[1]*x));
					} else {
						x = precision.Shape/(precision.Rate - dlogfss[0]);
					}
					if (x < 0) throw new Exception("x < 0");
				}
				if (double.IsNaN(x)) throw new Exception("x is nan");
				//Console.WriteLine("{0}: {1}", iter, x);
				if (MMath.AbsDiff(oldx, x, 1e-10) < 1e-10) {
					x = QReinitialize(sample, mean, precision, x);
					if (MMath.AbsDiff(oldx, x, 1e-10) < 1e-10) break;
				}
				if (iter == 1000-1) throw new Exception("not converging");
			}
			//x = r0;
			dlogfss = dlogfs(x, m, v);
			double dlogf = dlogfss[0];
			double ddlogf = dlogfss[1];
			a = precision.Shape - x*x*ddlogf;
			b = precision.Rate - (dlogf + x*ddlogf);
			if (a <= 0 || b <= 0) throw new Exception();
			if (double.IsNaN(a) || double.IsNaN(b)) throw new Exception("result is nan");
			return Gamma.FromShapeAndRate(a, b);
		}