Issue with the root-finding fzero function

2 views (last 30 days)
Amykelly
Amykelly on 2 Feb 2013
I made this below function to obtain a critical value, which is a root to another function involving an integral. I was supposed to get an answer to crit(5, 2, 0.2104, 60, 0.05) but Matlab gives me this error message(bounded by ~~~~~~~~). Could somebody figure out where this problem comes from, the integration or the root-finding method fzero? fzero looks more suspicious to me. I wonder if replacing fzero with another root-finding method would resolve the issue.
Your suggestions are highly appreciated.
Amy
function c=crit(r,s,q,nu,al)
a=q.*(1+q.^2).^(-1./2);
b=(1+q.^2).^(-1./2);
h=@intd;
function prob=intd(x)
E=@(t) fcdf((r-s)./s.*(sqrt(x).*(r.*t-x).^(1/2)- a.*b.*r.*t).^2./(a.^2.*r.*t-x).^2,s,r-s).*fpdf(t,r,nu) ;
up=x./(b.^2)./r;
down=x./r;
prob=fcdf(x/r,r,nu)+quadv(E,down,up)-1+al;
end
c=fzero(h,r*finv(1-al,r,nu),1e-6);
end
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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 crit>@(t)fcdf((r-s)./s.*(sqrt(x).*(r.*t-x).^(1/2)-a.*b.*r.*t).^2./(a.^2.*r.*t-x).^2,s,r-s).*fpdf(t,r,nu) (line 7) E=@(t) fcdf((r-s)./s.*(sqrt(x).*(r.*t-x).^(1/2)-a.*b.*r.*t).^2./(a.^2.*r.*t-x).^2,s,r-s).*fpdf(t,r,nu) ;
Error in quadv (line 61) y{j} = feval(f, x(j), varargin{:}); %#ok<AGROW>
Error in crit/intd (line 10) prob=fcdf(x/r,r,nu)+quadv(E,down,up)-1+al;
Error in fzero (line 343) a = x - dx; fa = FunFcn(a,varargin{:});
Error in crit (line 13) c=fzero(h,r*finv(1-al,r,nu),1e-6); ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  2 Comments
Walter Roberson
Walter Roberson on 2 Feb 2013
What is class() of each of your inputs ? Are any of them sparse?

Sign in to comment.

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 3 Feb 2013
If you "dbstop if error" and look at the class() of the parameters to betainc(), then are they all double? If they are, check to see if one of them has an imaginary component. Is your intd() well defined if x is negative? You might want to put a range for your x0, to prevent x from going negative.
Why are you passing three arguments to fzero? Your third argument is not an options structure.
  1 Comment
Amykelly
Amykelly on 6 Feb 2013
You are right. When the error occured, one of the inputs to betainc() was a complex number although all the inputs were double.
I attempted restricting x to a positive range inside fzero(). The same error message again! The debugging result confirmed that it was the complex input that caused the trouble. I looked the function intd() over. If x is straight positive, there is really no chance for the intergrand or domain ranges to have imaginary units.
Then I write intd() into an independent function, as follows.
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.*(sqrt(x).*(r.*t-x).^(1/2)- a.*b.*r.*t).^2./(a.^2.*r.*t-x).^2,s,r-s).*fpdf(t,r,nu) ;
up=x./(b.^2)./r;
down=x./r;
prob=fcdf(x/r,r,nu)+quadv(E,down,up)-1+al;
end
Use a function handle to turn intd() into univariate h().
h=@(x) intd(x,r,s,q,nu,al);
fplot(h,[1 10])
Suppose (r,s,q,nu,al)=(5, 2, 0.2104, 60, 0.05).Then h() gives me all reasonable answers when I attemp different values for x, e.g. the sequence (1:0.5:10). However, when I tried to plot h() by fplot(h,[1 10]), it fails on the same error message. Does this suggest that intd() or h() has no problems but fzero() or fplot() may impose some complex numbers on the inputs to intd() or h()? I am puzzled.

Sign in to comment.

Categories

Find more on Software Development Tools in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!