How can I integrate a function related to the modified bessel function of the first kind?

Thank you for reading my question!

The thing is that I'm tring to integrate a rice PDF, and it should be 1.

beta = 5 ;
sigma = 1.4/3;
distribution_func = @(x) x./sigma.^2.*exp(-(x.^2+beta.^2)./2./sigma.^2).*besseli(0,beta.*x./sigma.^2);
integral(distribution_func,0,100);

Warning: Inf or NaN value encountered.

As you can see, due to can go to Inf in a short range, besseli function can't be used.

And I think might be big alone, but it would be small with the scaling in distribution_func. So I try to constuct my own bessel function and it's like:

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) integral(@(theta) integ_func(theta,x),0,pi);
integral(distribution_func,0,1);

Error using integralCalc/finalInputChecks
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

The equation is right actually, but it seems that Integral can't be used in this way.

So I try trapz:

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) integral(@(theta) integ_func(theta,x),0,pi);
x = 0:0.01:1e3;
y = zeros(1,length(x));
for i = 1:length(x)
    y(i) = distribution_func(x(i));
end
disp(trapz(x,y));

It' ok, but I don't like it, I think it's clumsy and inaccurate under certain situations.

The last method I may try is double integral since there is two integration in my equation.

Could you give me some advice about some easy and accurate ways to reach my goal?

Thanks in advance!!

Answer 1

[Bruno Luong]

You can also "vectorize" distribution_func by applying arrayfun

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) arrayfun(@(x)integral(@(theta) integ_func(theta,x),0,pi), x);
integral(distribution_func,0,1)

[祥宇 崔]

Thanks!

I'll try yo figure out what is Arrayfun.

Answer 2

[祥宇 崔]

It seems that this works:

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) integral(@(theta) integ_func(theta,x),0,pi,'ArrayValued',true);
integral(distribution_func,0,10);

[祥宇 崔]

hh, it works!

Answer 3

[Sam Chak]

Hi @祥宇 崔

Have you checked if this besseli() function for Modified Bessel function of first kind is applicable for your case?

https://www.mathworks.com/help/matlab/ref/besseli.html

[祥宇 崔]

Hi,Chak!

Maybe not

But I think I can use the scale term in besseli to make it applicable