How to fastly compute a double integral of a self-defined function?

Hi everyone,

I need to compute , where is my self-defined function, M is infinity and in actual calculation I set it to be a large number and is the probability density function of standard bivariate normal. My code for a simple example looks as follows (here ):

clear;
tic
step=0.01; 
grid=-100:step:100; lg=length(grid);
x=repmat(grid',1,lg);
y=repmat(grid,lg,1);
tempz=[x(:) y(:)];
z=tempz(:,1).*tempz(:,2);
fz=mvnpdf(tempz);
t=sum(z.*fz);
toc

This simple example already takes me 30 seconds on a computing node of our university. In the real application is much more complex and defined by a separate function file. I'm wondering if there are any more efficient ways of doing this?

Answer 1

[Steven Lord]

Use the integral2 function.

[Tmat]

Thanks, Steven. I've tried it early on, it didn't work because of the matlab function mvnpdf(). See the simple example below:

@(x,y)mvnpdf(x,y)

>> integral2(fun,0,1,0,1)
Error using integral2Calc>integral2t/tensor (line 241)
Integrand output size does not match the input size.

[Steven Lord]

You're not calling mvnpdf correctly. The two-input form of mvnpdf treats the first input as coordinates and the second as the means. Use the x and y arrays integral2 passes into your integrand function to assemble the X matrix for a call to the one-input form of mvnpdf.

[Tmat]

Thanks, Steven. I also tried the following:

vfun=@(x)mvnpdf(x);
intvfun=integral2(vfun,-inf,0,-inf,0);

Again, it didn't work. What is a proper way of using it in my situation?

[Steven Lord]

fh = @(x, y) reshape(mvnpdf([x(:), y(:)]), size(x));
integral2(fh, -Inf, Inf, -Inf, Inf)

fh combines the two coordinate vectors integral2 will pass into it into the coordinate matrix mvnpdf requires, evaluates the PDF, and uses reshape to make the PDF values vector the same shape as the input.

[Tmat]

Thanks a lot, Steven. This worked!