View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.rng.examples.sampling;
19  
20  import java.io.PrintWriter;
21  import java.io.IOException;
22  import org.apache.commons.rng.UniformRandomProvider;
23  import org.apache.commons.rng.simple.RandomSource;
24  import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
25  import org.apache.commons.rng.sampling.distribution.MarsagliaNormalizedGaussianSampler;
26  import org.apache.commons.rng.sampling.distribution.BoxMullerNormalizedGaussianSampler;
27  import org.apache.commons.rng.sampling.distribution.ChengBetaSampler;
28  import org.apache.commons.rng.sampling.distribution.AhrensDieterExponentialSampler;
29  import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGammaSampler;
30  import org.apache.commons.rng.sampling.distribution.InverseTransformParetoSampler;
31  import org.apache.commons.rng.sampling.distribution.LogNormalSampler;
32  import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
33  import org.apache.commons.rng.sampling.distribution.GaussianSampler;
34  import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
35  
36  /**
37   * Approximation of the probability density by the histogram of the sampler output.
38   */
39  public class ProbabilityDensityApproximation {
40      /** Number of (equal-width) bins in the histogram. */
41      private final int numBins;
42      /** Number of samples to be generated. */
43      private final long numSamples;
44  
45      /**
46       * Application.
47       *
48       * @param numBins Number of "equal-width" bins.
49       * @param numSamples Number of samples.
50       */
51      private ProbabilityDensityApproximation(int numBins,
52                                              long numSamples) {
53          this.numBins = numBins;
54          this.numSamples = numSamples;
55      }
56  
57      /**
58       * @param sampler Sampler.
59       * @param min Right abscissa of the first bin: every sample smaller
60       * than that value will increment an additional bin (of infinite width)
61       * placed before the first "equal-width" bin.
62       * @param Left abscissa of the last bin: every sample larger than or
63       * equal to that value will increment an additional bin (of infinite
64       * width) placed after the last "equal-width" bin.
65       * @param output Filename.
66       */
67      private void createDensity(ContinuousSampler sampler,
68                                 double min,
69                                 double max,
70                                 String outputFile)
71          throws IOException {
72          final double binSize = (max - min) / numBins;
73          final long[] histogram = new long[numBins];
74  
75          long n = 0;
76          long belowMin = 0;
77          long aboveMax = 0;
78          while (++n < numSamples) {
79              final double r = sampler.sample();
80  
81              if (r < min) {
82                  ++belowMin;
83                  continue;
84              }
85  
86              if (r >= max) {
87                  ++aboveMax;
88                  continue;
89              }
90  
91              final int binIndex = (int) ((r - min) / binSize);
92              ++histogram[binIndex];
93          }
94  
95          final double binHalfSize = 0.5 * binSize;
96          final double norm = 1 / (binSize * numSamples);
97  
98          final PrintWriter out = new PrintWriter(outputFile);
99          out.println("# Sampler: " + sampler);
100         out.println("# Number of bins: " + numBins);
101         out.println("# Min: " + min + " (fraction of samples below: " + (belowMin / (double) numSamples) + ")");
102         out.println("# Max: " + max + " (fraction of samples above: " + (aboveMax / (double) numSamples) + ")");
103         out.println("# Bin width: " + binSize);
104         out.println("# Histogram normalization factor: " + norm);
105         out.println("#");
106         out.println("# " + (min - binHalfSize) + " " + (belowMin * norm));
107         for (int i = 0; i < numBins; i++) {
108             out.println((min + (i + 1) * binSize - binHalfSize) + " " + (histogram[i] * norm));
109         }
110         out.println("# " + (max + binHalfSize) + " " + (aboveMax * norm));
111         out.close();
112     }
113 
114     /**
115      * Program entry point.
116      *
117      * @param args Argument. They must be provided, in the following order:
118      * <ol>
119      *  <li>Number of "equal-width" bins.</li>
120      *  <li>Number of samples.</li>
121      * </ol>
122      * @throws IOException if failure occurred while writing to files.
123      */
124     public static void main(String[] args)
125         throws IOException {
126         final int numBins = Integer.valueOf(args[0]);
127         final long numSamples = Long.valueOf(args[1]);
128         final ProbabilityDensityApproximation app = new ProbabilityDensityApproximation(numBins, numSamples);
129 
130         final UniformRandomProvider rng = RandomSource.create(RandomSource.XOR_SHIFT_1024_S);
131 
132         final double gaussMean = 1;
133         final double gaussSigma = 2;
134         final double gaussMin = -9;
135         final double gaussMax = 11;
136         app.createDensity(new GaussianSampler(new ZigguratNormalizedGaussianSampler(rng),
137                                               gaussMean, gaussSigma),
138                           gaussMin, gaussMax, "gauss.ziggurat.txt");
139         app.createDensity(new GaussianSampler(new MarsagliaNormalizedGaussianSampler(rng),
140                                               gaussMean, gaussSigma),
141                           gaussMin, gaussMax, "gauss.marsaglia.txt");
142         app.createDensity(new GaussianSampler(new BoxMullerNormalizedGaussianSampler(rng),
143                                               gaussMean, gaussSigma),
144                           gaussMin, gaussMax, "gauss.boxmuller.txt");
145 
146         final double alphaBeta = 4.3;
147         final double betaBeta = 2.1;
148         final double betaMin = 0;
149         final double betaMax = 1;
150         app.createDensity(new ChengBetaSampler(rng, alphaBeta, betaBeta),
151                           betaMin, betaMax, "beta.case1.txt");
152         final double alphaBetaAlt = 0.5678;
153         final double betaBetaAlt = 0.1234;
154         app.createDensity(new ChengBetaSampler(rng, alphaBetaAlt, betaBetaAlt),
155                           betaMin, betaMax, "beta.case2.txt");
156 
157         final double meanExp = 3.45;
158         final double expMin = 0;
159         final double expMax = 60;
160         app.createDensity(new AhrensDieterExponentialSampler(rng, meanExp),
161                           expMin, expMax, "exp.txt");
162 
163         final double thetaGammaSmallerThanOne = 0.1234;
164         final double alphaGamma = 3.456;
165         final double gammaMin = 0;
166         final double gammaMax1 = 40;
167         app.createDensity(new AhrensDieterMarsagliaTsangGammaSampler(rng, alphaGamma, thetaGammaSmallerThanOne),
168                           gammaMin, gammaMax1, "gamma.case1.txt");
169         final double thetaGammaLargerThanOne = 2.345;
170         final double gammaMax2 = 70;
171         app.createDensity(new AhrensDieterMarsagliaTsangGammaSampler(rng, alphaGamma, thetaGammaLargerThanOne),
172                           gammaMin, gammaMax2, "gamma.case2.txt");
173 
174         final double scalePareto = 23.45;
175         final double shapePareto = 0.789;
176         final double paretoMin = 23;
177         final double paretoMax = 400;
178         app.createDensity(new InverseTransformParetoSampler(rng, scalePareto, shapePareto),
179                           paretoMin, paretoMax, "pareto.txt");
180 
181         final double loUniform = -9.876;
182         final double hiUniform = 5.432;
183         app.createDensity(new ContinuousUniformSampler(rng, loUniform, hiUniform),
184                           loUniform, hiUniform, "uniform.txt");
185 
186         final double scaleLogNormal = 2.345;
187         final double shapeLogNormal = 0.1234;
188         final double logNormalMin = 5;
189         final double logNormalMax = 25;
190         app.createDensity(new LogNormalSampler(new ZigguratNormalizedGaussianSampler(rng),
191                                                scaleLogNormal, shapeLogNormal),
192                           logNormalMin, logNormalMax, "lognormal.ziggurat.txt");
193         app.createDensity(new LogNormalSampler(new MarsagliaNormalizedGaussianSampler(rng),
194                                                scaleLogNormal, shapeLogNormal),
195                           logNormalMin, logNormalMax, "lognormal.marsaglia.txt");
196         app.createDensity(new LogNormalSampler(new BoxMullerNormalizedGaussianSampler(rng),
197                                                scaleLogNormal, shapeLogNormal),
198                           logNormalMin, logNormalMax, "lognormal.boxmuller.txt");
199     }
200 }