Relationships

Pearson Correlation

Normalize covariance by the product of the two standard deviations: r = cov(x,y) / (std_x × std_y), giving a unitless value in [−1, 1]. Two loops: one combined loop for means, one joint loop accumulating cov_sum, sq_x, and sq_y simultaneously. r = cov_sum / sqrt(sq_x × sq_y) — the n−1 denominators cancel, so ddof choice is irrelevant as long as it is consistent. With scipy, pearsonr(x, y).statistic; p-value in snapshot only.

By hand

Loop 1: accumulate sx and sy together → mx=3.0, my=4.0. Loop 2: for each index compute dx=x[i]−mx and dy=y[i]−my, then accumulate cov_sum+=dx×dy, sq_x+=dx², sq_y+=dy². cov_sum=8, sq_x=10, sq_y=8. r = 8/sqrt(80) = 2/√5 ≈ 0.8944. The n−1 factors in cov and std² both cancel in the ratio, so the result is independent of ddof choice.

naive.py
Replay: real traced execution (multi-file project)
import math
x = [1, 2, 3, 4, 5]
y = [2, 4, 4, 4, 6]
n = len(x)
sx = 0.0
sy = 0.0
for i in range(n):
    sx = sx + x[i]
    sy = sy + y[i]
mx = sx / n
my = sy / n
cov_sum = 0.0
sq_x = 0.0
sq_y = 0.0
for i in range(n):
    cov_sum = cov_sum + (x[i] - mx) * (y[i] - my)
    sq_x = sq_x + (x[i] - mx) ** 2
    sq_y = sq_y + (y[i] - my) ** 2
r = cov_sum / math.sqrt(sq_x * sq_y)
print('RESULT:', round(r, 4))

  1. import math

    1import math2x = [1, 2, 3, 4, 5]

  2. x ← [1, 2, 3, 4, 5]

    1import math2x = [1, 2, 3, 4, 5]3y = [2, 4, 4, 4, 6]
    values this step[1, 2, 3, 4, 5]x

  3. y ← [2, 4, 4, 4, 6]

    2x = [1, 2, 3, 4, 5]3y = [2, 4, 4, 4, 6]4n = len(x)
    values this step[2, 4, 4, 4, 6]y

  4. n ← 5

    3y = [2, 4, 4, 4, 6]4n = len(x)5sx = 0.0
    values this step5n

  5. sx ← 0.0

    4n = len(x)5sx = 0.06sy = 0.0
    values this step0.0sx

  6. sy ← 0.0

    5sx = 0.06sy = 0.07for i in range(n):
    values this step0.0sy

  7. i ← 0, sx ← 1.0, sy ← 2.0

    pass 1 of 5
    6sy = 0.07for i in range(n):8    sx = sx + x[i]9    sy = sy + y[i]10mx = sx / n
    values this step0i0.0 → 1.0sx0.0 → 2.0sy
    All 5 passes — pass 1 is the card above
    passisxsy
    100.0 → 1.00.0 → 2.0
    20 → 11.0 → 3.02.0 → 6.0
    31 → 23.0 → 6.06.0 → 10.0
    42 → 36.0 → 10.010.0 → 14.0
    53 → 410.0 → 15.014.0 → 20.0

  8. for i in range(n):

    6sy = 0.07for i in range(n):8    sx = sx + x[i]

  9. mx ← 3.0

    9    sy = sy + y[i]10mx = sx / n11my = sy / n
    values this step3.0mx

  10. my ← 4.0

    10mx = sx / n11my = sy / n12cov_sum = 0.0
    values this step4.0my

  11. cov_sum ← 0.0

    11my = sy / n12cov_sum = 0.013sq_x = 0.0
    values this step0.0cov_sum

  12. sq_x ← 0.0

    12cov_sum = 0.013sq_x = 0.014sq_y = 0.0
    values this step0.0sq_x

  13. sq_y ← 0.0

    13sq_x = 0.014sq_y = 0.015for i in range(n):
    values this step0.0sq_y

  14. i ← 0, cov_sum ← 4.0, sq_x ← 4.0, sq_y ← 4.0

    pass 1 of 5
    14sq_y = 0.015for i in range(n):16    cov_sum = cov_sum + (x[i] - mx) * (y[i] - my)17    sq_x = sq_x + (x[i] - mx) ** 218    sq_y = sq_y + (y[i] - my) ** 219r = cov_sum / math.sqrt(sq_x * sq_y)
    values this step4 → 0i0.0 → 4.0cov_sum0.0 → 4.0sq_x0.0 → 4.0sq_y
    All 5 passes — pass 1 is the card above
    passicov_sumsq_xsq_y
    14 → 00.0 → 4.00.0 → 4.00.0 → 4.0
    20 → 1—4.0 → 5.0—
    31 → 2———
    42 → 3—5.0 → 6.0—
    53 → 44.0 → 8.06.0 → 10.04.0 → 8.0

  15. for i in range(n):

    14sq_y = 0.015for i in range(n):16    cov_sum = cov_sum + (x[i] - mx) * (y[i] - my)

  16. r ← 0.8944271909999159

    18    sq_y = sq_y + (y[i] - my) ** 219r = cov_sum / math.sqrt(sq_x * sq_y)20print('RESULT:', round(r, 4))
    values this step0.8944271909999159r

  17. stdout ← RESULT: 0.8944

    19r = cov_sum / math.sqrt(sq_x * sq_y)20print('RESULT:', round(r, 4))
    values this stepRESULT: 0.8944stdout

With the library

scipy.stats.pearsonr(x, y) returns a result object; .statistic is the Pearson r, .pvalue is the two-tailed p-value for H₀: r=0. RESULT compares the correlation coefficient only — the p-value is shown in the snapshot for context but is not the comparable output.

library.py 
from scipy.stats import pearsonr
from dalib.display import set_display
set_display()

x = [1, 2, 3, 4, 5]
y = [2, 4, 4, 4, 6]
result = pearsonr(x, y)
r = result.statistic
p = result.pvalue
print('r:', round(r, 4))
print('p-value:', round(p, 4))
print('RESULT:', round(r, 4))

r: 0.8944
p-value: 0.0405
RESULT: 0.8944

Honesty

This lesson shows the computation exactly, on a tiny pinned sample. The arithmetic is correct and reproducible, but with a sample this small the result is not a valid statistical finding — it demonstrates the mechanism, not evidence. Real inference needs an adequate sample size and assumption checks (e.g. independence, linearity, and bivariate normality); the p-value / interval here should be read as "how the formula is computed," not as a conclusion about a population.

Implementation notes

  • r is bounded in [−1, 1]: 1 is perfect positive linear relationship, −1 perfect negative, 0 no linear relationship.
  • The ddof cancellation: cov = cov_sum/(n−1), std_x = sqrt(sq_x/(n−1)), std_y = sqrt(sq_y/(n−1)); r = [cov_sum/(n−1)] / [sqrt(sq_x/(n−1)) × sqrt(sq_y/(n−1))] = cov_sum / sqrt(sq_x × sq_y). The (n−1) fully cancels.
  • Cross-reference: covariance-by-hand (this chapter) for the unnormalized version and the n−1 convention.