Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: RE: Limits to df in invnchi2


From   "Steichen, Thomas" <STEICHT@rjrt.com>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: RE: Limits to df in invnchi2
Date   Tue, 16 Sep 2003 15:11:05 -0400

Jeff Pitblado, in reply to Jeffrey Simons, writes (in part):

> There are several approximations to this (inverse non-central
> chi-square) distribution, most 
> (if not all) of which are discussed in Johnson, Kotz, and 
> Balakrishnan (1995).  One approximation is based on a Central 
> Limit Argument.
> 
> 	Assume X2 is distributed as a non-central chi-square 
> random variable
> 	with n degrees of freedome and non-centrality parameter L, then
> 
> 		Z = (X2 - (n+L))/sqrt(2*(n+2*L))
> 
> 	will tend to follow the standard normal distribution 
> for fixed n and
> 	large L, or for fixed L and large n.
> 
> I compared the results between -invnchi2()- and this 
> approximation fixing n=200, using alpha = 0.05, and L ranging 
> between 5 and 99.  The relative difference between the two 
> ranged between 0.0040 and 0.0047.  The normal approximation 
> was consistently less than the result from -invnchi2()- 
> (minimum difference was 1.12, max was 1.38).

As Jeff says, there are a number of approximations. I repeated Jeff's
analysis with the 3 approximation formulas given in Abramowitz
& Stegun (1972) and show the results below. The first approximation, 
labeled "invnchi2A", is based on the central chi-square while 
approximations "invnchi2B" and "invnchi2C" are based on the central 
normal distribution. I report Stata's invnchi2 value, the result of 
each of these formulas and the relative difference of each from 
invnchi2, using invnchi2 as the reference.

invnchi2A, based on the central chi-square, is very poor, averaging 
5.037% relative difference and ranging from 0.064% to 11.831%. It was 
lower than the Stata result throughout and became progressively more 
divergent as lambda became larger.

invnchi2B is much better, at 0.139% average relative difference and 
range 0.108% to 0.167%. It was slightly lower throughout than the 
Stata value but became closer as lambda increased.

invnchi2C is quite good, at -0.020% average relative difference and 
range -0.001% to -0.040%. It was slightly higher throughout than the 
Stata value and became more divergent as lambda increased.

If I understood Jeff's relative difference notation, then the 
approximation he used ranged in relative difference from 0.40% to 0.47%.

Approximation invnchi2C seems quite acceptable (for this example) and 
could easily be programmed. It would be reasonable, though, to test
such an approximation formula over a wide range of values before 
offering it as a stand-alone program (or, better yet, a function).

Tom Steichen

------------------------------------

local p = .05
local v = 200

forvalues l = 5/99 {

  local a = `v' + `l'
  local b = `l'/`a'
  
  local stata = invnchi2(`v',`l',`p')
  local invnchi2A = (1+`b') * invchi2(`v', `p')
  local invnchi2B = ((1+`b')/2)*(invnorm(`p') + sqrt(2*`a'/(1+`b') - 1))^2
  local invnchi2C = `a'*(invnorm(`p')*sqrt(2/9*((1+`b')/`a')) + 1 - 2/9*((1+`b')/`a'))^3
   
  di  %5.0f `l' %10.3f `stata' _c
  di %10.3f `invnchi2A' %10.6f (`stata' - `invnchi2A')/`stata' _c
  di %10.3f `invnchi2B' %10.6f (`stata' - `invnchi2B')/`stata' _c
  di %10.3f `invnchi2C' %10.6f (`stata' - `invnchi2C')/`stata'
}

lambda invnchi2 invnchi2A  rel_difA invnchi2B  rel_difB invnchi2C  rel_difC 

    5   172.494   172.383  0.000642   172.206  0.001666   172.495 -0.000009
    6   173.340   173.180  0.000921   173.050  0.001668   173.340 -0.000005
    7   174.185   173.969  0.001242   173.895  0.001667   174.187 -0.000006
    8   175.031   174.751  0.001603   174.741  0.001661   175.033 -0.000011
    9   175.880   175.525  0.002019   175.587  0.001668   175.881 -0.000004
   10   176.726   176.292  0.002457   176.434  0.001654   176.729 -0.000017
   11   177.575   177.051  0.002947   177.281  0.001653   177.578 -0.000018
   12   178.423   177.804  0.003473   178.129  0.001648   178.427 -0.000022
   13   179.272   178.549  0.004033   178.978  0.001639   179.277 -0.000030
   14   180.124   179.287  0.004642   179.828  0.001642   180.128 -0.000026
   15   180.972   180.019  0.005268   180.678  0.001626   180.980 -0.000041
   16   181.824   180.744  0.005940   181.529  0.001622   181.832 -0.000043
   17   182.678   181.462  0.006658   182.380  0.001630   182.684 -0.000035
   18   183.529   182.173  0.007390   183.232  0.001619   183.537 -0.000044
   19   184.384   182.878  0.008165   184.085  0.001620   184.391 -0.000042
   20   185.235   183.577  0.008953   184.938  0.001603   185.246 -0.000058
   21   186.092   184.269  0.009797   185.792  0.001612   186.101 -0.000047
   22   186.943   184.955  0.010637   186.646  0.001589   186.956 -0.000069
   23   187.800   185.635  0.011532   187.501  0.001592   187.812 -0.000064
   24   188.657   186.308  0.012450   188.357  0.001592   188.669 -0.000062
   25   189.514   186.976  0.013392   189.213  0.001590   189.526 -0.000064
   26   190.371   187.638  0.014356   190.070  0.001584   190.384 -0.000068
   27   191.228   188.294  0.015343   190.927  0.001575   191.242 -0.000075
   28   192.085   188.944  0.016350   191.785  0.001564   192.101 -0.000084
   29   192.942   189.589  0.017378   192.643  0.001550   192.961 -0.000097
   30   193.802   190.228  0.018440   193.502  0.001547   193.820 -0.000097
   31   194.661   190.861  0.019520   194.361  0.001542   194.681 -0.000101
   32   195.521   191.489  0.020620   195.221  0.001534   195.542 -0.000107
   33   196.381   192.112  0.021736   196.081  0.001524   196.403 -0.000115
   34   197.243   192.729  0.022884   196.942  0.001525   197.265 -0.000112
   35   198.105   193.341  0.024049   197.804  0.001524   198.128 -0.000112
   36   198.965   193.948  0.025215   198.665  0.001506   198.990 -0.000127
   37   199.828   194.550  0.026411   199.528  0.001500   199.854 -0.000132
   38   200.693   195.147  0.027635   200.391  0.001505   200.718 -0.000125
   39   201.555   195.738  0.028860   201.254  0.001494   201.582 -0.000133
   40   202.418   196.325  0.030099   202.118  0.001481   202.447 -0.000145
   41   203.280   196.907  0.031352   202.982  0.001465   203.312 -0.000158
   42   204.148   197.484  0.032643   203.847  0.001474   204.178 -0.000147
   43   205.010   198.056  0.033921   204.712  0.001455   205.044 -0.000164
   44   205.878   198.624  0.035236   205.578  0.001459   205.911 -0.000157
   45   206.743   199.187  0.036551   206.444  0.001449   206.778 -0.000166
   46   207.609   199.745  0.037876   207.311  0.001436   207.645 -0.000176
   47   208.477   200.299  0.039225   208.178  0.001434   208.513 -0.000176
   48   209.344   200.849  0.040583   209.045  0.001431   209.382 -0.000177
   49   210.212   201.394  0.041952   209.913  0.001425   210.250 -0.000181
   50   211.080   201.934  0.043330   210.781  0.001417   211.120 -0.000186
   51   211.948   202.471  0.044717   211.650  0.001408   211.989 -0.000193
   52   212.816   203.003  0.046112   212.519  0.001396   212.859 -0.000203
   53   213.684   203.531  0.047516   213.389  0.001383   213.730 -0.000214
   54   214.555   204.054  0.048941   214.259  0.001380   214.601 -0.000214
   55   215.425   204.574  0.050372   215.129  0.001376   215.472 -0.000216
   56   216.296   205.089  0.051811   216.000  0.001370   216.344 -0.000219
   57   217.167   205.601  0.053257   216.871  0.001362   217.216 -0.000225
   58   218.037   206.109  0.054710   217.743  0.001352   218.088 -0.000232
   59   218.908   206.612  0.056169   218.615  0.001341   218.961 -0.000241
   60   219.781   207.112  0.057646   219.487  0.001340   219.834 -0.000239
   61   220.655   207.608  0.059128   220.360  0.001338   220.708 -0.000239
   62   221.528   208.100  0.060616   221.233  0.001334   221.582 -0.000241
   63   222.402   208.589  0.062109   222.106  0.001328   222.456 -0.000244
   64   223.275   209.073  0.063607   222.980  0.001321   223.331 -0.000249
   65   224.149   209.554  0.065109   223.855  0.001312   224.206 -0.000255
   66   225.022   210.032  0.066616   224.729  0.001301   225.081 -0.000263
   67   225.895   210.506  0.068127   225.604  0.001289   225.957 -0.000273
   68   226.769   210.976  0.069642   226.480  0.001275   226.833 -0.000284
   69   227.648   211.443  0.071184   227.355  0.001285   227.710 -0.000272
   70   228.521   211.906  0.072706   228.231  0.001268   228.587 -0.000286
   71   229.397   212.366  0.074242   229.108  0.001262   229.464 -0.000290
   72   230.273   212.823  0.075782   229.985  0.001254   230.341 -0.000295
   73   231.152   213.276  0.077335   230.862  0.001257   231.219 -0.000289
   74   232.026   213.726  0.078869   231.739  0.001235   232.097 -0.000309
   75   232.905   214.173  0.080428   232.617  0.001235   232.976 -0.000306
   76   233.784   214.616  0.081988   233.495  0.001234   233.855 -0.000305
   77   234.662   215.056  0.083550   234.374  0.001231   234.734 -0.000305
   78   235.539   215.493  0.085104   235.252  0.001215   235.614 -0.000319
   79   236.418   215.927  0.086670   236.132  0.001210   236.494 -0.000322
   80   237.296   216.358  0.088237   237.011  0.001203   237.374 -0.000326
   81   238.178   216.786  0.089816   237.891  0.001206   238.254 -0.000320
   82   239.057   217.211  0.091386   238.771  0.001196   239.135 -0.000327
   83   239.936   217.632  0.092956   239.651  0.001186   240.016 -0.000335
   84   240.815   218.051  0.094528   240.532  0.001174   240.898 -0.000344
   85   241.696   218.467  0.096111   241.413  0.001172   241.779 -0.000343
   86   242.578   218.880  0.097694   242.295  0.001169   242.662 -0.000344
   87   243.457   219.290  0.099267   243.176  0.001153   243.544 -0.000357
   88   244.341   219.697  0.100861   244.058  0.001159   244.427 -0.000349
   89   245.220   220.101  0.102434   244.941  0.001141   245.310 -0.000364
   90   246.105   220.503  0.104028   245.823  0.001144   246.193 -0.000358
   91   246.984   220.902  0.105602   246.706  0.001124   247.076 -0.000376
   92   247.868   221.298  0.107195   247.589  0.001125   247.960 -0.000372
   93   248.752   221.691  0.108788   248.473  0.001124   248.844 -0.000370
   94   249.637   222.082  0.110380   249.357  0.001123   249.729 -0.000369
   95   250.518   222.470  0.111962   250.241  0.001109   250.614 -0.000380
   96   251.400   222.855  0.113543   251.125  0.001095   251.499 -0.000392
   97   252.285   223.238  0.115133   252.010  0.001090   252.384 -0.000394
   98   253.169   223.618  0.116722   252.895  0.001084   253.269 -0.000397
   99   254.053   223.996  0.118310   253.780  0.001077   254.155 -0.000401
Averages:
        212.971   201.387    5.037%   212.678    0.139%   213.016   -0.020%

----------------------------------------



-----------------------------------------
CONFIDENTIALITY NOTE:  This e-mail message, including any
attachment(s), contains information that may be confidential,
protected by the attorney-client or other legal privileges, and/or
proprietary non-public information.  If you are not an intended
recipient of this message or an authorized assistant to an intended
recipient, please notify the sender by replying to this message and
then delete it from your system.  Use, dissemination, distribution,
or reproduction of this message and/or any of its attachments (if
any) by unintended recipients is not authorized and may be unlawful.


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2022 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index