Why, when I calculate the skewness of a series with constants, do I not get NaN?

Dear Matlab community,

I have a simulation of normally distributed numbers with a mean of 76 and a standard deviation of 0, which describes a series of constant numbers with the value 76. To this simulation, I apply a function that also yields constants and, therefore, its skewness should be also NaN. Why, when I calculate their skewness, do I not get NaN but 1? How do I correct this error?

Thanks in advance!

Np=20; 
cb=2; 
ab = 30*(Np/100)^(cb-1); 
n=100;
Qd = normrnd(76,0,n,1);
Pw=ab*((Qd/Np).^(cb-1)); 
skewness(Pw)
skewness(Qd)
var(Qd)

Answer 1

[Cris LaPierre]

I suspect this has to do with how floating point numbers are stored in a computer. 76 can be stored exactly as is, but 22.8000 is not exactly 22.8 in your computer. If you round Pw, you will see the skewness is also NaN.

format long
Np=20; 
cb=2; 
ab = 30*(Np/100)^(cb-1); 
Qd = normrnd(76,0,100,1);
Pw=ab*((Qd/Np).^(cb-1));
Pw(1)

ans = 
  22.799999999999997

skewness(Pw)

ans = 
    -1

var(Pw)

ans = 
     1.274926715508706e-29

skewness(round(Pw))

ans = 
   NaN

skewness(Qd)

ans = 
   NaN

var(Qd)

ans = 
     0

[Angelavtc]

@Cris LaPierre Thank you! It works! :)

[Angelavtc]

@Cris LaPierre However, this is creating another kind of problem for simulations with variance different than 0. For example when:

clear, clc;
Np=20; 
cb=1.1; 
ab = 30*(Np/100)^(cb-1); 
n=10000;
Qd = normrnd(108,3,n,1);
Pw=ab*((Qd/Np).^(cb-1)); 
Skew1=skewness(Pw)

Skew1 = -0.0787

Skew2=skewness(round(Pw))

Skew2 = 33.2883

Skew1 and Skew2 are completely different, being Skew1 the one close to 0 (my correct answer). How can I correct this?

Thank you!

[Cris LaPierre]

Consider taking a look at the algorithm the skewness function is using.