Function floating point issue
6 views (last 30 days)
Show older comments
I made a function intd(). Given all the arguments except x, I first want to plot the function intd() in terms of x. As is shown below, I use a function handle to turn intd() into univariate h().
Suppose (r,s,q,nu,al)=(5, 2, 0.2104, 60, 0.05).Then h() gives me the below error message when I plugged in 6.699999999999999 for x. By comparisons, when I let x be 6.7, h() returns a perfectly reasonable answer. So it looks like my function has some floating point issue. This issue particularly gets in the way when I need to find a root to function h(). I used fzero(h, initial value) to find the root and got the exact same error message.
How could I make my function insensitive to the floating issue? If I had to keep my function, would there be a way to modify fzero so that it could avoid the floating point issue?
function prob=intd(x,r,s,q,nu,al)
a=q.*(1+q.^2).^(-1./2);
b=(1+q.^2).^(-1./2);
E=@(t) fcdf((r-s)./s.*(x.*(r.*t-x.^2).^(1/2)-
a.*b.*r.*t).^2./(a.^2.*r.*t-x.^2).^2,s,r-s).*fpdf(t,r,nu) ;
up=x.^2./(b.^2)./r;
down=x.^2./r;
prob=fcdf(x.^2./r,r,nu)+quadv(E,down,up)-1+al;
end
h=@(x) intd(x,r,s,q,nu,al);
h(6.699999999999999)
h(6.7)
%________________ Error Message_______________
Error using betainc Inputs must be real, full, and double or single.
Error in fcdf (line 58) p(kk) = betainc(xx, v1(kk)/2, v2(kk)/2,'lower');
Error in intd>@(t)fcdf((r-s)./s.*(x.*(r.*t-x.^2).^(1/2)-a.*b.*r.*t).^2./(a.^2.*r.*t-x.^2).^2,s,r-s).*fpdf(t,r,nu) (line 4) E=@(t) fcdf((r-s)./s.*(x.*(r.*t-x.^2).^(1/2)-a.*b.*r.*t).^2./(a.^2.*r.*t-x.^2).^2,s,r-s).*fpdf(t,r,nu) ;
Error in quadv (line 61) y{j} = feval(f, x(j), varargin{:}); %#ok<AGROW>
Error in intd (line 7) prob=fcdf(x.^2./r,r,nu)+quadv(E,down,up)-1+al;
Error in @(x)intd(x,5,2,0.2104,60,0.05)
%______________________________________________________
>> h(6.7)
ans =
0.0500
0 Comments
Answers (2)
Sean de Wolski
on 12 Feb 2013
This error doesn't look like a floating point issue but rather a real issue. I.e. at some point it looks like your data became complex.
Run:
dbstop if error
This will stop with the debugger when the error occurs, you will be able to look at the stack information and identify which of your variables has an imaginary component.
2 Comments
Sean de Wolski
on 13 Feb 2013
This is why I encourage you to use dbstop if error. It will allow you to see exactly what everything is at the time of the error.
Image Analyst
on 13 Feb 2013
You call anonymous function h, which calls intd(). You pass in only x to h() but intd() requires r,s,q,nu,al in addition. Why are you not passing in those to your anonymous function h()?
See Also
Categories
Find more on Graph and Network Algorithms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!