function fzero with warning Imaginary parts of complex X and/or Y arguments ignored
4 ビュー (過去 30 日間)
古いコメントを表示
Hi everyone,
I am now encountering a problem that is a little complicated, hope you can read the detail below patiently and give me some advice, thanks.
My problem:
I am do some statistical problem and I write some code like this,
clc; clear
tic
U = 14; L = 7; d = (U-L)/2;
T = ( 3*U+L )/(3+1);
Du = U-T; Dl = T - L; d_star = min(Du,Dl);
du = d/Du; dl = d/Dl;
alpha = 0.05;
aLe = 0.05;
n = 100;
i = 1;
for ksi_given = 0:0.2:3
if ksi_given <= 0
sigma = fzero(@(x) aLe - (x^2)*( (dl^2*ksi_given^2+1)/d_star^2 ), [0 0.5] ) ;
mu = ksi_given*sigma + T;
B = n*( (dl*ksi_given)^2 + 1 )/aLe;
C(i,1) = (dl*ksi_given)^2/(d_star/sigma)^2 + 1/(d_star/sigma)^2;
else
sigma = fzero(@(x) aLe - (x^2)*( (du^2*ksi_given^2+1)/d_star^2 ), [0 0.5] ) ;
mu = ksi_given*sigma + T;
B = n*( (du*ksi_given)^2 + 1 )/aLe;
C(i,1) = (du*ksi_given)^2/(d_star/sigma)^2 + 1/(d_star/sigma)^2;
end
fun = @(ca,y) chi2cdf( B.*ca-y, n-1).* ( 1./ (2.*sqrt(y)) ).* ( (dl^-1).*normpdf( (dl^-1).*sqrt(y) + sqrt(n*ksi_given) ) + (du^-1).*normpdf( (du^-1).*sqrt(y) - sqrt(n*ksi_given) ) );
fun2 = @(ca) alpha - integral(@(y) fun(ca,y),0.0001,B*ca);
ca(i,1) = fzero(@(ca) alpha - integral(@(y) fun(ca,y),0.0001,B*ca), aLe);
i = i + 1;
end
toc
and my solving result comparing the standard answer is below:
My question are:
1. Why does I get the acceptable result only when ksi_given equal to 0 and 1 even though under the same code?
2. Actually, for this for loop, ksi_given must start from -3 to 3, but when I try it, MATLAB returns the error message:
Error using fzero (line 328) Function value at starting guess must be finite and real.
I am curious, so I try another code to check the case that ksi_given is negative:
clc; clear
tic
U = 14; L = 7; d = (U-L)/2;
T = ( 3*U+L )/(3+1);
Du = U-T; Dl = T - L; d_star = min(Du,Dl);
du = d/Du; dl = d/Dl;
alpha = 0.05;
aLe = 0.05;
n = 100;
ksi_given = -3;
if ksi_given <= 0
sigma = fzero(@(x) aLe - (x^2)*( (dl^2*ksi_given^2+1)/d_star^2 ), [0 0.5] ) ;
mu = ksi_given*sigma + T;
B = n*( (dl*ksi_given)^2 + 1 )/aLe;
C = (dl*ksi_given)^2/(d_star/sigma)^2 + 1/(d_star/sigma)^2;
b = d_star/sigma;
else
sigma = fzero(@(x) aLe - (x^2)*( (du^2*ksi_given^2+1)/d_star^2 ), [0 0.5] ) ;
mu = ksi_given*sigma + T;
B = ( n*( (du*ksi_given)^2 + 1 ) )/aLe;
C = (du*ksi_given)^2/(d_star/sigma)^2 + 1/(d_star/sigma)^2;
b = d_star/sigma;
end
fun = @(ca,y) chi2cdf( B.*ca-y, n-1).* ( 1./ (2.*sqrt(y)) ).* ( (dl^-1).*normpdf( (dl^-1).*sqrt(y) + sqrt(n*ksi_given) ) + (du^-1).*normpdf( (du^-1).*sqrt(y) - sqrt(n*ksi_given) ) );
fun2 =@(ca) alpha - integral(@(y) fun(ca,y),0.0001,B*ca);
c = linspace(1e-4,0.1,100);
for i=1:numel(c)
fc(i) = fun2(c(i));
end
plot(c,fc)
ca = fzero(@(ca) alpha - integral(@(y) fun(ca,y),0,B*ca), aLe);
toc
when ksi_given = -0.1, the plot and error message are
Warning: Imaginary parts of complex X and/or Y arguments ignored Error using fzero (line 328) Function value at starting guess must be finite and real.
when ksi_given = -3, the plot and error message are
Warning: Imaginary parts of complex X and/or Y arguments ignored Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.8e+57. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy. > In integralCalc/iterateScalarValued (line 372) In integralCalc/vadapt (line 132) In integralCalc (line 75) In integral (line 88) In @(ca)alpha-integral(@(y)fun(ca,y),0,B*ca) In fzero (line 303) Error using fzero (line 328) Function value at starting guess must be finite and real.
For my problem, I hope someone can give me some advice about the question below:
1. When ksi_given is positive, why I can only get the correct result when ksi_given is 0 and 1?
2. When ksi_given is negative, why does MATLAB couldn't solve the equation? Because of starting value(I have try many values...)? or some mathematical mistake about the integration ?
Thanks!
2 件のコメント
回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!