How to solve the common problem about the using fzero() to solve a integration expression?
Show older comments
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!
7 Comments
Walter Roberson
on 12 Oct 2018
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
on 12 Oct 2018
Edited: Chao-Zhen Liu
on 12 Oct 2018
Torsten
on 12 Oct 2018
Don't use a starting value for fzero, but a starting interval, e.g. [0.06,0.1].
Best wishes
Torsten.
Chao-Zhen Liu
on 12 Oct 2018
Chao-Zhen Liu
on 12 Oct 2018
Walter Roberson
on 12 Oct 2018
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
on 13 Oct 2018
Answers (0)
Categories
Find more on Numerical Integration and Differentiation in Help Center and File Exchange
on 11 Oct 2018
on 13 Oct 2018
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
