use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.
the class EMGainAnalysis method simulateFromPoissonGammaGaussian.
/**
* Randomly generate a histogram from poisson-gamma-gaussian samples
*
* @return The histogram
*/
private int[] simulateFromPoissonGammaGaussian() {
// Randomly sample
RandomGenerator random = new Well44497b(System.currentTimeMillis() + System.identityHashCode(this));
PoissonDistribution poisson = new PoissonDistribution(random, _photons, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
CustomGammaDistribution gamma = new CustomGammaDistribution(random, _photons, _gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
final int steps = simulationSize;
int[] sample = new int[steps];
for (int n = 0; n < steps; n++) {
if (n % 64 == 0)
IJ.showProgress(n, steps);
// Poisson
double d = poisson.sample();
// Gamma
if (d > 0) {
gamma.setShapeUnsafe(d);
d = gamma.sample();
}
// Gaussian
d += _noise * random.nextGaussian();
// Convert the sample to a count
sample[n] = (int) Math.round(d + _bias);
}
int max = Maths.max(sample);
int[] h = new int[max + 1];
for (int s : sample) h[s]++;
return h;
}
use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.
the class DiffusionRateTest method plotJumpDistances.
/**
* Plot a cumulative histogram and standard histogram of the jump distances.
*
* @param title
* the title
* @param jumpDistances
* the jump distances
* @param dimensions
* the number of dimensions for the jumps
* @param steps
* the steps
*/
private void plotJumpDistances(String title, DoubleData jumpDistances, int dimensions) {
// Cumulative histogram
// --------------------
double[] values = jumpDistances.values();
String title2 = title + " Cumulative Jump Distance " + dimensions + "D";
double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(values);
Plot2 jdPlot = new Plot2(title2, "Distance (um^2)", "Cumulative Probability", jdHistogram[0], jdHistogram[1]);
PlotWindow pw2 = Utils.display(title2, jdPlot);
if (Utils.isNewWindow())
idList[idCount++] = pw2.getImagePlus().getID();
// Plot the expected function
// This is the Chi-squared distribution: The sum of the squares of k independent
// standard normal random variables with k = dimensions. It is a special case of
// the gamma distribution. If the normals have non-unit variance the distribution
// is scaled.
// Chi ~ Gamma(k/2, 2) // using the scale parameterisation of the gamma
// s^2 * Chi ~ Gamma(k/2, 2*s^2)
// So if s^2 = 2D:
// 2D * Chi ~ Gamma(k/2, 4D)
double estimatedD = simpleD * simpleSteps;
double max = Maths.max(values);
double[] x = Utils.newArray(1000, 0, max / 1000);
double k = dimensions / 2.0;
double mean = 4 * estimatedD;
GammaDistribution dist = new GammaDistribution(k, mean);
double[] y = new double[x.length];
for (int i = 0; i < x.length; i++) y[i] = dist.cumulativeProbability(x[i]);
jdPlot.setColor(Color.red);
jdPlot.addPoints(x, y, Plot.LINE);
Utils.display(title2, jdPlot);
// Histogram
// ---------
title2 = title + " Jump " + dimensions + "D";
int plotId = Utils.showHistogram(title2, jumpDistances, "Distance (um^2)", 0, 0, Math.max(20, values.length / 1000));
if (Utils.isNewWindow())
idList[idCount++] = plotId;
// Recompute the expected function
for (int i = 0; i < x.length; i++) y[i] = dist.density(x[i]);
// Scale to have the same area
if (Utils.xValues.length > 1) {
final double area1 = jumpDistances.size() * (Utils.xValues[1] - Utils.xValues[0]);
final double area2 = dist.cumulativeProbability(x[x.length - 1]);
final double scaleFactor = area1 / area2;
for (int i = 0; i < y.length; i++) y[i] *= scaleFactor;
}
jdPlot = Utils.plot;
jdPlot.setColor(Color.red);
jdPlot.addPoints(x, y, Plot.LINE);
Utils.display(WindowManager.getImage(plotId).getTitle(), jdPlot);
}
use of org.apache.commons.math3.special.Gamma in project ffx by mjschnie.
the class Crystal method updateCrystal.
/**
* Update all Crystal variables that are a function of unit cell parameters.
*/
private void updateCrystal() {
switch(crystalSystem) {
case CUBIC:
case ORTHORHOMBIC:
case TETRAGONAL:
cos_alpha = 0.0;
sin_beta = 1.0;
cos_beta = 0.0;
sin_gamma = 1.0;
cos_gamma = 0.0;
beta_term = 0.0;
gamma_term = 1.0;
volume = a * b * c;
dVdA = b * c;
dVdB = a * c;
dVdC = a * b;
dVdAlpha = 0.0;
dVdBeta = 0.0;
dVdGamma = 0.0;
break;
case MONOCLINIC:
cos_alpha = 0.0;
sin_beta = sin(toRadians(beta));
cos_beta = cos(toRadians(beta));
sin_gamma = 1.0;
cos_gamma = 0.0;
beta_term = 0.0;
gamma_term = sin_beta;
volume = sin_beta * a * b * c;
dVdA = sin_beta * b * c;
dVdB = sin_beta * a * c;
dVdC = sin_beta * a * b;
dVdAlpha = 0.0;
dVdBeta = cos_beta * a * b * c;
dVdGamma = 0.0;
break;
case HEXAGONAL:
case TRICLINIC:
case TRIGONAL:
default:
double sin_alpha = sin(toRadians(alpha));
cos_alpha = cos(toRadians(alpha));
sin_beta = sin(toRadians(beta));
cos_beta = cos(toRadians(beta));
sin_gamma = sin(toRadians(gamma));
cos_gamma = cos(toRadians(gamma));
beta_term = (cos_alpha - cos_beta * cos_gamma) / sin_gamma;
gamma_term = sqrt(sin_beta * sin_beta - beta_term * beta_term);
volume = sin_gamma * gamma_term * a * b * c;
dVdA = sin_gamma * gamma_term * b * c;
dVdB = sin_gamma * gamma_term * a * c;
dVdC = sin_gamma * gamma_term * a * b;
double dbeta = 2.0 * sin_beta * cos_beta - (2.0 * (cos_alpha - cos_beta * cos_gamma) * sin_beta * cos_gamma) / (sin_gamma * sin_gamma);
double dgamma1 = -2.0 * (cos_alpha - cos_beta * cos_gamma) * cos_beta / sin_gamma;
double dgamma2 = cos_alpha - cos_beta * cos_gamma;
dgamma2 *= dgamma2 * 2.0 * cos_gamma / (sin_gamma * sin_gamma * sin_gamma);
dVdAlpha = (cos_alpha - cos_beta * cos_gamma) * sin_alpha / (sin_gamma * gamma_term) * a * b * c;
dVdBeta = 0.5 * sin_gamma * dbeta / gamma_term * a * b * c;
dVdGamma = (cos_gamma * gamma_term + 0.5 * sin_gamma * (dgamma1 + dgamma2) / gamma_term) * a * b * c;
break;
}
// a is the first row of A^(-1).
Ai[0][0] = a;
Ai[0][1] = 0.0;
Ai[0][2] = 0.0;
// b is the second row of A^(-1).
Ai[1][0] = b * cos_gamma;
Ai[1][1] = b * sin_gamma;
Ai[1][2] = 0.0;
// c is the third row of A^(-1).
Ai[2][0] = c * cos_beta;
Ai[2][1] = c * beta_term;
Ai[2][2] = c * gamma_term;
Ai00 = Ai[0][0];
Ai01 = Ai[0][1];
Ai02 = Ai[0][2];
Ai10 = Ai[1][0];
Ai11 = Ai[1][1];
Ai12 = Ai[1][2];
Ai20 = Ai[2][0];
Ai21 = Ai[2][1];
Ai22 = Ai[2][2];
// Invert A^-1 to get A
RealMatrix m = new Array2DRowRealMatrix(Ai, true);
m = new LUDecomposition(m).getSolver().getInverse();
A = m.getData();
// The columns of A are the reciprocal basis vectors
A00 = A[0][0];
A10 = A[1][0];
A20 = A[2][0];
A01 = A[0][1];
A11 = A[1][1];
A21 = A[2][1];
A02 = A[0][2];
A12 = A[1][2];
A22 = A[2][2];
// Reciprocal basis vector lengths
aStar = 1.0 / sqrt(A00 * A00 + A10 * A10 + A20 * A20);
bStar = 1.0 / sqrt(A01 * A01 + A11 * A11 + A21 * A21);
cStar = 1.0 / sqrt(A02 * A02 + A12 * A12 + A22 * A22);
if (logger.isLoggable(Level.FINEST)) {
logger.finest(format(" Reciprocal Lattice Lengths: (%8.3f, %8.3f, %8.3f)", aStar, bStar, cStar));
}
// Interfacial diameters from the dot product of the real and reciprocal vectors
interfacialRadiusA = (Ai00 * A00 + Ai01 * A10 + Ai02 * A20) * aStar;
interfacialRadiusB = (Ai10 * A01 + Ai11 * A11 + Ai12 * A21) * bStar;
interfacialRadiusC = (Ai20 * A02 + Ai21 * A12 + Ai22 * A22) * cStar;
// Divide by 2 to get radii.
interfacialRadiusA /= 2.0;
interfacialRadiusB /= 2.0;
interfacialRadiusC /= 2.0;
if (logger.isLoggable(Level.FINEST)) {
logger.finest(format(" Interfacial radii: (%8.3f, %8.3f, %8.3f)", interfacialRadiusA, interfacialRadiusB, interfacialRadiusC));
}
G[0][0] = a * a;
G[0][1] = a * b * cos_gamma;
G[0][2] = a * c * cos_beta;
G[1][0] = G[0][1];
G[1][1] = b * b;
G[1][2] = b * c * cos_alpha;
G[2][0] = G[0][2];
G[2][1] = G[1][2];
G[2][2] = c * c;
// invert G to yield Gstar
m = new Array2DRowRealMatrix(G, true);
m = new LUDecomposition(m).getSolver().getInverse();
Gstar = m.getData();
List<SymOp> symOps = spaceGroup.symOps;
int nSymm = symOps.size();
if (symOpsCartesian == null) {
symOpsCartesian = new ArrayList<>(nSymm);
} else {
symOpsCartesian.clear();
}
RealMatrix toFrac = new Array2DRowRealMatrix(A);
RealMatrix toCart = new Array2DRowRealMatrix(Ai);
for (int i = 0; i < nSymm; i++) {
SymOp symOp = symOps.get(i);
m = new Array2DRowRealMatrix(symOp.rot);
// rot_c = A^(-1).rot_f.A
RealMatrix rotMat = m.preMultiply(toCart.transpose()).multiply(toFrac.transpose());
// tr_c = tr_f.A^(-1)
double[] tr = toCart.preMultiply(symOp.tr);
symOpsCartesian.add(new SymOp(rotMat.getData(), tr));
}
}
use of org.apache.commons.math3.special.Gamma in project presto by prestodb.
the class MathFunctions method cauchyCdf.
@Description("Cauchy cdf for a given value, median, and scale (gamma)")
@ScalarFunction
@SqlType(StandardTypes.DOUBLE)
public static double cauchyCdf(@SqlType(StandardTypes.DOUBLE) double median, @SqlType(StandardTypes.DOUBLE) double scale, @SqlType(StandardTypes.DOUBLE) double value) {
checkCondition(scale > 0, INVALID_FUNCTION_ARGUMENT, "scale must be greater than 0");
CauchyDistribution distribution = new CauchyDistribution(null, median, scale, CauchyDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
return distribution.cumulativeProbability(value);
}
use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.
the class CameraModelAnalysis method convolveHistogram.
/**
* Convolve the histogram. The output is a discrete probability distribution.
*
* @param settings the settings
* @return The histogram
*/
private static double[][] convolveHistogram(CameraModelAnalysisSettings settings) {
// Find the range of the Poisson
final PoissonDistribution poisson = new PoissonDistribution(settings.getPhotons());
final int maxn = poisson.inverseCumulativeProbability(UPPER);
final double gain = getGain(settings);
final double noise = getReadNoise(settings);
final boolean debug = false;
// Build the Probability Mass/Density Function (PDF) of the distribution:
// either a Poisson (PMF) or Poisson-Gamma (PDF). The PDF is 0 at all
// values apart from the step interval.
// Note: The Poisson-Gamma is computed without the Dirac delta contribution
// at c=0. This allows correct convolution with the Gaussian of the dirac delta
// and the rest of the Poisson-Gamma (so matching the simulation).
final TDoubleArrayList list = new TDoubleArrayList();
double step;
String name;
int upsample = 100;
// Store the Dirac delta value at c=0. This must be convolved separately.
double dirac = 0;
// EM-CCD
if (settings.getMode() == MODE_EM_CCD) {
name = "Poisson-Gamma";
final double m = gain;
final double p = settings.getPhotons();
dirac = PoissonGammaFunction.dirac(p);
// Chose whether to compute a discrete PMF or a PDF using the approximation.
// Note: The delta function at c=0 is from the PMF of the Poisson. So it is
// a discrete contribution. This is omitted from the PDF and handled in
// a separate convolution.
// noise != 0;
final boolean discrete = false;
if (discrete) {
// Note: This is obsolete as the Poisson-Gamma function is continuous.
// Sampling it at integer intervals is not valid, especially for low gain.
// The Poisson-Gamma PDF should be integrated to form a discrete PMF.
step = 1.0;
double upper;
if (settings.getPhotons() < 20) {
upper = maxn;
} else {
// Approximate reasonable range of Poisson as a Gaussian
upper = settings.getPhotons() + 3 * Math.sqrt(settings.getPhotons());
}
final GammaDistribution gamma = new GammaDistribution(null, upper, gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
final int maxc = (int) gamma.inverseCumulativeProbability(0.999);
final int minn = Math.max(1, poisson.inverseCumulativeProbability(LOWER));
// See Ulbrich & Isacoff (2007). Nature Methods 4, 319-321, SI equation 3.
// Note this is not a convolution of a single Gamma distribution since the shape
// is modified not the count. So it is a convolution of a distribution made with
// a gamma of fixed count and variable shape.
// The count=0 is a special case.
list.add(PoissonGammaFunction.poissonGammaN(0, p, m));
final long total = (maxn - minn) * (long) maxc;
if (total < 1000000) {
// Full computation
// G(c) = sum n { (1 / n!) p^n e^-p (1 / ((n-1!)m^n)) c^n-1 e^-c/m }
// Compute as a log
// - log(n!) + n*log(p)-p -log((n-1)!) - n * log(m) + (n-1) * log(c) -c/m
// Note: Both methods work
LogFactorialCache lfc = LOG_FACTORIAL_CACHE.get().get();
if (lfc == null) {
lfc = new LogFactorialCache();
LOG_FACTORIAL_CACHE.set(new SoftReference<>(lfc));
}
lfc.ensureRange(minn - 1, maxn);
final double[] f = new double[maxn + 1];
final double logm = Math.log(m);
final double logp = Math.log(p);
for (int n = minn; n <= maxn; n++) {
f[n] = -lfc.getLogFactorialUnsafe(n) + n * logp - p - lfc.getLogFactorialUnsafe(n - 1) - n * logm;
}
// double total2 = total;
for (int c = 1; c <= maxc; c++) {
double sum = 0;
final double c_m = c / m;
final double logc = Math.log(c);
for (int n = minn; n <= maxn; n++) {
sum += StdMath.exp(f[n] + (n - 1) * logc - c_m);
}
// sum2 += pd[n] * gd[n].density(c);
list.add(sum);
// total += sum;
// This should match the approximation
// double approx = PoissonGammaFunction.poissonGamma(c, p, m);
// total2 += approx;
// System.out.printf("c=%d sum=%g approx=%g error=%g\n", c, sum2, approx,
// uk.ac.sussex.gdsc.core.utils.DoubleEquality.relativeError(sum2, approx));
}
// System.out.printf("sum=%g approx=%g error=%g\n", total, total2,
// uk.ac.sussex.gdsc.core.utils.DoubleEquality.relativeError(total, total2));
} else {
// Approximate
for (int c = 1; c <= maxc; c++) {
list.add(PoissonGammaFunction.poissonGammaN(c, p, m));
}
}
} else {
// This integrates the PDF using the approximation and up-samples together.
// Compute the sampling interval.
step = 1.0 / upsample;
// Reset
upsample = 1;
// Compute the integral of [-step/2:step/2] for each point.
// Use trapezoid integration.
final double step_2 = step / 2;
double prev = PoissonGammaFunction.poissonGammaN(0, p, m);
double next = PoissonGammaFunction.poissonGammaN(step_2, p, m);
list.add((prev + next) * 0.25);
double max = 0;
for (int i = 1; ; i++) {
// Each remaining point is modelling a PMF for the range [-step/2:step/2]
prev = next;
next = PoissonGammaFunction.poissonGammaN(i * step + step_2, p, m);
final double pp = (prev + next) * 0.5;
if (max < pp) {
max = pp;
}
if (pp / max < 1e-5) {
// Use this if non-zero since it has been calculated
if (pp != 0) {
list.add(pp);
}
break;
}
list.add(pp);
}
}
// Ensure the combined sum of PDF and Dirac is 1
final double expected = 1 - dirac;
// Require an odd number to get an even number (n) of sub-intervals:
if (list.size() % 2 == 0) {
list.add(0);
}
final double[] g = list.toArray();
// Number of sub intervals
final int n = g.length - 1;
// h = (a-b) / n = sub-interval width
final double h = 1;
double sum2 = 0;
double sum4 = 0;
for (int j = 1; j <= n / 2 - 1; j++) {
sum2 += g[2 * j];
}
for (int j = 1; j <= n / 2; j++) {
sum4 += g[2 * j - 1];
}
final double sum = (h / 3) * (g[0] + 2 * sum2 + 4 * sum4 + g[n]);
// Check
// System.out.printf("Sum=%g Expected=%g\n", sum * step, expected);
SimpleArrayUtils.multiply(g, expected / sum);
list.resetQuick();
list.add(g);
} else {
name = "Poisson";
// Apply fixed gain. Just change the step interval of the PMF.
step = gain;
for (int n = 0; n <= maxn; n++) {
list.add(poisson.probability(n));
}
final double p = poisson.probability(list.size());
if (p != 0) {
list.add(p);
}
}
// Debug
if (debug) {
final String title = name;
final Plot plot = new Plot(title, "x", "y");
plot.addPoints(SimpleArrayUtils.newArray(list.size(), 0, step), list.toArray(), Plot.LINE);
ImageJUtils.display(title, plot);
}
double zero = 0;
double[] pg = list.toArray();
// Sample Gaussian
if (noise > 0) {
step /= upsample;
pg = list.toArray();
// Convolve with Gaussian kernel
final double[] kernel = GaussianKernel.makeGaussianKernel(Math.abs(noise) / step, 6, true);
if (upsample != 1) {
// Use scaled convolution. This is faster that zero filling distribution g.
pg = Convolution.convolve(kernel, pg, upsample);
} else if (dirac > 0.01) {
// The Poisson-Gamma may be stepped at low mean causing wrap artifacts in the FFT.
// This is a problem if most of the probability is in the Dirac.
// Otherwise it can be ignored and the FFT version is OK.
pg = Convolution.convolve(kernel, pg);
} else {
pg = Convolution.convolveFast(kernel, pg);
}
// The convolution will have created a larger array so we must adjust the offset for this
final int radius = kernel.length / 2;
zero -= radius * step;
// Add convolution of the dirac delta function.
if (dirac != 0) {
// We only need to convolve the Gaussian at c=0
for (int i = 0; i < kernel.length; i++) {
pg[i] += kernel[i] * dirac;
}
}
// Debug
if (debug) {
String title = "Gaussian";
Plot plot = new Plot(title, "x", "y");
plot.addPoints(SimpleArrayUtils.newArray(kernel.length, -radius * step, step), kernel, Plot.LINE);
ImageJUtils.display(title, plot);
title = name + "-Gaussian";
plot = new Plot(title, "x", "y");
plot.addPoints(SimpleArrayUtils.newArray(pg.length, zero, step), pg, Plot.LINE);
ImageJUtils.display(title, plot);
}
zero = downSampleCdfToPmf(settings, list, step, zero, pg, 1.0);
pg = list.toArray();
zero = (int) Math.floor(zero);
step = 1.0;
// No convolution means we have the Poisson PMF/Poisson-Gamma PDF already
} else if (step != 1) {
// Sample to 1.0 pixel step interval.
if (settings.getMode() == MODE_EM_CCD) {
// Poisson-Gamma PDF
zero = downSampleCdfToPmf(settings, list, step, zero, pg, 1 - dirac);
pg = list.toArray();
zero = (int) Math.floor(zero);
// Add the dirac delta function.
if (dirac != 0) {
// Note: zero is the start of the x-axis. This value should be -1.
assert (int) zero == -1;
// Use as an offset to find the actual zero.
pg[-(int) zero] += dirac;
}
} else {
// Poisson PMF
// Simple non-interpolated expansion.
// This should be used when there is no Gaussian convolution.
final double[] pd = pg;
list.resetQuick();
// Account for rounding.
final Round round = getRound(settings);
final int maxc = round.round(pd.length * step + 1);
pg = new double[maxc];
for (int n = pd.length; n-- > 0; ) {
pg[round.round(n * step)] += pd[n];
}
if (pg[0] != 0) {
list.add(0);
list.add(pg);
pg = list.toArray();
zero--;
}
}
step = 1.0;
} else {
// Add the dirac delta function.
list.setQuick(0, list.getQuick(0) + dirac);
}
return new double[][] { SimpleArrayUtils.newArray(pg.length, zero, step), pg };
}
Aggregations