Yet another Normal quantile function for the 34s - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: Yet another Normal quantile function for the 34s (/thread-228275.html) Yet another Normal quantile function for the 34s - Dieter - 08-02-2012 As announced a few hours ago, here is my current (experimental ;-)) version of a 34s user code program for the Normal distribution's quantile function (QF); i.e. enter a probability p and get the corresponding z-value for which the Normal integral from -infinity to z equals p. Yes, the 34s (of course) provides such a function (Phi-1 on g-shift 4), which has undergone many changes since its first implementation. The last version I used (3.1 3208) still has some issues with arguments very close to (but less than) 0.5 and very close to 1. It does not provide the accuracy level that might be possible. Results are fine in SP mode, but in double precision there is some room for improvement. The following program works quite straightforward. First, it evaluates an initial estimate for the quantile: ```let q = min(p, 1-p) then, for q > 0.23 u = (0.5 - q) * sqrt(2 * pi) z = u + u3/6 else u = -2 * ln(q) z = sqrt(-2 * ln(q * sqrt(2 * pi * (u - 1)))) + 0.254/u ``` Then this estimate is improved by the method suggested in A&S, p. 954: ``` t = (upperCDF(z) - q) / PDF(z) z = z + t + t2*z/2 + t3*(2z+1)/6 ``` Finally the sign of z is adjusted: for p >= 0.5 return z, else -z. Now there is a pitfall in the evaluation of t. For z less than 1 and, more important, close to zero, the difference between CDF(z) and q may lose a more or less substantial number of digits. This problem can be solved if the Normal integral from 0 to z is evaluated, i.e. CDF(z) - 0.5. The 34s does not provide such a function, but there is a simple workaround: ``` CDF(z) - 0.5 = 1/2 Incomplete_regularized_Gamma_p(1/2, z2/2) ``` The desired difference then simply is (0.5 - q) minus this. The program uses this method for z =< 0.5. Switching at z = 1 should also be fine. When evaluated with sufficient precision (I did some tests in a 38-digit environment), two successive correction steps should provide a result with a max. error of approx. 5 units in the 35th significant digit (5D35), i.e. at the limit of what can be done in double precision. For z < 4E-9 the error of the initial estimate is less than 6D36, so no further correction is required, and z < 0.04 is so close that just one correction step keeps the error near 2D35 or below. At least when evaluated in 39-digit precision. ;-) If the correction method is expanded with a fourth term, the error can even be pushed down to somewhere in the 40th digit. Which may be an option if the algorithm is adopted for an internal function with 39-digit precision. In this case the threshold values mentioned in the previous paragraph should be changed to 2E-10 resp. 0.02. Provided the CDF resp. IGamma functions have all their 34 digits correct, the accuracy on a 34s in DP mode should be 33 digits. I did not do any extensive tests, so beware and see yourself what you get. According to my results, z > 1 usually is either exact or sometimes 1 ULP off. For z < 1 the last digit may be off by several units (I did not notice more than 3). And finally: here's the code. ```0001 *LBL'QF' 0002 CF 00 0003 SF 01 0004 ENTER[^] 0005 +/- 0006 INC X 0007 MIN 0008 STO 01 // q = min(p, 1-p) 0009 x[!=]? L 0010 SF 00 // Set Flag 0 if p < 0.5 (i.e. 1-p != q) 0011 . 0012 5 // pre-store 0.5 in R00 for later use 0013 STO 00 0014 # 023 0015 SDR 002 0016 RCL 01 0017 x>? Y 0018 GTO 01 0019 ENTER[^] // tail estimate 0020 LN 0021 STO+ X 0022 +/- 0023 STO Z 0024 DEC X 0025 # [pi] 0026 [times] 0027 STO+ X 0028 [sqrt] 0029 x 0030 LN 0031 STO+ X 0032 +/- 0033 [sqrt] 0034 # 254 0035 SDR 003 0036 RCL/ Z 0037 + 0038 STO 02 0039 GTO 02 0040 *LBL 01 // estimate for central region 0041 RCL 00 0042 RCL- Y 0043 # [pi] 0044 STO+ X 0045 [sqrt] 0046 [times] 0047 ENTER[^] 0048 x[^3] 0049 # 006 0050 / 0051 + 0052 STO 02 0053 # 004 0054 SDR 002 0055 x>? Y 0056 CF 01 // z < 0.04 requires just one correction step 0057 SDR 007 0058 x>? Y 0059 GTO 03 // z < 4E-9 needs no correction at all => exit 0060 *LBL 02 0061 RCL 02 0062 x>? 00 // "x>1?" should be fine either 0063 SKIP 009 0064 x[^2] 0065 RCL[times] 00 0066 RCL 00 0067 I[GAMMA][sub-p] 0068 RCL[times] 00 0069 RCL 00 0070 RCL- 01 // For q close to 0.5: 0071 RCL- Y // upperCDF(z) - q = (0.5-q) - IGamma_p(0.5, z2/2) 0072 SKIP 002 0073 [PHI][sub-u](x) // Otherwise evaluate Phi_u(z) - q directly 0074 RCL- 01 0075 RCL 02 0076 [phi](x) 0077 / 0078 FILL 0079 x[^3] // Apply a third-order correction 0080 RCL 02 // Evaluate all three terms individually 0081 x[^2] // and add them in order of their magnitude, 0082 STO+ X // smallest terms first, to preserve 0083 INC X // as much accuracy as possible 0084 [times] 0085 # 006 0086 / 0087 x[<->] Y 0088 x[^2] 0089 RCL[times] 02 0090 RCL[times] 00 0091 + 0092 + 0093 STO+ 02 0094 FS?C 01 0095 GTO 02 0096 *LBL 03 0097 RCL 02 // exit program 0098 FS?C 00 // set sign if required 0099 +/- 0100 END Registers: R00: 0.5 R01: q = min(p, 1-p) R02: |z| Flags: F00: Set for p < 0.5, else clear F01: Set if a second iteration is required, else clear ``` Still all I got is the 34s emulator. So I would appreciate some feedback from those who use a "real" hardware 34s. How fast does this program run? Switch to double precision (Mode menu, DBLON) and see how it performs. Do you find any errors that are larger than a few units in the last (34th) digit? Please note that this code has not been thoroughly tested - simply try it and see what you get. :-) I am sure there still is room for the one or other improvement. All suggestions are welcome. Dieter Edited: 2 Aug 2012, 6:11 p.m. Re: Yet another Normal quantile function for the 34s - Paul Dale - 08-03-2012 I had to rework things a little and might have introduced some problems but a new emulator is built that uses this code instead of the C helper routines. We've also gained # 1/2, [cmplx]# 1/2, IDIV and [cmplx]IDIV. - Pauli Re: Yet another Normal quantile function for the 34s - fhub - 08-03-2012 Quote: We've also gained # 1/2, [cmplx]# 1/2, IDIV and [cmplx]IDIV. Nice, thanks! :-) Franz Re: Yet another Normal quantile function for the 34s - Dieter - 08-03-2012 Pauli, Quote:I had to rework things a little and might have introduced some problems I tried the new emulator (3.1 3215) and indeed I found some problems. The Normal quantile now often has just 10 - 20 correct digits in DP mode. Sometimes even the SP result is off. Example: ``` SP 16-digit limit | | p = 0.1 => z = -1.281 5515 6554 4600 4460 5945 1107 2681 73 p = 0.9 => z = 1.281 5515 6554 4592 7643 6081 1042 0517 37 exact |z| = 1.281 5515 6554 4600 4669 6510 3329 4487 43 ``` The code I posted can not produce different results for p and 1-p (which is one of the improvements compared to the previous algorithm). It also returns the exact values for the given examples. So something else must be going on here. Dieter Edited: 3 Aug 2012, 10:07 a.m. Re: Yet another Normal quantile function for the 34s - Walter B - 08-03-2012 AFAIK you should try build 3225. HTH ;-) Using it, I get for ...```p = 0.1 z = -1.281 5515 6554 4600 4669 6510 3329 4487 42 , p = 0.9 z = 1.281 5515 6554 4600 4669 6510 3329 4487 42 . ``` Doesn't look too bad :-) Edited: 3 Aug 2012, 10:24 a.m. Re: Yet another Normal quantile function for the 34s - Walter B - 08-03-2012 More results:```p = 0.505 -> z = 0.012 533 469 508 069 263 161 142 848 523 149 61 p = 0.495 -> z = -0.012 533 469 508 069 263 161 142 848 523 149 61 p = 0.5001 -> z = 0.000 250 662 830 088 035 098 920 650 105 407 834 4 p = 0.4999 -> z = -0.000 250 662 830 088 035 098 920 650 105 407 834 4 p = 0.500 000 01 -> z = 0.000 250 662 827 463 100 076 490 926 438 018 476 5 p = 0.499 999 99 -> z = -0.000 250 662 827 463 100 076 490 926 438 018 476 5 ``` Can you check the exact solutions, please? Edited: 3 Aug 2012, 11:06 a.m. after one or more responses were posted Re: Yet another Normal quantile function for the 34s - Dieter - 08-03-2012 The emulator download includes version 3.1 3215, so this is what I tried since Pauli said the new algorithm is included there. I now manually updated to wp34sgui.exe 3.1 3225 and now I get... ``` p = 0.1 => z = -1.281 5515 6554 4600 4669 6510 3329 4487 43 p = 0.9 => z = 1.281 5515 6554 4600 4669 6510 3329 4487 43 ``` Yes, here the final digit is a 3, just as it is supposed to be. I wonder why you get a 2 instead. What is your setting for the rounding mode? Also, for 0.4999 and 0.5001 I get ...8343, which again is the exact result. Are you sure your version is 3.1 3225? BTW, 0.4999 is not what I would call "close to 0.5". ;-) The 34s should now return results like ``` p = 0.5 - 1E-8 => z = -2.506 6282 7463 1000 7649 0926 4380 1847 66 E-8 p = 0.5 + 1E-10 => z = 2.506 6282 7463 1000 5024 4201 4634 7205 82 E-10 p = 0.5 - 1E-30 => z = -2.506 6282 7463 1000 5024 1576 5284 8110 45 E-30 ``` ...which is exact within +/- 1 ULP. The values you posted generally show the same accuracy, but they differ in the last digit. Sometimes "my" results are dead on while the ones you posted are 1 ULP off, and in other cases it's just the other way round. For instance, for p = 0.4999 4999 the correct mantissa is ...65 38... Your result rounds correctly down to ...65, while I get it rounded up to ...66. For more results visit the Casio Keisan calculator, set it to 50 digits and check "Accuracy". Now enter  normalicd(0.5 - 10^(-n))  into the "expression" window and press "Execute". You will be prompted for n, so enter  1, 1, 32 there. This will create a table. But beware: Keisan does not always provide the same accuracy level as the 34s. At least it knows this, so it will only return as many digits as it thinks you can trust. ;-) Dieter Edited: 3 Aug 2012, 11:31 a.m. Re: Yet another Normal quantile function for the 34s - Walter B - 08-03-2012 Quote: Yes, here the final digit is a 3, just as it is supposed to be. I wonder why you get a 2 instead. Also, for 0.4999 and 0.5001 I get ...8343, which again is the exact result. Are you sure your version is 3.1 3225? Yes I am: VERS tells me 34S 3.1⎙ 3225 - nevertheless I get ...8742 for p = 0.9 and ...8344 for p = 0.5001 using g-shifted 4. Strange :-? Re: Yet another Normal quantile function for the 34s - fhub - 08-03-2012 Quote: Strange :-? Just cosmic noise! ;-) Re: Yet another Normal quantile function for the 34s - Dieter - 08-03-2012 Here's a short test routine that checks the results for p close to zero, i.e. |p - 0.5| < 1E-5. It uses the first terms of the series expansion for the normal quantile which, in this range, has 38 correct digits. Yes, the code is not exactly elegant, but it works. ;-) The output shows the magnitude of 0.5-p and the error in significant (!) digits, so -1E-34 means the internal routine returns a result that is 1 unit in the 34th digit low. ```0001 **LBL A 0002 # 034 0003 SDR 003 0004 # 005 0005 + 0006 STO 00 0007 **LBL 00 0008 RAN# 0009 . 0010 9 0011 [times] 0012 . 0013 1 0014 + 0015 SDR[->]00 0016 # 1/2 0017 + 0018 STO 01 0019 # 1/2 0020 - 0021 # [pi] 0022 STO+ X 0023 [sqrt] 0024 [times] 0025 FILL 0026 # 007 0027 y[^x] 0028 RCL L 0029 x! 0030 / 0031 # 127 0032 [times] 0033 x[<->] Y 0034 # 005 0035 y[^x] 0036 RCL L 0037 x! 0038 / 0039 # 007 0040 [times] 0041 + 0042 x[<->] Y 0043 x[^3] 0044 # 006 0045 / 0046 + 0047 + 0048 RCL 01 0049 [PHI][^-1](p) 0050 x[<->] Y 0051 - 0052 RCL L 0053 EXPT 0054 +/- 0055 DEC X 0056 x[<->] Y 0057 SDL[->]Y 0058 CL[alpha] 0059 [alpha]'1E-' 0060 [alpha]IP 00 0061 VW[alpha]+ X 0062 PSE 05 0063 ISG 00 0064 GTO 00 0065 CL[alpha] 0066 CLSTK 0067 END ``` Dieter Re: Yet another Normal quantile function for the 34s - Walter B - 08-03-2012 Quote: Just cosmic noise! ;-) Looks more like quantum noise to me at 10^(-34) ;-) Re: Yet another Normal quantile function for the 34s - fhub - 08-03-2012 Quote: Looks more like quantum noise to me at 10^(-34) ;-) WOW, so the WP34s is the first quantum calculator. :-))) Re: Yet another Normal quantile function for the 34s - Dieter - 08-05-2012 In the meantime I took a closer look at the way the quantile routine behaves at the very end of the working range, i.e. at p below 1E-6143 where numbers start getting sub-normal. Here they lose one significant digit for every power of ten, until at 1E-6176 just one digit is left. For this reason the result of the quantile routine will lose accuracy as well. Down to p = 1E-6147 the usual 33-34 digits can be expected, but below accuracy slowly fades away, just as the internal precision does. At the lower limit still 7 valid digits are left. I think the 34s should return only results as far as they can be trusted. Yes, there is no claim that all 30+ digits are correct, but maybe the following code can be added just for the user's convenience: ```.... ... 0094 FS?C 01 0095 GTO 02 0096 6 // new code starts here 0097 1 0098 4 0099 7 0100 RCL 01 0101 EXPT 0102 +/- 0103 x[<=]? Y 0104 GTO 03 // down to 1E-6147 everything is fine 0105 - 0106 # 034 0107 + // else the result is expected to have this number of digits 0108 # 007 0109 MAX // but not less than 7 0110 RCL 02 0111 RSD[->]Y // round the result accordingly 0112 STO 02 // new code ends here 0113 *LBL 03 0114 RCL 02 0115 FS?C 00 0116 +/- 0117 END ``` This will round the result to the number of digits that should be fine (within 1 ULP). I also have been thinking about possible underflow. Pauli, please read this carefully and correct me if I'm wrong. :-) Underflow may appear if the initial estimate is higher than the true quantile so that Phi_u(z) will become less than 1E-6176. The code I posted here provides an estimate that is less than 1E-5 higher than the true quantile, so that Phi_u(z) is 9.9...E-6177. While most classic HPs in this case would underflow to zero, the 34s will round this to the remaining one significiant digit, i.e. to 1E-6176, which means that simply no further correction term is applied and the result has the accuracy of the initial estimate, i.e. 7 significant digits. If this behaviour (round up instead of underflow) may eventually change, or if you just want to be 120% sure, simply insert a X[!=]0? test between line 73 and 74 (and change the preceding SKIP 002 to SKIP 003). But I think currently this is not required. Otherwise I am quite satisfied with the Normal quantile function now. In most cases the result is exact within 1 ULP, here and there I found errors like +/- 2 ULP. Please let me know if you find any problems of if you have suggestions for further improvements. Dieter Re: Yet another Normal quantile function for the 34s - Paul Dale - 08-06-2012 I've always been very reluctant to round results artifically. We're already paying at least a double rounding for many operations which isn't ideal & is less so in double precision where the number of guard digits is greatly reduced. - Pauli