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

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.

naive.py
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))
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.