Strange error using quad2d

I am trying to compute an integral using quad2d. The function that I'm integrating is also computed using quad2d. I suspect this might be the problem. I get the following error message

??? Error using ==> quad2d at 124 D must be a finite, scalar, floating point constant or a function handle.

Error in ==> H at 5 F=quad2d(@(x,z) fxw(x,z,lam,mu,sig),lb,ub,lb,w);

Error in ==> @(x,w)fxw(x,w,lam,mu,sig).*H(w,mu,lam,sig,V).*w

Error in ==> quad2d>tensor at 355 Z = FUN(X,Y); NFE = NFE + 1;

Error in ==> quad2d at 169 [Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in ==> QtleMax at 94 axm = quad2d(@(x,w) fxw(x,w,lam,mu,sig).*H(w,mu,lam,sig,V).*w , lb, ub, lb, ub)

Here is my code

close all;
clc;
clear all;

lam=[0.5 0.2 0.3];
 mu=[-2 1 2];
 sig=[1 0.5 1.5];
 avmu=lam*mu';
 mu = mu -avmu;

lb = min(mu)-6*max(max(sig),2);
ub = max(mu)+6*max(max(sig),2);
V = 38;
axm = quad2d(@(x,w) fxw(x,w,lam,mu,sig).*Hx(V,Fw(w,mu,lam,sig)).*w , lb, ub, lb, ub)

%%%%%%

function y = fx(x,lam,mu,sig) 
    k1=max(size(lam));
    k2=max(size(mu));
    k3=max(size(sig));
    if k1==k2 && k2==k3
        y=0;
        for i=1:k1
            y=y + lam(i)*normpdf(x,mu(i),sig(i));
        end
    end

function y = fxw(x,w,lam,mu,sig) 
y=normpdf(w-x,0,1).*fx(x,lam,mu,sig);
end

function y=Fw(w,lam,mu,sig) 
lb = min(mu)-6*max(max(sig),2);
ub = max(mu)+6*max(max(sig),2);
    y=quad2d(@(x,z) fxw(x,z,lam,mu,sig),lb,ub,lb,w);
end

function pr = Hx(V,F) 
odd=mod(V,2);
if odd==1
    pr= exp(sum(log(1:1:V-1)) - 2*sum(log(1:1:(V-1)/2)) + ((V-1)/2).*(log(F)+log(1-F)));
else
    pr =exp(sum(log(1:1:V-1)) - sum(log(1:1:(V/2-1))) - sum(log(1:1:(V/2)))-log(2) + (V/2-1).*(log(F) + log(1-F)));
end

Please help me to figure out what the problem is

Thanks in advance, B.

[Star Strider]

I don’t understand what you’re doing, so I won’t attempt to run your code.

For troubleshooting purposes, I would create an indpendent anonymous function from your integrand:

ig = @(x,w) fxw(x,w,lam,mu,sig).*Hx(V,Fw(w,mu,lam,sig)).*w;

and comment-out the quad2d call, until I got ‘ig’ (or whatever you choose to call it) running independently and producing output appropriate to what quad2d wants. (It’s always a good idea to test a function outside of functions that call it to be sure it returns what they want.) You can then pass ‘ig’ by name to quad2d, knowing that it works.

[B]

Thanks for your tip, I am computing the expected value of complicated distribution (median of a convolution of a Gaussian mixture and a Gaussian random variable).

The ig function does work, but oddly. It returns complex numbers. If I compute fxw() and Hx() separately, I get reasonably answers. I do not understand why this is happening. If you have any clue, please let me know.

Answer 1

[B]

Thank you Mike for your reply. It's good to know a proper vectorization using quad2d. However, after setting parameters for the model, I still run the code and the answer is an imaginary number -0.0266 - 0.0000i. I have played around changing parameters and I always get an imaginary component equal to zero, so it seems to me more like a number format issue rather than programming problem. I'm thinking of just taking the real part of the answer MATLAB gives me, although I would like to understand why this is happening.

Best

[Mike Hosea]

This comment moved to my answer.

[B]

Dear Mike, You nailed it! Now it works F is never supposed to be over 1 because it is a cumulative distribution function. I guess a slight rounding error returned a log of a negative number. Thank you very much!

Answer 2

[Mike Hosea]

I don't know when it will finish, but I made this adjustment to properly vectorize the Fw function (QUAD2D can't accept an array as a limit), and it is cranking away...

y=arrayfun(@(wscalar)quad2d(@(x,z)fxw(x,z,lam,mu,sig),lb,ub,lb,wscalar),w);

[Mike Hosea]

Instead of entering a new answer, it's best to use the comment field below the answer you're commenting on.

The complex results are coming from a calculation in Hx of the form n*(log(F) + log(1-F)), where n is a positive even number.

>> exp(18*(log(1.6) + log(1 - 1.6)))
ans =
  0.479603335372623 - 0.000000000000001i

You can also get F < 0 but very small because of roundoff error in the QUAD2D result when the integrand in the inner integration close to zero throughout the entire region. A simple fix that will speed things up slightly is to tack the following line on to the Hx function:

pr = real(pr);

Incidentally, MATLAB has some corrective code built-in that makes the following real:

>> (1.6*(1-1.6)).^18
ans =
   0.479603335372623

which you could exploit by factoring out that term, but I would just zero out the imaginary part there.