2 views (last 30 days)
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);
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
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.

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
Mike Hosea on 29 Sep 2014

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...
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.