1
2
3
4 package org.lcsim.math.probability;
5
6 import junit.framework.TestCase;
7
8
9
10
11
12
13 public class BivariateDistributionTest extends TestCase {
14
15
16 public void testBivariateDistribution() {
17
18
19 BivariateDistribution b = new BivariateDistribution();
20
21
22 int nx = 140;
23 double dx = 0.1;
24 double xmin = -0.5 * nx * dx;
25 b.xBins(nx, xmin, dx);
26
27
28 int ny = 140;
29 double dy = 0.1;
30 double ymin = -0.5 * ny * dy;
31 b.yBins(ny, ymin, dy);
32
33
34 double x0 = 0.526;
35 double y0 = -0.317;
36 double sigx0 = 0.642;
37 double sigy0 = 0.784;
38 double rho0 = 0.231;
39
40
41 double[][] bi = b.Calculate(x0, y0, sigx0, sigy0, rho0);
42
43
44 double xave = 0.;
45 double yave = 0.;
46 double xysum = 0.;
47 double xxsum = 0.;
48 double yysum = 0.;
49 double psum = 0.;
50 for (int i = 0; i < nx; i++) {
51 for (int j = 0; j < ny; j++) {
52 double x = xmin + dx * (i + 0.5);
53 double y = ymin + dy * (j + 0.5);
54 double prob = bi[i][j];
55 xave += prob * x;
56 yave += prob * y;
57 xxsum += prob * x * x;
58 yysum += prob * y * y;
59 xysum += prob * x * y;
60 psum += prob;
61 }
62 }
63
64
65 double sigx = Math.sqrt(xxsum - xave * xave);
66 double sigy = Math.sqrt(yysum - yave * yave);
67 double rho = (xysum - xave * yave) / (sigx * sigy);
68
69 System.out.println(" x ave: " + xave + " y ave: " + yave);
70 System.out.println(" x sd: " + sigx + " y sd: " + sigy + " rho: " + rho);
71 System.out.println("PSum: " + psum);
72
73
74
75
76 assertEquals("Probability sum failure", 1.0, psum, 5.0e-16*Math.sqrt(nx*ny));
77
78
79
80 assertEquals("x ave failure", x0, xave, 1e-10);
81 assertEquals("y ave failure", y0, yave, 1e-10);
82 assertEquals("sig x failure", sigx0, sigx, 2e-3);
83 assertEquals("sig y failure", sigy0, sigy, 2e-3);
84 assertEquals("rho failure", rho0, rho, 2e-3);
85
86 }
87 }