org.apache.commons.rng.examples.sampling.ProbabilityDensityApproximation Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of commons-rng-examples-sampling Show documentation
Show all versions of commons-rng-examples-sampling Show documentation
Contains examples of output from the samplers.
Code in this module is not part of the public API.
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.rng.examples.sampling;
import java.io.PrintWriter;
import java.io.IOException;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
import org.apache.commons.rng.sampling.distribution.MarsagliaNormalizedGaussianSampler;
import org.apache.commons.rng.sampling.distribution.BoxMullerNormalizedGaussianSampler;
import org.apache.commons.rng.sampling.distribution.ChengBetaSampler;
import org.apache.commons.rng.sampling.distribution.AhrensDieterExponentialSampler;
import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGammaSampler;
import org.apache.commons.rng.sampling.distribution.InverseTransformParetoSampler;
import org.apache.commons.rng.sampling.distribution.LogNormalSampler;
import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
import org.apache.commons.rng.sampling.distribution.GaussianSampler;
import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
/**
* Approximation of the probability density by the histogram of the sampler output.
*/
public class ProbabilityDensityApproximation {
/** Number of (equal-width) bins in the histogram. */
private final int numBins;
/** Number of samples to be generated. */
private final long numSamples;
/**
* Application.
*
* @param numBins Number of "equal-width" bins.
* @param numSamples Number of samples.
*/
private ProbabilityDensityApproximation(int numBins,
long numSamples) {
this.numBins = numBins;
this.numSamples = numSamples;
}
/**
* @param sampler Sampler.
* @param min Right abscissa of the first bin: every sample smaller
* than that value will increment an additional bin (of infinite width)
* placed before the first "equal-width" bin.
* @param Left abscissa of the last bin: every sample larger than or
* equal to that value will increment an additional bin (of infinite
* width) placed after the last "equal-width" bin.
* @param output Filename.
*/
private void createDensity(ContinuousSampler sampler,
double min,
double max,
String outputFile)
throws IOException {
final double binSize = (max - min) / numBins;
final long[] histogram = new long[numBins];
long n = 0;
long belowMin = 0;
long aboveMax = 0;
while (++n < numSamples) {
final double r = sampler.sample();
if (r < min) {
++belowMin;
continue;
}
if (r >= max) {
++aboveMax;
continue;
}
final int binIndex = (int) ((r - min) / binSize);
++histogram[binIndex];
}
final double binHalfSize = 0.5 * binSize;
final double norm = 1 / (binSize * numSamples);
final PrintWriter out = new PrintWriter(outputFile);
out.println("# Sampler: " + sampler);
out.println("# Number of bins: " + numBins);
out.println("# Min: " + min + " (fraction of samples below: " + (belowMin / (double) numSamples) + ")");
out.println("# Max: " + max + " (fraction of samples above: " + (aboveMax / (double) numSamples) + ")");
out.println("# Bin width: " + binSize);
out.println("# Histogram normalization factor: " + norm);
out.println("#");
out.println("# " + (min - binHalfSize) + " " + (belowMin * norm));
for (int i = 0; i < numBins; i++) {
out.println((min + (i + 1) * binSize - binHalfSize) + " " + (histogram[i] * norm));
}
out.println("# " + (max + binHalfSize) + " " + (aboveMax * norm));
out.close();
}
/**
* Program entry point.
*
* @param args Argument. They must be provided, in the following order:
*
* - Number of "equal-width" bins.
* - Number of samples.
*
* @throws IOException if failure occurred while writing to files.
*/
public static void main(String[] args)
throws IOException {
final int numBins = Integer.valueOf(args[0]);
final long numSamples = Long.valueOf(args[1]);
final ProbabilityDensityApproximation app = new ProbabilityDensityApproximation(numBins, numSamples);
final UniformRandomProvider rng = RandomSource.create(RandomSource.XOR_SHIFT_1024_S);
final double gaussMean = 1;
final double gaussSigma = 2;
final double gaussMin = -9;
final double gaussMax = 11;
app.createDensity(new GaussianSampler(new ZigguratNormalizedGaussianSampler(rng),
gaussMean, gaussSigma),
gaussMin, gaussMax, "gauss.ziggurat.txt");
app.createDensity(new GaussianSampler(new MarsagliaNormalizedGaussianSampler(rng),
gaussMean, gaussSigma),
gaussMin, gaussMax, "gauss.marsaglia.txt");
app.createDensity(new GaussianSampler(new BoxMullerNormalizedGaussianSampler(rng),
gaussMean, gaussSigma),
gaussMin, gaussMax, "gauss.boxmuller.txt");
final double alphaBeta = 4.3;
final double betaBeta = 2.1;
final double betaMin = 0;
final double betaMax = 1;
app.createDensity(new ChengBetaSampler(rng, alphaBeta, betaBeta),
betaMin, betaMax, "beta.case1.txt");
final double alphaBetaAlt = 0.5678;
final double betaBetaAlt = 0.1234;
app.createDensity(new ChengBetaSampler(rng, alphaBetaAlt, betaBetaAlt),
betaMin, betaMax, "beta.case2.txt");
final double meanExp = 3.45;
final double expMin = 0;
final double expMax = 60;
app.createDensity(new AhrensDieterExponentialSampler(rng, meanExp),
expMin, expMax, "exp.txt");
final double thetaGammaSmallerThanOne = 0.1234;
final double alphaGamma = 3.456;
final double gammaMin = 0;
final double gammaMax1 = 40;
app.createDensity(new AhrensDieterMarsagliaTsangGammaSampler(rng, alphaGamma, thetaGammaSmallerThanOne),
gammaMin, gammaMax1, "gamma.case1.txt");
final double thetaGammaLargerThanOne = 2.345;
final double gammaMax2 = 70;
app.createDensity(new AhrensDieterMarsagliaTsangGammaSampler(rng, alphaGamma, thetaGammaLargerThanOne),
gammaMin, gammaMax2, "gamma.case2.txt");
final double scalePareto = 23.45;
final double shapePareto = 0.789;
final double paretoMin = 23;
final double paretoMax = 400;
app.createDensity(new InverseTransformParetoSampler(rng, scalePareto, shapePareto),
paretoMin, paretoMax, "pareto.txt");
final double loUniform = -9.876;
final double hiUniform = 5.432;
app.createDensity(new ContinuousUniformSampler(rng, loUniform, hiUniform),
loUniform, hiUniform, "uniform.txt");
final double scaleLogNormal = 2.345;
final double shapeLogNormal = 0.1234;
final double logNormalMin = 5;
final double logNormalMax = 25;
app.createDensity(new LogNormalSampler(new ZigguratNormalizedGaussianSampler(rng),
scaleLogNormal, shapeLogNormal),
logNormalMin, logNormalMax, "lognormal.ziggurat.txt");
app.createDensity(new LogNormalSampler(new MarsagliaNormalizedGaussianSampler(rng),
scaleLogNormal, shapeLogNormal),
logNormalMin, logNormalMax, "lognormal.marsaglia.txt");
app.createDensity(new LogNormalSampler(new BoxMullerNormalizedGaussianSampler(rng),
scaleLogNormal, shapeLogNormal),
logNormalMin, logNormalMax, "lognormal.boxmuller.txt");
}
}
© 2015 - 2024 Weber Informatics LLC | Privacy Policy