Issue with the root-finding fzero function

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); ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[Walter Roberson]

What is class() of each of your inputs ? Are any of them sparse?

[Amykelly]

They are all simple scalars.

Answer 1

[Walter Roberson]

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.

[Amykelly]

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.