Nery Neto
Nery Neto

Reputation: 59

Precision Matlab and Python (numpy)

I'm converting a Matlab script to Python and I am getting different results in the 10**-4 order.

In matlab:

f_mean=f_mean+nanmean(f); 
f = f - nanmean(f);

f_t  = gradient(f);
f_tt = gradient(f_t);

if n_loop==1
  theta = atan2( sum(f.*f_tt), sum(f.^2) );
end

theta = -2.2011167e+03

In Python:

f_mean = f_mean + np.nanmean(vel)

vel = vel - np.nanmean(vel)

firstDerivative = np.gradient(vel)
secondDerivative = np.gradient(firstDerivative)


if numberLoop == 1:
    theta = np.arctan2(np.sum(vel * secondDerivative),
                       np.sum([vel**2]))

Although first and secondDerivative give the same results in Python and Matlab, f_mean is slightly different: -0.0066412 (Matlab) and -0.0066414 (Python); and so theta: -0.4126186 (M) and -0.4124718 (P). It is a small difference, but in the end leads to different results in my scripts.

I know some people asked about this difference, but always regarding std, which I get, but not regarding mean values. I wonder why it is.

Upvotes: 0

Views: 852

Answers (1)

Paul Panzer
Paul Panzer

Reputation: 53029

One possible source of the initial difference you describe (between means) could be numpy's use of pairwise summation which on large arrays will typically be appreciably more accurate than the naive method:

a = np.random.uniform(-1, 1, (10**6,))
a = np.r_[-a, a]
# so the sum should be zero

a.sum()
# 7.815970093361102e-14

# use cumsum to get naive summation:
a.cumsum()[-1]
# -1.3716805469243809e-11

Edit (thanks @sascha): for the last word and as a "provably exact" reference you could use math.fsum:

import math
math.fsum(a)
# 0.0

Don't have matlab, so can't check what they are doing.

Upvotes: 1

Related Questions