C# Gaussian.SetMeanAndPrecision方法代码示例

		public static Gaussian SampleAverageConditional(Gaussian sample, [SkipIfUniform] Gaussian mean, [Proper] Gamma precision, [Fresh] Gamma q)
		{
			if (mean.IsUniform() || sample.IsPointMass) return Gaussian.Uniform();
			if (precision.IsPointMass) return GaussianOp.SampleAverageConditional(mean, precision.Point);
			if (q.IsPointMass) throw new Exception();
			double mm, vm;
			mean.GetMeanAndVariance(out mm, out vm);
			if (sample.IsUniform()) {
				if (precision.Shape <= 1.0) {
					sample = Gaussian.FromNatural(1e-20, 1e-20);  // try to proceed instead of throwing an exception
					//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.");
				} else return Gaussian.FromMeanAndVariance(mm, vm + precision.GetMeanInverse());
			}
			double mx, vx;
			sample.GetMeanAndVariance(out mx, out vx);
			double m = mx-mm;
			double v = vx+vm;
			Gaussian sampleMarginal;
			if (true) {
				Gaussian qx = Qx(Gaussian.FromMeanAndVariance(m, v), precision, Gaussian.FromMeanAndPrecision(0, 1));
				double y = qx.GetMean();
				double[] dlogfx = dlogfxs(y, precision);
				double delta = 0.5*dlogfx[2]/(qx.Precision*qx.Precision);
				// sampleMean can be computed directly as:
				//double dlogz = dlogfx[0] + delta/v;
				//double sampleMean = mx + vx*dlogz;
				double my = y + delta;
				double vy = 1/qx.Precision + 4*delta*delta + 0.5*dlogfx[3]/(qx.Precision*qx.Precision*qx.Precision);
				if (vy < 0) throw new Exception();
				Gaussian yMsg = new Gaussian(my, vy)/(new Gaussian(m, v));
				Gaussian sampleMsg = DoublePlusOp.SumAverageConditional(yMsg, mean);
				sampleMarginal = sampleMsg*sample;
				if (!sampleMarginal.IsProper()) throw new Exception();
				if (sampleMarginal.IsPointMass) throw new Exception();
				//Console.WriteLine("{0} {1}", sampleMean, sampleMarginal);
			} else {
				double a = q.Shape;
				double b = q.Rate;
				double x = a/b;
				if (sample.IsPointMass || mean.IsPointMass) {
					double denom = 1 + x*v;
					double denom2 = denom*denom;
					Gaussian y = sample*mean;
					double my, vy;
					y.GetMeanAndVariance(out my, out vy);
					// 1-1/denom = x*v/denom
					// sampleMean = E[ (x*(mx*vm + mm*vx) + mx)/denom ]
					//            = E[ (mx*vm + mm*vx)/v*(1-1/denom) + mx/denom ]
					//            = E[ (mx/vx + mm/vm)/(1/vx + 1/vm)*(1-1/denom) + mx/denom ]
					// sampleVariance = var((my*(1-1/denom) + mx/denom)) + E[ (r*vx*vm + vx)/denom ]
					//                = (mx-my)^2*var(1/denom) + E[ vx*vm/v*(1-1/denom) + vx/denom ]
					double[] g = new double[] { 1/denom, -v/denom2, 2*v*v/(denom2*denom), -6*v*v*v/(denom2*denom2) };
					double [] dlogf = dlogfs(x, m, v);
					double edenom, vdenom;
					LaplaceMoments(q, g, dlogf, out edenom, out vdenom);
					double sampleMean = mx*edenom + my*(1 - edenom);
					double diff = mx-my;
					double sampleVariance = vx*edenom + vy*(1 - edenom) + diff*diff*vdenom;
					sampleMarginal = Gaussian.FromMeanAndVariance(sampleMean, sampleVariance);
				} else {
					// 1 - samplePrec*mPrec/denom = x*yprec/denom
					// sampleMean = E[ (x*ymprec + sampleMP*mPrec)/denom ]
					//            = E[ (1 - samplePrec*mPrec/denom)*ymprec/yprec + sampleMP*mPrec/denom ]
					// sampleVariance = var((1 - samplePrec*mPrec/denom)*ymprec/yprec + sampleMP*mPrec/denom) + E[ (x+mPrec)/denom ]
					//                = (sampleMP*mPrec - samplePrec*mPrec*ymprec/yprec)^2*var(1/denom) + E[ (1-samplePrec*mPrec/denom)/yprec + mPrec/denom ]
					double yprec = sample.Precision + mean.Precision;
					double ymprec = sample.MeanTimesPrecision + mean.MeanTimesPrecision;
					double denom = sample.Precision*mean.Precision + x*yprec;
					double denom2 = denom*denom;
					double[] g = new double[] { 1/denom, -yprec/denom2, 2*yprec*yprec/(denom2*denom), -6*yprec*yprec*yprec/(denom2*denom2) };
					double [] dlogf = dlogfs(x, m, v);
					double edenom, vdenom;
					LaplaceMoments(q, g, dlogf, out edenom, out vdenom);
					double sampleMean = sample.MeanTimesPrecision*mean.Precision*edenom + ymprec/yprec*(1-sample.Precision*mean.Precision*edenom);
					double diff = sample.MeanTimesPrecision*mean.Precision - sample.Precision*mean.Precision*ymprec/yprec;
					double sampleVariance = mean.Precision*edenom + (1 - sample.Precision*mean.Precision*edenom)/yprec + diff*diff*vdenom;
					sampleMarginal = Gaussian.FromMeanAndVariance(sampleMean, sampleVariance);
				}
			}
			Gaussian result = new Gaussian();
			result.SetToRatio(sampleMarginal, sample, GaussianOp.ForceProper);
			if (modified2) {
				if(!mean.IsPointMass) throw new Exception();
				result.SetMeanAndPrecision(mm, q.GetMean());
			} else if (modified && !sample.IsUniform()) {
				// heuristic method to avoid improper messages:
				// the message's mean must be E[mean] (regardless of context) and the variance is chosen to match the posterior mean when multiplied by context
				double sampleMean = sampleMarginal.GetMean();
				if (sampleMean != mm) {
					result.Precision = (sample.MeanTimesPrecision-sampleMean*sample.Precision)/(sampleMean - mm);
					if (result.Precision < 0) throw new Exception("internal: sampleMean is not between sample.Mean and mean.Mean");
					result.MeanTimesPrecision = result.Precision*mm;
				}
			}
			//if (!result.IsProper()) throw new Exception();
			if (result.Precision < -0.001) throw new Exception();
			if (double.IsNaN(result.MeanTimesPrecision) || double.IsNaN(result.Precision)) throw new Exception("result is nan");
			return result;
		}