Re: st: non-central chi-square?

[Dirk Enzmann <dirk.enzmann@uni-hamburg.de>]

About 2.5 years ago, Stas Kolenikov wrote:

http://www.stata.com/statalist/archive/2010-11/msg00781.html

I could not find an answer or solution, therefore I tried the following 
(it is not the solution Stas asked for, but - if correct, which I hope 
so - it is a quick and not too dirty way around the problem):

* --- start Stata code --------------------------

program define chi2den, rclass
  syntax, df(integer) x(real) [ ncp(integer 0) ]
  tempname chi2d
  sca `chi2d' = (nchi2(`df',`ncp',(1000*`x'+0.05)/1000) - ///
                 nchi2(`df',`ncp',(1000*`x'-0.05)/1000))*10000
  if `x'==0 & `df'==2 sca `chi2d' = 2*`chi2d'
  di `chi2d'
  return scalar density = `chi2d'
  return scalar ncp = `ncp'
  return scalar x = `x'
  return scalar df = `df'
end

* Example (chi² = 4.06, df = 2, ncp = 3):
chi2den, df(2) x(4.06) ncp(3)

* --- end Stata code ----------------------------

If you are looking for the density of a Chi² value of the central Chi² 
distribution simply omit the non-centrality parameter option -ncp()-.

Surely, this can be impoved (but not by me).

You can use a similar strategy to create density plots of the 
(noncentral) Chi² distribution:

* --- start Stata code --------------------------

cap program drop plotchi2den
program define plotchi2den
  syntax, df(integer) range(numlist min=2 max=2) ///
          [ ncp(integer 0) n(integer 300) * ]
  tempname nchi
  if `n' > 399 set matsize `=`n'+1'
  local from : word 1 of `range'
  local to : word 2 of `range'
  local interv = (`to'-`from')/`n'
  matrix `nchi' = J(`=`n'+1',2,.)
  forval i = 0/`n' {
    matrix `nchi'[`=`i'+1',1] = ///
           (nchi2(`df',`ncp',(`n'*`i'*`interv'+0.5)/`n') - ///
            nchi2(`df',`ncp',(`n'*`i'*`interv'-0.5)/`n'))*`n'
    if `i'==0 & `df' == 2 matrix `nchi'[1,1] = 2*`nchi'[1,1]
    matrix `nchi'[`=`i'+1',2] = `i'*`interv'
  }
  _matplot `nchi', recast(line) `options'
end

* Example df = 2, ncp = 3, range(chi²) = 0 - 25:
plotchi2den, df(2) ncp(3) range(0 25) ///
   title("Non-Central Chi{superscript: 2} Distribution") ///
   ytitle("Density") ylab(, angle(0)) ///
   xtitle("Chi{superscript: 2}") ///
   note("{it:df} = 2, {it:ncp} = 3")

* --- end Stata code ----------------------------

Of coure, an official (Stata) function -nchi2den- would be better.

Dirk


========================================
Dr. Dirk Enzmann
Institute of Criminal Sciences
Dept. of Criminology
Rothenbaumchaussee 33
D-20148 Hamburg
Germany

phone: +49-(0)40-42838.7498 (office)
       +49-(0)40-42838.4591 (Mrs Billon)
fax:   +49-(0)40-42838.2344
email: dirk.enzmann@uni-hamburg.de
http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Enzmann.html
========================================
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/