@@ -329,12 +329,8 @@ def nancorr(const float64_t[:, :] mat, bint cov=False, minp=None):
329
329
Py_ssize_t i, j, xi, yi, N, K
330
330
bint minpv
331
331
float64_t[:, ::1 ] result
332
- # Initialize to None since we only use in the no missing value case
333
- float64_t[::1 ] means= None , ssqds= None
334
332
ndarray[uint8_t, ndim= 2 ] mask
335
- bint no_nans
336
333
int64_t nobs = 0
337
- float64_t mean, ssqd, val
338
334
float64_t vx, vy, dx, dy, meanx, meany, divisor, ssqdmx, ssqdmy, covxy
339
335
340
336
N, K = (< object > mat).shape
@@ -346,57 +342,25 @@ def nancorr(const float64_t[:, :] mat, bint cov=False, minp=None):
346
342
347
343
result = np.empty((K, K), dtype = np.float64)
348
344
mask = np.isfinite(mat).view(np.uint8)
349
- no_nans = mask.all()
350
-
351
- # Computing the online means and variances is expensive - so if possible we can
352
- # precompute these and avoid repeating the computations each time we handle
353
- # an (xi, yi) pair
354
- if no_nans:
355
- means = np.empty(K, dtype = np.float64)
356
- ssqds = np.empty(K, dtype = np.float64)
357
-
358
- with nogil:
359
- for j in range (K):
360
- ssqd = mean = 0
361
- for i in range (N):
362
- val = mat[i, j]
363
- dx = val - mean
364
- mean += 1 / (i + 1 ) * dx
365
- ssqd += (val - mean) * dx
366
-
367
- means[j] = mean
368
- ssqds[j] = ssqd
369
345
370
346
with nogil:
371
347
for xi in range (K):
372
348
for yi in range (xi + 1 ):
373
- covxy = 0
374
- if no_nans:
375
- for i in range (N):
349
+ # Welford's method for the variance-calculation
350
+ # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
351
+ nobs = ssqdmx = ssqdmy = covxy = meanx = meany = 0
352
+ for i in range (N):
353
+ if mask[i, xi] and mask[i, yi]:
376
354
vx = mat[i, xi]
377
355
vy = mat[i, yi]
378
- covxy += (vx - means[xi]) * (vy - means[yi])
379
-
380
- ssqdmx = ssqds[xi]
381
- ssqdmy = ssqds[yi]
382
- nobs = N
383
-
384
- else :
385
- nobs = ssqdmx = ssqdmy = covxy = meanx = meany = 0
386
- for i in range (N):
387
- # Welford's method for the variance-calculation
388
- # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
389
- if mask[i, xi] and mask[i, yi]:
390
- vx = mat[i, xi]
391
- vy = mat[i, yi]
392
- nobs += 1
393
- dx = vx - meanx
394
- dy = vy - meany
395
- meanx += 1 / nobs * dx
396
- meany += 1 / nobs * dy
397
- ssqdmx += (vx - meanx) * dx
398
- ssqdmy += (vy - meany) * dy
399
- covxy += (vx - meanx) * dy
356
+ nobs += 1
357
+ dx = vx - meanx
358
+ dy = vy - meany
359
+ meanx += 1 / nobs * dx
360
+ meany += 1 / nobs * dy
361
+ ssqdmx += (vx - meanx) * dx
362
+ ssqdmy += (vy - meany) * dy
363
+ covxy += (vx - meanx) * dy
400
364
401
365
if nobs < minpv:
402
366
result[xi, yi] = result[yi, xi] = NaN
0 commit comments