Source Code Cross Referenced for Statistics.java in  » Science » weka » weka » core » Java Source Code / Java DocumentationJava Source Code and Java Documentation

Java Source Code / Java Documentation
1. 6.0 JDK Core
2. 6.0 JDK Modules
3. 6.0 JDK Modules com.sun
4. 6.0 JDK Modules com.sun.java
5. 6.0 JDK Modules sun
6. 6.0 JDK Platform
7. Ajax
8. Apache Harmony Java SE
9. Aspect oriented
10. Authentication Authorization
11. Blogger System
12. Build
13. Byte Code
14. Cache
15. Chart
16. Chat
17. Code Analyzer
18. Collaboration
19. Content Management System
20. Database Client
21. Database DBMS
22. Database JDBC Connection Pool
23. Database ORM
24. Development
25. EJB Server geronimo
26. EJB Server GlassFish
27. EJB Server JBoss 4.2.1
28. EJB Server resin 3.1.5
29. ERP CRM Financial
30. ESB
31. Forum
32. GIS
33. Graphic Library
34. Groupware
35. HTML Parser
36. IDE
37. IDE Eclipse
38. IDE Netbeans
39. Installer
40. Internationalization Localization
41. Inversion of Control
42. Issue Tracking
43. J2EE
44. JBoss
45. JMS
46. JMX
47. Library
48. Mail Clients
49. Net
50. Parser
51. PDF
52. Portal
53. Profiler
54. Project Management
55. Report
56. RSS RDF
57. Rule Engine
58. Science
59. Scripting
60. Search Engine
61. Security
62. Sevlet Container
63. Source Control
64. Swing Library
65. Template Engine
66. Test Coverage
67. Testing
68. UML
69. Web Crawler
70. Web Framework
71. Web Mail
72. Web Server
73. Web Services
74. Web Services apache cxf 2.0.1
75. Web Services AXIS2
76. Wiki Engine
77. Workflow Engines
78. XML
79. XML UI
Java
Java Tutorial
Java Open Source
Jar File Download
Java Articles
Java Products
Java by API
Photoshop Tutorials
Maya Tutorials
Flash Tutorials
3ds-Max Tutorials
Illustrator Tutorials
GIMP Tutorials
C# / C Sharp
C# / CSharp Tutorial
C# / CSharp Open Source
ASP.Net
ASP.NET Tutorial
JavaScript DHTML
JavaScript Tutorial
JavaScript Reference
HTML / CSS
HTML CSS Reference
C / ANSI-C
C Tutorial
C++
C++ Tutorial
Ruby
PHP
Python
Python Tutorial
Python Open Source
SQL Server / T-SQL
SQL Server / T-SQL Tutorial
Oracle PL / SQL
Oracle PL/SQL Tutorial
PostgreSQL
SQL / MySQL
MySQL Tutorial
VB.Net
VB.Net Tutorial
Flash / Flex / ActionScript
VBA / Excel / Access / Word
XML
XML Tutorial
Microsoft Office PowerPoint 2007 Tutorial
Microsoft Office Excel 2007 Tutorial
Microsoft Office Word 2007 Tutorial
Java Source Code / Java Documentation » Science » weka » weka.core 
Source Cross Referenced  Class Diagram Java Document (Java Doc) 


0001:        package weka.core;
0002:
0003:        /**
0004:         * Class implementing some distributions, tests, etc. The code is mostly adapted from the CERN
0005:         * Jet Java libraries:
0006:         * 
0007:         * Copyright 2001 University of Waikato
0008:         * Copyright 1999 CERN - European Organization for Nuclear Research.
0009:         * Permission to use, copy, modify, distribute and sell this software and its documentation for
0010:         * any purpose is hereby granted without fee, provided that the above copyright notice appear
0011:         * in all copies and that both that copyright notice and this permission notice appear in
0012:         * supporting documentation. 
0013:         * CERN and the University of Waikato make no representations about the suitability of this 
0014:         * software for any purpose. It is provided "as is" without expressed or implied warranty.
0015:         *
0016:         * @author peter.gedeck@pharma.Novartis.com
0017:         * @author wolfgang.hoschek@cern.ch
0018:         * @author Eibe Frank (eibe@cs.waikato.ac.nz)
0019:         * @author Richard Kirkby (rkirkby@cs.waikato.ac.nz)
0020:         * @version $Revision: 1.9 $
0021:         */
0022:        public class Statistics {
0023:
0024:            /** Some constants */
0025:            protected static final double MACHEP = 1.11022302462515654042E-16;
0026:            protected static final double MAXLOG = 7.09782712893383996732E2;
0027:            protected static final double MINLOG = -7.451332191019412076235E2;
0028:            protected static final double MAXGAM = 171.624376956302725;
0029:            protected static final double SQTPI = 2.50662827463100050242E0;
0030:            protected static final double SQRTH = 7.07106781186547524401E-1;
0031:            protected static final double LOGPI = 1.14472988584940017414;
0032:
0033:            protected static final double big = 4.503599627370496e15;
0034:            protected static final double biginv = 2.22044604925031308085e-16;
0035:
0036:            /*************************************************
0037:             *    COEFFICIENTS FOR METHOD  normalInverse()   *
0038:             *************************************************/
0039:            /* approximation for 0 <= |y - 0.5| <= 3/8 */
0040:            protected static final double P0[] = { -5.99633501014107895267E1,
0041:                    9.80010754185999661536E1, -5.66762857469070293439E1,
0042:                    1.39312609387279679503E1, -1.23916583867381258016E0, };
0043:            protected static final double Q0[] = {
0044:            /* 1.00000000000000000000E0,*/
0045:            1.95448858338141759834E0, 4.67627912898881538453E0,
0046:                    8.63602421390890590575E1, -2.25462687854119370527E2,
0047:                    2.00260212380060660359E2, -8.20372256168333339912E1,
0048:                    1.59056225126211695515E1, -1.18331621121330003142E0, };
0049:
0050:            /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
0051:             * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
0052:             */
0053:            protected static final double P1[] = { 4.05544892305962419923E0,
0054:                    3.15251094599893866154E1, 5.71628192246421288162E1,
0055:                    4.40805073893200834700E1, 1.46849561928858024014E1,
0056:                    2.18663306850790267539E0, -1.40256079171354495875E-1,
0057:                    -3.50424626827848203418E-2, -8.57456785154685413611E-4, };
0058:            protected static final double Q1[] = {
0059:            /*  1.00000000000000000000E0,*/
0060:            1.57799883256466749731E1, 4.53907635128879210584E1,
0061:                    4.13172038254672030440E1, 1.50425385692907503408E1,
0062:                    2.50464946208309415979E0, -1.42182922854787788574E-1,
0063:                    -3.80806407691578277194E-2, -9.33259480895457427372E-4, };
0064:
0065:            /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
0066:             * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
0067:             */
0068:            protected static final double P2[] = { 3.23774891776946035970E0,
0069:                    6.91522889068984211695E0, 3.93881025292474443415E0,
0070:                    1.33303460815807542389E0, 2.01485389549179081538E-1,
0071:                    1.23716634817820021358E-2, 3.01581553508235416007E-4,
0072:                    2.65806974686737550832E-6, 6.23974539184983293730E-9, };
0073:            protected static final double Q2[] = {
0074:            /*  1.00000000000000000000E0,*/
0075:            6.02427039364742014255E0, 3.67983563856160859403E0,
0076:                    1.37702099489081330271E0, 2.16236993594496635890E-1,
0077:                    1.34204006088543189037E-2, 3.28014464682127739104E-4,
0078:                    2.89247864745380683936E-6, 6.79019408009981274425E-9, };
0079:
0080:            /**
0081:             * Computes standard error for observed values of a binomial
0082:             * random variable.
0083:             *
0084:             * @param p the probability of success
0085:             * @param n the size of the sample
0086:             * @return the standard error
0087:             */
0088:            public static double binomialStandardError(double p, int n) {
0089:
0090:                if (n == 0) {
0091:                    return 0;
0092:                }
0093:                return Math.sqrt((p * (1 - p)) / (double) n);
0094:            }
0095:
0096:            /**
0097:             * Returns chi-squared probability for given value and degrees
0098:             * of freedom. (The probability that the chi-squared variate
0099:             * will be greater than x for the given degrees of freedom.)
0100:             *
0101:             * @param x the value
0102:             * @param v the number of degrees of freedom
0103:             * @return the chi-squared probability
0104:             */
0105:            public static double chiSquaredProbability(double x, double v) {
0106:
0107:                if (x < 0.0 || v < 1.0)
0108:                    return 0.0;
0109:                return incompleteGammaComplement(v / 2.0, x / 2.0);
0110:            }
0111:
0112:            /**
0113:             * Computes probability of F-ratio.
0114:             *
0115:             * @param F the F-ratio
0116:             * @param df1 the first number of degrees of freedom
0117:             * @param df2 the second number of degrees of freedom
0118:             * @return the probability of the F-ratio.
0119:             */
0120:            public static double FProbability(double F, int df1, int df2) {
0121:
0122:                return incompleteBeta(df2 / 2.0, df1 / 2.0, df2
0123:                        / (df2 + df1 * F));
0124:            }
0125:
0126:            /**
0127:             * Returns the area under the Normal (Gaussian) probability density
0128:             * function, integrated from minus infinity to <tt>x</tt>
0129:             * (assumes mean is zero, variance is one).
0130:             * <pre>
0131:             *                            x
0132:             *                             -
0133:             *                   1        | |          2
0134:             *  normal(x)  = ---------    |    exp( - t /2 ) dt
0135:             *               sqrt(2pi)  | |
0136:             *                           -
0137:             *                          -inf.
0138:             *
0139:             *             =  ( 1 + erf(z) ) / 2
0140:             *             =  erfc(z) / 2
0141:             * </pre>
0142:             * where <tt>z = x/sqrt(2)</tt>.
0143:             * Computation is via the functions <tt>errorFunction</tt> and <tt>errorFunctionComplement</tt>.
0144:             *
0145:             * @param a the z-value
0146:             * @return the probability of the z value according to the normal pdf
0147:             */
0148:            public static double normalProbability(double a) {
0149:
0150:                double x, y, z;
0151:
0152:                x = a * SQRTH;
0153:                z = Math.abs(x);
0154:
0155:                if (z < SQRTH)
0156:                    y = 0.5 + 0.5 * errorFunction(x);
0157:                else {
0158:                    y = 0.5 * errorFunctionComplemented(z);
0159:                    if (x > 0)
0160:                        y = 1.0 - y;
0161:                }
0162:                return y;
0163:            }
0164:
0165:            /**
0166:             * Returns the value, <tt>x</tt>, for which the area under the
0167:             * Normal (Gaussian) probability density function (integrated from
0168:             * minus infinity to <tt>x</tt>) is equal to the argument <tt>y</tt>
0169:             * (assumes mean is zero, variance is one).
0170:             * <p>
0171:             * For small arguments <tt>0 < y < exp(-2)</tt>, the program computes
0172:             * <tt>z = sqrt( -2.0 * log(y) )</tt>;  then the approximation is
0173:             * <tt>x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z)</tt>.
0174:             * There are two rational functions P/Q, one for <tt>0 < y < exp(-32)</tt>
0175:             * and the other for <tt>y</tt> up to <tt>exp(-2)</tt>. 
0176:             * For larger arguments,
0177:             * <tt>w = y - 0.5</tt>, and  <tt>x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2))</tt>.
0178:             *
0179:             * @param y0 the area under the normal pdf
0180:             * @return the z-value
0181:             */
0182:            public static double normalInverse(double y0) {
0183:
0184:                double x, y, z, y2, x0, x1;
0185:                int code;
0186:
0187:                final double s2pi = Math.sqrt(2.0 * Math.PI);
0188:
0189:                if (y0 <= 0.0)
0190:                    throw new IllegalArgumentException();
0191:                if (y0 >= 1.0)
0192:                    throw new IllegalArgumentException();
0193:                code = 1;
0194:                y = y0;
0195:                if (y > (1.0 - 0.13533528323661269189)) { /* 0.135... = exp(-2) */
0196:                    y = 1.0 - y;
0197:                    code = 0;
0198:                }
0199:
0200:                if (y > 0.13533528323661269189) {
0201:                    y = y - 0.5;
0202:                    y2 = y * y;
0203:                    x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
0204:                    x = x * s2pi;
0205:                    return (x);
0206:                }
0207:
0208:                x = Math.sqrt(-2.0 * Math.log(y));
0209:                x0 = x - Math.log(x) / x;
0210:
0211:                z = 1.0 / x;
0212:                if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */
0213:                    x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
0214:                else
0215:                    x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
0216:                x = x0 - x1;
0217:                if (code != 0)
0218:                    x = -x;
0219:                return (x);
0220:            }
0221:
0222:            /**
0223:             * Returns natural logarithm of gamma function.
0224:             *
0225:             * @param x the value
0226:             * @return natural logarithm of gamma function
0227:             */
0228:            public static double lnGamma(double x) {
0229:
0230:                double p, q, w, z;
0231:
0232:                double A[] = { 8.11614167470508450300E-4,
0233:                        -5.95061904284301438324E-4, 7.93650340457716943945E-4,
0234:                        -2.77777777730099687205E-3, 8.33333333333331927722E-2 };
0235:                double B[] = { -1.37825152569120859100E3,
0236:                        -3.88016315134637840924E4, -3.31612992738871184744E5,
0237:                        -1.16237097492762307383E6, -1.72173700820839662146E6,
0238:                        -8.53555664245765465627E5 };
0239:                double C[] = {
0240:                /* 1.00000000000000000000E0, */
0241:                -3.51815701436523470549E2, -1.70642106651881159223E4,
0242:                        -2.20528590553854454839E5, -1.13933444367982507207E6,
0243:                        -2.53252307177582951285E6, -2.01889141433532773231E6 };
0244:
0245:                if (x < -34.0) {
0246:                    q = -x;
0247:                    w = lnGamma(q);
0248:                    p = Math.floor(q);
0249:                    if (p == q)
0250:                        throw new ArithmeticException("lnGamma: Overflow");
0251:                    z = q - p;
0252:                    if (z > 0.5) {
0253:                        p += 1.0;
0254:                        z = p - q;
0255:                    }
0256:                    z = q * Math.sin(Math.PI * z);
0257:                    if (z == 0.0)
0258:                        throw new ArithmeticException("lnGamma: Overflow");
0259:                    z = LOGPI - Math.log(z) - w;
0260:                    return z;
0261:                }
0262:
0263:                if (x < 13.0) {
0264:                    z = 1.0;
0265:                    while (x >= 3.0) {
0266:                        x -= 1.0;
0267:                        z *= x;
0268:                    }
0269:                    while (x < 2.0) {
0270:                        if (x == 0.0)
0271:                            throw new ArithmeticException("lnGamma: Overflow");
0272:                        z /= x;
0273:                        x += 1.0;
0274:                    }
0275:                    if (z < 0.0)
0276:                        z = -z;
0277:                    if (x == 2.0)
0278:                        return Math.log(z);
0279:                    x -= 2.0;
0280:                    p = x * polevl(x, B, 5) / p1evl(x, C, 6);
0281:                    return (Math.log(z) + p);
0282:                }
0283:
0284:                if (x > 2.556348e305)
0285:                    throw new ArithmeticException("lnGamma: Overflow");
0286:
0287:                q = (x - 0.5) * Math.log(x) - x + 0.91893853320467274178;
0288:
0289:                if (x > 1.0e8)
0290:                    return (q);
0291:
0292:                p = 1.0 / (x * x);
0293:                if (x >= 1000.0)
0294:                    q += ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3)
0295:                            * p + 0.0833333333333333333333)
0296:                            / x;
0297:                else
0298:                    q += polevl(p, A, 4) / x;
0299:                return q;
0300:            }
0301:
0302:            /**
0303:             * Returns the error function of the normal distribution.
0304:             * The integral is
0305:             * <pre>
0306:             *                           x 
0307:             *                            -
0308:             *                 2         | |          2
0309:             *   erf(x)  =  --------     |    exp( - t  ) dt.
0310:             *              sqrt(pi)   | |
0311:             *                          -
0312:             *                           0
0313:             * </pre>
0314:             * <b>Implementation:</b>
0315:             * For <tt>0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2)</tt>; otherwise
0316:             * <tt>erf(x) = 1 - erfc(x)</tt>.
0317:             * <p>
0318:             * Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
0319:             * Java 2D Graph Package 2.4</A>,
0320:             * which in turn is a port from the
0321:             * <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A>
0322:             * Math Library (C).
0323:             *
0324:             * @param a the argument to the function.
0325:             */
0326:            static double errorFunction(double x) {
0327:                double y, z;
0328:                final double T[] = { 9.60497373987051638749E0,
0329:                        9.00260197203842689217E1, 2.23200534594684319226E3,
0330:                        7.00332514112805075473E3, 5.55923013010394962768E4 };
0331:                final double U[] = {
0332:                        //1.00000000000000000000E0,
0333:                        3.35617141647503099647E1, 5.21357949780152679795E2,
0334:                        4.59432382970980127987E3, 2.26290000613890934246E4,
0335:                        4.92673942608635921086E4 };
0336:
0337:                if (Math.abs(x) > 1.0)
0338:                    return (1.0 - errorFunctionComplemented(x));
0339:                z = x * x;
0340:                y = x * polevl(z, T, 4) / p1evl(z, U, 5);
0341:                return y;
0342:            }
0343:
0344:            /**
0345:             * Returns the complementary Error function of the normal distribution.
0346:             * <pre>
0347:             *  1 - erf(x) =
0348:             *
0349:             *                           inf. 
0350:             *                             -
0351:             *                  2         | |          2
0352:             *   erfc(x)  =  --------     |    exp( - t  ) dt
0353:             *               sqrt(pi)   | |
0354:             *                           -
0355:             *                            x
0356:             * </pre>
0357:             * <b>Implementation:</b>
0358:             * For small x, <tt>erfc(x) = 1 - erf(x)</tt>; otherwise rational
0359:             * approximations are computed.
0360:             * <p>
0361:             * Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
0362:             * Java 2D Graph Package 2.4</A>,
0363:             * which in turn is a port from the
0364:             * <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A>
0365:             * Math Library (C).
0366:             *
0367:             * @param a the argument to the function.
0368:             */
0369:            static double errorFunctionComplemented(double a) {
0370:                double x, y, z, p, q;
0371:
0372:                double P[] = { 2.46196981473530512524E-10,
0373:                        5.64189564831068821977E-1, 7.46321056442269912687E0,
0374:                        4.86371970985681366614E1, 1.96520832956077098242E2,
0375:                        5.26445194995477358631E2, 9.34528527171957607540E2,
0376:                        1.02755188689515710272E3, 5.57535335369399327526E2 };
0377:                double Q[] = {
0378:                        //1.0
0379:                        1.32281951154744992508E1, 8.67072140885989742329E1,
0380:                        3.54937778887819891062E2, 9.75708501743205489753E2,
0381:                        1.82390916687909736289E3, 2.24633760818710981792E3,
0382:                        1.65666309194161350182E3, 5.57535340817727675546E2 };
0383:
0384:                double R[] = { 5.64189583547755073984E-1,
0385:                        1.27536670759978104416E0, 5.01905042251180477414E0,
0386:                        6.16021097993053585195E0, 7.40974269950448939160E0,
0387:                        2.97886665372100240670E0 };
0388:                double S[] = {
0389:                        //1.00000000000000000000E0, 
0390:                        2.26052863220117276590E0, 9.39603524938001434673E0,
0391:                        1.20489539808096656605E1, 1.70814450747565897222E1,
0392:                        9.60896809063285878198E0, 3.36907645100081516050E0 };
0393:
0394:                if (a < 0.0)
0395:                    x = -a;
0396:                else
0397:                    x = a;
0398:
0399:                if (x < 1.0)
0400:                    return 1.0 - errorFunction(a);
0401:
0402:                z = -a * a;
0403:
0404:                if (z < -MAXLOG) {
0405:                    if (a < 0)
0406:                        return (2.0);
0407:                    else
0408:                        return (0.0);
0409:                }
0410:
0411:                z = Math.exp(z);
0412:
0413:                if (x < 8.0) {
0414:                    p = polevl(x, P, 8);
0415:                    q = p1evl(x, Q, 8);
0416:                } else {
0417:                    p = polevl(x, R, 5);
0418:                    q = p1evl(x, S, 6);
0419:                }
0420:
0421:                y = (z * p) / q;
0422:
0423:                if (a < 0)
0424:                    y = 2.0 - y;
0425:
0426:                if (y == 0.0) {
0427:                    if (a < 0)
0428:                        return 2.0;
0429:                    else
0430:                        return (0.0);
0431:                }
0432:                return y;
0433:            }
0434:
0435:            /**
0436:             * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
0437:             * Evaluates polynomial when coefficient of N is 1.0.
0438:             * Otherwise same as <tt>polevl()</tt>.
0439:             * <pre>
0440:             *                     2          N
0441:             * y  =  C  + C x + C x  +...+ C x
0442:             *        0    1     2          N
0443:             *
0444:             * Coefficients are stored in reverse order:
0445:             *
0446:             * coef[0] = C  , ..., coef[N] = C  .
0447:             *            N                   0
0448:             * </pre>
0449:             * The function <tt>p1evl()</tt> assumes that <tt>coef[N] = 1.0</tt> and is
0450:             * omitted from the array.  Its calling arguments are
0451:             * otherwise the same as <tt>polevl()</tt>.
0452:             * <p>
0453:             * In the interest of speed, there are no checks for out of bounds arithmetic.
0454:             *
0455:             * @param x argument to the polynomial.
0456:             * @param coef the coefficients of the polynomial.
0457:             * @param N the degree of the polynomial.
0458:             */
0459:            static double p1evl(double x, double coef[], int N) {
0460:
0461:                double ans;
0462:                ans = x + coef[0];
0463:
0464:                for (int i = 1; i < N; i++)
0465:                    ans = ans * x + coef[i];
0466:
0467:                return ans;
0468:            }
0469:
0470:            /**
0471:             * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
0472:             * <pre>
0473:             *                     2          N
0474:             * y  =  C  + C x + C x  +...+ C x
0475:             *        0    1     2          N
0476:             *
0477:             * Coefficients are stored in reverse order:
0478:             *
0479:             * coef[0] = C  , ..., coef[N] = C  .
0480:             *            N                   0
0481:             * </pre>
0482:             * In the interest of speed, there are no checks for out of bounds arithmetic.
0483:             *
0484:             * @param x argument to the polynomial.
0485:             * @param coef the coefficients of the polynomial.
0486:             * @param N the degree of the polynomial.
0487:             */
0488:            static double polevl(double x, double coef[], int N) {
0489:
0490:                double ans;
0491:                ans = coef[0];
0492:
0493:                for (int i = 1; i <= N; i++)
0494:                    ans = ans * x + coef[i];
0495:
0496:                return ans;
0497:            }
0498:
0499:            /**
0500:             * Returns the Incomplete Gamma function.
0501:             * @param a the parameter of the gamma distribution.
0502:             * @param x the integration end point.
0503:             */
0504:            static double incompleteGamma(double a, double x) {
0505:
0506:                double ans, ax, c, r;
0507:
0508:                if (x <= 0 || a <= 0)
0509:                    return 0.0;
0510:
0511:                if (x > 1.0 && x > a)
0512:                    return 1.0 - incompleteGammaComplement(a, x);
0513:
0514:                /* Compute  x**a * exp(-x) / gamma(a)  */
0515:                ax = a * Math.log(x) - x - lnGamma(a);
0516:                if (ax < -MAXLOG)
0517:                    return (0.0);
0518:
0519:                ax = Math.exp(ax);
0520:
0521:                /* power series */
0522:                r = a;
0523:                c = 1.0;
0524:                ans = 1.0;
0525:
0526:                do {
0527:                    r += 1.0;
0528:                    c *= x / r;
0529:                    ans += c;
0530:                } while (c / ans > MACHEP);
0531:
0532:                return (ans * ax / a);
0533:            }
0534:
0535:            /**
0536:             * Returns the Complemented Incomplete Gamma function.
0537:             * @param a the parameter of the gamma distribution.
0538:             * @param x the integration start point.
0539:             */
0540:            static double incompleteGammaComplement(double a, double x) {
0541:
0542:                double ans, ax, c, yc, r, t, y, z;
0543:                double pk, pkm1, pkm2, qk, qkm1, qkm2;
0544:
0545:                if (x <= 0 || a <= 0)
0546:                    return 1.0;
0547:
0548:                if (x < 1.0 || x < a)
0549:                    return 1.0 - incompleteGamma(a, x);
0550:
0551:                ax = a * Math.log(x) - x - lnGamma(a);
0552:                if (ax < -MAXLOG)
0553:                    return 0.0;
0554:
0555:                ax = Math.exp(ax);
0556:
0557:                /* continued fraction */
0558:                y = 1.0 - a;
0559:                z = x + y + 1.0;
0560:                c = 0.0;
0561:                pkm2 = 1.0;
0562:                qkm2 = x;
0563:                pkm1 = x + 1.0;
0564:                qkm1 = z * x;
0565:                ans = pkm1 / qkm1;
0566:
0567:                do {
0568:                    c += 1.0;
0569:                    y += 1.0;
0570:                    z += 2.0;
0571:                    yc = y * c;
0572:                    pk = pkm1 * z - pkm2 * yc;
0573:                    qk = qkm1 * z - qkm2 * yc;
0574:                    if (qk != 0) {
0575:                        r = pk / qk;
0576:                        t = Math.abs((ans - r) / r);
0577:                        ans = r;
0578:                    } else
0579:                        t = 1.0;
0580:
0581:                    pkm2 = pkm1;
0582:                    pkm1 = pk;
0583:                    qkm2 = qkm1;
0584:                    qkm1 = qk;
0585:                    if (Math.abs(pk) > big) {
0586:                        pkm2 *= biginv;
0587:                        pkm1 *= biginv;
0588:                        qkm2 *= biginv;
0589:                        qkm1 *= biginv;
0590:                    }
0591:                } while (t > MACHEP);
0592:
0593:                return ans * ax;
0594:            }
0595:
0596:            /**
0597:             * Returns the Gamma function of the argument.
0598:             */
0599:            static double gamma(double x) {
0600:
0601:                double P[] = { 1.60119522476751861407E-4,
0602:                        1.19135147006586384913E-3, 1.04213797561761569935E-2,
0603:                        4.76367800457137231464E-2, 2.07448227648435975150E-1,
0604:                        4.94214826801497100753E-1, 9.99999999999999996796E-1 };
0605:                double Q[] = { -2.31581873324120129819E-5,
0606:                        5.39605580493303397842E-4, -4.45641913851797240494E-3,
0607:                        1.18139785222060435552E-2, 3.58236398605498653373E-2,
0608:                        -2.34591795718243348568E-1, 7.14304917030273074085E-2,
0609:                        1.00000000000000000320E0 };
0610:
0611:                double p, z;
0612:                int i;
0613:
0614:                double q = Math.abs(x);
0615:
0616:                if (q > 33.0) {
0617:                    if (x < 0.0) {
0618:                        p = Math.floor(q);
0619:                        if (p == q)
0620:                            throw new ArithmeticException("gamma: overflow");
0621:                        i = (int) p;
0622:                        z = q - p;
0623:                        if (z > 0.5) {
0624:                            p += 1.0;
0625:                            z = q - p;
0626:                        }
0627:                        z = q * Math.sin(Math.PI * z);
0628:                        if (z == 0.0)
0629:                            throw new ArithmeticException("gamma: overflow");
0630:                        z = Math.abs(z);
0631:                        z = Math.PI / (z * stirlingFormula(q));
0632:
0633:                        return -z;
0634:                    } else {
0635:                        return stirlingFormula(x);
0636:                    }
0637:                }
0638:
0639:                z = 1.0;
0640:                while (x >= 3.0) {
0641:                    x -= 1.0;
0642:                    z *= x;
0643:                }
0644:
0645:                while (x < 0.0) {
0646:                    if (x == 0.0) {
0647:                        throw new ArithmeticException("gamma: singular");
0648:                    } else if (x > -1.E-9) {
0649:                        return (z / ((1.0 + 0.5772156649015329 * x) * x));
0650:                    }
0651:                    z /= x;
0652:                    x += 1.0;
0653:                }
0654:
0655:                while (x < 2.0) {
0656:                    if (x == 0.0) {
0657:                        throw new ArithmeticException("gamma: singular");
0658:                    } else if (x < 1.e-9) {
0659:                        return (z / ((1.0 + 0.5772156649015329 * x) * x));
0660:                    }
0661:                    z /= x;
0662:                    x += 1.0;
0663:                }
0664:
0665:                if ((x == 2.0) || (x == 3.0))
0666:                    return z;
0667:
0668:                x -= 2.0;
0669:                p = polevl(x, P, 6);
0670:                q = polevl(x, Q, 7);
0671:                return z * p / q;
0672:            }
0673:
0674:            /**
0675:             * Returns the Gamma function computed by Stirling's formula.
0676:             * The polynomial STIR is valid for 33 <= x <= 172.
0677:             */
0678:            static double stirlingFormula(double x) {
0679:
0680:                double STIR[] = { 7.87311395793093628397E-4,
0681:                        -2.29549961613378126380E-4, -2.68132617805781232825E-3,
0682:                        3.47222221605458667310E-3, 8.33333333333482257126E-2, };
0683:                double MAXSTIR = 143.01608;
0684:
0685:                double w = 1.0 / x;
0686:                double y = Math.exp(x);
0687:
0688:                w = 1.0 + w * polevl(w, STIR, 4);
0689:
0690:                if (x > MAXSTIR) {
0691:                    /* Avoid overflow in Math.pow() */
0692:                    double v = Math.pow(x, 0.5 * x - 0.25);
0693:                    y = v * (v / y);
0694:                } else {
0695:                    y = Math.pow(x, x - 0.5) / y;
0696:                }
0697:                y = SQTPI * y * w;
0698:                return y;
0699:            }
0700:
0701:            /**
0702:             * Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>.
0703:             *
0704:             * @param aa the alpha parameter of the beta distribution.
0705:             * @param bb the beta parameter of the beta distribution.
0706:             * @param xx the integration end point.
0707:             */
0708:            public static double incompleteBeta(double aa, double bb, double xx) {
0709:
0710:                double a, b, t, x, xc, w, y;
0711:                boolean flag;
0712:
0713:                if (aa <= 0.0 || bb <= 0.0)
0714:                    throw new ArithmeticException("ibeta: Domain error!");
0715:
0716:                if ((xx <= 0.0) || (xx >= 1.0)) {
0717:                    if (xx == 0.0)
0718:                        return 0.0;
0719:                    if (xx == 1.0)
0720:                        return 1.0;
0721:                    throw new ArithmeticException("ibeta: Domain error!");
0722:                }
0723:
0724:                flag = false;
0725:                if ((bb * xx) <= 1.0 && xx <= 0.95) {
0726:                    t = powerSeries(aa, bb, xx);
0727:                    return t;
0728:                }
0729:
0730:                w = 1.0 - xx;
0731:
0732:                /* Reverse a and b if x is greater than the mean. */
0733:                if (xx > (aa / (aa + bb))) {
0734:                    flag = true;
0735:                    a = bb;
0736:                    b = aa;
0737:                    xc = xx;
0738:                    x = w;
0739:                } else {
0740:                    a = aa;
0741:                    b = bb;
0742:                    xc = w;
0743:                    x = xx;
0744:                }
0745:
0746:                if (flag && (b * x) <= 1.0 && x <= 0.95) {
0747:                    t = powerSeries(a, b, x);
0748:                    if (t <= MACHEP)
0749:                        t = 1.0 - MACHEP;
0750:                    else
0751:                        t = 1.0 - t;
0752:                    return t;
0753:                }
0754:
0755:                /* Choose expansion for better convergence. */
0756:                y = x * (a + b - 2.0) - (a - 1.0);
0757:                if (y < 0.0)
0758:                    w = incompleteBetaFraction1(a, b, x);
0759:                else
0760:                    w = incompleteBetaFraction2(a, b, x) / xc;
0761:
0762:                /* Multiply w by the factor
0763:                   a      b   _             _     _
0764:                   x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */
0765:
0766:                y = a * Math.log(x);
0767:                t = b * Math.log(xc);
0768:                if ((a + b) < MAXGAM && Math.abs(y) < MAXLOG
0769:                        && Math.abs(t) < MAXLOG) {
0770:                    t = Math.pow(xc, b);
0771:                    t *= Math.pow(x, a);
0772:                    t /= a;
0773:                    t *= w;
0774:                    t *= gamma(a + b) / (gamma(a) * gamma(b));
0775:                    if (flag) {
0776:                        if (t <= MACHEP)
0777:                            t = 1.0 - MACHEP;
0778:                        else
0779:                            t = 1.0 - t;
0780:                    }
0781:                    return t;
0782:                }
0783:                /* Resort to logarithms.  */
0784:                y += t + lnGamma(a + b) - lnGamma(a) - lnGamma(b);
0785:                y += Math.log(w / a);
0786:                if (y < MINLOG)
0787:                    t = 0.0;
0788:                else
0789:                    t = Math.exp(y);
0790:
0791:                if (flag) {
0792:                    if (t <= MACHEP)
0793:                        t = 1.0 - MACHEP;
0794:                    else
0795:                        t = 1.0 - t;
0796:                }
0797:                return t;
0798:            }
0799:
0800:            /**
0801:             * Continued fraction expansion #1 for incomplete beta integral.
0802:             */
0803:            static double incompleteBetaFraction1(double a, double b, double x) {
0804:
0805:                double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
0806:                double k1, k2, k3, k4, k5, k6, k7, k8;
0807:                double r, t, ans, thresh;
0808:                int n;
0809:
0810:                k1 = a;
0811:                k2 = a + b;
0812:                k3 = a;
0813:                k4 = a + 1.0;
0814:                k5 = 1.0;
0815:                k6 = b - 1.0;
0816:                k7 = k4;
0817:                k8 = a + 2.0;
0818:
0819:                pkm2 = 0.0;
0820:                qkm2 = 1.0;
0821:                pkm1 = 1.0;
0822:                qkm1 = 1.0;
0823:                ans = 1.0;
0824:                r = 1.0;
0825:                n = 0;
0826:                thresh = 3.0 * MACHEP;
0827:                do {
0828:                    xk = -(x * k1 * k2) / (k3 * k4);
0829:                    pk = pkm1 + pkm2 * xk;
0830:                    qk = qkm1 + qkm2 * xk;
0831:                    pkm2 = pkm1;
0832:                    pkm1 = pk;
0833:                    qkm2 = qkm1;
0834:                    qkm1 = qk;
0835:
0836:                    xk = (x * k5 * k6) / (k7 * k8);
0837:                    pk = pkm1 + pkm2 * xk;
0838:                    qk = qkm1 + qkm2 * xk;
0839:                    pkm2 = pkm1;
0840:                    pkm1 = pk;
0841:                    qkm2 = qkm1;
0842:                    qkm1 = qk;
0843:
0844:                    if (qk != 0)
0845:                        r = pk / qk;
0846:                    if (r != 0) {
0847:                        t = Math.abs((ans - r) / r);
0848:                        ans = r;
0849:                    } else
0850:                        t = 1.0;
0851:
0852:                    if (t < thresh)
0853:                        return ans;
0854:
0855:                    k1 += 1.0;
0856:                    k2 += 1.0;
0857:                    k3 += 2.0;
0858:                    k4 += 2.0;
0859:                    k5 += 1.0;
0860:                    k6 -= 1.0;
0861:                    k7 += 2.0;
0862:                    k8 += 2.0;
0863:
0864:                    if ((Math.abs(qk) + Math.abs(pk)) > big) {
0865:                        pkm2 *= biginv;
0866:                        pkm1 *= biginv;
0867:                        qkm2 *= biginv;
0868:                        qkm1 *= biginv;
0869:                    }
0870:                    if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
0871:                        pkm2 *= big;
0872:                        pkm1 *= big;
0873:                        qkm2 *= big;
0874:                        qkm1 *= big;
0875:                    }
0876:                } while (++n < 300);
0877:
0878:                return ans;
0879:            }
0880:
0881:            /**
0882:             * Continued fraction expansion #2 for incomplete beta integral.
0883:             */
0884:            static double incompleteBetaFraction2(double a, double b, double x) {
0885:
0886:                double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
0887:                double k1, k2, k3, k4, k5, k6, k7, k8;
0888:                double r, t, ans, z, thresh;
0889:                int n;
0890:
0891:                k1 = a;
0892:                k2 = b - 1.0;
0893:                k3 = a;
0894:                k4 = a + 1.0;
0895:                k5 = 1.0;
0896:                k6 = a + b;
0897:                k7 = a + 1.0;
0898:                ;
0899:                k8 = a + 2.0;
0900:
0901:                pkm2 = 0.0;
0902:                qkm2 = 1.0;
0903:                pkm1 = 1.0;
0904:                qkm1 = 1.0;
0905:                z = x / (1.0 - x);
0906:                ans = 1.0;
0907:                r = 1.0;
0908:                n = 0;
0909:                thresh = 3.0 * MACHEP;
0910:                do {
0911:                    xk = -(z * k1 * k2) / (k3 * k4);
0912:                    pk = pkm1 + pkm2 * xk;
0913:                    qk = qkm1 + qkm2 * xk;
0914:                    pkm2 = pkm1;
0915:                    pkm1 = pk;
0916:                    qkm2 = qkm1;
0917:                    qkm1 = qk;
0918:
0919:                    xk = (z * k5 * k6) / (k7 * k8);
0920:                    pk = pkm1 + pkm2 * xk;
0921:                    qk = qkm1 + qkm2 * xk;
0922:                    pkm2 = pkm1;
0923:                    pkm1 = pk;
0924:                    qkm2 = qkm1;
0925:                    qkm1 = qk;
0926:
0927:                    if (qk != 0)
0928:                        r = pk / qk;
0929:                    if (r != 0) {
0930:                        t = Math.abs((ans - r) / r);
0931:                        ans = r;
0932:                    } else
0933:                        t = 1.0;
0934:
0935:                    if (t < thresh)
0936:                        return ans;
0937:
0938:                    k1 += 1.0;
0939:                    k2 -= 1.0;
0940:                    k3 += 2.0;
0941:                    k4 += 2.0;
0942:                    k5 += 1.0;
0943:                    k6 += 1.0;
0944:                    k7 += 2.0;
0945:                    k8 += 2.0;
0946:
0947:                    if ((Math.abs(qk) + Math.abs(pk)) > big) {
0948:                        pkm2 *= biginv;
0949:                        pkm1 *= biginv;
0950:                        qkm2 *= biginv;
0951:                        qkm1 *= biginv;
0952:                    }
0953:                    if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
0954:                        pkm2 *= big;
0955:                        pkm1 *= big;
0956:                        qkm2 *= big;
0957:                        qkm1 *= big;
0958:                    }
0959:                } while (++n < 300);
0960:
0961:                return ans;
0962:            }
0963:
0964:            /**
0965:             * Power series for incomplete beta integral.
0966:             * Use when b*x is small and x not too close to 1.  
0967:             */
0968:            static double powerSeries(double a, double b, double x) {
0969:
0970:                double s, t, u, v, n, t1, z, ai;
0971:
0972:                ai = 1.0 / a;
0973:                u = (1.0 - b) * x;
0974:                v = u / (a + 1.0);
0975:                t1 = v;
0976:                t = u;
0977:                n = 2.0;
0978:                s = 0.0;
0979:                z = MACHEP * ai;
0980:                while (Math.abs(v) > z) {
0981:                    u = (n - b) * x / n;
0982:                    t *= u;
0983:                    v = t / (a + n);
0984:                    s += v;
0985:                    n += 1.0;
0986:                }
0987:                s += t1;
0988:                s += ai;
0989:
0990:                u = a * Math.log(x);
0991:                if ((a + b) < MAXGAM && Math.abs(u) < MAXLOG) {
0992:                    t = gamma(a + b) / (gamma(a) * gamma(b));
0993:                    s = s * t * Math.pow(x, a);
0994:                } else {
0995:                    t = lnGamma(a + b) - lnGamma(a) - lnGamma(b) + u
0996:                            + Math.log(s);
0997:                    if (t < MINLOG)
0998:                        s = 0.0;
0999:                    else
1000:                        s = Math.exp(t);
1001:                }
1002:                return s;
1003:            }
1004:
1005:            /**
1006:             * Main method for testing this class.
1007:             */
1008:            public static void main(String[] ops) {
1009:
1010:                System.out.println("Binomial standard error (0.5, 100): "
1011:                        + Statistics.binomialStandardError(0.5, 100));
1012:                System.out.println("Chi-squared probability (2.558, 10): "
1013:                        + Statistics.chiSquaredProbability(2.558, 10));
1014:                System.out.println("Normal probability (0.2): "
1015:                        + Statistics.normalProbability(0.2));
1016:                System.out.println("F probability (5.1922, 4, 5): "
1017:                        + Statistics.FProbability(5.1922, 4, 5));
1018:                System.out.println("lnGamma(6): " + Statistics.lnGamma(6));
1019:            }
1020:        }
www.java2java.com | Contact Us
Copyright 2009 - 12 Demo Source and Support. All rights reserved.
All other trademarks are property of their respective owners.