1
2
3
4 package org.lcsim.util.probability;
5
6 import org.lcsim.util.probability.Erf;
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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
35 private double[] _w6 = {0.1713244923791705, 0.3607615730481384, 0.4679139345726904};
36 private double[] _x6 = {0.9324695142031522, 0.6612093864662647, 0.2386191860831970};
37
38
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
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
56
57
58
59
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
70
71
72
73
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
84
85
86
87
88
89
90
91
92
93
94
95
96
97 public double[][] Calculate(double x0, double y0, double sigx, double sigy,
98 double rho) {
99
100
101 for (int i = 0; i < _nx + 1; i++) {
102 _h[i] = (_xmin + i * _dx - x0) / sigx;
103 }
104
105
106 for (int j = 0; j < _ny + 1; j++) {
107 _k[j] = (_ymin + j * _dy - y0) / sigy;
108 }
109
110
111 double[][] bi = new double[_nx][_ny];
112
113
114 for (int i = 0; i < _nx + 1; i++) {
115 for (int j = 0; j < _ny + 1; j++) {
116
117
118 double prob = GenzCalc(_h[i], _k[j], rho);
119
120
121
122
123
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
147 int ng;
148 double[] w;
149 double[] x;
150
151 if (Math.abs(rho) < 0.3) {
152
153 ng = 3;
154 w = _w6;
155 x = _x6;
156 } else if (Math.abs(rho) < 0.75) {
157
158 ng = 6;
159 w = _w12;
160 x = _x12;
161 } else {
162
163 ng = 10;
164 w = _w20;
165 x = _x20;
166 }
167
168
169 double bvn = 0.;
170 double h = dh;
171 double k = dk;
172 double hk = h * k;
173
174
175 if (Math.abs(rho) < 0.925) {
176
177
178
179
180
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
191
192 bvn = bvn * asr / (2. * twopi) + Erf.phi(-h) * Erf.phi(-k);
193
194 } else {
195
196
197
198
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 }