[SciPy-User] scipy.interpolate.rbf sensitive to input noise ?

[denis]

A followup to the 15feb thread "creating a 3D surface plot from
collected data":
the doc for scipy.interpolate.rbf says that 3 of the 7 radial
functions
use a scale paramater `epsilon`, 4 do not.
Odd -- how can some be scale-free ?
Running rbf on 100 1d random.uniform input points
(after changing the linalg.solve in rbf.py to lstsq)
shows that all but linear and gaussian are noisy,
giving interpolants way outside the input range:

    N: 100  ngrid: 50  sigma: 0.1
    # min  max  max |delta|  Rbf
     -1.0  1.0   0.5  gaussian
     -1.0  1.2   0.3  linear
     -1.4  2.1   2.5  thin-plate
     -1.0  2.2   3.0  inverse multiquadric
     -1.0  2.4   3.2  multiquadric
     -2.1 12.8   7.4  quintic
    -10.6  7.0  11.7  cubic

Should Rbf be moved off to noisymethods/...  until experts can look it
over,
or have I done something stupid ?
(A back-of-the-envelope method for choosing epsilon aka r0 would be
nice,
and a rough noise amplification model would be nice too.)

Also, the doc should have a link to
http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
and should say that the initial rbf() is N^2.

cheers
  -- denis

<pre><code>
    # from http://scipy.org/Cookbook/RadialBasisFunctions
    import sys
    import numpy as np
    from rbf import Rbf  # lstsq
    # from scipy.interpolate import Rbf
    Rbfuncs = ( 'gaussian', 'linear', 'thin-plate', 'inverse
multiquadric',
        'multiquadric', 'quintic', 'cubic', )

    N = 100
    ngrid = 50
    sigma = .1
    exec "\n".join( sys.argv[1:] )  # N=
    np.random.seed(1)
    np.set_printoptions( 2, threshold=50, edgeitems=5,
suppress=True )  # .2f

    x = np.sort( np.random.uniform( 0, 10, N+1 ))
    y = np.sin(x) + np.random.normal( 0, sigma, x.shape )
    xi = np.linspace( 0, 10, ngrid+1 )

    print "N: %d  ngrid: %d  sigma: %.2g" % (N, ngrid, sigma)
    print "# min  max  max |delta|  Rbf"
    for func in Rbfuncs:
        rbf = Rbf( x, y, function=func )
        fi = rbf(xi)
        maxdiff = np.amax( np.abs( np.diff( fi )))
        print "%5.1f %4.1f  %4.1f  %s" % (
            fi.min(), fi.max(), maxdiff, func )
</pre></code>