annotate src/samer/maths/random/GeneralisedExponential.java @ 8:5e3cbbf173aa tip

Reorganise some more
author samer
date Fri, 05 Apr 2019 22:41:58 +0100
parents bf79fb79ee13
children
rev   line source
samer@0 1 /*
samer@0 2 * Copyright (c) 2000, Samer Abdallah, King's College London.
samer@0 3 * All rights reserved.
samer@0 4 *
samer@0 5 * This software is provided AS iS and WITHOUT ANY WARRANTY;
samer@0 6 * without even the implied warranty of MERCHANTABILITY or
samer@0 7 * FITNESS FOR A PARTICULAR PURPOSE.
samer@0 8 */
samer@0 9
samer@0 10 package samer.maths.random;
samer@0 11 import samer.core.*;
samer@0 12 import samer.core.types.*;
samer@0 13 import samer.maths.*;
samer@0 14
samer@0 15 // This is going to work by the rejection method
samer@0 16 // (see Numerical Recipes in C)
samer@0 17 // We'll use the cauchy distribution as our
samer@0 18 // bounding function
samer@0 19
samer@0 20 // NB ONLY WORKS FOR alpha<=2!
samer@0 21
samer@0 22 public class GeneralisedExponential implements Generator
samer@0 23 {
samer@0 24 BaseRandom cauchy = new RectifiedCauchy();
samer@0 25 VDouble alpha = new VDouble("alpha", 1);
samer@0 26
samer@0 27 public void dispose() { cauchy.dispose(); alpha.dispose(); }
samer@0 28 public double next()
samer@0 29 {
samer@0 30 double x, f, g, a=alpha.value;
samer@0 31 do {
samer@0 32 // first, get a Cauchy distributed RV
samer@0 33 x = 2*cauchy.next();
samer@0 34 f = 4/(4+x*x); // cauchy
samer@0 35 g = Math.exp(-Math.pow(x,a)); // exponential
samer@0 36 } while (f*cauchy.rnd.nextDouble()>g);
samer@0 37 return x;
samer@0 38 }
samer@0 39
samer@0 40 public void next(double [] x)
samer@0 41 {
samer@0 42 for (int i=0; i<x.length; i++) x[i]=next();
samer@0 43 }
samer@0 44 }