<p dir="auto">The implementation of the non-central chi-squared distribution in scipy (scipy.stats.ncx2, scipy.special.chndtr) is not accurate for large values of the non-centrality parameter. E.g. consider</p>
<h3 dir="auto">Reproducing Code example:</h3>
<p dir="auto"><code class="notranslate">from scipy.stats import ncx2</code></p>
<p dir="auto"><code class="notranslate">print(ncx2.cdf(1200,2,1000))</code><br>
<code class="notranslate">print(ncx2.cdf(1500,2,1000))</code><br>
<code class="notranslate">print(ncx2.cdf(2000,2,1000))</code><br>
<code class="notranslate">print(ncx2.cdf(5000,2,1000))</code></p>
<p dir="auto">which produces</p>
<p dir="auto">0.998604279948<br>
0.999933004449<br>
0.999933004449<br>
0.999933004449</p>
<p dir="auto">instead of converging to one.</p>
<p dir="auto">The <a href="https://github.com/scipy/scipy/blob/master/scipy/special/cdflib/cdfchn.f">implementation</a> in scipy is based on Formula 26.4.25 of Abramowitz and Stegun, Handbook of Mathematical Functions (1966).</p>
<p dir="auto">There are faster and more accurate algorithms for implementing this function, especially for large non-centrality parameter (ncp). For example, <a href="http://stat.ethz.ch/R-manual/R-devel/library/stats/html/Chisquare.html" rel="nofollow">R</a> apparently uses Ding (1992) for ncp&gt;80. However <a href="https://www.boost.org/doc/libs/1_38_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/nc_chi_squared_dist.html" rel="nofollow">boost</a> states that Ding can take a long time to converge for large ncp, and recommends a different method. I am wondering if it is worth upgrading the implementation in scipy to work for large non-centrality?</p>
<h3 dir="auto">Scipy/Numpy/Python version information:</h3>
<blockquote>
<blockquote>
<blockquote>
<p dir="auto">import sys, scipy, numpy; print(scipy.<strong>version</strong>, numpy.<strong>version</strong>, sys.version_info)<br>
0.18.0 1.11.1 sys.version_info(major=3, minor=5, micro=2, releaselevel='final', serial=0)</p>
</blockquote>
</blockquote>
</blockquote>
<p dir="auto">I'm getting unexpected behavior from <code class="notranslate">scipy.ncx2.sf</code>.  The survival function should eventually go to zero for large enough inputs, but instead I'm seeing a kind of "foot", i.e., it is bottoming out at a positive value.  The following interactive I/O demonstrates the issue.</p>
<h3 dir="auto">Reproducing code <strong>example:</strong></h3>
<p dir="auto">(Pdb) v_T_norm[::50]<br>
array([  4.99515382e+00,   1.07617327e+01,   2.31854502e+01,<br>
4.99515382e+01,   1.07617327e+02,   2.31854502e+02,<br>
4.99515382e+02,   1.07617327e+03,   2.31854502e+03,<br>
4.99515382e+03,   1.07617327e+04,   2.31854502e+04,<br>
4.99515382e+04])<br>
(Pdb) nu<br>
20.0<br>
(Pdb) lam<br>
499.51538166556196<br>
(Pdb) ncx2.sf(v_T_norm[::50], df=nu, nc=lam)<br>
array([  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,<br>
1.00000000e+00,   1.00000000e+00,   1.00000000e+00,<br>
6.64666099e-01,   4.16325934e-05,   4.16325934e-05,<br>
4.16325934e-05,   4.16325934e-05,   4.16325934e-05,<br>
4.16325934e-05])</p>