Compute all pairwise Pearson r values for three columns and assemble the 3×3 symmetric correlation matrix (diagonal = 1). Use sum() for column means, then one joint loop accumulating six sums (three cross-products and three squared-deviation sums). Three distinct off-diagonal r values reveal that each pair has a different relationship: strong positive (a,b), moderate negative (a,c), weak negative (b,c). With numpy, np.corrcoef([a, b, c]).

By hand

Data: a=[1,2,3,4], b=[1,2,4,3], c=[4,1,3,2]. All means are 2.5 (sum(col)/n). One joint loop for all six accumulators: three cross-products (ss_ab, ss_ac, ss_bc) and three squared-deviation sums (ss_a, ss_b, ss_c). Results: ss_ab=4, ss_ac=−2, ss_bc=−1, ss_a=ss_b=ss_c=5. Three distinct r values: r_ab=4/√25=0.8, r_ac=−2/√25=−0.4, r_bc=−1/√25=−0.2.

naive.py
Replay: real traced execution (multi-file project)
import math
a = [1, 2, 3, 4]
b = [1, 2, 4, 3]
c = [4, 1, 3, 2]
n = len(a)
ma = sum(a) / n
mb = sum(b) / n
mc = sum(c) / n
ss_ab = 0.0
ss_ac = 0.0
ss_bc = 0.0
ss_a = 0.0
ss_b = 0.0
ss_c = 0.0
for i in range(n):
    ss_ab = ss_ab + (a[i] - ma) * (b[i] - mb)
    ss_ac = ss_ac + (a[i] - ma) * (c[i] - mc)
    ss_bc = ss_bc + (b[i] - mb) * (c[i] - mc)
    ss_a = ss_a + (a[i] - ma) ** 2
    ss_b = ss_b + (b[i] - mb) ** 2
    ss_c = ss_c + (c[i] - mc) ** 2
r_ab = ss_ab / math.sqrt(ss_a * ss_b)
r_ac = ss_ac / math.sqrt(ss_a * ss_c)
r_bc = ss_bc / math.sqrt(ss_b * ss_c)
print('RESULT:', (round(r_ab, 4), round(r_ac, 4), round(r_bc, 4)))
  1. import math

    1import math2a = [1, 2, 3, 4]
  2. a ← [1, 2, 3, 4]

    1import math2a = [1, 2, 3, 4]3b = [1, 2, 4, 3]
    values this step[1, 2, 3, 4]a
  3. b ← [1, 2, 4, 3]

    2a = [1, 2, 3, 4]3b = [1, 2, 4, 3]4c = [4, 1, 3, 2]
    values this step[1, 2, 4, 3]b
  4. c ← [4, 1, 3, 2]

    3b = [1, 2, 4, 3]4c = [4, 1, 3, 2]5n = len(a)
    values this step[4, 1, 3, 2]c
  5. n ← 4

    4c = [4, 1, 3, 2]5n = len(a)6ma = sum(a) / n
    values this step4n
  6. ma ← 2.5

    5n = len(a)6ma = sum(a) / n7mb = sum(b) / n
    values this step2.5ma
  7. mb ← 2.5

    6ma = sum(a) / n7mb = sum(b) / n8mc = sum(c) / n
    values this step2.5mb
  8. mc ← 2.5

    7mb = sum(b) / n8mc = sum(c) / n9ss_ab = 0.0
    values this step2.5mc
  9. ss_ab ← 0.0

    8mc = sum(c) / n9ss_ab = 0.010ss_ac = 0.0
    values this step0.0ss_ab
  10. ss_ac ← 0.0

    9ss_ab = 0.010ss_ac = 0.011ss_bc = 0.0
    values this step0.0ss_ac
  11. ss_bc ← 0.0

    10ss_ac = 0.011ss_bc = 0.012ss_a = 0.0
    values this step0.0ss_bc
  12. ss_a ← 0.0

    11ss_bc = 0.012ss_a = 0.013ss_b = 0.0
    values this step0.0ss_a
  13. ss_b ← 0.0

    12ss_a = 0.013ss_b = 0.014ss_c = 0.0
    values this step0.0ss_b
  14. ss_c ← 0.0

    13ss_b = 0.014ss_c = 0.015for i in range(n):
    values this step0.0ss_c
  15. i ← 0, ss_ab ← 2.25, ss_ac ← -2.25, ss_bc ← -2.25, ss_a ← 2.25

    pass 1 of 4
    14ss_c = 0.015for i in range(n):16    ss_ab = ss_ab + (a[i] - ma) * (b[i] - mb)17    ss_ac = ss_ac + (a[i] - ma) * (c[i] - mc)18    ss_bc = ss_bc + (b[i] - mb) * (c[i] - mc)19    ss_a = ss_a + (a[i] - ma) ** 220    ss_b = ss_b + (b[i] - mb) ** 221    ss_c = ss_c + (c[i] - mc) ** 222r_ab = ss_ab / math.sqrt(ss_a * ss_b)
    values this step0i0.0 2.25ss_ab0.0 -2.25ss_ac0.0 -2.25ss_bc0.0 2.25ss_a0.0 2.25ss_b0.0 2.25ss_c
    All 4 passes — pass 1 is the card above
    passiss_abss_acss_bcss_ass_bss_c
    100.0 2.250.0 -2.250.0 -2.250.0 2.250.0 2.250.0 2.25
    20 12.25 2.5-2.25 -1.5-2.25 -1.52.25 2.52.25 2.52.25 4.5
    31 22.5 3.25-1.5 -1.25-1.5 -0.752.5 2.752.5 4.754.5 4.75
    42 33.25 4.0-1.25 -2.0-0.75 -1.02.75 5.04.75 5.04.75 5.0
  16. for i in range(n):

    14ss_c = 0.015for i in range(n):16    ss_ab = ss_ab + (a[i] - ma) * (b[i] - mb)
  17. r_ab ← 0.8

    21    ss_c = ss_c + (c[i] - mc) ** 222r_ab = ss_ab / math.sqrt(ss_a * ss_b)23r_ac = ss_ac / math.sqrt(ss_a * ss_c)
    values this step0.8r_ab
  18. r_ac ← -0.4

    22r_ab = ss_ab / math.sqrt(ss_a * ss_b)23r_ac = ss_ac / math.sqrt(ss_a * ss_c)24r_bc = ss_bc / math.sqrt(ss_b * ss_c)
    values this step-0.4r_ac
  19. r_bc ← -0.2

    23r_ac = ss_ac / math.sqrt(ss_a * ss_c)24r_bc = ss_bc / math.sqrt(ss_b * ss_c)25print('RESULT:', (round(r_ab, 4), round(r_ac, 4), round(r_bc, 4)))
    values this step-0.2r_bc
  20. stdout ← RESULT: (0.8, -0.4, -0.2)

    24r_bc = ss_bc / math.sqrt(ss_b * ss_c)25print('RESULT:', (round(r_ab, 4), round(r_ac, 4), round(r_bc, 4)))
    values this stepRESULT: (0.8, -0.4, -0.2)stdout

With the library

np.corrcoef([a, b, c]) returns the full 3×3 matrix in one call. Diagonals are always 1.0; the three unique off-diagonal values are m[0,1], m[0,2], m[1,2]. Note the three relationships are visibly distinct — not all equal.

library.py
import numpy as np
from dalib.display import set_display
set_display()

a = [1, 2, 3, 4]
b = [1, 2, 4, 3]
c = [4, 1, 3, 2]
m = np.corrcoef([a, b, c])
print('r_ab:', round(float(m[0, 1]), 4))
print('r_ac:', round(float(m[0, 2]), 4))
print('r_bc:', round(float(m[1, 2]), 4))
print('RESULT:', (round(float(m[0, 1]), 4), round(float(m[0, 2]), 4), round(float(m[1, 2]), 4)))
r_ab: 0.8
r_ac: -0.4
r_bc: -0.2
RESULT: (0.8, -0.4, -0.2)

Honesty

This lesson shows the computation exactly, on a tiny pinned sample. The pairwise correlations are correct and reproducible, but with this few rows they demonstrate the mechanism, not evidence — a correlation matrix on a handful of observations is not a valid statistical finding. Real inference needs an adequate sample size and assumption checks (independence, linearity, bivariate normality); read each entry as "how the coefficient is computed," not as a conclusion about the population.

Implementation notes

  • sum(col)/n treats sum() as a primitive (mean computation was shown in pearson-correlation; here the focus is the cross-product accumulation).
  • The six-accumulator joint loop computes all three pairwise r values in a single O(n) pass, avoiding three separate paired loops.
  • Interpreting the matrix: (a,b)=0.8 — b tracks a closely (they differ only in positions 2 and 3); (a,c)=−0.4 — c partially inverts a; (b,c)=−0.2 — weak negative, little linear structure.
  • For k columns the number of unique off-diagonal pairs is k(k−1)/2 (here k=3 gives 3 pairs). np.corrcoef handles arbitrary k in one call.
  • Cross-reference: pearson-correlation (this chapter) for the single-pair formula that each matrix cell applies.