← Back to branch summary

~frgomes/jquantlib/trunk

« back to all changes in this revision

Viewing changes to jquantlib-experimental/src/main/java/org/jquantlib/experimental/math/distributions/TESTICN.java

  • Committer: Richard Gomes
  • Date: 2015-03-11 00:40:46 UTC
  • Revision ID: rgomes.info@gmail.com-20150311004046-7pck4lxbmc2iy4gj
Remove jquantlib-experimental.

Show diffs side-by-side

added added

removed removed

Lines of Context:
jquantlib-experimental/src/main/java/org/jquantlib/experimental/math/distributions/TESTICN.java
1
 
/*
2
 
 Copyright (C) 2008 Richard Gomes
3
 
 
4
 
 This source code is release under the BSD License.
5
 
 
6
 
 This file is part of JQuantLib, a free-software/open-source library
7
 
 for financial quantitative analysts and developers - http://jquantlib.org/
8
 
 
9
 
 JQuantLib is free software: you can redistribute it and/or modify it
10
 
 under the terms of the JQuantLib license.  You should have received a
11
 
 copy of the license along with this program; if not, please email
12
 
 <jquant-devel@lists.sourceforge.net>. The license is also available online at
13
 
 <http://www.jquantlib.org/index.php/LICENSE.TXT>.
14
 
 
15
 
 This program is distributed in the hope that it will be useful, but WITHOUT
16
 
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17
 
 FOR A PARTICULAR PURPOSE.  See the license for more details.
18
 
 
19
 
 JQuantLib is based on QuantLib. http://quantlib.org/
20
 
 When applicable, the original copyright notice follows this notice.
21
 
 */
22
 
 
23
 
package org.jquantlib.experimental.math.distributions;
24
 
 
25
 
/**
26
 
 *  CODE FROM  http://home.online.no/~pjacklam/notes/invnorm/impl/karimov/StatUtil.java
27
 
 * @author Q. Boiler
28
 
 */
29
 
public class TESTICN {
30
 
/**
31
 
 * Class contains the implementation of:
32
 
 * - Inverse Normal Cummulative Distribution Function Algorythm
33
 
 * - Error Function Algorythm
34
 
 * - Complimentary Error Function Algorythm
35
 
 *
36
 
 * @author Sherali Karimov (sherali.karimov@proxima-tech.com)
37
 
 */
38
 
/* ********************************************
39
 
 * Original algorythm and Perl implementation can
40
 
 * be found at:
41
 
 * http://www.math.uio.no/~jacklam/notes/invnorm/index.html
42
 
 * Author:
43
 
 *  Peter J. Acklam
44
 
 *  jacklam@math.uio.no
45
 
 * ****************************************** */
46
 
  private static final double P_LOW  = 0.02425D;
47
 
  private static final double P_HIGH = 1.0D - P_LOW;
48
 
 
49
 
  // Coefficients in rational approximations.
50
 
  private static final double ICDF_A[] =
51
 
  { -3.969683028665376e+01,  2.209460984245205e+02,
52
 
    -2.759285104469687e+02,  1.383577518672690e+02,
53
 
    -3.066479806614716e+01,  2.506628277459239e+00 };
54
 
 
55
 
  private static final double ICDF_B[] =
56
 
  { -5.447609879822406e+01,  1.615858368580409e+02,
57
 
    -1.556989798598866e+02,  6.680131188771972e+01,
58
 
    -1.328068155288572e+01 };
59
 
 
60
 
  private static final double ICDF_C[] =
61
 
  { -7.784894002430293e-03, -3.223964580411365e-01,
62
 
    -2.400758277161838e+00, -2.549732539343734e+00,
63
 
    4.374664141464968e+00,  2.938163982698783e+00 };
64
 
 
65
 
  private static final double ICDF_D[] =
66
 
  { 7.784695709041462e-03,  3.224671290700398e-01,
67
 
    2.445134137142996e+00,  3.754408661907416e+00 };
68
 
 
69
 
  public static double getInvCDF(double d, boolean highPrecision)
70
 
  {
71
 
    // Define break-points.
72
 
    // variable for result
73
 
    double z = 0;
74
 
 
75
 
    if(d == 0) z = Double.NEGATIVE_INFINITY;
76
 
    else if(d == 1) z = Double.POSITIVE_INFINITY;
77
 
    else if(Double.isNaN(d) || d < 0 || d > 1) z = Double.NaN;
78
 
 
79
 
    // Rational approximation for lower region:
80
 
    else if( d < P_LOW )
81
 
    {
82
 
      double q  = Math.sqrt(-2*Math.log(d));
83
 
      z = (((((ICDF_C[0]*q+ICDF_C[1])*q+ICDF_C[2])*q+ICDF_C[3])*q+ICDF_C[4])*q+ICDF_C[5]) / ((((ICDF_D[0]*q+ICDF_D[1])*q+ICDF_D[2])*q+ICDF_D[3])*q+1);
84
 
    }
85
 
 
86
 
    // Rational approximation for upper region:
87
 
    else if ( P_HIGH < d )
88
 
    {
89
 
      double q  = Math.sqrt(-2*Math.log(1-d));
90
 
      z = -(((((ICDF_C[0]*q+ICDF_C[1])*q+ICDF_C[2])*q+ICDF_C[3])*q+ICDF_C[4])*q+ICDF_C[5]) / ((((ICDF_D[0]*q+ICDF_D[1])*q+ICDF_D[2])*q+ICDF_D[3])*q+1);
91
 
    }
92
 
   // Rational approximation for central region:
93
 
    else
94
 
    {
95
 
      double q = d - 0.5D;
96
 
      double r = q * q;
97
 
      z = (((((ICDF_A[0]*r+ICDF_A[1])*r+ICDF_A[2])*r+ICDF_A[3])*r+ICDF_A[4])*r+ICDF_A[5])*q / (((((ICDF_B[0]*r+ICDF_B[1])*r+ICDF_B[2])*r+ICDF_B[3])*r+ICDF_B[4])*r+1);
98
 
    }
99
 
    if(highPrecision) z = refine(z, d);
100
 
    return z;
101
 
  }
102
 
 
103
 
//C------------------------------------------------------------------
104
 
//C  Coefficients for approximation to  erf  in first interval
105
 
//C------------------------------------------------------------------
106
 
  private static final double ERF_A[] =
107
 
  { 3.16112374387056560E00, 1.13864154151050156E02,
108
 
    3.77485237685302021E02, 3.20937758913846947E03,
109
 
    1.85777706184603153E-1 };
110
 
 
111
 
  private static final double ERF_B[] =
112
 
  { 2.36012909523441209E01, 2.44024637934444173E02,
113
 
    1.28261652607737228E03, 2.84423683343917062E03 };
114
 
 
115
 
//C------------------------------------------------------------------
116
 
//C  Coefficients for approximation to  erfc  in second interval
117
 
//C------------------------------------------------------------------
118
 
  private static final double ERF_C[] =
119
 
  { 5.64188496988670089E-1, 8.88314979438837594E0,
120
 
    6.61191906371416295E01, 2.98635138197400131E02,
121
 
    8.81952221241769090E02, 1.71204761263407058E03,
122
 
    2.05107837782607147E03, 1.23033935479799725E03,
123
 
    2.15311535474403846E-8 };
124
 
 
125
 
  private static final double ERF_D[] =
126
 
  { 1.57449261107098347E01,1.17693950891312499E02,
127
 
    5.37181101862009858E02,1.62138957456669019E03,
128
 
    3.29079923573345963E03,4.36261909014324716E03,
129
 
    3.43936767414372164E03,1.23033935480374942E03 };
130
 
 
131
 
//C------------------------------------------------------------------
132
 
//C  Coefficients for approximation to  erfc  in third interval
133
 
//C------------------------------------------------------------------
134
 
  private static final double ERF_P[] =
135
 
  { 3.05326634961232344E-1,3.60344899949804439E-1,
136
 
    1.25781726111229246E-1,1.60837851487422766E-2,
137
 
    6.58749161529837803E-4,1.63153871373020978E-2 };
138
 
 
139
 
  private static final double ERF_Q[] =
140
 
  { 2.56852019228982242E00,1.87295284992346047E00,
141
 
    5.27905102951428412E-1,6.05183413124413191E-2,
142
 
    2.33520497626869185E-3 };
143
 
 
144
 
  private static final double PI_SQRT = Math.sqrt(Math.PI);
145
 
  private static final double THRESHOLD = 0.46875D;
146
 
 
147
 
/* **************************************
148
 
 * Hardware dependant constants were calculated
149
 
 * on Dell "Dimension 4100":
150
 
 * - Pentium III 800 MHz
151
 
 * running Microsoft Windows 2000
152
 
 * ************************************* */
153
 
  private static final double X_MIN = Double.MIN_VALUE;
154
 
  private static final double X_INF = Double.MAX_VALUE;
155
 
  private static final double X_NEG = -9.38241396824444;
156
 
  private static final double X_SMALL = 1.110223024625156663E-16;
157
 
  private static final double X_BIG = 9.194E0;
158
 
  private static final double X_HUGE = 1.0D / (2.0D * Math.sqrt(X_SMALL));
159
 
  private static final double X_MAX = Math.min(X_INF, (1 / (Math.sqrt(Math.PI) * X_MIN)));
160
 
 
161
 
  private static double calerf(double X, int jint)
162
 
  {
163
 
/* ******************************************
164
 
 * ORIGINAL FORTRAN version can be found at:
165
 
 * http://www.netlib.org/specfun/erf
166
 
 ********************************************
167
 
C------------------------------------------------------------------
168
 
C
169
 
C THIS PACKET COMPUTES THE ERROR AND COMPLEMENTARY ERROR FUNCTIONS
170
 
C   FOR REAL ARGUMENTS  ARG.  IT CONTAINS TWO FUNCTION TYPE
171
 
C   SUBPROGRAMS,  ERF  AND  ERFC  (OR  DERF  AND  DERFC),  AND ONE
172
 
C   SUBROUTINE TYPE SUBPROGRAM,  CALERF.  THE CALLING STATEMENTS
173
 
C   FOR THE PRIMARY ENTRIES ARE
174
 
C
175
 
C                   Y=ERF(X)     (OR   Y=DERF(X) )
176
 
C   AND
177
 
C                   Y=ERFC(X)    (OR   Y=DERFC(X) ).
178
 
C
179
 
C   THE ROUTINE  CALERF  IS INTENDED FOR INTERNAL PACKET USE ONLY,
180
 
C   ALL COMPUTATIONS WITHIN THE PACKET BEING CONCENTRATED IN THIS
181
 
C   ROUTINE.  THE FUNCTION SUBPROGRAMS INVOKE  CALERF  WITH THE
182
 
C   STATEMENT
183
 
C          CALL CALERF(ARG,RESULT,JINT)
184
 
C   WHERE THE PARAMETER USAGE IS AS FOLLOWS
185
 
C
186
 
C      FUNCTION                     PARAMETERS FOR CALERF
187
 
C       CALL              ARG                  RESULT          JINT
188
 
C     ERF(ARG)      ANY REAL ARGUMENT         ERF(ARG)          0
189
 
C     ERFC(ARG)     ABS(ARG) .LT. XMAX        ERFC(ARG)         1
190
 
C
191
 
C   THE MAIN COMPUTATION EVALUATES NEAR MINIMAX APPROXIMATIONS
192
 
C   FROM "RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE ERROR FUNCTION"
193
 
C   BY W. J. CODY, MATH. COMP., 1969, PP. 631-638.  THIS
194
 
C   TRANSPORTABLE PROGRAM USES RATIONAL FUNCTIONS THAT THEORETICALLY
195
 
C       APPROXIMATE  ERF(X)  AND  ERFC(X)  TO AT LEAST 18 SIGNIFICANT
196
 
C   DECIMAL DIGITS.  THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC
197
 
C   SYSTEM, THE COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER
198
 
C   SELECTION OF THE MACHINE-DEPENDENT CONSTANTS.
199
 
C
200
 
C  AUTHOR: W. J. CODY
201
 
C          MATHEMATICS AND COMPUTER SCIENCE DIVISION
202
 
C          ARGONNE NATIONAL LABORATORY
203
 
C          ARGONNE, IL 60439
204
 
C
205
 
C  LATEST MODIFICATION: JANUARY 8, 1985
206
 
C
207
 
C------------------------------------------------------------------
208
 
*/
209
 
    double result = 0;
210
 
    double Y = Math.abs(X);
211
 
    double Y_SQ, X_NUM, X_DEN;
212
 
 
213
 
    if(Y <= THRESHOLD)
214
 
    {
215
 
      Y_SQ = 0.0D;
216
 
      if (Y > X_SMALL) Y_SQ = Y * Y;
217
 
      X_NUM = ERF_A[4] * Y_SQ;
218
 
      X_DEN = Y_SQ;
219
 
      for (int i=0;i<3;i++)
220
 
      {
221
 
        X_NUM = (X_NUM + ERF_A[i]) * Y_SQ;
222
 
        X_DEN = (X_DEN + ERF_B[i]) * Y_SQ;
223
 
      }
224
 
      result = X * (X_NUM + ERF_A[3]) / (X_DEN + ERF_B[3]);
225
 
      if (jint != 0)  result = 1 - result;
226
 
      if(jint == 2) result = Math.exp(Y_SQ) * result;
227
 
      return result;
228
 
    }
229
 
    else if (Y <= 4.0D)
230
 
    {
231
 
      X_NUM = ERF_C[8] * Y;
232
 
      X_DEN = Y;
233
 
      for (int i=0;i<7;i++)
234
 
      {
235
 
        X_NUM = (X_NUM + ERF_C[i]) * Y;
236
 
        X_DEN = (X_DEN + ERF_D[i]) * Y;
237
 
      }
238
 
      result = (X_NUM + ERF_C[7]) / (X_DEN + ERF_D[7]);
239
 
      if(jint != 2)
240
 
      {
241
 
        Y_SQ = Math.round(Y*16.0D)/16.0D;
242
 
        double del = (Y - Y_SQ) * (Y + Y_SQ);
243
 
        result = Math.exp(-Y_SQ * Y_SQ) * Math.exp(-del) * result;
244
 
      }
245
 
    }
246
 
    else
247
 
    {
248
 
      result = 0.0D;
249
 
      if( Y >= X_BIG && (jint != 2 || Y >= X_MAX));
250
 
      else if(Y >= X_BIG && Y >= X_HUGE) result = PI_SQRT / Y;
251
 
      else
252
 
      {
253
 
        Y_SQ = 1.0D / (Y * Y);
254
 
        X_NUM = ERF_P[5] * Y_SQ;
255
 
        X_DEN = Y_SQ;
256
 
        for (int i=0;i<4; i++)
257
 
        {
258
 
          X_NUM = (X_NUM + ERF_P[i]) * Y_SQ;
259
 
          X_DEN = (X_DEN + ERF_Q[i]) * Y_SQ;
260
 
        }
261
 
        result = Y_SQ * (X_NUM + ERF_P[4]) / (X_DEN + ERF_Q[4]);
262
 
        result = (PI_SQRT - result) / Y;
263
 
        if(jint != 2)
264
 
        {
265
 
          Y_SQ = Math.round(Y*16.0D)/16.0D;
266
 
          double del = (Y - Y_SQ) * (Y + Y_SQ);
267
 
          result = Math.exp(-Y_SQ * Y_SQ) * Math.exp(-del) * result;
268
 
        }
269
 
      }
270
 
    }
271
 
 
272
 
    if(jint == 0)
273
 
    {
274
 
      result = (0.5D - result) + 0.5D;
275
 
      if(X < 0) result = -result;
276
 
    }
277
 
    else if(jint == 1)
278
 
    {
279
 
      if(X < 0) result = 2.0D - result;
280
 
    }
281
 
    else
282
 
    {
283
 
      if(X < 0)
284
 
      {
285
 
        if(X < X_NEG) result = X_INF;
286
 
        else
287
 
        {
288
 
          Y_SQ = Math.round(X*16.0D)/16.0D;
289
 
          double del = (X - Y_SQ) * (X + Y_SQ);
290
 
          Y = Math.exp(Y_SQ * Y_SQ) * Math.exp(del);
291
 
          result = (Y + Y) - result;
292
 
        }
293
 
      }
294
 
    }
295
 
    return result;
296
 
  }
297
 
 
298
 
  public static double erf(double d){ return calerf(d, 0); }
299
 
  public static double erfc(double d){ return calerf(d, 1); }
300
 
  public static double erfcx(double d){ return calerf(d, 2); }
301
 
 
302
 
/* ****************************************************
303
 
 * Refining algorytm is based on Halley rational method
304
 
 * for finding roots of equations as described at:
305
 
 * http://www.math.uio.no/~jacklam/notes/invnorm/index.html
306
 
 * by:
307
 
 *  Peter J. Acklam
308
 
 *  jacklam@math.uio.no
309
 
 * ************************************************** */
310
 
  public static double refine(double x, double d)
311
 
  {
312
 
    if( d > 0 && d < 1)
313
 
    {
314
 
      double e = 0.5D * erfc(-x/Math.sqrt(2.0D)) - d;
315
 
      double u = e * Math.sqrt(2.0D*Math.PI) * Math.exp((x*x)/2.0D);
316
 
      x = x - u/(1.0D + x*u/2.0D);
317
 
    }
318
 
    return x;
319
 
  }
320
 
}