lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/digitization
diff -N BivariateDistribution.java
--- BivariateDistribution.java 29 Mar 2009 02:59:10 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,224 +0,0 @@
-/*
- * Calculate the probability integral for a set of bins in the x-y plane
- * of a bivariate normal distribution (i.e., a 2D Gaussian probability).
- *
- * The evaluation of the probability integrals is described in:
- *
- * Alan Genz, "Numerical Computation of Rectangular Bivariate and Trivariate
- * Normal and t Probabilities" in Statistics and Computing 14, 151 (2004).
- *
- * The integration code is adapted from the MATLAB source at:
- *
- * http://www.math.wsu.edu/faculty/genz/homepage
- *
- */
-package org.lcsim.contrib.Partridge.digitization;
-
-import org.apache.commons.math.MathException;
-import org.apache.commons.math.special.Erf;
-
-/**
- *
- * @author Richard Partridge
- */
-public class BivariateDistribution {
-
- private int _nx;
- private int _ny;
- private double _xmin;
- private double _ymin;
- private double _dx;
- private double _dy;
- private double[] _h;
- private double[] _k;
- private double[] _w6 = {0.1713244923791705, 0.3607615730481384, 0.4679139345726904};
- private double[] _x6 = {0.9324695142031522, 0.6612093864662647, 0.2386191860831970};
- private double[] _w12 = {.04717533638651177, 0.1069393259953183, 0.1600783285433464,
- 0.2031674267230659, 0.2334925365383547, 0.2491470458134029};
- private double[] _x12 = {0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
- 0.5873179542866171, 0.3678314989981802, 0.1252334085114692};
- private double[] _w20 = {.01761400713915212, .04060142980038694, .06267204833410906,
- .08327674157670475, 0.1019301198172404, 0.1181945319615184,
- 0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
- 0.1527533871307259};
- private double[] _x20 = {0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
- 0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
- 0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
- 0.07652652113349733};
-
- public void xBins(int nx, double xmin, double dx) {
- _nx = nx;
- _xmin = xmin;
- _dx = dx;
- _h = new double[_nx];
- }
-
- public void yBins(int ny, double ymin, double dy) {
- _ny = ny;
- _ymin = ymin;
- _dy = dy;
- _k = new double[_ny];
- }
-
- public double[][] Calculate(double x0, double y0, double sigx, double sigy,
- double rho) throws MathException {
-
- // Calculate the scaled x coordinate for each bin edge
- for (int i = 0; i < _nx + 1; i++) {
- _h[i] = (_xmin + i * _dx - x0) / sigx;
- }
-
- // Calculate the scaled y coordinate for each bin edge
- for (int j = 0; j < _ny + 1; j++) {
- _k[j] = (_ymin + j * _dx - y0) / sigy;
- }
-
- // Create the array that will hold the binned probabilities
- double[][] bi = new double[_nx][_ny];
-
- // Loop over the bin vertices
- for (int i = 0; i < _nx + 1; i++) {
- for (int j = 0; j < _ny + 1; j++) {
-
- // Calculate the probability for x>h and y>k for this vertex
- double prob = GenzCalc(_h[i], _k[j], rho);
-
- // Add or subtract this probability from the affected bins.
- // The bin probability for bin (0,0) is the sum of the Genz
- // probabilities for the (0,0) and (1,1) vertices MINUS the
- // sum of the probabilities for the (0,1) and (1,0) vertices
- if (i < _nx && j < _ny) {
- bi[i][j] += prob;
- }
- if (i > 0 && j > 0) {
- bi[i - 1][j - 1] += prob;
- }
- if (i > 0 && j < _ny) {
- bi[i - 1][j] -= prob;
- }
- if (i < _nx && j > 0) {
- bi[i][j - 1] -= prob;
- }
- }
- }
-
- return bi;
- }
-
- private double GenzCalc(double dh, double dk, double rho) throws MathException {
-
- double twopi = 2. * Math.PI;
-
- // Declare the Gauss-Legendre constants
- int ng;
- double[] w;
- double[] x;
-
- if (Math.abs(rho) < 0.3) {
- // for rho < 0.3 use 6 point Gauss-Legendre integration
- ng = 3;
- w = _w6;
- x = _x6;
- } else if (Math.abs(rho) < 0.75) {
- // for 0.3 < rho < 0.75 use 12 point Gauss-Legendre integration
- ng = 6;
- w = _w12;
- x = _x12;
- } else {
- // for rho > 0.75 use 20 point Gauss-Legendre integration
- ng = 10;
- w = _w20;
- x = _x20;
- }
-
- // Initialize the probability and some local variables
- double bvn = 0.;
- double h = dh;
- double k = dk;
- double hk = h * k;
-
- // For rho < 0.925, integrate equation 3 in Genz paper
- if (Math.abs(rho) < 0.925) {
- double hs = (h * h + k * k) / 2.;
- double asr = Math.asin(rho);
- double sn;
- for (int i = 0; i < ng; i++) {
- sn = Math.sin(asr * (1 - x[i]) / 2.);
- bvn += w[i] * Math.exp((sn * hk - hs) / (1 - sn * sn));
- sn = Math.sin(asr * (1 + x[i]) / 2.);
- bvn += w[i] * Math.exp((sn * hk - hs) / (1 - sn * sn));
- }
- bvn = bvn * asr / (2. * twopi) + phid(-h) * phid(-k);
-
- } else {
- // rho > 0.925 - integrate equation 6 in Genz paper
- if (rho < 0.) {
- k = -k;
- hk = -hk;
- }
- if (Math.abs(rho) < 1.) {
- double as = (1 - rho) * (1 + rho);
- double a = Math.sqrt(as);
- double bs = (h - k) * (h - k);
- double c = (4. - hk) / 8.;
- double d = (12. - hk) / 16.;
- double asr = -(bs / as + hk) / 2.;
- if (asr > -100.) {
- bvn = a * Math.exp(asr) *
- (1. - c * (bs - as) * (1. - d * bs / 5.) / 3. +
- c * d * as * as / 5.);
- }
-
- if (-hk < 100.) {
- double b = Math.sqrt(bs);
- bvn -= Math.exp(-hk / 2.) * Math.sqrt(twopi) * phid(-b / a) *
- b * (1 - c * bs * (1 - d * bs / 5.) / 3.);
- }
- a = a / 2.;
- for (int i = 0; i < ng; i++) {
- for (int j = 0; j < 2; j++) {
- int is = -1;
- if (j > 0) {
- is = 1;
- }
- double xs = Math.pow(a * (is * x[i] + 1), 2);
- double rs = Math.sqrt(1 - xs);
- asr = -(bs / xs + hk) / 2;
- if (asr > -100) {
- double sp = (1 + c * xs * (1 + d * xs));
- double ep = Math.exp(-hk * (1 - rs) / (2 * (1 + rs))) / rs;
- bvn += a * w[i] * Math.exp(asr) * (ep - sp);
- }
- }
- }
- bvn = -bvn / twopi;
- }
- if (rho > 0) {
- bvn = bvn + phid(-Math.max(h, k));
- } else {
- bvn = -bvn;
- if (k > h) {
- bvn += phid(k) - phid(h);
- }
- }
- }
- //
- return Math.max(0, Math.min(1, bvn));
- }
-
- private double phid(double arg) {
-
- // Use the Apache error function routine for now...catch it's inability
- // to handle all cases and potentially iterate excessively
- double erfc;
- try {
- erfc = (1. - Erf.erf(-arg / Math.sqrt(2)));
- } catch (MathException e) {
- // If we ever get here, we should find a better error function
- System.out.println("Apache error function appears to be poorly coded");
- System.out.println("Non-convergent error function for arg = " + arg);
- erfc = 0.;
- }
- return 0.5 * erfc;
- }
-}