View Javadoc

1   /*
2    *  Class BivariateDistribution
3    */
4   package org.lcsim.util.probability;
5   
6   import org.lcsim.util.probability.Erf;
7   
8   /**
9    * Calculate the probability integral for a set of bins in the x-y plane
10   * of a bivariate normal distribution (i.e., a 2D Gaussian probability).
11   *<p>
12   * The evaluation of the probability integrals is described in:
13   *<p>
14   * Alan Genz, "Numerical Computation of Rectangular Bivariate and Trivariate
15   * Normal and t Probabilities" in Statistics and Computing 14, 151 (2004).
16   *<p>
17   * The integration code is adapted from the FORTRAN source at:
18   *<p>
19   *   http://www.math.wsu.edu/faculty/genz/homepage
20   *<p>
21   * @author Richard Partridge
22   */
23  public class BivariateDistribution {
24  
25      private int _nx;
26      private int _ny;
27      private double _xmin;
28      private double _ymin;
29      private double _dx;
30      private double _dy;
31      private double[] _h;
32      private double[] _k;
33  
34      //  Weights and coordinates for 6 point Gauss-Legendre integration
35      private double[] _w6 = {0.1713244923791705, 0.3607615730481384, 0.4679139345726904};
36      private double[] _x6 = {0.9324695142031522, 0.6612093864662647, 0.2386191860831970};
37  
38      //  Weights and coordinates for 12 point Gauss-Legendre integration
39      private double[] _w12 = {.04717533638651177, 0.1069393259953183, 0.1600783285433464,
40          0.2031674267230659, 0.2334925365383547, 0.2491470458134029};
41      private double[] _x12 = {0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
42          0.5873179542866171, 0.3678314989981802, 0.1252334085114692};
43  
44      //  Weights and coordinates for 20 point Gauss-Legendre integration
45      private double[] _w20 = {.01761400713915212, .04060142980038694, .06267204833410906,
46          .08327674157670475, 0.1019301198172404, 0.1181945319615184,
47          0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
48          0.1527533871307259};
49      private double[] _x20 = {0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
50          0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
51          0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
52          0.07652652113349733};
53  
54      /**
55       * Set the locations of the x-coordinate bins
56       *
57       * @param nx number of x coordinate bins
58       * @param xmin minimum x coordinate
59       * @param dx width of x coordinate bins
60       */
61      public void xBins(int nx, double xmin, double dx) {
62          _nx = nx;
63          _xmin = xmin;
64          _dx = dx;
65          _h = new double[_nx + 1];
66      }
67  
68      /**
69       * Set the locations of the y-coordinate bins
70       *
71       * @param ny number of y coordinate bins
72       * @param ymin minimum y coordinate
73       * @param dy width of y coordinate bins
74       */
75      public void yBins(int ny, double ymin, double dy) {
76          _ny = ny;
77          _ymin = ymin;
78          _dy = dy;
79          _k = new double[_ny + 1];
80      }
81  
82      /**
83       * Integrate the Gaussian probability distribution over each x-y bins,
84       * which must be defined before calling this method.
85       * <p>
86       * The output is a double array that gives the binned probability
87       * distribution.  The first array index is used to indicate the bin in x
88       * and the second array index is used to indicate the bin in y.
89       * <p>
90       * @param x0 mean x coordinate of Gaussian distribution
91       * @param y0 mean y coordinate of Gaussian distribution
92       * @param sigx x coordinate standard deviation
93       * @param sigy y coordinate standard deviation
94       * @param rho x-y correlation coefficient
95       * @return probability distribution
96       */
97      public double[][] Calculate(double x0, double y0, double sigx, double sigy,
98              double rho) {
99  
100         //  Calculate the scaled x coordinate for each bin edge
101         for (int i = 0; i < _nx + 1; i++) {
102             _h[i] = (_xmin + i * _dx - x0) / sigx;
103         }
104 
105         //  Calculate the scaled y coordinate for each bin edge
106         for (int j = 0; j < _ny + 1; j++) {
107             _k[j] = (_ymin + j * _dy - y0) / sigy;
108         }
109 
110         //  Create the array that will hold the binned probabilities
111         double[][] bi = new double[_nx][_ny];
112 
113         //  Loop over the bin vertices
114         for (int i = 0; i < _nx + 1; i++) {
115             for (int j = 0; j < _ny + 1; j++) {
116 
117                 //  Calculate the probability for x>h and y>k for this vertex
118                 double prob = GenzCalc(_h[i], _k[j], rho);
119 
120                 //  Add or subtract this probability from the affected bins.
121                 //  The bin probability for bin (0,0) is the sum of the Genz
122                 //  probabilities for the (0,0) and (1,1) vertices MINUS the
123                 //  sum of the probabilities for the (0,1) and (1,0) vertices
124                 if (i > 0 && j > 0) {
125                     bi[i - 1][j - 1] += prob;
126                 }
127                 if (i > 0 && j < _ny) {
128                     bi[i - 1][j] -= prob;
129                 }
130                 if (i < _nx && j > 0) {
131                     bi[i][j - 1] -= prob;
132                 }
133                 if (i < _nx && j < _ny) {
134                     bi[i][j] += prob;
135                 }
136             }
137         }
138 
139         return bi;
140     }
141 
142     private double GenzCalc(double dh, double dk, double rho) {
143 
144         double twopi = 2. * Math.PI;
145 
146         //  Declare the Gauss-Legendre constants
147         int ng;
148         double[] w;
149         double[] x;
150 
151         if (Math.abs(rho) < 0.3) {
152             //  for rho < 0.3 use 6 point Gauss-Legendre integration
153             ng = 3;
154             w = _w6;
155             x = _x6;
156         } else if (Math.abs(rho) < 0.75) {
157             //  for 0.3 < rho < 0.75 use 12 point Gauss-Legendre integration
158             ng = 6;
159             w = _w12;
160             x = _x12;
161         } else {
162             //  for rho > 0.75 use 20 point Gauss-Legendre integration
163             ng = 10;
164             w = _w20;
165             x = _x20;
166         }
167 
168         //  Initialize the probability and some local variables
169         double bvn = 0.;
170         double h = dh;
171         double k = dk;
172         double hk = h * k;
173 
174         //  For rho < 0.925, integrate equation 3 in the Genz paper
175         if (Math.abs(rho) < 0.925) {
176 
177             //  More or less direct port of Genz code follows
178             //  It is fairly easy to match this calculation against equation 3 of
179             //  Genz's paper if you take into account that you need to change
180             //  variables so the integration argument spans the range -1 to 1
181             double hs = (h * h + k * k) / 2.;
182             double asr = Math.asin(rho);
183             double sn;
184             for (int i = 0; i < ng; i++) {
185                 sn = Math.sin(asr * (1 - x[i]) / 2.);
186                 bvn += w[i] * Math.exp((sn * hk - hs) / (1 - sn * sn));
187                 sn = Math.sin(asr * (1 + x[i]) / 2.);
188                 bvn += w[i] * Math.exp((sn * hk - hs) / (1 - sn * sn));
189             }
190             //  The factor of asr/2 comes from changing variables so the
191             //  integration is over the range -1 to 1 instead of 0 - asin(rho)
192             bvn = bvn * asr / (2. * twopi) + Erf.phi(-h) * Erf.phi(-k);
193 
194         } else {
195             //  rho > 0.925 - integrate equation 6 in Genz paper with the
196             //  extra term in the Taylor expansion given in equation 7.
197             //  The rest of this code is pretty dense and is a pretty direct
198             //  port of Genz's code.
199 
200             if (rho < 0.) {
201                 k = -k;
202                 hk = -hk;
203             }
204 
205             if (Math.abs(rho) < 1.) {
206 
207                 double as = (1 - rho) * (1 + rho);
208                 double a = Math.sqrt(as);
209                 double bs = (h - k) * (h - k);
210                 double c = (4. - hk) / 8.;
211                 double d = (12. - hk) / 16.;
212                 double asr = -(bs / as + hk) / 2.;
213 
214                 if (asr > -100.) {
215                     bvn = a * Math.exp(asr) *
216                             (1. - c * (bs - as) * (1. - d * bs / 5.) / 3. +
217                             c * d * as * as / 5.);
218                 }
219 
220                 if (-hk < 100.) {
221                     double b = Math.sqrt(bs);
222                     bvn -= Math.exp(-hk / 2.) * Math.sqrt(twopi) * Erf.phi(-b / a) *
223                             b * (1 - c * bs * (1 - d * bs / 5.) / 3.);
224                 }
225 
226                 a = a / 2.;
227                 for (int i = 0; i < ng; i++) {
228                     for (int j = 0; j < 2; j++) {
229                         int is = -1;
230                         if (j > 0) {
231                             is = 1;
232                         }
233                         double xs = Math.pow(a * (is * x[i] + 1), 2);
234                         double rs = Math.sqrt(1 - xs);
235                         asr = -(bs / xs + hk) / 2;
236 
237                         if (asr > -100) {
238                             double sp = (1 + c * xs * (1 + d * xs));
239                             double ep = Math.exp(-hk * (1 - rs) / (2 * (1 + rs))) / rs;
240                             bvn += a * w[i] * Math.exp(asr) * (ep - sp);
241                         }
242                     }
243                 }
244 
245                 bvn = -bvn / twopi;
246             }
247 
248             if (rho > 0) {
249                 bvn = bvn + Erf.phi(-Math.max(h, k));
250             } else {
251                 bvn = -bvn;
252                 if (k > h) {
253                     bvn += Erf.phi(k) - Erf.phi(h);
254                 }
255             }
256         }
257          
258         return Math.max(0, Math.min(1, bvn));
259     }
260 }