Marcum Q

[Karol Kowalik]

Dear All

I would like to implement Marcum Q function of a and b defined as:

Q(a,b) = integral from b to inf of [x * exp(-(x^2+a^2)/2) * I_0(ax)] dx

where: I_0(ax) is is the modified Bessel function of the first kind of 
zero order.


Thus I have implemented Marcum Q in Octave as:

function mar = marcumq (a, b)       #
  f = @(x) x .* exp(-((x.^2) .+ (a.^2))./2) .* besseli(0, a .* x);
  mar=quadl(f,b,Inf);
endfunction

However, when I try to compute it for example for a=2 and b=4, then I 
the result is NaN:

octave:2> marcumq(2,4)
a 2.000000 b 4.000000
warning: quadl: interval contains no more machine number
warning: quadl: required tolerance may not be met
marcumq NaN
ans = NaN + NaNi

Could you please help me to find out what is wrong with my implementation?

Many thanks
Karol