Strange error using quad2d

2 views (last 30 days)
B
B on 9 Sep 2014
Commented: Mike Hosea on 29 Sep 2014
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.
  2 Comments
B
B on 9 Sep 2014
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.

Sign in to comment.

Accepted Answer

B
B on 28 Sep 2014
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
  3 Comments
Mike Hosea
Mike Hosea on 29 Sep 2014
You accepted your own answer instead of mine.

Sign in to comment.

More Answers (1)

Mike Hosea
Mike Hosea on 26 Sep 2014
Edited: Mike Hosea on 26 Sep 2014
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);
  1 Comment
Mike Hosea
Mike Hosea on 29 Sep 2014
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!