samer@0: package samer.maths.random;
samer@0: import samer.maths.*;
samer@0: import samer.core.*;
samer@0: 
samer@0: /**
samer@0:  sample from posterior
samer@0:  instead of supplying weighting factor,
samer@0:  reject a fraction of the samples 
samer@0:  make take a long time if A<<sigma
samer@0: */
samer@0: 
samer@0: class PosteriorSampler extends BaseRandom
samer@0: {
samer@0: 	double A, sigma;
samer@0: 	double std, var, mean;
samer@0: 
samer@0: 	public PosteriorSampler() {}
samer@0: 
samer@0: 	public void setA(double a) { A=a; std=sigma/A; var=std*std; }
samer@0: 	public void setSigma(double s) { sigma=s; std=sigma/A; var=std*std; }
samer@0: 	public void setx(double x) { mean=x/A; }
samer@0: 
samer@0: 	public double next() 
samer@0: 	{
samer@0: 		double v;
samer@0: 		boolean ok=false;
samer@0: 
samer@0: 		do {
samer@0: 			v=mean + std*rnd.nextGaussian();
samer@0: 
samer@0: 			if (v<=-var) {
samer@0: 				v=v+var;
samer@0: 				if (mean<0) {
samer@0: 					double k=Math.exp(2*mean);
samer@0: 					double u=rnd.nextDouble();
samer@0: 					if (u<k) ok=true;
samer@0: 				} else ok=true;
samer@0: 			} else if (v>=var) {
samer@0: 				v=v-var;
samer@0: 				if (mean>0) {
samer@0: 					double k=Math.exp(-2*mean);
samer@0: 					double u=rnd.nextDouble();
samer@0: 					if (u<k) ok=true;
samer@0: 				} else ok=true;
samer@0: 			}
samer@0: 		} while (!ok);
samer@0: 
samer@0: 		return v;
samer@0: 	}
samer@0: }