Yet another Normal quantile function for the 34s



#2

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.


#3

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


#4

Quote:
We've also gained # 1/2, [cmplx]# 1/2, IDIV and [cmplx]IDIV.

Nice, thanks! :-)

Franz

#5

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.


#6

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.


#7

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


#8

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.


#9

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&#9113; 3225 - nevertheless I get ...8742 for p = 0.9 and ...8344 for p = 0.5001 using g-shifted 4. Strange :-?

#10

Quote:
Strange :-?

Just cosmic noise! ;-)

#11

Quote:

Just cosmic noise! ;-)


Looks more like quantum noise to me at 10^(-34) ;-)

#12

Quote:
Looks more like quantum noise to me at 10^(-34) ;-)

WOW, so the WP34s is the first quantum calculator. :-)))
#13

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

#14

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


#15

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP50g: Writing a function that returns a function Chris de Castro 2 294 12-10-2013, 06:49 PM
Last Post: Han
  WP-34S function execution speed ? Gene Wright 4 248 09-04-2013, 05:40 PM
Last Post: Paul Dale
  SINC function (WP 34S) Gerson W. Barbosa 22 822 08-24-2013, 03:05 AM
Last Post: Didier Lachieze
  [WP 34s] Summation Function (Sigma) Paul C 11 456 01-29-2013, 07:42 AM
Last Post: C.Ret
  Summary: Normal CDF and quantile function Dieter 3 189 05-13-2012, 02:20 AM
Last Post: Paul Dale
  [WP34S] Curious Bug in Inverse Normal Function Les Wright 61 1,212 05-01-2012, 02:44 AM
Last Post: Paul Dale
  SUM function on WP-34S Thomas Scott 11 446 04-18-2012, 06:10 PM
Last Post: fhub
  32E's Normal & Inverse Normal Distribution Matt Agajanian 27 930 04-14-2012, 06:07 PM
Last Post: Dieter
  HP 32sII Integration Error of Standard Normal Curve Anthony (USA) 4 243 03-14-2012, 03:18 AM
Last Post: Nick_S
  [wp 34s] wp 34s picture and scan Jeroen Van Nieuwenhove 2 156 10-27-2011, 09:02 PM
Last Post: Les Wright

Forum Jump: