001 /*
002 * Created on Aug 30, 2007
003 */
004 package com.x8ing.mc.distribution;
005
006 import cern.jet.random.Normal;
007 import cern.jet.random.engine.MersenneTwister;
008
009 /**
010 * Uses an open library by cern.ch
011 * <p>
012 * http://spi.cern.ch/extsoft/packages.php?pkg=Colt
013 * <p>
014 * http://dsd.lbl.gov/~hoschek/colt/
015 *
016 * @author Patrick Heusser
017 */
018 public class RandomDistributionGauss implements RandomDistribution {
019
020 private Normal normal = null;
021
022 private static final int MAX_SEARCH_ITERATIONS = 5000;
023
024 private double standardDeviation = 0;
025
026 private double min = -1;
027
028 private double max = -1;
029
030 /**
031 * <p>
032 *
033 * nStdDeviations = percentage <br>
034 * 1 = 0.682689492137 <br>
035 * 2 = 0.954499736104 <br>
036 * 3 = 0.997300203937 <br>
037 * 4 = 0.999936657516 <br>
038 * 5 = 0.999999426697 <br>
039 * 6 = 0.999999998027
040 *
041 * @param min
042 * the minimal return value
043 * @param max
044 * the maximum return value
045 * @param nStdDeviations:
046 * e.g. 1,2,3 etc... as described above..
047 * @return a value within range, even if it's hard bound!
048 */
049 public RandomDistributionGauss(double min, double max, int nStdDeviations) {
050
051 super();
052
053 normal = new Normal(-1, 1, new MersenneTwister());
054
055 init(min, max, nStdDeviations);
056
057 }
058
059 public void init(double min, double max, int nStdDeviations) {
060
061 if (nStdDeviations < 0.5 || nStdDeviations > 20) {
062
063 throw new IllegalStateException("error while initialization. 'stdDeviation must be within 0.5-6"
064 + ". value where: min=" + min + " max=" + max + " nStdDeviations=" + nStdDeviations);
065
066 }
067
068 double mean = (((double) (max + min)) / 2);
069
070 this.standardDeviation = (max - mean) / nStdDeviations;
071 this.max = max;
072 this.min = min;
073
074 if (min >= max) {
075 throw new IllegalStateException("error while initialization. 'min' value must be smaller "
076 + "then 'max' value. but values where: min=" + min + " max=" + max + " nStdDeviations=" + nStdDeviations);
077 }
078
079 normal.setState(mean, standardDeviation);
080
081 }
082
083 /**
084 * @see com.x8ing.mc.distribution.RandomDistribution#getNextRandomNumber()
085 */
086 public double getNextRandomNumber() {
087
088 double ret = 0;
089
090 int searchloop = 0;
091
092 do {
093 searchloop++;
094 ret = normal.nextDouble();
095
096 // limit the endless loop
097 if (searchloop >= MAX_SEARCH_ITERATIONS) {
098 ret = min;
099 System.err.println("warning, searched to long for gaussian number in hard limited range."
100 + " please adjust numbers. min=" + min + " max=" + max + " standardDeviation=" + standardDeviation);
101 break;
102 }
103
104 } while ((ret > max || ret < min));
105
106 return ret;
107
108 }
109
110 }
|