View Javadoc

1   /*
2    * BivariateDistributionTest class
3    */
4   package org.lcsim.math.probability;
5   
6   import junit.framework.TestCase;
7   
8   /**
9    * Test case for BivariateDistribution class
10   *
11   * @author Richard Partridge
12   */
13  public class BivariateDistributionTest extends TestCase {
14  
15      /** Creates a new instance of HelicalTrackFitterTest */
16      public void testBivariateDistribution() {
17  
18          //  Instantiate the BivariateDistribution class
19          BivariateDistribution b = new BivariateDistribution();
20  
21          //  Set up the x coordinate binning
22          int nx = 140;
23          double dx = 0.1;
24          double xmin = -0.5 * nx * dx;
25          b.xBins(nx, xmin, dx);
26  
27          //  Set up the y coordinate binning
28          int ny = 140;
29          double dy = 0.1;
30          double ymin = -0.5 * ny * dy;
31          b.yBins(ny, ymin, dy);
32  
33          //  Set the bivariate Gaussian parameters to some semi-random values
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          //  Calculate the bivariate probabilities for our x-y bins
41          double[][] bi = b.Calculate(x0, y0, sigx0, sigy0, rho0);
42  
43          //  Now calculate our parameter estimates from the binned data
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          //  Calculate the measured error matrix
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          //  Test that probability is conserved - this is the key test that
74          //  the method is working.  Estimate a 5 sigma round-off error from
75          //  summing nx*ny bins assuming 1e-16 precision per bin.
76          assertEquals("Probability sum failure", 1.0, psum, 5.0e-16*Math.sqrt(nx*ny));
77  
78          //  Make some crude tests that the Gaussian parameters are reasonable
79          //  These are not precisely measured due to the coarse binning
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  }