samer@0: /*
samer@0:  *	Copyright (c) 2000, Samer Abdallah, King's College London.
samer@0:  *	All rights reserved.
samer@0:  *
samer@0:  *	This software is provided AS iS and WITHOUT ANY WARRANTY; 
samer@0:  *	without even the implied warranty of MERCHANTABILITY or 
samer@0:  *	FITNESS FOR A PARTICULAR PURPOSE.
samer@0:  */
samer@0: 
samer@0: package samer.maths.random;
samer@0: import  samer.core.*;
samer@0: import  samer.core.types.*;
samer@0: import  samer.maths.*;
samer@0: 
samer@0: // This is going to work by the rejection method
samer@0: // (see Numerical Recipes in C)
samer@0: // We'll use the cauchy distribution as our
samer@0: // bounding function
samer@0: 
samer@0: // NB ONLY WORKS FOR alpha<=2!
samer@0: 
samer@0: public class GeneralisedExponential implements Generator
samer@0: {
samer@0: 	BaseRandom	cauchy = new RectifiedCauchy();
samer@0: 	VDouble		alpha = new VDouble("alpha", 1);
samer@0: 
samer@0: 	public void dispose() { cauchy.dispose(); alpha.dispose(); }
samer@0: 	public double next()
samer@0: 	{ 
samer@0: 		double x, f, g, a=alpha.value;
samer@0: 		do {
samer@0: 			// first, get a Cauchy distributed RV
samer@0: 			x = 2*cauchy.next();
samer@0: 			f = 4/(4+x*x); // cauchy
samer@0: 			g = Math.exp(-Math.pow(x,a)); // exponential
samer@0: 		} while (f*cauchy.rnd.nextDouble()>g);
samer@0: 		return x; 
samer@0: 	}
samer@0: 
samer@0: 	public void next(double [] x)
samer@0: 	{
samer@0: 		for (int i=0; i<x.length; i++) x[i]=next();
samer@0: 	}
samer@0: }