How to solve the common problem about the using fzero() to solve a integration expression?

Hi everyone,

I have encountered the problem and get this feedback "Warning: Infinite or Not-a-Number value encountered. "from MATLAB when I use integral() and fzero()...

About this problem, I have tried many advice from my past post on this website but I still could not solve it, hope you can give me some advice.

My code:

    tic
  U = 14; L = 7; d = (U-L)/2;
  T = ( 2*U+L )/(2+1);
  Du = U-T; Dl = T - L; d_star = min(Du,Dl);
  du = d/Du; dl = d/Dl;
  alpha = 0.10;
  aLe = 0.05;
  N = 2000; % simmulation times
  k = 1;
  for  ksi_given =  [-2,-1,-0.5 ,0, 0.5,1,2]
      if  ksi_given <= 0
          sigma = fzero(@(x) aLe - (x^2)*(  (dl^2*ksi_given^2+1)/d_star^2  ), 0.4 ) ;
          mu = ksi_given*sigma + T;
      else
          sigma = fzero(@(x) aLe - (x^2)*(  (du^2*ksi_given^2+1)/d_star^2  ), 0.4 ) ;
          mu = ksi_given*sigma + T;
      end
      every_sigma(k) = sigma;
      every_mu(k) = mu;
      j = 1;
      for n = [25,50,100,150,200]
          i = 1;
          for m = 1:1:N
              rdata = normrnd(mu,sigma,1,n);
              xbar = mean(rdata);
              sd = std(rdata);
              aLehat_rdata = ( max( (xbar-T)*d/Du,(T-xbar)*d/Dl )/d_star )^2 + var(rdata)/d_star^2;
              kursi_max = 0.5; % method1
              kursi_hat = (xbar-T)/sd; %method2
              fun_max = @(y,C) chi2cdf( (n*( (du.*kursi_max).^2+1 )./C).*aLehat_rdata-y, n-1).* ( 1./ (2.*sqrt(y)) ).* (   (dl^-1).*normpdf( (dl^-1).*sqrt(y) + sqrt(n)*(kursi_max) ) + (du^-1).*normpdf( (du^-1).*sqrt(y) - sqrt(n)*(kursi_max) )    );
              exact_UCB_max(i) = fzero(@(C) alpha - integral(@(y) fun_max(y,C),0.00001,(n.*( (du.*kursi_max).^2+1 )./C).*aLehat_rdata), 0.06 );
              if kursi_hat <= 0
                  fun2 = @(y,C) chi2cdf( (n*( (dl.*kursi_hat).^2+1 )./C).*aLehat_rdata-y, n-1).* ( 1./ (2.*sqrt(y)) ).* (   (dl^-1).*normpdf( (dl^-1).*sqrt(y) + sqrt(n)*(kursi_hat) ) + (du^-1).*normpdf( (du^-1).*sqrt(y) - sqrt(n)*(kursi_hat) )    );
                  %                 plotfun = @(C) alpha - integral(@(y) fun2(y,C),0.00001,(n.*( (dl.*kursi_hat).^2+1 )./C).*aLehat_rdata);
                  %                 c = linspace(1e-2,0.2,100);
                  %                 for i=1:numel(c)
                  %                     fc(i) = plotfun(c(i));
                  %                 end
                  %                 plot(c,fc)
                  exact_UCB_hat(i) = fzero(@(C) alpha - integral(@(y) fun2(y,C),0.00001,(n.*( (dl.*kursi_hat).^2+1 )./C).*aLehat_rdata), 0.06 ); 
              else
                  fun2 = @(y,C) chi2cdf( (n*( (du.*kursi_hat).^2+1 )./C).*aLehat_rdata-y, n-1).* ( 1./ (2.*sqrt(y)) ).* (   (dl^-1).*normpdf( (dl^-1).*sqrt(y) + sqrt(n)*(kursi_hat) ) + (du^-1).*normpdf( (du^-1).*sqrt(y) - sqrt(n)*(kursi_hat) )    );
                  exact_UCB_hat(i) = fzero(@(C) alpha - integral(@(y) fun2(y,C),0.00001,(n.*( (du.*kursi_hat).^2+1 )./C).*aLehat_rdata), 0.06 );
              end
              each_aLehat(i) = aLehat_rdata;
              i = i + 1;
          end
          CR_1(j,k) = mean(aLe<exact_UCB_max); %method1
          CR_2(j,k) = mean(aLe<exact_UCB_hat);  %method2
          average_aLehat(j,k) = mean(each_aLehat);
          sd_aLehat(j,k) = std(each_aLehat);
          average_exact_UCB_max(j,k) = mean(exact_UCB_max);
          sd_exact_UCB_max(j,k) = std(exact_UCB_max);
          average_exact_UB_hat(j,k) = mean(exact_UCB_hat);
          sd_exact_UB_hat(j,k) = std(exact_UCB_hat);
          j = j + 1;
      end
      k = k + 1;
  end
  toc

Error message:

Warning: Infinite or Not-a-Number value encountered. 

> In integralCalc/iterateScalarValued (line 349)

In integralCalc/vadapt (line 132)

In integralCalc (line 75)

In integral (line 88)

In @(C)alpha-integral(@(y)fun_max(y,C),0.00001,(n.*((du.*kursi_max).^2+1)./C).*aLehat_rdata)

In fzero (line 522)

About this problem, I am very sure that the problem will be caused by the two lines code below,

exact_UCB_max(i) = fzero(@(C) alpha - integral(@(y) fun_max(y,C),0.00001,(n.*( (du.*kursi_max).^2+1 )./C).*aLehat_rdata), 0.06 );

exact_UCB_max(i) = fzero(@(C) alpha - integral(@(y) fun_max(y,C),0.00001,(n.*( (du.*kursi_max).^2+1 )./C).*aLehat_rdata), 0.06 );

because I have encountered the similar problem before, so I try to modify the lower bound of integral from 0 to 0.00001 and use the plot technique like below to decide the starting value of fzero().

plotfun = @(C) alpha - integral(@(y) fun2(y,C),0.00001,(n.*( (dl.*kursi_hat).^2+1 )./C).*aLehat_rdata);
c = linspace(1e-2,0.2,100);
for i=1:numel(c)
    fc(i) = plotfun(c(i));
end
plot(c,fc)

Finally, I use 0.06 as the starting value for fzero() but it still return the same error message ...

Does anyone have any idea of this problem to let me fix it?

Very thank you!

[Walter Roberson]

I am still chasing this down, but I did manage to establish that one of the places that fun_max is evaluated at, the value would be about - 1.1804876878895667635963349937353e359 + 4.9524983241287364847887259189447e357i which overflows to -inf + inf*i if converted to double.

[Chao-Zhen Liu]

Hi Walter,

Can you tell me how to evaluate the fun_max value like what you do? And according to your comment, it means that the value will contain imaginary part?

Thanks!

[Torsten]

Don't use a starting value for fzero, but a starting interval, e.g. [0.06,0.1].

Best wishes

Torsten.

[Chao-Zhen Liu]

Hi Torsten,

I use the [0.06 0.1], [0.05 0.1],...until [1e-4 0.1] as my starting interval, but MATLAB still return the error message below,

Error using fzero (line 290)

The function values at the interval endpoints must differ in sign.

Did you have any other idea?

Thanks!

[Chao-Zhen Liu]

Hi Walter,

Can you tell me how to evaluate the fun_max value like what you do? And according to your comment, it means that the value will contain imaginary part?

Thanks!

[Walter Roberson]

Can you tell me how to evaluate the fun_max value like what you do?

It involves my adding symbolic replacements for the functions such as chi2cdf that you use, so that I can call the functions with symbolic numbers and get a more precise result than the nan or inf that you are getting returned. Unfortunately there are parts of the existing numeric statistical functions that do not accept symbolic inputs, because they try to do things like

a(a < 0) = nan;

to try to deal with values out of range. This is a problem if the inputs involve symbolic variables, because symbolic_variable relation constant is not decideable without a particular value for the variable -- so I have to rewrite in terms of detecting symbolic variables and using piecewise() to do the tests.

... it is a bit slow going to make the changes.

[Chao-Zhen Liu]

Hi Walter,

Thank for your help and explanation!

But it seems like my integration expression could not work by numerical integration and it is also have to do many pretreatment if we want to use symbolic integration, right?