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

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!

[Torsten]

In fun, you divide by t, but your integration starts at t=0 ...

[Chao-Zhen Liu]

Hi Torsten,

I have changed the lower bound of integration from 0 to 0.001, but it still return the same warnings... How can I do some remedy for it?

Thanks!

[Torsten]

I suggest you try to plot "fun" before integration to see what the problem is.

[Chao-Zhen Liu]

Hi Torsten,

I have try it like this:

    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(0.01,0.1,10);
    for q=1:numel(c)
        fc(q) = plotfun(c(q));
    end
    plot(c,fc)

but I receive the message;

Subscripted assignment dimension mismatch.

Why is the reason? Please tell me...

Thanks!

[Chao-Zhen Liu]

Hi Toresten,

I have try to plot the fun like this:

    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)

Does it means that we could not do integration by using fun because it's a y = 0 line?

Thanks!

[Torsten]

Please show your complete code to generate the graphics (i.e. how you defined the constants within "plotfun").

[Chao-Zhen Liu]

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]

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]

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]

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)

[Chao-Zhen Liu]

Hi Walter,

You means that I have to use Symbolic solution and can not use numerical integration to solve my question?

Thanks!

[Walter Roberson]

No, it means that you will need to use a true function instead of an anonymous function for your formula, and your true function will need to detect the 0 / 0 case and substituted 0 instead of NaN in order to proceed. You would then get a numeric result. Which will just happen to be completely 0 for all locations, because all of the numeric results are far far far smaller than can be represented by IEEE 754 binary double precision floating point.

[Walter Roberson]

Note that you have a for loop for kursi, and inside that loop you reset j and i to 1, so for each i and j combination, pdf_Lehat(j,i) is assigned to 6 times (one time for each different kursi value.) The result is going to reflect only the last kursi value, with all of the work for the previous kursi values thrown away.

[Walter Roberson]

I did some tests with int() and vpa() of the result. MATLAB gave back the same expression for int(), which means it does not know a closed-form expression for the integral. But it also gave back the int() in response to vpa() of the int(), which means that it could not find a convergent answer to within 32 decimal places (somehow I accidentally had 1000 decimal places active, so indeed it was not able to find a convergent solution within 1000 decimal places.)

There appears to be an issue right near 1. At 1, the function becomes 0, but at 1-10^-delta the slope is about 10^-(22*delta) -- e.g., at 1-1E-199, the value of the function is about 1E-19677 but at 1-1E-200 the value of the function is about 1E-19699 . Makes numeric integration difficult.

It isn't just MATLAB: Maple would want hours and hours to try to find a single closed form solution.

I would suggest you investigate the possibility of an error in your expressions.

[Chao-Zhen Liu]

Hi Walter,

Thank for your advice. I will double check the possible error in my expressions and modify the for loop to collect all values for each kursi.

Sincerely thank the help from you and Torsten!!

[Walter Roberson]

One of the functions to be integrated comes out (with particular random sample) as

 -(2*2^(1/2)*14^(1/2)*exp(17500*t - 17500)*(35000*t - 35000)^197*((2251799813685248*exp(-(350*3^(1/2) - 200^(1/2)/2)^2/2))/5644425081792261 + (2251799813685248*exp(-(350*3^(1/2) + 200^(1/2)/2)^2/2))/5644425081792261)*(1/t)^(1/2))/(2143938466757130852473958595130865258556619703211863714207255462832859419843826015972543007930252357592916363484036789873155201291342207545698304901117569170096411970882415771484375*Pi^(1/2))

This can be generalized to

A*exp(b*(t-1))*(c*(t-1))^197*(1/t)^(1/2)

with t not exceeding 1, the t-1 part is negative, so you have a component like exp(-17500) involved. This is overwhelmingly close to 0.

I see some hints that with sufficient effort it might be possible to find a closed form integral for this function, involving a bunch of exp() and erf() or erfi() and some increasingly large coefficients (probably combinatorially large). But the exp(-b) are really small....

[Chao-Zhen Liu]

Hi Walter,

Can you briefly tell me why you use erf() in here or what's the meaning of erf()?

Thanks!

[Walter Roberson]

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]

Hi Walter,

I think I have do additional study about erf()...

Because I still do not know why should I use erf() in the integration about my problem and which part of the integration should I add the erf()

My original idea is that use log() first and then use exp() to handle the notably bigger or smaller value for what I want to integrate and change the lower bound of t to avoid t = 0 for integration, but it seem to be not work...

Can you tell me more about the usage of erf() or erfi() for my problem?

Thanks!

[Walter Roberson]

What I said was that it looks like potentially there might be a closed form solution for the integral given in the code that could be expressed with erf() or erfi() . However, finding the coefficients for the closed form might be difficult.

But even if you had a closed form solution, the values would still be in the approximate range of 1E-81000, which is about the 250th power of the smallest representable double precision number.

Effectively you cannot compute this with double precision arithmetic, and need on the order of 100 digits of software floating point to get a meaningful answer (that will be in the range of 1E-81000 or so.)

You really need to recheck your equations before you get into potential mechanisms to accurately calculate values that small.

If you did have code that could return (say) log10 of the values, so you had a 3D array of values near -81000, then what would you do with those values?

[Chao-Zhen Liu]

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]

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)