How small is Zero?

I run some statistical analyses which result in p-values=0. Are they < the minimum positive value of a double precision (<2e-308), which is odd to me?

I am worried if another precision is used under the hood in the Statistics and Machine Learning Toolbox (e.g., corr).

Answer 1

[John D'Errico]

It is easy for that to happen. Of course, you don't actually tell us WHAT test you did. But this just means you have an event out in the tails. And it need not be even that far out.

For example, what is the probability of a standard Normally distributed number being as far out as -40?

format long g
Z = (-40:-35)';
[Z,normcdf(Z)]

ans = 6×2

-40 0 -39 0 -38 2.88542835100396e-316 -37 5.72557122252514e-300 -36 4.18262406579739e-284 -35 1.12491070647255e-268

So below 38 sigma, the probability just underflows. Actually, in context, -39*sigma really is a long way out. But outliers exist. It may mean your assumptions of normality may be in question.

But it is not at all impossible for a statistical test to yield an underflow as you have reported. You just need some data that justifies a result out in the tails.

[Ahmed Ramadan]

Great example and answer, John. My tests include corr and fitlme.

To follow up, do you know why the cdf values go smaller than the min real of a double's positive value (2.22507e-308)? There are e-316 in your example. So, what's the actual minimum (how small is zero)?

If you like to see some data, Figure 5 (last row) of this open-access paper is my analysis. I thought the min was machine epsilon (see the 2e-16 values in the figure) but then figured out it's even smaller

https://www.frontiersin.org/articles/10.3389/fpain.2023.1072786/full

[Les Beckham]

format long e
eps(0)

ans = 
    4.940656458412465e-324

[Walter Roberson]

smallest_normal_number = realmin

smallest_normal_number = 2.2251e-308

smallest_representable_number = eps(smallest_normal_number)

smallest_representable_number = 4.9000e-324

IEEE 754 Double Precision has two ranges.

In the range that is used nearly all of the time, the representation is sign * 2^(exponent+ bias) * (2^53 + mantissa) . So for example, 1.0 exactly is stored as [0, -53+bias, 0] meaning (-1)^0*2^(-53) * (2^53 + 0) -- a mantissa of 0 and an exponent that cancels out the 2^53 giving 1. And 2.0 exactly is stored as [0, -52+bias, 0] meaning (-1)^0*2^(-52)*(2^53 + 0) canceling out to 2^(53-52) giving 2^1 == 2.0 . In this "normalized" range, in each case, the mantissa is an integer at most 2^53-1 and the representation is always such that doubling the number just increases the exponent by 1 without changing the mantissa.

The second range applies only to numbers that are less than 2^(-1022) . For such numbers, the representation is instead a fixed 2^(-1074) * mantissa where mantissa is at most 2^53-1 . In that range, doubling the number leaves the exponent fixed and requires doubling the mantissa itself. The increment between adjcent representable numbers in this "denormalized" range is fixed at 2^(-1074), whereas the increment numbers in the "normalized" range varies with floor(log2(abs()) of the number.

[Ahmed Ramadan]

Thank you, Les and Walter!

I do appreciate it.

Answer 2

[Walter Roberson]

format long g
A = sym(floor(randn(5,2) * 16)) / 16

A = 

cord = corr(double(A))

cord = 2×2

1 -0.864423766518445 -0.864423766518445 1

cors = corr(A)

cors = 

corsv = vpa(cors, 16)

corsv = 

corsd = double(cors)

corsd = 2×2

1 -0.864423766518445 -0.864423766518445 1

That is, the corr() function happens to be able to run on symbolic numbers, and will provide exact results over a rather wide range -- 10^-10000 not being a problem for example. So you could test your "exact 0"