Mercurial > hg > jslab
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 } |