Fitting data on a Bessel function (of 2nd kind) with scipy

Question

I have data that looks like a Bessel function of second kind (and arguably should be described by one of these functions).

I've been trying to do so using the scipy optimize toolbox following the documentation example but without success so far: I get the following error message

ValueError: array must not contain infs or NaNs

I'd say that the divergence in 0 is the cause of the problem.

By the way, I'm trying to fit two parameter, the index of the Bessel function and a scale factor in the varialbe, a and b in K_a(bx). Is it a problem to try to fit on a discrete space (a takes value in the natural integers).

My code looks like this at the moment:

from scipy.special import yn #importing the Bessel functions
from scipy.optimize import curve_fit

def func(var, a, b):
    return yn(b*var,a)
popt, pcov = curve_fit(func, x, y) # x and y are my data points

Answer 1

First of all, you're passing to yn the arguments in the wrong order, it should be yn(a,b*var) instead of yn(b*var,a). Likely it was this mistake leading to the function yn blowing up to inf.

As a second point, as you suspected, scipy will truncate your a to a float when you call yn, raising a RuntimeWarning. You're better off optimizing only with respect to your scaling variable b, and then investigate different values for the integer order a. You can do that either by hand, or in a loop.

I'll discuss some converge issues starting with an example, fitting sin(x)/2 with yn(1,x) on [1,2*pi].

from scipy.special import yn
from scipy.optimize import curve_fit
from numpy import sin,linspace,pi

a=1#choose the order a here!
func = lambda var,b : yn(a,b*var)

x=linspace(1,2*pi,100)
y=sin(x)/2.#we fit this as an example
[b], pcov = curve_fit(func, x, y) # x and y are my data points
print b

If you now change your domain to x=linspace(0,2*pi,100)[1:], it happens that curve_fit will not converge. This is because in the part of the domain close to zero the optimization algorithm will try to squeeze yn towards the axis. This leads to large values of b, that in turn result in a strongly oscillatory behaviour of the function (try plot(x,yn(a,10*x))) up to a discretization limit (try plot(x,yn(a,10*x))) where things get even worse.

The moral is that if your data is close to zero for low x, you should start fitting a bit far away from zero to get decent convergence.

As a side note, usually Ka(x) refers to the MODIFIED Bessel function of the second kind, while yn is the Bessel function of the 2nd kind, more often referred to as Ia(x).