[WP-34S] Chi-Square Distribution

[Sanjeev Visvanatha]

Hello,

I was comparing output between HP48sx MATHLIB [Ref. 1] and WP-34S (v 3.2 3375) for some of the probability distributions. For this purpose, I was replicating the examples from the HP-21S user's manual (pg. 42-52)

I noticed that the Inverse of the Chi-Square distribution for the WP-34S does not seem (to me) to produce the same result as what is listed in the 21S manual pg. 50.

Degrees of freedom = 19
Upper tail probability = 0.001

Keystrokes on WP-34S:

19 [STO] J
0.001 [h] [PROB] [ϰ²INV]
This results in 5.4068 (in FIX4)

The 21S example on pg. 50 of its user manual gives a result of 43.8202 (function name is ϰ²p). This is also consistent with HP-48sx MATHLIB function IUTPC(19, 0.001).

For reference, the IOP for the v3.2 manual of the WP-34S indicates that ϰ²INV should be equal to ϰ²p of the 21S. Is my understanding of how to use this distribution on the WP-34S incorrect?

Thank you for any insights you may have!

[1] HP 48SX Engineering Mathematics Library, John F. Holland

[walter b]

(01-11-2014 01:37 PM)Sanjeev Visvanatha Wrote:  Hello,

I was comparing output between HP48sx MATHLIB [Ref. 1] and WP-34S (v 3.2 3375) for some of the probability distributions. For this purpose, I was replicating the examples from the HP-21S user's manual (pg. 42-52)

I noticed that the Inverse of the Chi-Square distribution for the WP-34S does not seem (to me) to produce the same result as what is listed in the 21S manual pg. 50.

Degrees of freedom = 19
Upper tail probability = 0.001

Keystrokes on WP-34S:

19 [STO] J
0.001 [h] [PROB] [ϰ²INV]
This results in 5.4068 (in FIX4)

The 21S example on pg. 50 of its user manual gives a result of 43.8202 (function name is ϰ²p). This is also consistent with HP-48sx MATHLIB function IUTPC(19, 0.001).

For reference, the IOP for the v3.2 manual of the WP-34S indicates that ϰ²INV should be equal to ϰ²p of the 21S. Is my understanding of how to use this distribution on the WP-34S incorrect?

Thank you for any insights you may have!

[1] HP 48SX Engineering Mathematics Library, John F. Holland

5.4068 is correct for 0.1%, while 43.8202 is correct for 99.9% Integrating form zero to p, the value returned must increase for growing p (see pp. 54f of the manual quoted). The HP-21S calculates with the error probability instead.

d:-)

[Dieter]

(01-11-2014 03:02 PM)walter b Wrote:  

(01-11-2014 01:37 PM)Sanjeev Visvanatha Wrote:  Degrees of freedom = 19
Upper tail probability = 0.001

5.4068 is correct for 0.1%, while 43.8202 is correct for 99.9% Integrating form zero to p, the value returned must increase for growing p (see pp. 54f of the manual quoted). The HP-21S calculates with the error probability instead.

I tried this example on my 34s (v. 3.2 3405) and was puzzled about the execution time: about 25 seconds. So I wrote a short quick-and-dirty program based on my original contributions for the chi-square quantile (providing a decent initial guess, followed by a slightly modified Newton iteration), and even this simple user-code program got the result in not much more than four (!) seconds (SP), resp. 5-6 seconds in DP mode. The initial guess is approx. 5,6 (resp. 42,8 for p=0,999) and requires four or five iterations for SP accuracy and one or two more for DP.

Another example:
n = 20, p = 0,05 requires 17 s with the internal function and about 4 resp. 5 s in user code.

Which leads to the question: what's going on there?

Dieter

[Paul Dale]

The distribution code hasn't changed since we did all that work on it. Looking at the code, it boils down to doing the initial estimate and then running the Newton solver. The internal code will be running in double precision, although the convergence criteria is precision mode dependent.


Is it possible your 34S is in one of the slow clock states?


- Pauli

Code:

                XLBL"QF_CHI2"
                        xIN MONADIC
                        GSB chi2_param
                        GSB qf_check_probability
                        GSB qf_chi2_est
                        _INT DIST_CHI2
                        XEQ qf_newton
                        xOUT xOUT_NORMAL

#define R_P     .00
#define R_T     .01
qf_chi2_est::           LocR 2
                        STO R_P
                        _INT 15
                        x<? J
                                JMP qf_chi2_est_ll
                        _INT 1
                        JMP qf_chi2_est_ld
qf_chi2_est_ll::        _INT 85
                        SDR 2
qf_chi2_est_ld::        +/-
                        RCL[times] J
                        e[^x]
                        x[<=]? R_P
                                JMP qf_chi2_large
                        Num 1/2
                        RCL[times] J
                        ENTER[^]        // n/2 n/2 ? ?
                        LN[GAMMA]
                        RCL/ Y
                        e[^x]
                        RCL Y           // n/2 exp n/2 ?
                        RCL[times] R_P  // pn/2 exp n/2 ?
                        RCL Z
                        [^x][sqrt]y     // xroot exp n/2
                        [times]
                        RCL+ X
                        JMP qf_chi2_est_exit

qf_chi2_large::         RCL R_P
                        GSB qf_q_est
                        _INT 97
                        SDR 2
                        [times]         // .97q
                        Num Chi2        // 0.2214
                        RCL/ J          // .2214/n .97q
                        ENTER[^]        // .2214/N .2214/N .97q
                        [sqrt]
                        RCL[times] Z    // sqrt(.2214/N)*.97q .2214/N .97q
                        RCL- Y          // sqrt(.2214/N)*.97q-.2214/N .2214/N .97q
                        INC X
                        x[^3]
                        RCL[times] J    // ch3 .2214/N .97q
                        _INT 6
                        RCL[times] J
                        _INT 16
                        +
                        SWAP
                        x[<=]? Y
                                JMP qf_chi2_est_exit
                        RCL R_P
                        GSB log1m
                        _INT 150
                        RCL[times] J
                        /
                        INC X
                        [times]

qf_chi2_est_exit::      RCL R_P
                        RTN

[Dieter]

(01-11-2014 11:11 PM)Paul Dale Wrote:  The distribution code hasn't changed since we did all that work on it.
...
Is it possible your 34S is in one of the slow clock states?

Just to be sure I set it to FAST mode. But I do not think this is the problem - the point is that the iteration in user code runs so much faster than the XROM routine.

Maybe there's a problem with the initial guess. Could you provide the first guesses for df=20 and p = 1E-10, 0,05, 0,95 and 0,99999? Or is there a way to single-step through the routine so that I can do this myself? The initial guess should be within approx. ±10% of the final result.

Back then we also discussed a modification of the Newton iteration in that the correction term applied in each iteration should be limited to, say, ±1 or ±5%. This prevents divergence and oscillation that may otherwise occur in some cases.

Dieter

[Sanjeev Visvanatha]

(01-11-2014 03:02 PM)walter b Wrote:  5.4068 is correct for 0.1%, while 43.8202 is correct for 99.9% Integrating form zero to p, the value returned must increase for growing p (see pp. 54f of the manual quoted). The HP-21S calculates with the error probability instead.

Thank you for clarifying my misunderstanding. I wonder if you consider the statement on pg. 126 of the WP-34S manual correct then? :

"... ϰ²INV equals ϰ²p in the HP-21S ..."

[Sanjeev Visvanatha]

(01-11-2014 06:18 PM)Dieter Wrote:  Another example:
n = 20, p = 0,05 requires 17 s with the internal function and about 4 resp. 5 s in user code.

In SLOW mode, this takes about 18s for me, and about 12s in FAST mode.

[walter b]

(01-12-2014 03:13 PM)Sanjeev Visvanatha Wrote:  

(01-11-2014 03:02 PM)walter b Wrote:  5.4068 is correct for 0.1%, while 43.8202 is correct for 99.9% Integrating form zero to p, the value returned must increase for growing p (see pp. 54f of the manual quoted). The HP-21S calculates with the error probability instead.

Thank you for clarifying my misunderstanding. I wonder if you consider the statement on pg. 126 of the WP-34S manual correct then? :

"... ϰ²INV equals ϰ²p in the HP-21S ..."

I replaced it by "... χ²INV corresponds to χ²p in the HP-21S ..." and put an extra footnote on p. 55 stating:

The HP-21S returns Q instead of P. Thus, also the inverse works differently there.

Also modified the respective IOP entries for F, t, and Φ this afternoon. Thank you for your head up. Anyway, however, those modifications will show up no earlier than with the second (revised) edition of the printed manual. Please be patient.

d:-)

[Sanjeev Visvanatha]

(01-12-2014 05:57 PM)walter b Wrote:  I replaced it by "... χ²INV corresponds to χ²p in the HP-21S ..." and put an extra footnote on p. 55 stating:

The HP-21S returns Q instead of P. Thus, also the inverse works differently there.

Also modified the respective IOP entries for F, t, and Φ this afternoon. Thank you for your head up. Anyway, however, those modifications will show up no earlier than with the second (revised) edition of the printed manual. Please be patient.

d:-)

Perhaps you would consider an errata sheet for those of us that have recently purchased the current printed manual

[Paul Dale]

(01-12-2014 02:37 PM)Dieter Wrote:  Maybe there's a problem with the initial guess. Could you provide the first guesses for df=20 and p = 1E-10, 0,05, 0,95 and 0,99999? Or is there a way to single-step through the routine so that I can do this myself? The initial guess should be within approx. ±10% of the final result.

It should be possible. Turning trace mode on lets you single step through xrom. I'm not sure how to do this in any emulator except the console one (which you really don't want to use These programs could be typed in in user mode, which often is the best way to debug them.

Yeah, the initial guess is bad for some input values:

Code:

 1E-10    0.9454441852127046    0.9057457376233529524406618674391
 .05    10.85081139418259    30.903880572067197129947058754473
 .95    31.41043284423093    30.903880572067197129947058754473
 .99999    59.04455038680165    57.7753881925718842706491638346524
 .9    28.41198058430563    28.0176453217379586052786281495866
 .4    16.26585648501278    19.343125586513683013556720398377658

Seems good for larger values and really small ones. Pretty bad at .05 however.

Time for some digging into code.


- Pauli

[Dieter]

(01-12-2014 09:38 PM)Paul Dale Wrote:  Turning trace mode on lets you single step through xrom.
...
Yeah, the initial guess is bad for some input values:
...
Seems good for larger values and really small ones. Pretty bad at .05 however.

Seems we found the first error. Take a look at the initial estimates for 0,05 and 0,95. They are the same. In the center, the guess is based on the Normal estimate, so that I assume there might be a problem there. You may want to check if the Normal estimate returns the same absolute value with opposite signs for 0,05 and 0,95.

Back then we discussed different possible ways to provide a good estimate for the Chi^2 quantile. So the version I tried yesterday may be different from the one that actually got implemented. But anyway, here are the results from a slightly simplified method and df=20:

Code:

  p        quantile     estimate
1E-10      0,94544      0,905...
 0,05      10,8508      10,93...
 0,95      31,4104      31,21...
 0,99999   59,0446      58,46...

Not bad for a first estimate, I think. SP accuracy usually does not need more than four or maybe five iterations. If (!) the guess is implemented correctly.

Dieter

[Paul Dale]

Think I found the problem. The normal quantile estimate wasn't dealing with small input probabilities properly. I don't know if this fix breaks anything else


- Pauli

[walter b]

(01-12-2014 07:23 PM)Sanjeev Visvanatha Wrote:  Perhaps you would consider an errata sheet for those of us that have recently purchased the current printed manual

OK, for you and all the others: Here are the errors found and reported in almost one year since the manual was printed. One sheet of paper will do.

d:-)

[Sanjeev Visvanatha]

(01-13-2014 07:55 PM)walter b Wrote:  OK, for you and all the others: Here are the errors found and reported in almost one year since the manual was printed. One sheet of paper will do.

Thank you!

[Dieter]

(01-13-2014 01:08 AM)Paul Dale Wrote:  Think I found the problem. The normal quantile estimate wasn't dealing with small input probabilities properly. I don't know if this fix breaks anything else

The Normal estimate is crucial not only for the Normal quantile. A problem with small probabilities would cause inaccurate results with the Normal quantile as well - which I have not noticed yet.

The Normal estimate distinguishes between two cases: 0,5 >= p > 0,23 and p < 0,23. So an error with small probabilities would affect every p < 0,23. That's why I wonder what exactly the problem was. Could you post the Normal estimate code before and after your last correction?

Dieter

[walter b]

(01-14-2014 02:33 PM)Dieter Wrote:  Could you post the Normal estimate code before and after your last correction?

You find it at sourceforge

[Dieter]

(01-14-2014 02:46 PM)walter b Wrote:  You find it at sourceforge

Yes, you will find virtually anything "on the internet". #-)

Sorry, but I am simply not able to reliably find a certain piece of information on the sourceforge website. How could I, since I do not even know which piece of code is included in which file and in which directory? You seem to know better, so may I ask you for two direct links to the old and new code?

Dieter

[Paul Dale]

The problem I fixed added a +/- at the end if p<.5
The number coming out was correct, just with the wrong sign.

Code below. The change is to add FS? F_SMALL +/- right down at the bottom.

- Pauli

Code:

qf_q_est::    LocR 003

qf_q_int_est::    ENTER[^]
        +/-
        INC X
        MIN
        STO R_MINP        // q = min(p, 1-p)
        x[!=]? L
            SF F_SMALL    // Set Flag if p < 0.5 (i.e. 1-p != q)
        Num 1/2            // pre-store 0.5 in a register for later use
        STO R_HALF
        _INT 023
        SDR 002
        RCL R_MINP
        x>? Y
            JMP qf_q_mid
        ENTER[^]        // tail estimate
        LN
        STO+ X
        +/-
        STO Z
        DEC X
        Num [pi]
        [times]
        STO+ X
        [sqrt]
        [times]
        LN
        STO+ X
        +/-
        [sqrt]
        _INT 254
        SDR 003
        RCL/ Z
        +
        STO R_ABSZ
        JMP qf_q_signfix

qf_q_mid::    Num 1/2
        RCL- Y
        Num [pi]
        STO+ X
        [sqrt]
        [times]
        ENTER[^]
        x[^3]
        _INT 006
        /
        +
        STO R_ABSZ
        _INT 004
        SDR 002
        x>? Y
            CF F_ITER        // z < 0.04 requires just one correction step
        SDR 007
        x>? Y
            SF F_EXACT        // z < 4E-9 needs no correction at all, flag for exit
qf_q_signfix::    FS? F_SMALL
            +/-
        RTN

[Dieter]

(01-14-2014 10:07 PM)Paul Dale Wrote:  The problem I fixed added a +/- at the end if p<.5
The number coming out was correct, just with the wrong sign.

Code below. The change is to add FS? F_SMALL +/- right down at the bottom.

Thank you. But as far as I can tell the new code contains a substantial error: For p in the center (i.e. between 0,23 and 0,77) the estimate uses the code following label qf_q_mid. Then the estimate is compared with two threshold values (0,04 and 4E-9), and finally label qf_q_signfix is reached where the sign gets adjusted. But at this point X does not hold the absolute value of the estimate (R_ABSZ), but 4E-9 (!). The estimate is in Y. As far as I can tell there is a X<>Y, Roll down or RCL R_ABSZ missing.

This error must be relatively new since otherwise the Normal quantile would not be able to return exact results - the quantile algorithm uses at most two iterations, and this only works if the initial estimate is sufficiently accurate.

Or do I miss something here?

Dieter

[Paul Dale]

Something I'll have to look into. Not going to get a good lump of time in the near future unfortunately, so this might take a while. Unless someone else wants to investigate.


- Pauli