Re: Test Failure

[Ben Pfaff]

I've now investigated the problem more thoroughly.  There is at
least one red herring here, but it may be masking a real
problem.  I'll explain further.

Jason Stover <address@hidden> writes:

> I get a different failure, but not the one John reported:
>
> 2c2
> <      .01     2.00      .02      .01      .49      .99 
> ---
>>      .01     2.00      .02      .01      .50      .99 
> compare output for chisq distribution
> FAILED
> FAIL: tests/expressions/randist.sh
>
> Which is the same type of error as John's failed test: The second line
> is rounded up as it should be, but the value stored in
> tests/expressions/randist/chisq.out (shown on the first line) is
> not. pspp seems to have the correct value stored, as can be seen by
> telling the PRINT command to print more decimal places. The correct
> value to many digits is 4.950249168745841e-01. The reported value in
> this case should be rounded to 0.50, if the user wants only two decimal
> places displayed.

There seems to be the assumption here that the fifth column is
supposed to be the value of pdf.chisq(0.02, 2).  If I do ask for
that, then I get it, correctly rounded regardless of the number
of decimal places requested.  That is, the following input:

DATA LIST LIST/x df.
COMPUTE pdf = pdf.chisq (x, df).
PRINT OUTFILE='-'/df x pdf * df x pdf (3f9.6).
BEGIN DATA.
0.02 2 /* Should show .495025 in column 5 when x=0.02. */
9.21 2 /* Should show .005001 in column 5 when x=9.21.*/
END DATA.

produces the following output, which exactly matches your values
above[*]:

2.00000000000000000000
0.02000000000000000042
0.49502491687458405334
2.00000000000000000000
0.02000000000000000042
0.49502491687458405334
    2.00      .02      .50  2.000000  .020000  .495025

2.00000000000000000000
9.21000000000000085265
0.00500085100235273926
2.00000000000000000000
9.21000000000000085265
0.00500085100235273926
    2.00     9.21      .01  2.000000 9.210000  .005001


The catch is that this isn't what chisq is actually computing.
It is actually printing pdf.chisq(idf.chisq(0.01, 2), 2).
Here is code that illustrates that[+]:

DATA LIST LIST/P df.
COMPUTE x = idf.chisq (P, df).
COMPUTE pdf = pdf.chisq (x, df).
PRINT OUTFILE='-'/P df x pdf * df x pdf (3f20.16).
BEGIN DATA.
0.01 2 /* Should show .495025 in column 5 when x=0.02. */
0.99 2 /* Should show .005001 in column 5 when x=9.21.*/
END DATA.

and the corresponding output:

0.01000000000000000021
2.00000000000000000000
0.02010067170700288383
0.49499999999999999556
2.00000000000000000000
0.02010067170700288383
0.49499999999999999556
     .01     2.00      .02      .49   2.0000000000000000   .0201006717070029   
.4950000000000000

0.98999999999999999112
2.00000000000000000000
9.21034037197620492066
0.00499999999999994459
2.00000000000000000000
9.21034037197620492066
0.00499999999999994459
     .99     2.00     9.21      .00   2.0000000000000000  9.2103403719762049   
.0049999999999999

As you can see, the position where we're rounding is very, very
close to .5.  With my naive grasp of numerics, I can imagine that
running on a different processor or even using a different
compiler could make a difference.  (Different compilers might
choose to convert between 80-bit FPU doubles and 64-bit memory
doubles at different times, or one might use SSE floating point
and the other 80x87 floating point.)

So the real question may be, how should we adjust the chi-square
test to avoid this edge case?  Perhaps we should just pick a
different initial value of P that does not come so close to a
rounding error.


> The doc says that COMPUTE assigns a format of F8.2 to any
> numeric target variable. Is there no way to specify a different
> format?

I looked at the SPSS manuals and this follows SPSS.  The SPSS
manuals suggest using FORMATS to change the numeric format, and
of course PRINT FORMATS is similar.

[*] I applied the following patch to print.c to output data to 20
decimal places, one per line, in addition to the format actually
requested:

@@ -946,6 +946,7 @@ print_trns_proc (void *trns_, struct cca
 
        memset (buf, ' ', t->max_width);
        len = 0;
+        putchar ('\n');
        break;
 
       case PRT_CONST:
@@ -956,6 +957,7 @@ print_trns_proc (void *trns_, struct cca
        break;
 
       case PRT_VAR:
+        printf ("%.20f\n", case_num (c, i->u.v.v->fv));
         data_out (&buf[i->fc], &i->u.v.f, case_data (c, i->u.v.v->fv));
        len = i->fc + i->u.v.f.w;
        break;

[+] Why is it done that way?  Because then I don't have to know
the proper range of values of x for the various pdf, cdf,
etc. functions.  I know the range of P and feeding it through idf
gives me the proper range for x.  This is lazy, I admit.
-- 
"Sanity is not statistical."
--George Orwell