EXPLANATION FILE OF PROGRAM INVNORM =================================== Rational Polynomials -------------------- Our attention th us far has been largely directed at polynomial approxima- tions of the form N i f(x) = P(x) = Sum a x (2.7.1) i=0 i This the most popular type of approximation partly because it usuaIly works weIl (after some change of variables perhaps), and partly because very simple and effective techniques exist for determining the ai. However, there are ma- ny cases in which this structure is not appropriate. For example, a function that has an infinity at some value of x, say xo. We might approach this pro- blem by replacing the variable x in equation (2.7.1) with 1/(x - xo), thus creating a Laurent series. However, in some cases, this may also not be appropriate. There is another polynomial approximation form that is often very useful for dealing with such curve shapes, the rational polynomial: M i R(x) = P(x)/ Q(x) where Q(x) = Sum b x (2.7.2) i=0 i We will not go into the theory behind the rational polynomial or the ways in which it can be derived. Instead, a few simple rules of thumb along with some examples will be presented. For a good discussion regarding rational polyno- mials, see Ref. 15. We will define the foIlowing four parameters: 1) M(P): The lowest degree appearing in P(x). For example, for P(x) = 6x^3 + 7x^4 + 2x^5, M(P) = 3. 2) N(P): The highest degree appearing in P(x). For the example above, N(P) = 5. 3) M(Q): The lowest degree appearing in Q(x). 4) N(Q): The highest degree appearing in Q(x). The relative values of these parameters define the asymptotic behavior of P(x)/Q(x). By examining the behavior of the function to be approximated at x = 0 and x --> inf., strong clues as to the form of P(x) and Q(x) can be obtained. For example, if f(x) = constant at x = 0, then M(P) = M(Q). If instead, f(x) --> inf. near x = 0, then M(P) < M(Q). If f(x) = 0 at x = 0, then M(P) > M(Q). There are also three general types of behavior that exist as x becomes infinite: f(x) --> 0 [N(P) < N(Q)]; f(x) = constant [N(P) = N(Q)]; and f(x) --> inf. [N(P) > N(Q)]. By examining the shape of f(x), we can usual- ly get an idea of wether or not a rational polynomial is caIled for, and we can even determine some of the relative properties of P(x) and Q(x). As a numerical example, we will consider an approximation related to the inverse of the complementary error function. We will define Q(x) as inf. Q(x) = (1/sqrt(2pi) Sum exp(-t^2/2) dt (2.7.3) x The normal distribution probability density function is p(x) = (1/sqrt(2pi) exp(-t^2/2 (2.7.4) Therefore, Q(x) is one minus the cumulative normal distribution function: Q(x) = 1 - P(x). Also. Q(x) is related to the complementary error function: Q(x) = 1/2 erfc(x). An approximation to the complementary error function was presented earlier. We will now in effect approximate its inverse. The object is to de termine the value of x that corresponds with given Q. Abramowitz and Stegun give the following rational polynomial approximation for Xo (0 < Q <= 0.5): c0 + c1 t + c2 t^2 Xo = t - ------------------------ 1 + d1 t + d2 t^2 + d3 t where t = sqrt(ln l/Q^2) c0 = 2.515517 c1 = 0.802853 c2 = 0.010328 d1 = 1.432788 d2 = 0.189269 d3 = 0.001308 The error in this approximation is |E(Q)| < 0.0005 (Note that the error is referenced to Q, not Xo. The error in Xo is relative to how much Q would have to be altered to give a corresponding change in X). A little ingenuity was applied in choosing the functional form for t. The choice was based on the asymptotic form relating Xo and Q as Q --> O. A program for applying this approximation is shown in INVNORM. The approxi- mation associated with that function also involved a logarithmic change of variables. Other rational polynomial approximations are available From the literature. A few are given below: Type I cos x for -1 <= x <= 1 2 4 6 1 + a2 x + a4 x + a6 x cos x = ---------------------- (Ref. 15) 2 4 6 1 + b2 x + b4 x + b6 x |error| <= 2 X 10-11 a2 = - 0.470 595 788 392 a4 = 0.027 388 289 676 a6 = - 0.000 372 342 269 b2 = 0.029 404 211 608 b4 = 0.000 423 728 814 b6 = 0.000 003 235 543 Type II (Ref. 6) erf(x) for 0 <= x < inf. 1 erf(x) = 1 - ----------------------------------------------- 1 + b1x + b2x^2 + b3x^3 + b4x^4 + b5x^5 + b6x^6 |error| <= 3 X 10-7 b1 = 0.07052 30784 b2 = 0.00927 05272 b3 = 0.00027 65672 b4 = 0.04228 20123 b5 = 0.00015 20143 b6 = 0.00004 30638 Type III tanh µx for -1 <= x <= 1 µ = 1/2 ln 3 3 a1 x + a3 x tanh µx = ----------------- (Ref. 15) 2 4 1 + b2 x + b4 x |error| <= 6 X 10-9 al = 0.549306 14401 a3 = 0.01574 011995 b2 = 0.12923360954 b4 = 0.00085 891904 Although rational polynomials offer an interesting flexibility in functional aproximation, they do not appear to be very popular. This is probably due to the dificulty that exists in generating values for the coefficients. [15] Fike, C.T. Computer Evaluation of Mathematical Functions,Englewood Cliffs, NJ: Prentice-Hall, 1968. From [BIBLMI 01] -------------------------------------------------------- End of file invnorm.txt