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