Re: [OctDev] regress and rcoplot (partially written, please help)

[William P. Y. H. <wil...@gm...>]

On 2/20/06, Paul Kienzle <pki...@us...> wrote:
>
> Okay, I added a test of the Longley data set from the NIST Statistical
> Reference Database.
>
> The performance is pretty bad (norm of answer-certified value is about
> 10^-7), but that's a bit better than the competition.
>
> Someone with more R experience than me will have to see if it performs
> any better.
>
> Also, if you can convert bint to standard error values that would be
> convenient.
>
>         - Paul
>
> ----
>
> % Longley data from the NIST Statistical Reference Dataset
> Z =3D [  60323    83.0   234289   2356     1590    107608  1947
>         61122    88.5   259426   2325     1456    108632  1948
>         60171    88.2   258054   3682     1616    109773  1949
>         61187    89.5   284599   3351     1650    110929  1950
>         63221    96.2   328975   2099     3099    112075  1951
>         63639    98.1   346999   1932     3594    113270  1952
>         64989    99.0   365385   1870     3547    115094  1953
>         63761   100.0   363112   3578     3350    116219  1954
>         66019   101.2   397469   2904     3048    117388  1955
>         67857   104.6   419180   2822     2857    118734  1956
>         68169   108.4   442769   2936     2798    120445  1957
>         66513   110.8   444546   4681     2637    121950  1958
>         68655   112.6   482704   3813     2552    123366  1959
>         69564   114.2   502601   3931     2514    125368  1960
>         69331   115.7   518173   4806     2572    127852  1961
>         70551   116.9   554894   4007     2827    130081  1962 ];
> % Results certified by NIST using 100 digit arithmetic
> % b and standard error in b
> V =3D [  -3482258.63459582         890420.383607373
>          15.0618722713733         84.9149257747669
>         -0.358191792925910E-01    0.334910077722432E-01
>         -2.02022980381683         0.488399681651699
>         -1.03322686717359         0.214274163161675
>         -0.511041056535807E-01    0.226073200069370
>          1829.15146461355         455.478499142212 ];
> Rsq =3D 0.995479004577296;
> F =3D 330.285339234588;
> y =3D Z(:,1); X =3D [ones(rows(Z),1), Z(:,2:end)];
> [b, bint, r, rint, stats] =3D regress (y, X, 0.05);
> assert(b, V(:,1), 6e-7)
> assert(stats(1), Rsq, 7e-14)
> assert(stats(2), F, 5e-9)
>
>

Yes, you can convert bint to the standard error values with (0.05 is
alpha, 9 is the degree of freedom):
((bint(:, 1) - bint(:, 2)) / 2) / tinv (0.05 / 2, 9)
but the accuracy will be affected. (It's not so accurate anyway...)

--
William Poetra Yoga Hadisoeseno