This document contains a fit to AiGauss(z) = erfc(z)/2 with an error ε < 10-15 × AiGauss(z). #BestFit
The code is in ..\..\zips\FittingBrack.zip. #ZipFile is an annotated directory of this file.
The function is defined by
(1.1)
This is approximated as
Where
(1.3)
The form in (1.2) has 15 digit accuracy for the exponential and also for the Pade approximate. This is better than the form Pade = log( AiGau(z)) which has 15 digit accuracy in the log, but uses two of these digits at z = -6 where the lot is approximately -36 to simply represent the order of magnitude of the function.
The first fit showed errors in the integral version of AiGauss itself ../aigau15/Fitting%20Brack.doc .htm
The code for producing AiGauss data is saved here as ..\..\zips\genint.zip Compiling was done with
Wfl386 tgleg67
This produced a welcome extra digit of accuracy beyond that available with Watfor.
Aigau2.dir
AIGAU.DIR, OUTPUT FILE NAME
..\TestData\BRACK.OUT, DATA SET NAME
1, DIMENSION OF X IN F(X)
8,8 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS
0, 0 fits to exp(-x**2) times a Pade, 0 to a Pade
1E-25,1E-15, err(i)= max(ep1,ep2*dat(i)) this fit is to a part in 100000, NUMBER OF MINIMIZATION STEPS 1015 real*8 limit.
.9, INITIAL FR DECREASE DESIRED as we lowe precision small jumps are harder than big ones
1 VARY O FIX FITTED CONS ERROR IN CONS
0 0.5000000000000000 +- 0.000E+00
1 -0.7780697362924373 +- 0.916E-04
1 0.5990638419338665 +- 0.129E-03
1 -0.2771270371606460 +- 0.891E-04
1 0.8027282341142808E-01 +- 0.354E-04
1 -0.1379770398066123E-01 +- 0.805E-05
1 0.1113061469828807E-02 +- 0.857E-06
1 0
1 -2.684518639716602 +- 0.183E-03
1 3.227282589432518 +- 0.465E-03
1 -2.263586666841055 +- 0.520E-03
1 1.006883612355409 +- 0.330E-03
1 -0.2865296587487073 +- 0.127E-03
1 0.4891172945745505E-01 +- 0.285E-04
1 -0.3945696943660320E-02 +- 0.304E-05
1 0
CHIS=0.3550E+08 FOR 5001 DATA POINTS
CHI USED IN CALCULATING ERRORS IS 0.3550E+08
Figure 1 fi-fA(xi)
In this case fi is exp(z2) × AiGauss(z), and fA(x) = Brack(x) with the 8 numerator and denominator constants shown above. The quantitiy minimized emphasizes the left side of figure 1. The numerical thickness in the data to the left indicates how far the fit has to go to get to the noise level.
I:\PUBLIC~1\nlfit\EXPTIM~1>mpnlfit aigau2.dir the complete output is below, but this is
NADIR=aigau2.dir a good time for coffee.
NAOUT=AIGAU.DIR
NADAT=..\TestData\BRACK.OUT
ND= 1
VORF CONS
0 0.500000
1 -0.778070
1 0.599064
1 -0.277127
1 0.802728E-01
1 -0.137977E-01
1 0.111306E-02
1 0.000000E+00
1 -2.68452
1 3.22728
1 -2.26359
1 1.00688
1 -0.286530
1 0.489117E-01
1 -0.394570E-02
1 0.000000E+00
FIT IS BEING MADE TO DATA IN FILE
..\TestData\BRACK.OUT
FR,CHI,CHB,CHL 0.900000000 35499382.3 0.100000000E+67 0.500000000E+67
FR,CHI,CHB,CHL 0.729000000 25968696.2 25879049.7 35499382.3
FR,CHI,CHB,CHL 0.387420489 10184521.6 10060805.0 25968696.2
FR,CHI,CHB,CHL 0.581497370E-01 556637.656 592227.255 10184521.6
from bnraph going low cannot get to .01 * chi with any lambda exp(-175)
FR,CHI,CHB,CHL 0.100000000E-01 16867.8334 6673.35984 556637.656
from bnraph going low
FR,CHI,CHB,CHL 0.100000000E-01 9421.78133 9466.63002 16867.8334
from bnraph going low
FR,CHI,CHB,CHL 0.100000000E-01 9411.07375 9407.58683 9421.78133
from bnraph going low
FR,CHI,CHB,CHL 0.109000000 9423.80202 9409.60583 9411.07375
IREV= 1 IEND=99993
from bnraph going low
FR,CHI,CHB,CHL 0.910900000 9423.80202 9409.60583 9411.07375
IREV= 2 IEND=99993
from bnraph going low
FR,CHI,CHB,CHL 0.991090000 9423.80202 9409.60583 9411.07375
IREV= 20 IEND=99993 I modified the code, so that if it keeps getting the same value it jumps to the last IREV within 20-100, the chi-square at this level is dominated by truncation errors.
AT END FR, CHI 0.999109000 9411.07375
TEMP2= 0.000000000E+00
NEW SECTION FR= 0.9991090000000000 nlfit is set to automatically restart
FR,CHI,CHB,CHL 0.999999990 9411.07375 0.100000000E+67 0.500000000E+67
FR,CHI,CHB,CHL 0.999999970 9411.07375 9411.07347 9411.07375
IREV= 1 IEND=99992
FR,CHI,CHB,CHL 0.999999997 9411.07375 9411.07372 9411.07375
IREV= 2 IEND=99992
FR,CHI,CHB,CHL 0.999999999 9411.07375 9411.07374 9411.07375
AT END FR, CHI 0.999999999 9411.07375
TEMP2= 0.999999901E-09
NEW SECTION FR= 0.9999999990000000
CHI AT END OF FIT 9411.0737524199120000
0 0.5000000000000000 +- 0.0000000000000000
1 -0.8591435477130216 +- 1.6272308649864600D-004
1 0.7341466366171062 +- 2.5735696428003250D-004
1 -0.3871217792735996 +- 2.0090444822930290D-004
1 0.1342351250044111 +- 9.4349896125745370D-005
1 -0.0304744298576865 +- 2.7806336792956480D-005
1 0.0042054633161981 +- 4.8818201328104450D-006
1 -2.7497424824452200D-004 +- 4.0523682417858010D-007
1 -2.8466662625211940 +- 3.2544617298786200D-004
1 3.6804121795519390 +- 0.0008819329591858
1 -2.8327305035912660 +- 0.0010715148057411
1 1.4258647627189870 +- 0.0007606586755338
1 -0.4833027070743471 +- 3.4312244331430270D-004
1 0.1085165756053468 +- 9.9288237674976140D-005
1 -0.0149079729638434 +- 1.7305634415493940D-005
1 0.0009747584530339 +- 1.4365264327455150D-006
FINAL CHIT= 9411.0737524199120000
CHI/DEG OF FREE 1.8874997497833760
RESULTS IN FILE AIGAU.DIR
The entire run takes about hour on a 1 Giga-Hertz laptop.
Aigau.dir
aigau2.dir, OUTPUT FILE NAME
..\TestData\BRACK.OUT, DATA SET NAME
1, DIMENSION OF X IN F(X)
8,8 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS
0, 0 fits to exp(-x**2) times a Pade, 0 to a Pade
1E-25,1E-15, err(i)= max(ep1,ep2*dat(i))
100000, NUMBER OF MINIMIZATION STEPS
.9, INITIAL FR DECREASE DESIRED
1 VARY O FIX FITTED CONS ERROR IN CONS
0 0.5000000000000000 +- 0.000E+00
1 -0.8591435477130216 +- 0.163E-03
1 0.7341466366171062 +- 0.257E-03
1 -0.3871217792735996 +- 0.201E-03
1 0.1342351250044111 +- 0.943E-04
1 -0.3047442985768653E-01 +- 0.278E-04
1 0.4205463316198079E-02 +- 0.488E-05
1 -0.2749742482445220E-03 +- 0.405E-06
1 -2.846666262521194 +- 0.325E-03
1 3.680412179551939 +- 0.882E-03
1 -2.832730503591266 +- 0.107E-02
1 1.425864762718987 +- 0.761E-03
1 -0.4833027070743471 +- 0.343E-03
1 0.1085165756053468 +- 0.993E-04
1 -0.1490797296384342E-01 +- 0.173E-04
1 0.9747584530339248E-03 +- 0.144E-05
CHIS= 9411. FOR 5001 DATA POINTS
CHI USED IN CALCULATING ERRORS IS 9411.
Figure 2 (f-fA)/(1e-15*f)
The region -6<z<4 is approaching its accuracy limit. This is probably caused by truncation errors in evaluating the integrals. To the extent that these effects are “random” the underlying curve is more accurate by a factor of √(Nd/Ncons) = 17. Both regions have systematic errors on the order of 2 implying that this fit does not quite reach the goal.
The pass with 9,9 unlike the earlier ones, does not drop chi-square by orders of magnitude. The numerical noise ~ 100 also stops a systematic approach to the minimum, once its vicinity has been reached.
FIT IS BEING MADE TO DATA IN FILE
..\TestData\BRACK.OUT
FR,CHI,CHB,CHL 0.900000000 9612.87997 0.100000000E+67 0.500000000E+67
FR,CHI,CHB,CHL 0.729000000 7079.62515 7007.78950 9612.87997
from bnraph going low
FR,CHI,CHB,CHL 0.387420489 6713.67300 6599.35985 7079.62515
from bnraph going low
FR,CHI,CHB,CHL 0.152334763 6711.56567 6691.10203 6713.67300
from bnraph going low
FR,CHI,CHB,CHL 0.120885292 6689.95567 6705.41629 6711.56567
from bnraph going low
FR,CHI,CHB,CHL 0.100000000E-01 6677.24298 6685.54279 6689.95567
from bnraph going low
FR,CHI,CHB,CHL 0.100000000E-01 6639.01696 6674.19657 6677.24298
from bnraph going low
FR,CHI,CHB,CHL 0.100000000E-01 6694.62463 6635.20438 6639.01696
IREV= 1 IEND=99993
from bnraph going low
FR,CHI,CHB,CHL 0.901000000 6694.62463 6635.20438 6639.01696
IREV= 2 IEND=99993
from bnraph going low
FR,CHI,CHB,CHL 0.990100000 6694.62463 6635.20438 6639.01696
IREV= 20 IEND=99993
FR,CHI,CHB,CHL 0.999010000 6639.01696 6639.01696 6639.01696
AT END FR, CHI 0.999010000 6639.01696
TEMP2= 0.000000000E+00
NEW SECTION FR= 0.9990100000000001
FR,CHI,CHB,CHL 0.999999990 6639.01696 0.100000000E+67 0.500000000E+67
FR,CHI,CHB,CHL 0.999999970 6639.01696 6639.01676 6639.01696 At this IREV= 1 IEND=99992 level there are no small predictions rather than
FR,CHI,CHB,CHL 0.999999997 6639.01696 6639.01694 6639.01696 the
IREV= 2 IEND=99992 traditional fr 1, would do better with fr .9 so
FR,CHI,CHB,CHL 0.999999999 6639.01696 6639.01695 6639.01696 that a larger
AT END FR, CHI 0.999999999 6639.01696 distance in parameter space is
TEMP2= 0.999999953E-09 attempted.
NEW SECTION FR= 0.9999999990000000
Figure 3 (fi-fA(xi))/(fi × 1e-15)
The systematic errors appear to be on the order of 5 × 10-16. The “random” errors result in function errors that are reduced from the fluctuationa above by a factor approximately equal to √5000/18 = 16 so that these are also less than 1 × 10-15.
Following is last constant file with all updates
aigau.dir, OUTPUT FILE NAME
..\TestData\brack.OUT, DATA SET NAME
1, DIMENSION OF X IN F(X)
9,9 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS
0, 0 fits to exp(-x**2) times a Pade, 0 to a Pade
1E-30,1E-15, err(i)= max(ep1,ep2*dat(i))
100000, NUMBER OF MINIMIZATION STEPS
.9, INITIAL FR DECREASE DESIRED
1 VARY O FIX FITTED CONS ERROR IN CONS
0 0.5000000000000000 +- 0.000E+00
1 -1.065781843119121 +- 0.254E-01
1 1.098849939447941 +- 0.441E-01
1 -0.7059835832256011 +- 0.380E-01
1 0.3064328550407249 +- 0.202E-01
1 -0.9175866513292681E-01 +- 0.709E-02
1 0.1853701251033871E-01 +- 0.163E-02
1 -0.2323566331513835E-02 +- 0.228E-03
1 0.1400748062111785E-03 +- 0.152E-04
1 -3.259942853333778 +- 0.507E-01
1 4.876151480518089 +- 0.145
1 -4.406424837428250 +- 0.189
1 2.661133284964304 +- 0.147
1 -1.118888036108821 +- 0.746E-01
1 0.3293940293195529 +- 0.256E-01
1 -0.6596029956005756E-01 +- 0.581E-02
1 0.8236827162542478E-02 +- 0.810E-03
1 -0.4965522790119376E-03 +- 0.540E-04
CHIS= 6368. FOR 5001 DATA POINTS
CHI USED IN CALCULATING ERRORS IS 6368.
These constants above give
Note that the signs of the constants are such that for negative z, every term adds. There is no cancellation and the truncation error is minimized. There is a slight gain in fitting the results of the integral directly. The constants below are slightly better.
The code was designed to fit AiGauss directly. Starting with the constants above, a number of formats needed changing in ..\TestData\genint\TGLEG67.FOR to get the full range of the data The data above with much smaller derivatives seems unaffected by these changes. The final fit begins with the constants from fitting brack above and a chi-square of 4647. This is lower than that for brack. The values used now are direct integration values while brack’s fit multiplied these values by exp(z2). It appears that this introduces an extra 2000 into the chi-square.
aigau2.dir, OUTPUT FILE NAME
..\TestData\aigauss.OUT, DATA SET NAME
1, DIMENSION OF X IN F(X)
9,9 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS
1, 1 fits to exp(-x**2) times a Pade, 0 to a Pade
1E-30,1E-15, err(i)= max(ep1,ep2*dat(i))
100000, NUMBER OF MINIMIZATION STEPS
.9, INITIAL FR DECREASE DESIRED
1 VARY O FIX FITTED CONS ERROR IN CONS
0 0.5000000000000000 +- 0.000E+00
1 -1.502885941371286 +- 0.652E-01
1 1.851590712127221 +- 0.113
1 -1.350498684163228 +- 0.964E-01
1 0.6470300251477153 +- 0.510E-01
1 -0.2101536907945629 +- 0.178E-01
1 0.4549469611716167E-01 +- 0.405E-02
1 -0.6057434494797435E-02 +- 0.562E-03
1 0.3854291140973252E-03 +- 0.371E-04
1 -4.134151049838160 +- 0.130
1 7.368071342513580 +- 0.372
1 -7.633077301200458 +- 0.482
1 5.148920726637210 +- 0.372
1 -2.373642689043555 +- 0.188
1 0.7557095175197603 +- 0.640E-01
1 -0.1619578118794217 +- 0.144E-01
1 0.2147303973187547E-01 +- 0.199E-02
1 -0.1366310760644399E-02 +- 0.131E-03
CHIS= 4580. FOR 5001 DATA POINTS
CHI USED IN CALCULATING ERRORS IS 4984.
The minimum without the truncation correction was 4519. Simply re-reading the constants goes to 4619 4566 (This is the random/truncation zone). The constants above are the best for representing AiGauss.
Figure 4 Final fit to AiGauss. -- Black is data, blue is fit, red is |data-fit|
Figure 5 (data-fit)/(data ×10-15)
There is a truncation like error possibly in the Gauss Quadrature for -5.5<z<-3.5. There is a small systematic effect at about -2.0 of a size < 1. In general the “random” error in the fit is less than that in the data by a factor of √(5000/18) = 16.7 that the error is < 10-15.
The values below -6 are quickly generated to this accuracy using the asymptotic expansion ..\..\TestData\genint\Expanding the integrand.htm. doc. -- ..\..\TestData\genint\genint.zip Those above 0 are 1-the value below zero. Thus the code to find AiGauss everywhere is ready to write.
The code for making this fit is in ..\..\zips\FittingBrack.zip
aigau2.dir -- sample direction file for fitting brack.out
AIGAUSS.OUT -- data of the form exp(-z2)*Brack(z) generated by direct integration of Gaussian
BRACK.OUT -- data of the form Brack(z) Integral has been multiplied by exp(z2) with slight loss of accuracy
transmp.exe -- Bailey’s code for converting to multi-precision -- ..\..\..\..\MultiplePrecision\bailey\Welcome.htm
mpfun.for --- Bailey’s set of multi-precision routines
nlfit.wpj Watcom ide’s for making nlfit
nlfit.tgt
cnlfit.bat -- bat file for using wfl386 to make nlfit
mpnlfit.wpj Watcom ide’s for making mpnlfit
mpnlfit.tgt
cmpnlfit.bat bat file for using wfl386 to make mpnlfit
MAOPEN.FOR -- used by Watfor
MSOPEN.FOR used by watcom
WSYSTEM.FOR -- used by Watfor
openwsys.for used by watcom
nlfpade.for - nlfit slightly modified to have Nnum, Nden, eps1, and eps2. Passes flag for fitting Brack or exp*Brack to poly
MPnlfpade.for Generated from nlfpade by running transmp.exe
ROBMIN.FOR --- minimization code. Run transmp.exe and enter robmin.for to produce
MProbmin.for --- MProbmin
polyexpa.for --- The poly routine for making exp(x2)×Brack(x) or simply Brack(x)