Skip to content
Explore Automated Analysis

How to calculate a Pearson correlation and confidence interval in SciPy (with Code)

Editorial · Published

The Short Answer (TL;DR)

Call scipy.stats.pearsonr on matched numeric observations, then use its confidence_interval method for a Fisher-transformed interval. This example reports the default two-sided normal-model p-value and an approximate 95% interval, with explicit checks for finite, nonconstant input.

from scipy import stats
result = stats.pearsonr([1,2,3,4,5,6,7,8,9,10,11,12], [2,4,1,5,3,7,4,8,6,10,7,9], alternative="two-sided", method=None)
print(result.statistic, result.pvalue, result.confidence_interval(confidence_level=0.95, method=None))

Detailed Implementation

This example uses synthetic data and does not need a network connection. Install the recorded package versions in a separate Python environment before running it.

python -m pip install numpy==2.5.4 scipy==1.18.1
from pathlib import Path
import numpy as np
from scipy import stats

# One pair per independent unit; these fixed synthetic values need no download.
x = np.arange(1.0, 13.0)
y = np.array([2, 4, 1, 5, 3, 7, 4, 8, 6, 10, 7, 9], dtype=float)
if x.ndim != 1 or y.ndim != 1 or x.size != y.size or x.size <= 3:
    raise ValueError("Use matched one-dimensional arrays with more than three pairs.")
if not (np.isfinite(x).all() and np.isfinite(y).all()):
    raise ValueError("Resolve non-finite values jointly without breaking the pairs.")
if np.ptp(x) == 0 or np.ptp(y) == 0:
    raise ValueError("Correlation is undefined for a constant input.")
result = stats.pearsonr(x, y, alternative="two-sided", method=None)
interval = result.confidence_interval(confidence_level=0.95, method=None)
report = (
    f"Independent pairs: {x.size}\n"
    f"Pearson r: {result.statistic:.8f}\n"
    f"Default two-sided p-value: {result.pvalue:.8g}\n"
    f"Approximate 95% Fisher interval: [{interval.low:.8f}, {interval.high:.8f}]\n"
)
print(report)
Path("pearson-correlation-output.txt").write_text(report, encoding="utf-8")

The report contains twelve matched pairs, r approximately 0.82482190, a two-sided p-value approximately 0.00096039945, and a 95% Fisher interval from approximately 0.47646103 to 0.94933840. Both the short answer and full implementation use the same observations and inferential settings.

Which observations belong together?

Each x value and the y value at the same position describe the same independent unit. Equal array lengths do not establish correct matching. Resolve missing or non-finite observations as pairs, record exclusions, and inspect a scatterplot for outliers, nonlinear patterns and restricted ranges. Repeated measurements from one unit need an analysis that retains that dependence.

What assumptions support the reported inference?

The coefficient describes linear association. The default p-value tests zero population correlation using the normal-model null distribution; independent pairs from a bivariate normal population support this calculation. The Fisher interval is an approximation, especially with few pairs or correlations near the boundaries. Non-normal or influential data can make either inference misleading even when the code runs.

What is the interval estimating?

The interval concerns the population correlation, not the slope, a future response, or the fraction of points lying inside a band. A 95% confidence procedure describes repeated-sampling coverage under its assumptions; it is not a 95% probability that this fixed population parameter lies in these realized endpoints. Positive correlation does not establish a causal relationship.

Mathematical Context

Center both variables and divide their cross-product sum by the square root of the product of their centered sum of squares. For these twelve pairs the result is approximately 0.824822. For the Fisher interval transform r with arctanh, use the approximate standard error 1/sqrt(n-3), add and subtract the 0.975 standard-normal quantile times that standard error, and transform the endpoints back with tanh. The resulting interval is approximately [0.476461, 0.949338]. This requires n greater than three for the stated standard error. The checks separately reconstruct r, the transformed interval, and the default two-sided probability using both the beta null distribution and its t-test equivalence.

\[ z=\operatorname{arctanh}(r),\qquad \rho_{L,U}=\tanh\left(z\mp\frac{\Phi^{-1}(0.975)}{\sqrt{n-3}}\right) \]

Reproducibility

Executed with Python 3.12.14; numpy 2.5.4, scipy 1.18.1. Automated checks confirm the stated output properties.

Primary documentation

Related guides