Could not integral: Infinite or Not-a-Number value encountered

Question

Chao-Zhen Liu 2018-9-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/420487-could-not-integral-infinite-or-not-a-number-value-encountered

评论： Walter Roberson 2018-10-2

采纳的回答： Walter Roberson

在 MATLAB Online 中打开

Hi everyone,

Does anyone can tell me what's wrong with my code? I always receives the warnings:

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

U = 14; L = 7; d = (U-L)/2;
d_star = d;
T = ( U+L )/(1+1);  % symmeric
alpha = 0.10;
Le = 0.05;
for kursi = [0 0.25 0.5 1 2 3]
  sigma = fzero( @(sigma) Le - (sigma./d)^2 - ( kursi.*sigma./d).^2, 1 ) ;
  mu = kursi*sigma + T;
    j = 1;
    for n = [25 50 100 150 200]
        delta = ( n.^(1/2) ).*kursi ; 
        B = (n*d_star^2)/sigma^2;
        i= 1;
        for x = 7:0.01:14
            sample = normrnd(mu,sigma,1,n);
%             fK = (  2^(-(n-1)/2)/gamma((n-1)/2)  ).*((B.*x.*(1-t)).^(n-3)/2  ).*exp(-B.*x.*(1-t)/2);     
            fun = @(t) (  sqrt(  (B^3).*x./t  )./2  ).*(   (  2^(-(n-1)/2)  / gamma((n-1)/2)  ).*(  (B.*x.*(1-t)).^((n-3)/2)  ).*exp(-B.*x.*(1-t)/2)   ).*(   normpdf(sqrt(B*x*T)+delta,0,1) + normpdf(sqrt(B*x*T)-delta,0,1)   );
            pdf_Lehat(j,i) = integral(@(t) fun(t),0,1);
            i = i + 1;
        end       
        j = j + 1;
    end 
end
x = 7:0.01:14;
plot(x, pdf_Lehat(1,:)); hold on
plot(x, pdf_Lehat(2,:)); hold on
plot(x, pdf_Lehat(3,:)); hold on
plot(x, pdf_Lehat(4,:)); hold on
plot(x, pdf_Lehat(5,:)); hold on
xlabel('X')

I guess the problem may be the handle ,fun, especially the mid part of the code (i.e. the above code, fK). Hope you can give me some advice, thanks!

9 个评论
显示 7更早的评论隐藏 7更早的评论

Chao-Zhen Liu 2018-9-27

编辑：Chao-Zhen Liu 2018-9-27

在 MATLAB Online 中打开

Hi Torsten,

It's the whole complete code, hope you can give me some advice,

Thanks!

 U = 14; L = 7; d = (U-L)/2;
d_star = d;
T = ( U+L )/(1+1);  % symmeric
alpha = 0.10;
Le = 0.05;
for kursi = [0 0.25 0.5 1 2 3]
  sigma = fzero( @(sigma) Le - (sigma./d)^2 - ( kursi.*sigma./d).^2, 1 ) ;
  mu = kursi*sigma + T;
    j = 1;
    for n = [25 50 100 150 200]
        delta = ( n.^(1/2) ).*kursi ; 
        B = (n*d_star^2)/sigma^2;
        i= 1;
        for x = 7:0.01:14
            sample = normrnd(mu,sigma,1,n);
%             fK = (  2^(-(n-1)/2)/gamma((n-1)/2)  ).*((B.*x.*(1-t)).^(n-3)/2  ).*exp(-B.*x.*(1-t)/2);     
            fun = @(t) (  sqrt(  (B^3).*x./t  )./2  ).*(   (  2^(-(n-1)/2)  / gamma((n-1)/2)  ).*(  (B.*x.*(1-t)).^((n-3)/2)  ).*exp(-B.*x.*(1-t)/2)   ).*(   normpdf(sqrt(B*x*T)+delta,0,1) + normpdf(sqrt(B*x*T)-delta,0,1)   );
            i = i + 1;
        end       
        j = j + 1;
    end 
end
plotfun = @(t) (  sqrt(  (B^3).*x./t  )./2  ).*(   (  2^(-(n-1)/2)  / gamma((n-1)/2)  ).*(  (B.*x.*(1-t)).^((n-3)/2)  ).*exp(-B.*x.*(1-t)/2)   ).*(   normpdf(sqrt(B*x*T)+delta,0,1) + normpdf(sqrt(B*x*T)-delta,0,1)   );
c = linspace(1,1.001,10);
    for i=1:numel(c)
        fc(i) = plotfun(c(i));
    end
    plot(c,fc)

Torsten 2018-9-27

编辑：Torsten 2018-9-27

In the evaluation of plotfun, you use B=4.0e4, x=14, n=200 and T=10.5.

Now specify a value for t and evaluate all parts of "plotfun" separately for these parameter values for B,x,n and T. See where there might be problems in the evaluation (e.g. gamma((n-1)/2)= gamma(199/2) seems too huge, 2^(-(n-1)/2)=2^(-199/2) seems too small).

Best wishes

Torsten.

Chao-Zhen Liu 2018-9-29

编辑：Chao-Zhen Liu 2018-9-30

在 MATLAB Online 中打开

Hi Torsten,

Thanks!

I have try to use log() or gammaln() and then use exp() to make some adjustment for the part that are too small or too big but it still could not work... could you have any idea?

If you want to see my code, please tell me and I will attach it right way when I am back to my computer.

Thanks!

========================================================================

    U = 14; L = 7; d = (U-L)/2;
d_star = d;
T = ( U+L )/(1+1);  % symmeric
alpha = 0.10;
Le = 0.05;
for kursi = [0 0.25 0.5 1 2 3]
  sigma = fzero( @(sigma) Le - (sigma./d)^2 - ( kursi.*sigma./d).^2, 1 ) ;
  mu = kursi*sigma + T;
    j = 1;
    for n = [25 50 100 150 200]
        delta = ( n.^(1/2) ).*kursi ; 
        B = (n*d_star^2)/sigma^2;
        i= 1;
        for x = 7:0.01:14
            sample = normrnd(mu,sigma,1,n);
%             fK = (  exp(log(2^(-(n-1)/2)))/exp(gammaln((n-1)/2))  ).*((B.*x.*(1-t)).^(n-3)/2  ).*exp(-B.*x.*(1-t)/2); 
            fun = @(t) (  sqrt(  exp(log(B^3)).*x./t  )./2  ).*(  exp(log(2^(-(n-1)/2)))/exp(gammaln((n-1)/2))  ).*exp(log(((B.*x.*(1-t)).^(n-3)/2  ))).*exp(-B.*x.*(1-t)/2).*(   normpdf(sqrt(exp(log(B*x*T)))+delta,0,1) + normpdf(exp(log(sqrt(B*x*T)))-delta,0,1)   );
            pdf_Lehat(j,i) = integral(@(t) fun(t),0.001,1);
            i = i + 1;
        end       
        j = j + 1;
    end 
end

Above is my code, what's wrong with my code? Thanks!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2018-9-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/420487-could-not-integral-infinite-or-not-a-number-value-encountered#answer_339080

The values of your integral are so small that they cannot be represented in double precision, and can barely be represented in the Symbolic Toolbox either. Values like 2*10^(-87012)

12 个评论
显示 10更早的评论隐藏 10更早的评论

Walter Roberson 2018-10-2

https://en.wikipedia.org/wiki/Error_function

I am not sure exactly how the erf or erfi arises in the integration, but apparently it does.

> simplify(int(A*exp(17500*(t-1))*(c*(t-1))^52*(1/t)^(1/2), t = 0 .. 1));

(139238746316160036479077520768376966254588477973825110187492738499994582454165856795118814677464953434308587460556690007028481114413130852159054721140000420116240988347245595841918273275908261304560369620917187872032236889/48661106313192264386250963245828545408130594296380877494812011718750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)*A*7^(1/2)*Pi^(1/2)*erfi(50*7^(1/2))*exp(-17500)*c^52-(3978363568774181572545017756089532100541695675193619366051081131334796786307621726077051976481550595871688523386229024885870479122006858993684260857373900206027743249831191634128209034915690641375057323875500203531889/69515866161703234837501376065469350583043706137686967849731445312500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)*A*c^52

The corresponding output for ^51 is

-(3978135901541453747049118581136237546045251906077855081775711463189655152723121851998258542138524704941831534393219115835604350312928461399865092952957204190293256172887192524997904004108271420294656493457339267251863/1390317323234064696750027521309387011660874122753739356994628906250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)*A*exp(-17500)*Pi^(1/2)*7^(1/2)*erfi(50*7^(1/2))*c^51+(113664273494599364480416077804955902043809572097757321783043375422587033965465269170557543293328209869842890195548358149037271199623233169014039664868418090109386228996164916838652382493969546574570189065140516863/1986167604620092423928610744727695730944105889648199081420898437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)*A*c^51

The coefficients almost but not quite balance out between the two sides. To 100 decimal places, the one for ^51 is -0.104845e-98*A*c^51 and the one for ^52 is 0.104847e-98*A*c^52 . The ratio between the two is just a hair different than -c.

However, testing with numeric integration near ^157 gives a notably smaller value than would be implied by extrapolation of these results; perhaps the estimate breaks down at higher values, or perhaps I did not use sufficient intermediate digits.

Chao-Zhen Liu 2018-10-2

Hi Walter,

About what you question me, I do not understand because I do not know why I have a 3D array of values near -81000... but I will recheck my equations as the top priority thing!

Thanks!

Walter Roberson 2018-10-2

Your term exp(-B.*x.*(1-t)/2) is responsible. The -B*x/2 is coming out at about 35000 and the 1-t flips that to about exp(-35000 *t)

请先登录，再进行评论。

Could not integral: Infinite or Not-a-Number value encountered

9 个评论
显示 7更早的评论隐藏 7更早的评论

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Could not integral: Infinite or Not-a-Number value encountered

9 个评论 显示 7更早的评论隐藏 7更早的评论

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

9 个评论
显示 7更早的评论隐藏 7更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论