 |
 |
| Commit in lcsim/src/org/lcsim/math on MAIN | |||
| chisq/ChisqProb.java | -247 | 1.4 removed | |
| coordinatetransform/CovarianceMatrixTransformer.java | -102 | 1.1 removed | |
| distribution/LandauDistribution.java | -137 | 1.1 removed | |
| /MoyalDistribution.java | -25 | 1.2 removed | |
| /RandLandau.java | -582 | 1.1 removed | |
| interpolation/BilinearInterpolator.java | -195 | 1.2 removed | |
| moments/CentralMomentsCalculator.java | -174 | 1.1 removed | |
| probability/BivariateDistribution.java | -258 | 1.1 removed | |
| /Erf.java | -283 | 1.1 removed | |
| -2003 | |||
remove classes from org.lcsim.math; moved to lcsim-math pkg
diff -N ChisqProb.java --- ChisqProb.java 6 Jul 2007 22:06:28 -0000 1.4 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,247 +0,0 @@
-package org.lcsim.math.chisq;
-/** A utility class to return the probability that a value of chi-squared <b> x </b> measured
- * for a system with <b> n </b>degrees of freedom would be exceeded by chance. Based on routines
- * in the book <em> "Numerical Recipes: The Art of Scientific Computing"</em>.
- * @author Norman Graf
- * @version $Id: ChisqProb.java,v 1.4 2007/07/06 22:06:28 ngraf Exp $
- *
- */
-
-public class ChisqProb
-{
-
- /**Returns the incomplete gamma function P(a, x).
- *
- * <br clear="all" /><table border="0" width="100%"><tr><td>
- * <table align="center"><tr><td nowrap="nowrap" align="center">
- * <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
- * > </td><td nowrap="nowrap" align="center">
- * <font face="symbol">g</font
- * >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
- * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
- * <font face="symbol">�</font
- * > </td><td nowrap="nowrap" align="center">
- * 1<hr noshade="noshade" /><font face="symbol">G</font
- * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * <font size="-1"><i>x</i></font><!--sup
- * --><br /><font face="symbol">�<br />�<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
- * <i>e</i><sup> <font face="symbol">-</font
- * > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
- * > 1</sup> <i>dt</i> </td></tr></table>
- * </td></tr></table>
- *
- * This is the probability that the observed chi-square for a correct
- * model should be less than a value of chi-square.
- * @param a the number of degrees of freedom
- * @param x the value of chi-squared
- * @return The incomplete gamma function P(a, x).
- */
- public static double gammp(double a, double x)
- {
- if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine gammp");
- if (x < (a+1.0))
- { //Use the series representation.
- return gser(a/2.,x/2.);
- }
- else
- { //Use the continued fraction representation and take its complement.
- return 1.0-gcf(a/2.,x/2.);
- }
- }
-
- /**Returns the incomplete gamma function Q(a, x)= 1 - P(a, x).
- *
- * <br clear="all" /><table border="0" width="100%"><tr><td>
- * <table align="center"><tr><td nowrap="nowrap" align="center">
- * <i>Q</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
- * > 1 <font face="symbol">-</font
- * > <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
- * > </td><td nowrap="nowrap" align="center">
- * <font face="symbol">G</font
- * >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
- * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * <font size="-1"><font face="symbol">�</font
- * ></font><!--sup
- * --><br /><font face="symbol">�<br />�<br /></font><font size="-1"><i>x</i></font> <br /></td><td nowrap="nowrap" align="center">
- * <i>e</i><sup> <font face="symbol">-</font
- * > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
- * > 1</sup> <i>dt</i> </td></tr></table>
- * </td></tr></table>
- *
- * This is the probability that the observed chi-square
- * will exceed the value of chi-square by chance even for a correct model.
- * @param a the number of degrees of freedom
- * @param x the value of chi-squared
- * @return The incomplete gamma function Q(a, x)= 1 - P(a, x).
- */
- public static double gammq(double a, double x)
- {
-
- if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine gammq");
- if (x < (a+1.0))
- { //Use the series representation and take its complement.
- return 1.0-gser(a/2.,x/2.);
- }
- else
- { //Use the continued fraction representation.
- return gcf(a/2.,x/2.);
- }
- }
-
- /**Returns the incomplete gamma function P(a, x) evaluated by its series representation.
- *
- * <br clear="all" /><table border="0" width="100%"><tr><td>
- * <table align="center"><tr><td nowrap="nowrap" align="center">
- * <font face="symbol">g</font
- * >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
- * > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td nowrap="nowrap" align="center">
- * <font size="-1"><font face="symbol">�</font
- * ></font><!--sup
- * --><br /><font size="+3"><font face="symbol">�<br />
- * </font></font><font size="-1"><i>n</i> = 0</font> <br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * <font face="symbol">G</font
- * >(<i>a</i>)
- * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><font face="symbol">G</font
- * >(<i>a</i> + 1 + <i>n</i>)<br /></td><td nowrap="nowrap" align="center">
- * <i>x</i><sup><i>n</i></sup></td></tr></table>
- * </td></tr></table>
- *
- * @param a the number of degrees of freedom/2
- * @param x the value of chi-squared/2
- * @return The incomplete gamma P(a, x) evaluated by its series representation.
- */
- public static double gser(double a, double x)
- {
- int ITMAX=100;
- double EPS=3.0e-7;
-
- int n;
- double sum,del,ap;
- double gln=gammln(a);
- if (x <= 0.0)
- {
- if (x < 0.0) System.out.println("x less than 0 in routine gser");
- return 0.0;
- }
- else
- {
- ap=a;
- del=sum=1.0/a;
- for (n=1;n<=ITMAX;n++)
- {
- ++ap;
- del *= x/ap;
- sum += del;
- if (Math.abs(del) < Math.abs(sum)*EPS)
- {
- return sum*Math.exp(-x+a*Math.log(x)-(gln));
- }
- }
- System.out.println("a too large, ITMAX too small in routine gser");
- return 0.;
- }
- }
-
- /**Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation.
- *
- * <br clear="all" /><table border="0" width="100%"><tr><td>
- * <table align="center"><tr><td nowrap="nowrap" align="center">
- * <font face="symbol">G</font
- * >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
- * > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td align="left" class="cl"><font face="symbol">
- * �<br />�
- * </font> </td><td nowrap="nowrap" align="center">
- * 1
- * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * 1 <font face="symbol">-</font
- * > <i>a</i>
- * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>1 + <br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * 1
- * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * 2 <font face="symbol">-</font
- * > <i>a</i>
- * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>1 + <br /></td><td nowrap="nowrap" align="center">
- * </td><td nowrap="nowrap" align="center">
- * 2
- * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
- * <font face="symbol">�</font
- * > </td><td align="left" class="cl"><font face="symbol">
- * �<br />�
- * </font></td><td nowrap="nowrap" align="center">
- * (<i>x</i> > 0)</td></tr></table>
- * </td></tr></table>
- *
- * @param a the number of degrees of freedom/2
- * @param x the value of chi-squared/2
- * @return The incomplete gamma Q(a, x) evaluated by its continued fraction representation.
- */
- public static double gcf(double a, double x)
- {
- int ITMAX=100; //Maximum allowed number of iteration
- double EPS=3.0e-7; //Relative accuracy.
- double FPMIN=1.0e-30; //Number near the smallest representable doubleing-point number.
- {
- int i;
- double an,b,c,d,del,h;
- double gln=gammln(a);
- b=x+1.0-a; //Set up for evaluating continued fraction by modified Lentz's method with b0 = 0.
- c=1.0/FPMIN;
- d=1.0/b;
- h=d;
- for (i=1;i<=ITMAX;i++)
- { //Iterate to convergence.
- an = -i*(i-a);
- b += 2.0;
- d=an*d+b;
- if (Math.abs(d) < FPMIN) d=FPMIN;
- c=b+an/c;
- if (Math.abs(c) < FPMIN) c=FPMIN;
- d=1.0/d;
- del=d*c;
- h *= del;
- if (Math.abs(del-1.0) < EPS) break;
- }
- if (i > ITMAX) System.out.println("a too large, ITMAX too small in gcf");
- return Math.exp(-x+a*Math.log(x)-(gln))*h; //Put factors in front.
- }
- }
-
- /**Returns the logarithm of the Gamma function of x, for x>0.
- *
- * <br clear="all" /><table border="0" width="100%"><tr><td>
- * <table align="center"><tr><td nowrap="nowrap" align="center">
- * <font face="symbol">G</font
- * >(<i>z</i>) = </td><td nowrap="nowrap" align="center">
- * <font size="-1"><font face="symbol">�</font
- * ></font><!--sup
- * --><br /><font face="symbol">�<br />�<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
- * <i>t</i><sup><i>z</i> <font face="symbol">-</font
- * > 1</sup> <i>e</i><sup> <font face="symbol">-</font
- * > <i>t</i></sup> <i>dt</i> </td></tr></table>
- * </td></tr></table>
- *
- *@param x
- *@return The value ln(Gamma(x))
- */
- public static double gammln(double x)
- {
- double y,tmp,ser;
- double[] cof=
- {76.18009172947146,-86.50532032941677,
- 24.01409824083091,-1.231739572450155,
- 0.1208650973866179e-2,-0.5395239384953e-5};
- int j;
- y=x;
- tmp=x+5.5;
- tmp -= (x+0.5)*Math.log(tmp);
- ser=1.000000000190015;
- for (j=0;j<=5;j++) ser += cof[j]/++y;
- return -tmp+Math.log(2.5066282746310005*ser/x);
- }
-}
\ No newline at end of file
diff -N CovarianceMatrixTransformer.java --- CovarianceMatrixTransformer.java 31 Mar 2006 01:08:36 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,102 +0,0 @@
-/*
- * CovarianceMatrixTransformer.java
- *
- * Created on March 30, 2006, 3:35 PM
- *
- * $Id: CovarianceMatrixTransformer.java,v 1.1 2006/03/31 01:08:36 ngraf Exp $
- */
-
-package org.lcsim.math.coordinatetransform;
-import static java.lang.Math.sqrt;
-import static java.lang.Math.cos;
-import static java.lang.Math.sin;
-
-/**
- * Utility class to handle transformations of covariance matrices
- * Inspired by Nick Sinev's CMTransform.java.
- * @author Norman Graf
- */
-public class CovarianceMatrixTransformer
-{
-
- /** No need to create, all methods static*/
- private CovarianceMatrixTransformer()
- {
- }
-
- /**
- * Convert covariance matrix elements in cartesian coordinates (x,y) to cylindrical (r,phi).
- * @param x The cartesian x coordinate.
- * @param y The cartesian y coordinate.
- * @param sxx The covariance matrix term for x at (x,y).
- * @param syy The covariance matrix term for y at (x,y).
- * @param sxy The covariance matrix term for xy at (x,y).
- * @return The packed r-phi covariance matrix terms:
- * <table>
- * <tr><td> cov[0] </td><td> r-r </td><tr>
- * <tr><td> cov[1] </td><td> phi-phi </td><tr>
- * <tr><td> cov[2] </td><td> r-phi</td><tr>
- * </table>
- */
- public static double[] xy2rphi(double x, double y, double sxx, double syy, double sxy)
- {
- double[] cov = new double[3];
- double xx = x*x;
- double xy = x*y;
- double yy = y*y;
- double rr = xx+yy;
- double r = sqrt(rr);
- double oneOverR2 = 1/rr;
-
- // srr
- cov[0] = oneOverR2*(sxx*xx + syy*yy + 2.*sxy*xy);
-
- // sphiphi
- cov[1] = oneOverR2*oneOverR2*(sxx*yy - 2.*sxy*xy + syy*xx);
-
- // srphi
- cov[2] = oneOverR2*(-sxx*xy + sxy*(xx-yy) + syy*xy)/r;
-
- return cov;
- }
-
- /**
- *Convert covariance matrix elements in cylindrical (r,phi) to cartesian coordinates (x,y).
- * @param r The cylindrical radius.
- * @param phi The cylindrical angle.
- * @param srr The covariance matrix term for r at (r, phi).
- * @param sff The covariance matrix term for phi at (r, phi).
- * @param srf The covariance matrix term for r-phi at (r, phi).
- * @return The packed x-y covariance matrix terms:
- * <table>
- * <tr><td> cov[0] </td><td> x-x </td><tr>
- * <tr><td> cov[1] </td><td> y-y </td><tr>
- * <tr><td> cov[2] </td><td> x-y</td><tr>
- * </table>
- */
-
- public static double[] rphi2xy(double r, double phi, double srr, double sff, double srf)
- {
- double[] cov = new double[3];
- // cosine^2(phi)
- double cf = cos(phi);
- double cc = cf*cf;
- //sine^2(phi)
- double sf = sin(phi);
- double ss = sf*sf;
-
- // cosine(phi)*sine(phi)
- double cs = cf*sf;
-
- // sxx
- cov[0] = srr*cc - 2.*srf*r*cs + sff*r*r*ss;
-
- // syy
- cov[1] = srr*ss + 2.*srf*r*cs + sff*r*r*cc;
-
- // sxy
- cov[2] = srr*cs + srf*r*(cc - ss) - sff*r*r*cs;
-
- return cov;
- }
-}
diff -N LandauDistribution.java --- LandauDistribution.java 18 Jul 2005 07:12:46 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,137 +0,0 @@
-/*
- * LandauDistribution.java
- * This corresponds to the CLHEP version, which
- * is based on the CERNLIB routine denlan (G110).
- * Note that there is some variation in whether
- * to divide by sigma, and whether the peak
- * coincides with the most probable value.
- *
- * Created on May 23, 2005, 4:21 PM
- */
-
-package org.lcsim.math.distribution;
-
-import static java.lang.Math.sqrt;
-import static java.lang.Math.exp;
-import static java.lang.Math.log;
-/**
- *
- * @author ngraf
- */
-public class LandauDistribution
-{
- private double _peak;
- private double _width;
-
- /** Creates a new instance of LandauDistribution */
- public LandauDistribution(double peak, double width)
- {
- _peak = peak;
- _width = width;
- }
-
- public double peak()
- {
- return _peak;
- }
-
- public double width()
- {
- return _width;
- }
-
- public double value(double x)
- {
- double xs = _peak + 0.222782*_width;
- return denlan((x-xs)/_width)/_width;
- }
-
- private double denlan(double x)
- {
- /* Initialized data */
-
- /* System generated locals */
- double ret_val, r__1;
-
- /* Local variables */
- double u, v;
- v = x;
- if (v < -5.5)
- {
- u = exp(v + 1.);
- ret_val = exp(-1 / u) / sqrt(u) * .3989422803 * ((a1[0] + (a1[
- 1] + a1[2] * u) * u) * u + 1);
- }
- else if (v < -1.)
- {
- u = exp(-v - 1);
- ret_val = exp(-u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[
- 4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] +
- q1[4] * v) * v) * v) * v);
- }
- else if (v < 1.)
- {
- ret_val = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) *
- v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v)
- * v);
- }
- else if (v < 5.)
- {
- ret_val = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) *
- v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v)
- * v);
- }
- else if (v < 12.)
- {
- u = 1 / v;
- /* Computing 2nd power */
- r__1 = u;
- ret_val = r__1 * r__1 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u)
- * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] *
- u) * u) * u) * u);
- }
- else if (v < 50.)
- {
- u = 1 / v;
- /* Computing 2nd power */
- r__1 = u;
- ret_val = r__1 * r__1 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u)
- * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] *
- u) * u) * u) * u);
- }
- else if (v < 300.)
- {
- u = 1 / v;
- /* Computing 2nd power */
- r__1 = u;
- ret_val = r__1 * r__1 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u)
- * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] *
- u) * u) * u) * u);
- }
- else
- {
- u = 1 / (v - v * log(v) / (v + 1));
- /* Computing 2nd power */
- r__1 = u;
- ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + 1);
- }
- return ret_val;
-
- }
-
- static final double[] p1 = { .4259894875,-.124976255,.039842437,-.006298287635,.001511162253 };
- static final double[] q5 = { 1.,156.9424537,3745.310488,9834.698876,66924.28357 };
- static final double[] p6 = { 1.000827619,664.9143136,62972.92665,475554.6998,-5743609.109 };
- static final double[] q6 = { 1.,651.4101098,56974.73333,165917.4725,-2815759.939 };
- static final double[] a1 = { .04166666667,-.01996527778,.02709538966 };
- static final double[] a2 = { -1.84556867,-4.284640743 };
- static final double[] q1 = { 1.,-.3388260629,.09594393323,-.01608042283,.003778942063 };
- static final double[] p2 = { .1788541609,.1173957403,.01488850518,-.001394989411,1.283617211e-4 };
- static final double[] q2 = { 1.,.7428795082,.3153932961,.06694219548,.008790609714 };
- static final double[] p3 = { .1788544503,.09359161662,.006325387654,6.611667319e-5,-2.031049101e-6 };
- static final double[] q3 = { 1.,.6097809921,.2560616665,.04746722384,.006957301675 };
- static final double[] p4 = { .9874054407,118.6723273,849.279436,-743.7792444,427.0262186 };
- static final double[] q4 = { 1.,106.8615961,337.6496214,2016.712389,1597.063511 };
- static final double[] p5 = { 1.003675074,167.5702434,4789.711289,21217.86767,-22324.9491 };
-
-}
diff -N MoyalDistribution.java --- MoyalDistribution.java 28 Jun 2006 01:54:50 -0000 1.2 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,25 +0,0 @@
-package org.lcsim.math.distribution;
-
-import java.util.Random;
-
-public class MoyalDistribution
-{
-
- private String[] parameters = {"N", "x0", "sigma"};
- private double[] parValues = {1, 1, 1};
- private String title = "MoyalDistribution";
- private String[] variables = {"x"};
-
- public double value(double x)
- {
- double var = x;
- double n = parValues[0];
- double x0 = parValues[1];
- double s = parValues[2];
- double lambda = (var-x0)/s;
- double arg = lambda + Math.exp( -1.*lambda );
- double num = Math.exp( -1.*arg );
- return n*Math.sqrt( num/(2*Math.PI) );
- }
-
-}
diff -N RandLandau.java --- RandLandau.java 18 Jul 2005 07:12:46 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,582 +0,0 @@
-/*
- * RandLandau.java
- *
- * Created on June 6, 2005, 4:05 PM
- *
- * To change this template, choose Tools | Options and locate the template under
- * the Source Creation and Management node. Right-click the template and choose
- * Open. You can then make changes to the template in the Source Editor.
- */
-
-package org.lcsim.math.distribution;
-
-
-import static java.lang.Math.log;
-import java.util.Random;
-/**
- *
- * @author ngraf
- */
-public class RandLandau
-{
- private Random _ran = new Random();
- private double _sigma;
- private double _mean;
-
- /** Creates a new instance of RandLandau */
- public RandLandau()
- {
- _mean = 0.;
- _sigma = 1.0;
- }
-
- public RandLandau(double mean, double sigma)
- {
- _mean = mean;
- _sigma = sigma;
- }
-
-
- public double shoot()
- {
- return _mean+_sigma*transform(_ran.nextDouble());
- }
-
- private double transform(double r)
- {
-
- double u = r * TABLE_MULTIPLIER;
- int index = (int) u;
- double du = u - index;
-
- // du is scaled such that the we dont have to multiply by TABLE_INTERVAL
- // when interpolating.
-
- // Five cases:
- // A) Between .070 and .800 the function is so smooth, straight
- // linear interpolation is adequate.
- // B) Between .007 and .070, and between .800 and .980, quadratic
- // interpolation is used. This requires the same 4 points as
- // a cubic spline (thus we need .006 and .981 and .982) but
- // the quadratic interpolation is accurate enough and quicker.
- // C) Below .007 an asymptotic expansion for low negative lambda
- // (involving two logs) is used; there is a pade-style correction
- // factor.
- // D) Above .980, a simple pade approximation is made (asymptotic to
- // 1/(1-r)), but...
- // E) the coefficients in that pade are different above r=.999.
-
- if ( index >= 70 && index <= 800 )
- { // (A)
-
- double f0 = inverseLandau [index];
- double f1 = inverseLandau [index+1];
- return f0 + du * (f1 - f0);
-
- }
- else if ( index >= 7 && index <= 980 )
- { // (B)
-
- double f_1 = inverseLandau [index-1];
- double f0 = inverseLandau [index];
- double f1 = inverseLandau [index+1];
- double f2 = inverseLandau [index+2];
-
- return f0 + du * (f1 - f0 - .25*(1-du)* (f2 -f1 - f0 + f_1) );
-
- }
- else if ( index < 7 )
- { // (C)
-
- final double n0 = 0.99858950;
- final double n1 = 34.5213058; final double d1 = 34.1760202;
- final double n2 = 17.0854528; final double d2 = 4.01244582;
-
- double logr = log(r);
- double x = 1/logr;
- double x2 = x*x;
-
- double pade = (n0 + n1*x + n2*x2) / (1.0 + d1*x + d2*x2);
-
- return ( - log( -.91893853 - logr ) -1 ) * pade;
-
- }
- else if ( index <= 999 )
- { // (D)
-
- final double n0 = 1.00060006;
- final double n1 = 263.991156; final double d1 = 257.368075;
- final double n2 = 4373.20068; final double d2 = 3414.48018;
-
- double x = 1-r;
- double x2 = x*x;
-
- return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
-
- }
- else
- { // (E)
-
- final double n0 = 1.00001538;
- final double n1 = 6075.14119; final double d1 = 6065.11919;
- final double n2 = 734266.409; final double d2 = 694021.044;
-
- double x = 1-r;
- double x2 = x*x;
-
- return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
-
- }
-
- } // transform()
-
-
-
-//
-// Table of values of inverse Landau, from r = .060 to .982
-//
-
-// Since all these are this is static to this compilation unit only, the
-// info is establised a priori and not at each invocation.
-
- static final float TABLE_INTERVAL = .001f;
- static final int TABLE_END = 982;
- static final float TABLE_MULTIPLIER = 1.0f/TABLE_INTERVAL;
-
-// Here comes the big (4K bytes) table ---
-//
-// inverseLandau[ n ] = the inverse cdf at r = n*TABLE_INTERVAL = n/1000.
-//
-// Credit CERNLIB for these computations.
-//
-// This data is float because the algortihm does not benefit from further
-// data accuracy. The numbers below .006 or above .982 are moot since
-// non-table-based methods are used below r=.007 and above .980.
-
- static final float[] inverseLandau = {
-
- 0.0f, // .000
- 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, // .001 - .005
- -2.244733f, -2.204365f, -2.168163f, -2.135219f, -2.104898f, // .006 - .010
- -2.076740f, -2.050397f, -2.025605f, -2.002150f, -1.979866f,
- -1.958612f, -1.938275f, -1.918760f, -1.899984f, -1.881879f, // .020
- -1.864385f, -1.847451f, -1.831030f, -1.815083f, -1.799574f,
- -1.784473f, -1.769751f, -1.755383f, -1.741346f, -1.727620f, // .030
- -1.714187f, -1.701029f, -1.688130f, -1.675477f, -1.663057f,
- -1.650858f, -1.638868f, -1.627078f, -1.615477f, -1.604058f, // .040
- -1.592811f, -1.581729f, -1.570806f, -1.560034f, -1.549407f,
- -1.538919f, -1.528565f, -1.518339f, -1.508237f, -1.498254f, // .050
- -1.488386f, -1.478628f, -1.468976f, -1.459428f, -1.449979f,
- -1.440626f, -1.431365f, -1.422195f, -1.413111f, -1.404112f, // .060
- -1.395194f, -1.386356f, -1.377594f, -1.368906f, -1.360291f,
- -1.351746f, -1.343269f, -1.334859f, -1.326512f, -1.318229f, // .070
- -1.310006f, -1.301843f, -1.293737f, -1.285688f, -1.277693f,
- -1.269752f, -1.261863f, -1.254024f, -1.246235f, -1.238494f, // .080
- -1.230800f, -1.223153f, -1.215550f, -1.207990f, -1.200474f,
- -1.192999f, -1.185566f, -1.178172f, -1.170817f, -1.163500f, // .090
- -1.156220f, -1.148977f, -1.141770f, -1.134598f, -1.127459f,
- -1.120354f, -1.113282f, -1.106242f, -1.099233f, -1.092255f, // .100
-
- -1.085306f, -1.078388f, -1.071498f, -1.064636f, -1.057802f,
- -1.050996f, -1.044215f, -1.037461f, -1.030733f, -1.024029f,
- -1.017350f, -1.010695f, -1.004064f, -.997456f, -.990871f,
- -.984308f, -.977767f, -.971247f, -.964749f, -.958271f,
- -.951813f, -.945375f, -.938957f, -.932558f, -.926178f,
- -.919816f, -.913472f, -.907146f, -.900838f, -.894547f,
- -.888272f, -.882014f, -.875773f, -.869547f, -.863337f,
- -.857142f, -.850963f, -.844798f, -.838648f, -.832512f,
- -.826390f, -.820282f, -.814187f, -.808106f, -.802038f,
- -.795982f, -.789940f, -.783909f, -.777891f, -.771884f, // .150
- -.765889f, -.759906f, -.753934f, -.747973f, -.742023f,
- -.736084f, -.730155f, -.724237f, -.718328f, -.712429f,
- -.706541f, -.700661f, -.694791f, -.688931f, -.683079f,
- -.677236f, -.671402f, -.665576f, -.659759f, -.653950f,
- -.648149f, -.642356f, -.636570f, -.630793f, -.625022f,
- -.619259f, -.613503f, -.607754f, -.602012f, -.596276f,
- -.590548f, -.584825f, -.579109f, -.573399f, -.567695f,
- -.561997f, -.556305f, -.550618f, -.544937f, -.539262f,
- -.533592f, -.527926f, -.522266f, -.516611f, -.510961f,
- -.505315f, -.499674f, -.494037f, -.488405f, -.482777f, // .200
-
- -.477153f, -.471533f, -.465917f, -.460305f, -.454697f,
- -.449092f, -.443491f, -.437893f, -.432299f, -.426707f,
- -.421119f, -.415534f, -.409951f, -.404372f, -.398795f,
- -.393221f, -.387649f, -.382080f, -.376513f, -.370949f,
- -.365387f, -.359826f, -.354268f, -.348712f, -.343157f,
- -.337604f, -.332053f, -.326503f, -.320955f, -.315408f,
- -.309863f, -.304318f, -.298775f, -.293233f, -.287692f,
- -.282152f, -.276613f, -.271074f, -.265536f, -.259999f,
- -.254462f, -.248926f, -.243389f, -.237854f, -.232318f,
- -.226783f, -.221247f, -.215712f, -.210176f, -.204641f, // .250
- -.199105f, -.193568f, -.188032f, -.182495f, -.176957f,
- -.171419f, -.165880f, -.160341f, -.154800f, -.149259f,
- -.143717f, -.138173f, -.132629f, -.127083f, -.121537f,
- -.115989f, -.110439f, -.104889f, -.099336f, -.093782f,
- -.088227f, -.082670f, -.077111f, -.071550f, -.065987f,
- -.060423f, -.054856f, -.049288f, -.043717f, -.038144f,
- -.032569f, -.026991f, -.021411f, -.015828f, -.010243f,
- -.004656f, .000934f, .006527f, .012123f, .017722f,
- .023323f, .028928f, .034535f, .040146f, .045759f,
- .051376f, .056997f, .062620f, .068247f, .073877f, // .300
-
- .079511f, .085149f, .090790f, .096435f, .102083f,
- .107736f, .113392f, .119052f, .124716f, .130385f,
- .136057f, .141734f, .147414f, .153100f, .158789f,
- .164483f, .170181f, .175884f, .181592f, .187304f,
- .193021f, .198743f, .204469f, .210201f, .215937f,
- .221678f, .227425f, .233177f, .238933f, .244696f,
- .250463f, .256236f, .262014f, .267798f, .273587f,
- .279382f, .285183f, .290989f, .296801f, .302619f,
- .308443f, .314273f, .320109f, .325951f, .331799f,
- .337654f, .343515f, .349382f, .355255f, .361135f, // .350
- .367022f, .372915f, .378815f, .384721f, .390634f,
- .396554f, .402481f, .408415f, .414356f, .420304f,
- .426260f, .432222f, .438192f, .444169f, .450153f,
- .456145f, .462144f, .468151f, .474166f, .480188f,
- .486218f, .492256f, .498302f, .504356f, .510418f,
- .516488f, .522566f, .528653f, .534747f, .540850f,
- .546962f, .553082f, .559210f, .565347f, .571493f,
- .577648f, .583811f, .589983f, .596164f, .602355f,
- .608554f, .614762f, .620980f, .627207f, .633444f,
- .639689f, .645945f, .652210f, .658484f, .664768f, // .400
-
- .671062f, .677366f, .683680f, .690004f, .696338f,
- .702682f, .709036f, .715400f, .721775f, .728160f,
- .734556f, .740963f, .747379f, .753807f, .760246f,
- .766695f, .773155f, .779627f, .786109f, .792603f,
- .799107f, .805624f, .812151f, .818690f, .825241f,
- .831803f, .838377f, .844962f, .851560f, .858170f,
- .864791f, .871425f, .878071f, .884729f, .891399f,
- .898082f, .904778f, .911486f, .918206f, .924940f,
- .931686f, .938446f, .945218f, .952003f, .958802f,
- .965614f, .972439f, .979278f, .986130f, .992996f, // .450
- .999875f, 1.006769f, 1.013676f, 1.020597f, 1.027533f,
- 1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
- 1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
- 1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
- 1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
- 1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
- 1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
- 1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
- 1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
- 1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f, // .500
-
- 1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
- 1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
- 1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
- 1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
- 1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
- 1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
- 1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
- 1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
- 1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
- 1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f, // .550
- 1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
- 1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
- 1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
- 1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
- 1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
- 2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
- 2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
- 2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
- 2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
- 2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f, // .600
-
- 2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
- 2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
- 2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
- 2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
- 2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
- 2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
- 2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
- 2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
- 2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
- 2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f, // .650
- 2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
- 2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
- 2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
- 3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
- 3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
- 3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
- 3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
- 3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
- 3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
- 3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f, // .700
-
- 3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
- 3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
- 3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
- 3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
- 3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
- 3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
- 4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
- 4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
- 4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
- 4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
- 4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f, // .750
- 4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
- 4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
- 4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
- 4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
- 5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
- 5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
- 5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
- 5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
- 5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f, // .800
-
- 5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
- 5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
- 6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
- 6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
- 6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
- 6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
- 6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
- 7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
- 7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
- 7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f, // .850
- 7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
- 8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
- 8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
- 8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
- 9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
- 9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
- 9.840558f, 9.922186f, 10.005107f, 10.089353f, 10.174959f,
- 10.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
- 10.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
- 11.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f, // .900
-
- 11.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
- 12.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
- 13.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
- 13.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
- 14.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
- 15.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
- 16.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
- 17.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
- 19.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
- 20.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f, // .950
- 22.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
- 25.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f, // .960
- 28.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
- 32.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f, // .970
- 37.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
- 44.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f, // .980
- 56.123356f, 59.103894f, // .982
-
- }; // End of the inverseLandau table
-
-
-// root method
-
- static final double[] f = {
- 0 , 0 , 0 ,0 ,0 ,-2.244733,
- -2.204365,-2.168163,-2.135219,-2.104898,-2.076740,-2.050397,
- -2.025605,-2.002150,-1.979866,-1.958612,-1.938275,-1.918760,
- -1.899984,-1.881879,-1.864385,-1.847451,-1.831030,-1.815083,
- -1.799574,-1.784473,-1.769751,-1.755383,-1.741346,-1.727620,
- -1.714187,-1.701029,-1.688130,-1.675477,-1.663057,-1.650858,
- -1.638868,-1.627078,-1.615477,-1.604058,-1.592811,-1.581729,
- -1.570806,-1.560034,-1.549407,-1.538919,-1.528565,-1.518339,
- -1.508237,-1.498254,-1.488386,-1.478628,-1.468976,-1.459428,
- -1.449979,-1.440626,-1.431365,-1.422195,-1.413111,-1.404112,
- -1.395194,-1.386356,-1.377594,-1.368906,-1.360291,-1.351746,
- -1.343269,-1.334859,-1.326512,-1.318229,-1.310006,-1.301843,
- -1.293737,-1.285688,-1.277693,-1.269752,-1.261863,-1.254024,
- -1.246235,-1.238494,-1.230800,-1.223153,-1.215550,-1.207990,
- -1.200474,-1.192999,-1.185566,-1.178172,-1.170817,-1.163500,
- -1.156220,-1.148977,-1.141770,-1.134598,-1.127459,-1.120354,
- -1.113282,-1.106242,-1.099233,-1.092255,
- -1.085306,-1.078388,-1.071498,-1.064636,-1.057802,-1.050996,
- -1.044215,-1.037461,-1.030733,-1.024029,-1.017350,-1.010695,
- -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247,
- -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
- -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
- -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
- -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
- -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
- -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
- -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
- -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
- -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
- -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
- -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
- -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
- -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
- -.499674, -.494037, -.488405, -.482777,
- -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
- -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
- -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
- -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
- -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
- -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
- -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
- -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
- -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
- -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
- -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
- -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
- -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
- -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
- -.010243, -.004656, .000934, .006527, .012123, .017722,
- .023323, .028928, .034535, .040146, .045759, .051376,
- .056997, .062620, .068247, .073877,
- .079511, .085149, .090790, .096435, .102083, .107736,
- .113392, .119052, .124716, .130385, .136057, .141734,
- .147414, .153100, .158789, .164483, .170181, .175884,
- .181592, .187304, .193021, .198743, .204469, .210201,
- .215937, .221678, .227425, .233177, .238933, .244696,
- .250463, .256236, .262014, .267798, .273587, .279382,
- .285183, .290989, .296801, .302619, .308443, .314273,
- .320109, .325951, .331799, .337654, .343515, .349382,
- .355255, .361135, .367022, .372915, .378815, .384721,
- .390634, .396554, .402481, .408415, .414356, .420304,
- .426260, .432222, .438192, .444169, .450153, .456145,
- .462144, .468151, .474166, .480188, .486218, .492256,
- .498302, .504356, .510418, .516488, .522566, .528653,
- .534747, .540850, .546962, .553082, .559210, .565347,
- .571493, .577648, .583811, .589983, .596164, .602355,
- .608554, .614762, .620980, .627207, .633444, .639689,
- .645945, .652210, .658484, .664768,
- .671062, .677366, .683680, .690004, .696338, .702682,
- .709036, .715400, .721775, .728160, .734556, .740963,
- .747379, .753807, .760246, .766695, .773155, .779627,
- .786109, .792603, .799107, .805624, .812151, .818690,
- .825241, .831803, .838377, .844962, .851560, .858170,
- .864791, .871425, .878071, .884729, .891399, .898082,
- .904778, .911486, .918206, .924940, .931686, .938446,
- .945218, .952003, .958802, .965614, .972439, .979278,
- .986130, .992996, .999875, 1.006769, 1.013676, 1.020597,
- 1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
- 1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
- 1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
- 1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
- 1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
- 1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
- 1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
- 1.332819, 1.340454, 1.348108, 1.355780,
- 1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
- 1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
- 1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
- 1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
- 1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
- 1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
- 1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
- 1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
- 1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
- 1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
- 1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
- 1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
- 1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
- 2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
- 2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
- 2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
- 2.212257, 2.222533, 2.232845, 2.243195,
- 2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
- 2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
- 2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
- 2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
- 2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
- 2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
- 2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
- 2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
- 2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
- 2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
- 2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
- 3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
- 3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
- 3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
- 3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
- 3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
- 3.478500, 3.494372, 3.510328, 3.526370,
- 3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
- 3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
- 3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
- 3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
- 3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
- 4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
- 4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
- 4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
- 4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
- 4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
- 4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
- 4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
- 4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
- 5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
- 5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
- 5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
- 5.673552, 5.704634, 5.735986, 5.767610,
- 5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
- 5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
- 6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
- 6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
- 6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
- 6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
- 7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
- 7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
- 7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
- 8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
- 8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
- 8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
- 9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
- 9.681054, 9.760191, 9.840558, 9.922186,10.005107,10.089353,
- 10.174959,10.261958,10.350389,10.440287,10.531693,10.624646,
- 10.719188,10.815362,10.913214,11.012789,11.114137,11.217307,
- 11.322352,11.429325,11.538283,11.649285,
- 11.762390,11.877664,11.995170,12.114979,12.237161,12.361791,
- 12.488946,12.618708,12.751161,12.886394,13.024498,13.165570,
- 13.309711,13.457026,13.607625,13.761625,13.919145,14.080314,
- 14.245263,14.414134,14.587072,14.764233,14.945778,15.131877,
- 15.322712,15.518470,15.719353,15.925570,16.137345,16.354912,
- 16.578520,16.808433,17.044929,17.288305,17.538873,17.796967,
- 18.062943,18.337176,18.620068,18.912049,19.213574,19.525133,
- 19.847249,20.180480,20.525429,20.882738,21.253102,21.637266,
- 22.036036,22.450278,22.880933,23.329017,23.795634,24.281981,
- 24.789364,25.319207,25.873062,26.452634,27.059789,27.696581,
- 28.365274,29.068370,29.808638,30.589157,31.413354,32.285060,
- 33.208568,34.188705,35.230920,36.341388,37.527131,38.796172,
- 40.157721,41.622399,43.202525,44.912465,46.769077,48.792279,
- 51.005773,53.437996,56.123356,59.103894 };
-
- public double nextLandau()
- {
- if (_sigma <= 0) return 0;
- double ranlan, x, u, v;
- x = _ran.nextDouble();
- u = 1000.*x;
- int i = (int) u;
- u -= i;
- if (i >= 70 && i < 800)
- {
- ranlan = f[i-1] + u*(f[i] - f[i-1]);
- }
- else if (i >= 7 && i <= 980)
- {
- ranlan = f[i-1] + u*(f[i]-f[i-1]-0.25*(1-u)*(f[i+1]-f[i]-f[i-1]+f[i-2]));
- }
- else if (i < 7)
- {
- v = log(x);
- u = 1/v;
- ranlan = ((0.99858950+(3.45213058E1+1.70854528E1*u)*u)/
- (1 +(3.41760202E1+4.01244582 *u)*u))*
- (log(-0.91893853-v)-1);
- }
- else
- {
- u = 1-x;
- v = u*u;
- if (x <= 0.999)
- {
- ranlan = (1.00060006+2.63991156E2*u+4.37320068E3*v)/
- ((1 +2.57368075E2*u+3.41448018E3*v)*u);
- }
- else
- {
- ranlan = (1.00001538+6.07514119E3*u+7.34266409E5*v)/
- ((1 +6.06511919E3*u+6.94021044E5*v)*u);
- }
- }
- return _mean + _sigma*ranlan;
- }
-
-}
diff -N BilinearInterpolator.java --- BilinearInterpolator.java 6 Jun 2008 07:16:32 -0000 1.2 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,195 +0,0 @@
-/*
- * BilinearInterpolator.java
- *
- * Created on June 3, 2008, 4:04 PM
- *
- * $Id: BilinearInterpolator.java,v 1.2 2008/06/06 07:16:32 ngraf Exp $
- */
-//
-package org.lcsim.math.interpolation;
-
-import static java.lang.Math.abs;
-
-/**
- * A class to provide interpolated values for values determined at discrete points
- * on a 2D grid
- *
- * @author Norman Graf
- */
-public class BilinearInterpolator
-{
- private double[] _x;
- private double _xmin;
- private double _xmax;
- private int _xDim;
- private double[] _y;
- private double _ymin;
- private double _ymax;
- private int _yDim;
- private double[][] _val;
- /**
- * Creates a new instance of BilinearInterpolator
- * @param x Array of first independent variable at which values are known
- * @param y Array of second independent variable at which values are known
- * @param z Array of values at the (x,y) points
- */
- public BilinearInterpolator(double[] x, double[] y, double[][] z)
- {
- _x = new double[x.length];;
- System.arraycopy(x,0,_x,0, x.length);
- _xDim = _x.length;
- _xmin = _x[0];
- _xmax = _x[_xDim-1];
-
- _y = new double[y.length];
- System.arraycopy(y,0,_y,0, y.length);
- _yDim = _y.length;
- _ymin = _y[0];
- _ymax = _y[_yDim-1];
-
-
-
- _val = new double[z.length][z[0].length];
- for(int i=0; i< z.length; ++i)
- {
- System.arraycopy(z[i],0,_val[i],0,z[i].length);
- }
- }
-
- //TODO add protection for values out of range
- /**
- * Return the value at an arbitrary x,y point, using bilinear interpolation
- * @param x the first independent variable
- * @param y the second independent variable
- * @return the interpolated value at (x,y)
- */
- public double interpolateValueAt(double x, double y)
- {
- return trueBilinear( x, y );
-// return polin2(_x, _y, _val, x, y);
- }
-
- /*
- * Given arrays x1a[1..m] and x2a[1..n] of independent variables, and a submatrix of function
- * values ya[1..m][1..n], tabulated at the grid points defined by x1a and x2a; and given values
- * x1 and x2 of the independent variables; this routine returns an interpolated function value y.
- */
- double polin2(double[] x1a, double[] x2a, double[][] ya, double x1, double x2)
- {
- int m = x1a.length;
- int n = x2a.length;
- double[] ymtmp = new double[m];
- double[] yntmp = new double[n];
- for (int j=0; j<m; ++j) //Loop over rows.
- {
- // copy the row into temporary storage
- for(int k = 0; k<n; ++k)
- {
- yntmp[k] = ya[j][k];
- }
- ymtmp[j] = polint(x2a, yntmp, x2 ); //Interpolate answer into temporary storage.
- }
- return polint(x1a, ymtmp, x1); //Do the final interpolation.
- }
-
- /*
- * Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y.
- * If P(x) is the polynomial of degree N -1 such that P(xai) = ya, i ; i = 1... n, then
- * the returned value y = P(x).
- */
- double polint( double[] xa, double[] ya, double x)
- {
- int ns=0;
- int n = xa.length;
- double den,dif,dift,ho,hp,w;
- double dy;
- double[] c = new double[n];
- double[] d = new double[n];
- dif=abs(x-xa[0]);
-
- for (int i=0;i<n;++i)
- { //Here we find the index ns of the closest table entry,
- dift = abs(x-xa[i]);
- if ( dift < dif)
- {
- ns=i;
- dif=dift;
- }
- c[i]=ya[i]; // and initialize the table of c's and d's.
- d[i]=ya[i];
- }
- double y = ya[ns--]; // This is the initial approximation to y.
-// ns = ns-1;
- for (int m=1; m<n ; ++m)
- {
- //For each column of the table,
- for (int i=0; i<n-m;++i)
- { //we loop over the current c's and d's and update them.
- ho=xa[i]-x;
- hp=xa[i+m]-x;
- w=c[i+1]-d[i];
- den = ho-hp;
- if (den == 0.0) System.err.println("Error in routine polint");
- //This error can occur only if two input xa's are (to within roundo ) identical.
- den=w/den;
- d[i]=hp*den; //Here the c's and d's are updated.
- c[i]=ho*den;
- }
- if(2*ns < (n-m))
- {
- dy = c[ns+1];
- }
- else
- {
- dy = d[ns];
- ns = ns-1;
- }
- y += dy;
- //After each column in the tableau is completed, we decide which correction, c or d,
- //we want to add to our accumulating value of y, i.e., which path to take through the
- //tableau|forking up or down. We do this in such a way as to take the most \straight
- //line" route through the tableau to its apex, updating ns accordingly to keep track of
- //where we are. This route keeps the partial approximations centered (insofar as possible)
- //on the target x. The last dy added is thus the error indication.
- }
- return y;
- }
-
- double trueBilinear(double x, double y)
- {
- // find bin for x
- int ixlo = 0;
- int ixhi = 0;
- for(int i=0; i<_xDim-1; ++i)
- {
- if(x>=_x[i] && x<=_x[i+1])
- {
- ixlo = i;
- ixhi = i+1;
- break;
- }
- }
-
- // find bin for y
- int iylo = 0;
- int iyhi = 0;
- for(int i=0; i<_yDim-1; ++i)
- {
- if(y>=_y[i] && y<=_y[i+1])
- {
- iylo = i;
- iyhi = i+1;
- break;
- }
- }
- double v1 = _val[ixlo][iylo];
- double v2 = _val[ixhi][iylo];
- double v3 = _val[ixhi][iyhi];
- double v4 = _val[ixlo][iyhi];
-
- double t = (x - _x[ixlo])/(_x[ixhi] - _x[ixlo]);
- double u = (y - _y[iylo])/(_y[iyhi] - _y[iylo]);
-
- return (1-t)*(1-u)*v1 +t*(1-u)*v2+t*u*v3+(1-t)*u*v4;
- }
-}
\ No newline at end of file
diff -N CentralMomentsCalculator.java --- CentralMomentsCalculator.java 12 Jul 2006 06:45:32 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,174 +0,0 @@
-/*
- * CentralMomentsCalculator.java
- *
- * Created on April 5, 2006, 11:25 AM
- *
- * $Id: CentralMomentsCalculator.java,v 1.1 2006/07/12 06:45:32 ngraf Exp $
- */
-
-package org.lcsim.math.moments;
-
-import Jama.EigenvalueDecomposition;
-import Jama.Matrix;
-
-/**
- * calculates rotational and translational invariant moments of spatial distributions
- * @author Norman Graf
- */
-public class CentralMomentsCalculator
-{
- // 0
- private double m000;
-
- // 1
- private double m100, m010, m001;
-
- // 2
- private double m110, m101, m011;
- private double m200, m020, m002;
-
- // centroids
- private double xc, yc, zc;
-
- // invariants
- private double J1, J2, J3;
-
-
-
- private Matrix _tensor = new Matrix(3,3);
-
- // eigenvalues
- private double[] _eigenvalues = new double[3];
-
- // eigenvectors
- private Matrix _eigenvectors;
-
- // direction cosines for the cluster direction
- private double[] _dircos = new double[3];
-
-
- /**
- * Creates a new instance of CentralMomentsCalculator
- */
- public CentralMomentsCalculator()
- {
- }
-
- public void calculateMoments(double[] x, double[] y, double[] z, double[] w)
- {
- reset();
- //TODO make this more efficient.
- int n = x.length;
- // TODO check that all arrays are the same size.
-
- // calculate centroids
- for(int i=0; i<n; ++i)
- {
- m000 += w[i];
- m100 += x[i]*w[i];
- m010 += y[i]*w[i];
- m001 += z[i]*w[i];
- }
-
- xc = m100/m000;
- yc = m010/m000;
- zc = m001/m000;
-
- // on to the higher moments wrt centroid
- double xa, ya, za;
- for(int i=0; i<n; ++i)
- {
- xa = x[i]-xc;
- ya = y[i]-yc;
- za = z[i]-zc;
-
- m110 += xa*ya*w[i];
- m101 += xa*za*w[i];
- m011 += ya*za*w[i];
-
- m200 += xa*xa*w[i];
- m020 += ya*ya*w[i];
- m002 += za*za*w[i];
- }
-
- // normalize
- m110 /= m000;
- m101 /= m000;
- m011 /= m000;
-
- m200 /= m000;
- m020 /= m000;
- m002 /= m000;
-
- J1 = m200 + m020 + m002;
- J2 = m200*m020 + m200*m002 + m020*m002 - m110*m110 - m101*m101 - m011*m011;
- J3 = m200*m020*m002 + 2.*m110*m101*m011 - m002*m110*m110 - m020*m101*m101 - m200*m011*m011;
-
- // now for eigenvalues, eigenvectors
- _tensor.set(0, 0, m020 + m002);
- _tensor.set(1,0, - m110);
- _tensor.set(2,0, - m101);
- _tensor.set(0,1, - m110);
- _tensor.set(1,1, + m200 + m002);
- _tensor.set(2,1, - m011);
- _tensor.set(0,2, - m101);
- _tensor.set(1,2, - m011);
- _tensor.set(2,2, + m200 + m020);
-
- EigenvalueDecomposition eig = _tensor.eig();
- _eigenvalues = eig.getRealEigenvalues();
-// System.out.println("eigenvalues: "+_eigenvalues[0]+" "+_eigenvalues[1]+" "+_eigenvalues[2]);
- _eigenvectors = eig.getV();
-// System.out.println("eigenvectors:");
-// _eigenvectors.print(4,4);
-
- // direction cosines are the eigenvector elements corresponding to the lowest eigenvalue
- // note that eigenvalues are sorted in ascending order by EigenvalueDecomposition
- _dircos[0] = -_eigenvectors.get(0,0);
- _dircos[1] = -_eigenvectors.get(1,0);
- _dircos[2] = -_eigenvectors.get(2,0);
-// System.out.println("dircos: "+_dircos[0]+" "+_dircos[1]+" "+_dircos[2]);
-
- }
-
- public double[] eigenvalues()
- {
- return _eigenvalues;
- }
-
- public double[] directionCosines()
- {
- return _dircos;
- }
-
- public double[] centroid()
- {
- return new double[] {xc, yc, zc};
- }
-
- public double[] centroidVariance()
- {
- return new double[] {m200, m020, m002, m110, m101, m011};
- }
-
- public double[] invariants()
- {
-
- return new double[] {J1, J2, J3};
- }
-
- void reset()
- {
- m000 = 0.;
- m100 = 0.;
- m010 = 0.;
- m001 = 0.;
- m110 = 0.;
- m101 = 0.;
- m011 = 0.;
- m200 = 0.;
- m020 = 0.;
- m002 = 0.;
- }
-
-}
diff -N BivariateDistribution.java --- BivariateDistribution.java 31 Mar 2009 17:10:38 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,258 +0,0 @@
-/*
- * Class BivariateDistribution
- */
-package org.lcsim.math.probability;
-
-/**
- * Calculate the probability integral for a set of bins in the x-y plane
- * of a bivariate normal distribution (i.e., a 2D Gaussian probability).
- *<p>
- * The evaluation of the probability integrals is described in:
- *<p>
- * Alan Genz, "Numerical Computation of Rectangular Bivariate and Trivariate
- * Normal and t Probabilities" in Statistics and Computing 14, 151 (2004).
- *<p>
- * The integration code is adapted from the FORTRAN source at:
- *<p>
- * http://www.math.wsu.edu/faculty/genz/homepage
- *<p>
- * @author Richard Partridge
- */
-public class BivariateDistribution {
-
- private int _nx;
- private int _ny;
- private double _xmin;
- private double _ymin;
- private double _dx;
- private double _dy;
- private double[] _h;
- private double[] _k;
-
- // Weights and coordinates for 6 point Gauss-Legendre integration
- private double[] _w6 = {0.1713244923791705, 0.3607615730481384, 0.4679139345726904};
- private double[] _x6 = {0.9324695142031522, 0.6612093864662647, 0.2386191860831970};
-
- // Weights and coordinates for 12 point Gauss-Legendre integration
- private double[] _w12 = {.04717533638651177, 0.1069393259953183, 0.1600783285433464,
- 0.2031674267230659, 0.2334925365383547, 0.2491470458134029};
- private double[] _x12 = {0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
- 0.5873179542866171, 0.3678314989981802, 0.1252334085114692};
-
- // Weights and coordinates for 20 point Gauss-Legendre integration
- private double[] _w20 = {.01761400713915212, .04060142980038694, .06267204833410906,
- .08327674157670475, 0.1019301198172404, 0.1181945319615184,
- 0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
- 0.1527533871307259};
- private double[] _x20 = {0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
- 0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
- 0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
- 0.07652652113349733};
-
- /**
- * Set the locations of the x-coordinate bins
- *
- * @param nx number of x coordinate bins
- * @param xmin minimum x coordinate
- * @param dx width of x coordinate bins
- */
- public void xBins(int nx, double xmin, double dx) {
- _nx = nx;
- _xmin = xmin;
- _dx = dx;
- _h = new double[_nx + 1];
- }
-
- /**
- * Set the locations of the y-coordinate bins
- *
- * @param ny number of y coordinate bins
- * @param ymin minimum y coordinate
- * @param dy width of y coordinate bins
- */
- public void yBins(int ny, double ymin, double dy) {
- _ny = ny;
- _ymin = ymin;
- _dy = dy;
- _k = new double[_ny + 1];
- }
-
- /**
- * Integrate the Gaussian probability distribution over each x-y bins,
- * which must be defined before calling this method.
- * <p>
- * The output is a double array that gives the binned probability
- * distribution. The first array index is used to indicate the bin in x
- * and the second array index is used to indicate the bin in y.
- * <p>
- * @param x0 mean x coordinate of Gaussian distribution
- * @param y0 mean y coordinate of Gaussian distribution
- * @param sigx x coordinate standard deviation
- * @param sigy y coordinate standard deviation
- * @param rho x-y correlation coefficient
- * @return probability distribution
- */
- public double[][] Calculate(double x0, double y0, double sigx, double sigy,
- double rho) {
-
- // Calculate the scaled x coordinate for each bin edge
- for (int i = 0; i < _nx + 1; i++) {
- _h[i] = (_xmin + i * _dx - x0) / sigx;
- }
-
- // Calculate the scaled y coordinate for each bin edge
- for (int j = 0; j < _ny + 1; j++) {
- _k[j] = (_ymin + j * _dy - y0) / sigy;
- }
-
- // Create the array that will hold the binned probabilities
- double[][] bi = new double[_nx][_ny];
-
- // Loop over the bin vertices
- for (int i = 0; i < _nx + 1; i++) {
- for (int j = 0; j < _ny + 1; j++) {
-
- // Calculate the probability for x>h and y>k for this vertex
- double prob = GenzCalc(_h[i], _k[j], rho);
-
- // Add or subtract this probability from the affected bins.
- // The bin probability for bin (0,0) is the sum of the Genz
- // probabilities for the (0,0) and (1,1) vertices MINUS the
- // sum of the probabilities for the (0,1) and (1,0) vertices
- if (i > 0 && j > 0) {
- bi[i - 1][j - 1] += prob;
- }
- if (i > 0 && j < _ny) {
- bi[i - 1][j] -= prob;
- }
- if (i < _nx && j > 0) {
- bi[i][j - 1] -= prob;
- }
- if (i < _nx && j < _ny) {
- bi[i][j] += prob;
- }
- }
- }
-
- return bi;
- }
-
- private double GenzCalc(double dh, double dk, double rho) {
-
- double twopi = 2. * Math.PI;
-
- // Declare the Gauss-Legendre constants
- int ng;
- double[] w;
- double[] x;
-
- if (Math.abs(rho) < 0.3) {
- // for rho < 0.3 use 6 point Gauss-Legendre integration
- ng = 3;
- w = _w6;
- x = _x6;
- } else if (Math.abs(rho) < 0.75) {
- // for 0.3 < rho < 0.75 use 12 point Gauss-Legendre integration
- ng = 6;
- w = _w12;
- x = _x12;
- } else {
- // for rho > 0.75 use 20 point Gauss-Legendre integration
- ng = 10;
- w = _w20;
- x = _x20;
- }
-
- // Initialize the probability and some local variables
- double bvn = 0.;
- double h = dh;
- double k = dk;
- double hk = h * k;
-
- // For rho < 0.925, integrate equation 3 in the Genz paper
- if (Math.abs(rho) < 0.925) {
-
- // More or less direct port of Genz code follows
- // It is fairly easy to match this calculation against equation 3 of
- // Genz's paper if you take into account that you need to change
- // variables so the integration argument spans the range -1 to 1
- double hs = (h * h + k * k) / 2.;
- double asr = Math.asin(rho);
- double sn;
- for (int i = 0; i < ng; i++) {
- sn = Math.sin(asr * (1 - x[i]) / 2.);
- bvn += w[i] * Math.exp((sn * hk - hs) / (1 - sn * sn));
- sn = Math.sin(asr * (1 + x[i]) / 2.);
- bvn += w[i] * Math.exp((sn * hk - hs) / (1 - sn * sn));
- }
- // The factor of asr/2 comes from changing variables so the
- // integration is over the range -1 to 1 instead of 0 - asin(rho)
- bvn = bvn * asr / (2. * twopi) + Erf.phi(-h) * Erf.phi(-k);
-
- } else {
- // rho > 0.925 - integrate equation 6 in Genz paper with the
- // extra term in the Taylor expansion given in equation 7.
- // The rest of this code is pretty dense and is a pretty direct
- // port of Genz's code.
-
- if (rho < 0.) {
- k = -k;
- hk = -hk;
- }
-
- if (Math.abs(rho) < 1.) {
-
- double as = (1 - rho) * (1 + rho);
- double a = Math.sqrt(as);
- double bs = (h - k) * (h - k);
- double c = (4. - hk) / 8.;
- double d = (12. - hk) / 16.;
- double asr = -(bs / as + hk) / 2.;
-
- if (asr > -100.) {
- bvn = a * Math.exp(asr) *
- (1. - c * (bs - as) * (1. - d * bs / 5.) / 3. +
- c * d * as * as / 5.);
- }
-
- if (-hk < 100.) {
- double b = Math.sqrt(bs);
- bvn -= Math.exp(-hk / 2.) * Math.sqrt(twopi) * Erf.phi(-b / a) *
- b * (1 - c * bs * (1 - d * bs / 5.) / 3.);
- }
-
- a = a / 2.;
- for (int i = 0; i < ng; i++) {
- for (int j = 0; j < 2; j++) {
- int is = -1;
- if (j > 0) {
- is = 1;
- }
- double xs = Math.pow(a * (is * x[i] + 1), 2);
- double rs = Math.sqrt(1 - xs);
- asr = -(bs / xs + hk) / 2;
-
- if (asr > -100) {
- double sp = (1 + c * xs * (1 + d * xs));
- double ep = Math.exp(-hk * (1 - rs) / (2 * (1 + rs))) / rs;
- bvn += a * w[i] * Math.exp(asr) * (ep - sp);
- }
- }
- }
-
- bvn = -bvn / twopi;
- }
-
- if (rho > 0) {
- bvn = bvn + Erf.phi(-Math.max(h, k));
- } else {
- bvn = -bvn;
- if (k > h) {
- bvn += Erf.phi(k) - Erf.phi(h);
- }
- }
- }
-
- return Math.max(0, Math.min(1, bvn));
- }
-}
\ No newline at end of file
diff -N Erf.java --- Erf.java 31 Mar 2009 17:10:38 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,283 +0,0 @@
-/*
- * Class Erf
- *
- */
-package org.lcsim.math.probability;
-
-/**
- *
- * Calculates the following probability integrals:
- *<p>
- * erf(x) <br>
- * erfc(x) = 1 - erf(x) <br>
- * phi(x) = 0.5 * erfc(-x/sqrt(2)) <br>
- * phic(x) = 0.5 * erfc(x/sqrt(2))
- *<p>
- * Note that phi(x) gives the probability for an observation smaller than x for
- * a Gaussian probability distribution with zero mean and unit standard
- * deviation, while phic(x) gives the probability for an observation larger
- * than x.
- *<p>
- * The algorithms for erf(x) and erfc(x) are based on Schonfelder's work.
- * See J.L. Schonfelder, "Chebyshev Expansions for the Error and Related
- * Functions", Math. Comp. 32, 1232 (1978). The calculations of phi(x)
- * and phic(x) are trivially calculated using the erfc(x) algorithm.
- *<p>
- * Schonfelder's algorithms are advertised to provide "30 digit" accurracy.
- * Since this level of accuracy exceeds the machine precision for doubles,
- * summation terms whose relative weight is below machine precision are
- * dropped.
- *<p>
- * In this algorithm, we calculate
- *<p>
- * erf(x) = x* y(t) for |x| < 2 <br>
- * erf(x) = 1 - exp(-x*x) * y(t) / x for x >= 2 <br>
- * erfc(|x|) = exp(-x*x)*y(|x|) <br>
- *<p>
- * The functions y(x) are expanded in terms of Chebyshev polynomials, where
- * there is a different set of coefficients a[r] for each of the above 3 cases.
- *<p>
- * y(x) = Sum'( a[r] * T(r, t) )
- *<p>
- * The notation Sum' indicates that the r = 0 term is divided by 2.
- *<p>
- * The variable t is defined as
- *<p>
- * t = ( x*x - 2 ) / 2 for erf(x) with x < 2 <br>
- * t = ( 21 - 2*x*x ) / (5 + 2*x*x) for erf(x) with x >= 2 <br>
- * t = ( 4*|x| - 15 ) / ( 4*|x| + 15 ) for erfc(x)
- *<p>
- * The code and implementation are based on Alan Genz's FORTRAN source code
- * that can be found at http://www.math.wsu.edu/faculty/genz/homepage.
- *<p>
- * Genz's code was a bit tricky to "reverse engineer", so we go through the
- * way these calculations are performed in some detail. Rather than calculate
- * y(x) directly, he calculates
- *<p>
- * bm = Sum( a[r] * U(r, t) ) r = 0 : N <br>
- * bp = Sum( a[r] * U(r-2, t) ) r = 2 : N
- *<p>
- * where U(r, t) are Chebyshev polynomials of the second kind. The coefficients
- * a[r] decrease with r, and the value of N is chosen where a[N] / a[0] is
- * ~10^-16, reflecting the machine precision for doubles.
- *<p>
- * The Chebyshev polynomials of the second kind U(r, t) are calculated using the
- * recursion relation:
- *<p>
- * U(r, t) = 2 * t * U(r-1, t) - U(r-2, t)
- *<p>
- * Genz uses the identity
- *<p>
- * T(r, t) = ( U(r, t) - U(r-2, t) ) / 2
- *<p>
- * to calculate y(x)
- *<p>
- * y(x) = ( bm - bp ) / 2.
- *<p>
- * Note that we get the correct contributions for the r = 0 and r = 1 terms by
- * ignoring these terms in the bp sum, including getting the desired factor
- * of 1/2 in the contribution from the r = 0 term.
- *
- * @author Richard Partridge
- */
-public class Erf {
-
- private static double rtwo = 1.414213562373095048801688724209e0;
-
- // Coefficients for the erf(x) calculation with |x| < 2
- private static double[] a1 = {
- 1.483110564084803581889448079057e0,
- -3.01071073386594942470731046311e-1,
- 6.8994830689831566246603180718e-2,
- -1.3916271264722187682546525687e-2,
- 2.420799522433463662891678239e-3,
- -3.65863968584808644649382577e-4,
- 4.8620984432319048282887568e-5,
- -5.749256558035684835054215e-6,
- 6.11324357843476469706758e-7,
- -5.8991015312958434390846e-8,
- 5.207009092068648240455e-9,
- -4.23297587996554326810e-10,
- 3.1881135066491749748e-11,
- -2.236155018832684273e-12,
- 1.46732984799108492e-13,
- -9.044001985381747e-15,
- 5.25481371547092e-16};
-
- // Coefficients for the err(x) calculation with x > 2
- private static double[] a2 = {
- 1.077977852072383151168335910348e0,
- -2.6559890409148673372146500904e-2,
- -1.487073146698099509605046333e-3,
- -1.38040145414143859607708920e-4,
- -1.1280303332287491498507366e-5,
- -1.172869842743725224053739e-6,
- -1.03476150393304615537382e-7,
- -1.1899114085892438254447e-8,
- -1.016222544989498640476e-9,
- -1.37895716146965692169e-10,
- -9.369613033737303335e-12,
- -1.918809583959525349e-12,
- -3.7573017201993707e-14,
- -3.7053726026983357e-14,
- 2.627565423490371e-15,
- -1.121322876437933e-15,
- 1.84136028922538e-16};
-
- // Coefficients for the erfc(x) calculation
- private static double[] a3 = {
- 6.10143081923200417926465815756e-1,
- -4.34841272712577471828182820888e-1,
- 1.76351193643605501125840298123e-1,
- -6.0710795609249414860051215825e-2,
- 1.7712068995694114486147141191e-2,
- -4.321119385567293818599864968e-3,
- 8.54216676887098678819832055e-4,
- -1.27155090609162742628893940e-4,
- 1.1248167243671189468847072e-5,
- 3.13063885421820972630152e-7,
- -2.70988068537762022009086e-7,
- 3.0737622701407688440959e-8,
- 2.515620384817622937314e-9,
- -1.028929921320319127590e-9,
- 2.9944052119949939363e-11,
- 2.6051789687266936290e-11,
- -2.634839924171969386e-12,
- -6.43404509890636443e-13,
- 1.12457401801663447e-13,
- 1.7281533389986098e-14,
- -4.264101694942375e-15,
- -5.45371977880191e-16,
- 1.58697607761671e-16,
- 2.0899837844334e-17,
- -5.900526869409e-18};
-
- /**
- * Calculate the error function
- *
- * @param x argument
- * @return error function
- */
- public static double erf(double x) {
-
- // Initialize
- double xa = Math.abs(x);
- double erf;
-
- // Case 1: |x| < 2
- if (xa < 2.) {
-
- // First calculate 2*t
- double tt = x*x - 2.;
-
- // Initialize the recursion variables.
- double bm = 0.;
- double b = 0.;
- double bp = 0.;
-
- // Calculate bm and bp as defined above
- for (int i = 16; i >= 0; i--) {
- bp = b;
- b = bm;
- bm = tt * b - bp + a1[i];
- }
-
- // Finally, calculate erfc using the Chebyshev polynomial identity
- erf = x * (bm - bp) / 2.;
-
- } else {
-
- // Case 2: |x| >= 2
-
- // First calculate 2*t
- double tt = (42. - 4 * xa*xa) / (5. + 2 * xa*xa);
-
- // Initialize the recursion variables.
- double bm = 0.;
- double b = 0.;
- double bp = 0.;
-
- // Calculate bm and bp as defined above
- for (int i = 16; i >= 0; i--) {
- bp = b;
- b = bm;
- bm = tt * b - bp + a2[i];
- }
-
- // Finally, calculate erfc using the Chebyshev polynomial identity
- erf = 1. - Math.exp(-x * x) * (bm - bp) / (2. * xa);
-
- // Take care of negative argument for case 2
- if (x < 0.) erf = -erf;
- }
-
- // Finished both cases!
- return erf;
- }
-
- /**
- * Calculate the error function complement
- * @param x argument
- * @return error function complement
- */
- public static double erfc(double x) {
-
- // Initialize
- double xa = Math.abs(x);
- double erfc;
-
- // Set phi to 0 when the argument is too big
- if (xa > 100.) {
- erfc = 0.;
- } else {
-
- // First calculate 2*t
- double tt = (8. * xa - 30.) / (4. * xa + 15.);
-
- // Initialize the recursion variables.
- double bm = 0.;
- double b = 0.;
- double bp = 0.;
-
- // Calculate bm and bp as defined above
- for (int i = 24; i >= 0; i--) {
- bp = b;
- b = bm;
- bm = tt * b - bp + a3[i];
- }
-
- // Finally, calculate erfc using the Chebyshev polynomial identity
- erfc = Math.exp(-x * x) * (bm - bp) / 2.;
- }
-
- // Cacluate erfc for negative arguments
- if (x < 0.) erfc = 2. - erfc;
-
- return erfc;
- }
-
- /**
- * Calcualate the probability for an observation smaller than x for a
- * Gaussian probability distribution with zero mean and unit standard
- * deviation
- *
- * @param x argument
- * @return probability integral
- */
- public static double phi(double x) {
- return 0.5 * erfc( -x / rtwo);
- }
-
- /**
- * Calculate the probability for an observation larger than x for a
- * Gaussian probability distribution with zero mean and unit standard
- * deviation
- *
- * @param x argument
- * @return probability integral
- */
- public static double phic(double x) {
- return 0.5 * erfc(x / rtwo);
- }
-}
-