<br><br><div class="gmail_quote">On Mon, Oct 5, 2009 at 7:29 PM, Anne Archibald <span dir="ltr">&lt;<a href="mailto:peridot.faceted@gmail.com">peridot.faceted@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
2009/10/5 Charles R Harris &lt;<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>&gt;:<br>
<div><div></div><div class="h5">&gt;<br>
&gt;<br>
&gt; On Mon, Oct 5, 2009 at 3:20 PM, Anne Archibald &lt;<a href="mailto:peridot.faceted@gmail.com">peridot.faceted@gmail.com</a>&gt;<br>
&gt; wrote:<br>
&gt;&gt;<br>
&gt;&gt; 2009/10/5 David Goldsmith &lt;<a href="mailto:d.l.goldsmith@gmail.com">d.l.goldsmith@gmail.com</a>&gt;:<br>
&gt;&gt; &gt; Just curious, Anne: have you anything in particular in mind (i.e., are<br>
&gt;&gt; &gt; there<br>
&gt;&gt; &gt; some small - or gaping - holes in scipy (IYO, of course) which you know<br>
&gt;&gt; &gt; could be filled by a careful implementation of something(s) extant in<br>
&gt;&gt; &gt; the<br>
&gt;&gt; &gt; literature)?<br>
&gt;&gt;<br>
&gt;&gt; Well, not exactly - the examples I had in mind were minor and/or in<br>
&gt;&gt; the past. I ran into problems with scipy&#39;s hyp2f1, for example, so I<br>
&gt;&gt; went and looked up the best algorithm I could find for it (and I think<br>
&gt;&gt; I contributed that code). I wanted the Kuiper test as an alternative<br>
&gt;&gt; to the Kolmogorov-Smirnov test (it&#39;s invariant under cyclic<br>
&gt;&gt; permutations, and is sensitive to different features of the<br>
&gt;&gt; distribution) so I looked up the test and the special function needed<br>
&gt;&gt; to interpret its results. (I haven&#39;t contributed this to scipy yet,<br>
&gt;&gt; mostly because I chose an interface that&#39;s not particularly compatible<br>
&gt;&gt; with that for scipy&#39;s K-S test.) And on a larger scale, that&#39;s what<br>
&gt;&gt; scipy.spatial&#39;s kdtree implementation is.<br>
&gt;&gt;<br>
&gt;&gt; For examples where I think a bit of lit review plus implementation<br>
&gt;&gt; work might help, I&#39;d say that the orthogonal polynomials could use<br>
&gt;&gt; some work - the generic implementation in scipy.special falls apart<br>
&gt;&gt; rapidly as you go to higher orders. I always implement my own<br>
&gt;&gt; Chebyshev polynomials using the cos(n*arccos(x)) expression, for<br>
&gt;&gt; example, and special implementations for the others might be very<br>
&gt;&gt; useful.<br>
&gt;&gt;<br>
&gt;<br>
&gt; At what order does the scipy implementation of the Chebyshev polynomials<br>
&gt; fall apart? I looked briefly at that package a long time ago, but never used<br>
&gt; it. I ask so I can check the chebyshev module that is going into numpy.<br>
<br>
</div></div>By n=30 they are off by as much as 0.0018 on [-1,1]; n=31 they are off<br>
by 0.1, and by n=35 they are off by four - not great for values that<br>
should be in the interval [-1,1]. This may seem like an outrageously<br>
high degree for a polynomial, but there&#39;s no reason they need be this<br>
bad, and one could quite reasonably want to use an order this high,<br>
say for function approximation.<br>
<br></blockquote><div><br>Hmm, I get an maximum error of about 1e-13 for n=100 using my routine, which isn&#39;t great and can be improved a bit with a few tricks, but is probably good enough.  That&#39;s using simple Clenshaw recursion for the evaluation. There must be something seriously wrong with the scipy version. I routinely use fits with n &gt; 50 because truncating such a series gives a good approximation to the minmax fit and it&#39;s also nice to see how the coefficients fall off to estimate the size of the truncation error.<br>
<br></div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
I think this inaccuracy is probably inevitable in a scheme that<br>
computes values using a recurrence relation</blockquote><div><br>Not so, using the Cheybshev recurrence in either direction should be stable for |x| &lt;= 1. It&#39;s like multiplying with a complex number of modulus 1, i.e., a unitary matrix.<br>
 <br></div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"> and something like it probably occurs for all the orthogonal polynomials that don&#39;t have<br>

special-purpose evaluators.<br></blockquote><div><br>Depends on the recursion, the direction of the recursion, and the domain.<br><br>Chuck<br><br></div></div>