[SciPy-dev] Numerical Recipes (was tagging 0.7rc1 this weekend?)
Jarrod Millman
millman@berkeley....
Wed Dec 3 14:35:44 CST 2008
Given the discussion about Numerical Recipes derived-code, I did a
little grepping and found 7 instances where code is referencing
Numerical Recipes. Before branching for 0.7, I would like to see this
resolved. I am not sure what we should do, but the license for using
Numerical Recipes derived code clearly does not permit us to include
it in scipy:
http://www.nr.com/licenses/redistribute.html
Below I include all the instances I could easily find. Thoughts?
--------------------------------------------------
In scipy/optimize/zeros.py:
def ridder(f, a, b, args=(),
xtol=_xtol, rtol=_rtol, maxiter=_iter,
full_output=False, disp=False):
"""Find root of f in [a,b]
Ridder routine to find a zero of the function
f between the arguments a and b. f(a) and f(b) can not
have the same signs. Faster than bisection, but not
generaly as fast as the brent rountines. A description
may be found in a recent edition of Numerical Recipes.
The routine here is modified a bit to be more careful
of tolerance.
def brentq(f, a, b, args=(),
xtol=_xtol, rtol=_rtol, maxiter=_iter,
full_output=False, disp=False):
"""Find root of f in [a,b]
The classic Brent routine to find a zero of the function
f between the arguments a and b. f(a) and f(b) can not
have the same signs. Generally the best of the routines here.
It is a safe version of the secant method that uses inverse
quadratic extrapolation. The version here is a slight
modification that uses a different formula in the extrapolation
step. A description may be found in Numerical Recipes, but the
code here is probably easier to understand.
In scipy/signal/medianfilter.c:
/*
* This Quickselect routine is based on the algorithm described in
* "Numerical recipies in C", Second Edition,
* Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
*/
In scipy/special/orthogonal.py:
"""
A collection of functions to find the weights and abscissas for
Gaussian Quadrature.
These calculations are done by finding the eigenvalues of a
tridiagonal matrix whose entries are dependent on the coefficients
in the recursion formula for the orthogonal polynomials with the
corresponding weighting function over the interval.
Many recursion relations for orthogonal polynomials are given:
a1n f_n+1 (x) = (a2n + a3n x ) f_n (x) - a4n f_n-1 (x)
The recursion relation of interest is
P_n+1 (x) = (x - A_n) P_n (x) - B_n P_n-1 (x)
where P has a different normalization than f.
The coefficients can be found as:
A_n = -a2n / a3n
B_n = ( a4n / a3n sqrt(h_n-1 / h_n))**2
where
h_n = int_a^b w(x) f_n(x)^2
assume:
P_0(x) = 1
P_-1(x) == 0
See Numerical Recipies in C, page 156 and
Abramowitz and Stegun p. 774, 782
In scipy/stats/stats.py:
def ttest_ind(a, b, axis=0):
"""Calculates the t-obtained T-test on TWO INDEPENDENT samples of scores
a, and b. From Numerical Recipies, p.483. Axis can equal None (ravel
array first), or an integer (the axis over which to operate on a and b).
def ttest_rel(a,b,axis=None):
"""Calculates the t-obtained T-test on TWO RELATED samples of scores, a
and b. From Numerical Recipies, p.483. Axis can equal None (ravel array
first), or an integer (the axis over which to operate on a and b).
def ks_2samp(data1, data2):
""" Computes the Kolmogorov-Smirnof statistic on 2 samples. Modified
from Numerical Recipies in C, page 493. Returns KS D-value, prob. Not
ufunc- like.
More information about the Scipy-dev
mailing list