# --- Setup: load our helpers and the data --------------------------------
import sys, pathlib
for _p in [pathlib.Path.cwd(), pathlib.Path.cwd().parent]:
    if (_p / 'lib' / 'statslab.py').exists():
        sys.path.insert(0, str(_p / 'lib')); break
import statslab as sl

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sl.use_course_style()

AMES_URL = 'https://www.openintro.org/data/csv/ames.csv'   # same file + url as the core lesson
ames = sl.load_csv('ames.csv', url=AMES_URL)
print(f'Loaded a {type(ames).__name__} with shape {ames.shape}  (rows, columns)')

Loaded a DataFrame with shape (2930, 82)  (rows, columns)

print('shape  :', ames.shape)
print('the two columns we will study:', ['area', 'price'])
print('\ndtypes of those two columns:')
print(ames[['area', 'price']].dtypes)
ames[['area', 'price']].head()

shape  : (2930, 82)
the two columns we will study: ['area', 'price']

dtypes of those two columns:
area     int64
price    int64
dtype: object

print('rows missing area or price:', ames[['area', 'price']].isna().any(axis=1).sum())

ames = ames.dropna(subset=['area', 'price'])
print('rows kept after dropna     :', len(ames))

# Pull out the two columns we'll regress. Each is a Series of 2,930 numbers.
x = ames['area']      # predictor: above-grade living area (sq ft)
y = ames['price']     # response : sale price (dollars)
print('\nx is a', type(x).__name__, 'of', len(x), 'values; first 5 areas :', x.head().to_list())
print('y is a', type(y).__name__, 'of', len(y), 'values; first 5 prices:', y.head().to_list())

# A boolean mask in action: just the big homes.
big_mask = x > 3000
print(f'\nHomes over 3,000 sq ft: {big_mask.sum()} of {len(ames)}')

rows missing area or price: 0
rows kept after dropna     : 2930

x is a Series of 2930 values; first 5 areas : [1656, 896, 1329, 2110, 1629]
y is a Series of 2930 values; first 5 prices: [215000, 105000, 172000, 244000, 189900]

Homes over 3,000 sq ft: 26 of 2930

dx = x - x.mean()      # every home's area minus the mean area, all at once
dy = y - y.mean()      # every home's price minus the mean price

r_by_hand = (dx * dy).sum() / np.sqrt((dx**2).sum() * (dy**2).sum())
r_pandas  = x.corr(y)

print(f'r by hand     : {r_by_hand:.6f}')
print(f'x.corr(y)     : {r_pandas:.6f}')
print('match?', np.isclose(r_by_hand, r_pandas))

r = r_by_hand
print(f'\n-> r = {r:.2f}: a strong, positive linear trend — bigger homes really do go with bigger prices.')

r by hand     : 0.706780
x.corr(y)     : 0.706780
match? True

-> r = 0.71: a strong, positive linear trend — bigger homes really do go with bigger prices.

n = len(x)
sx, sy = x.std(ddof=1), y.std(ddof=1)         # sample SDs (divide by n-1)
cov_xy = (dx * dy).sum() / (n - 1)            # sample covariance, same ddof=1
var_x  = (dx**2).sum()  / (n - 1)             # sample variance of x,   same ddof=1

slope_a = r * (sy / sx)                       # way 1: r * (s_y / s_x)
slope_b = cov_xy / var_x                      # way 2: cov(x, y) / var(x)

print(f'slope via r*(sy/sx)     : {slope_a:.6f}')
print(f'slope via cov(x,y)/var(x): {slope_b:.6f}')
print('match?', np.isclose(slope_a, slope_b))

b1 = slope_a
print(f'\n-> slope b1 = {b1:.2f}: each extra square foot is associated with about ${b1:.0f} more in price.')

slope via r*(sy/sx)     : 111.694001
slope via cov(x,y)/var(x): 111.694001
match? True

-> slope b1 = 111.69: each extra square foot is associated with about $112 more in price.

a = y.mean() - b1 * x.mean()
print(f'intercept a = {a:,.2f} dollars  (the predicted price of a 0-sq-ft "home" — not physically meaningful)')

yhat = a + b1 * x            # 2,930 predicted prices, computed in one line
print(f'\nyhat is a {type(yhat).__name__} of {len(yhat)} predictions; first 5:')
print(np.round(yhat.head().to_numpy(), 0))

intercept a = 13,289.63 dollars  (the predicted price of a 0-sq-ft "home" — not physically meaningful)

yhat is a Series of 2930 predictions; first 5:
[198255. 113367. 161731. 248964. 195239.]

e = y - yhat                 # residuals for all 2,930 homes, one line
print(f'first 5 residuals (dollars): {np.round(e.head().to_numpy(), 0)}')
print(f'mean residual (should be ~0): {e.mean():.4f} dollars')

first 5 residuals (dollars): [16745. -8367. 10269. -4964. -5339.]
mean residual (should be ~0): -0.0000 dollars

ss_res = (e**2).sum()                 # leftover variation around the line
ss_tot = ((y - y.mean())**2).sum()    # total variation around the mean
R2_sums = 1 - ss_res / ss_tot         # way 1: from the sums of squares
R2_rsq  = r**2                        # way 2: just r, squared

print(f'R^2 from 1 - SS_res/SS_tot : {R2_sums:.6f}')
print(f'R^2 from r**2              : {R2_rsq:.6f}')
print('match?', np.isclose(R2_sums, R2_rsq))
print(f'\n-> About {R2_sums*100:.0f}% of the variation in sale price is explained by living area alone.')

R^2 from 1 - SS_res/SS_tot : 0.499538
R^2 from r**2              : 0.499538
match? True

-> About 50% of the variation in sale price is explained by living area alone.

xn = x.to_numpy()
yn = y.to_numpy()
b1_np, a_np = np.polyfit(xn, yn, 1)      # returns [slope, intercept]

print(f'our slope     : {b1:.6f}      np.polyfit slope    : {b1_np:.6f}')
print(f'our intercept : {a:,.4f}   np.polyfit intercept: {a_np:,.4f}')
print('slope match?    ', np.isclose(b1, b1_np))
print('intercept match?', np.isclose(a, a_np))

our slope     : 111.694001      np.polyfit slope    : 111.694001
our intercept : 13,289.6344   np.polyfit intercept: 13,289.6344
slope match?     True
intercept match? True

fig, (axL, axR) = plt.subplots(1, 2, figsize=(13, 5.0))

# LEFT: scatter + our fitted line.
line_x = np.array([xn.min(), xn.max()])
line_y = a + b1 * line_x                  # our own intercept & slope
axL.scatter(xn, yn, s=10, alpha=0.35, color='#4c72b0', label='homes')
axL.plot(line_x, line_y, color='#D55E00', lw=2.5,
         label=f'our fit: price = {a:,.0f} + {b1:.0f} x area')
axL.set_xlabel('living area (sq ft)'); axL.set_ylabel('sale price (dollars)')
axL.set_title('Our least-squares line through the Ames homes')
axL.legend(loc='upper left')

# RIGHT: residuals vs x, with a zero line.
axR.scatter(xn, e.to_numpy(), s=10, alpha=0.35, color='#4c72b0')
axR.axhline(0, color='#D55E00', lw=2, ls='--', label='residual = 0 (a perfect prediction)')
axR.set_xlabel('living area (sq ft)'); axR.set_ylabel('residual = actual − predicted (dollars)')
axR.set_title('Residual plot: a flat band that fans out for big homes')
axR.legend(loc='upper left')

plt.tight_layout(); plt.show()

# The loop way (how you might first think of it):
resid_loop = []
for xi, yi in zip(xn, yn):
    resid_loop.append(yi - (a + b1 * xi))
resid_loop = np.array(resid_loop)

# The vectorized way (what NumPy/pandas is for) — the whole column at once:
resid_vec = yn - (a + b1 * xn)

print('same residuals either way?', np.allclose(resid_loop, resid_vec))
print(f'first 3, loop      : {np.round(resid_loop[:3], 1)}')
print(f'first 3, vectorized: {np.round(resid_vec[:3], 1)}')
print('\nThe vectorized line reads cleaner and runs far faster on big data — and it is what the')
print('library does under the hood. (`np.polyfit`/`stats.linregress` solve the same least-squares')
print('problem for you; here we did it by hand to see exactly what they return.)')

same residuals either way? True
first 3, loop      : [16745.1 -8367.5 10269. ]
first 3, vectorized: [16745.1 -8367.5 10269. ]

The vectorized line reads cleaner and runs far faster on big data — and it is what the
library does under the hood. (`np.polyfit`/`stats.linregress` solve the same least-squares
problem for you; here we did it by hand to see exactly what they return.)

def fit_line(x, y):
    '''Least-squares fit of y on x. Returns dict(slope, intercept, r, r2).'''
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    dx = x - x.mean()
    dy = y - y.mean()
    r = (dx * dy).sum() / np.sqrt((dx**2).sum() * (dy**2).sum())
    slope = (dx * dy).sum() / (dx**2).sum()      # = cov(x,y)/var(x); the ddof cancels
    intercept = y.mean() - slope * x.mean()
    return {'slope': slope, 'intercept': intercept, 'r': r, 'r2': r**2}

# Reuse it on our two columns — and confirm it reproduces the by-hand numbers:
result = fit_line(ames['area'], ames['price'])
print('fit_line(area, price):')
for k, v in result.items():
    print(f'  {k:9}: {v:,.6f}')

print('\nmatches our step-by-step work?',
      np.isclose(result['slope'], b1),
      np.isclose(result['intercept'], a),
      np.isclose(result['r'], r),
      np.isclose(result['r2'], R2_sums))
print('\nCaution: a strong fit (high r, high R^2) shows association, NOT causation —')
print('bigger homes are *associated with* higher prices; the line alone never proves size *causes* price.')

fit_line(area, price):
  slope    : 111.694001
  intercept: 13,289.634365
  r        : 0.706780
  r2       : 0.499538

matches our step-by-step work? True True True True

Caution: a strong fit (high r, high R^2) shows association, NOT causation —
bigger homes are *associated with* higher prices; the line alone never proves size *causes* price.

Tool	What it is	What we use it for
pandas	tables (`DataFrame`) and columns (`Series`)	loading the CSV, selecting columns, masking rows
NumPy	fast math on whole arrays at once	the actual formulas (sum, sqrt, the fitted line)
matplotlib	plotting	the scatter-with-fit and the residual plot

The code behind it: fit a regression line by hand¶

1) The goal¶

2) The data in code¶

Selecting columns, masking rows, and missing values¶

3) Build it step by step¶

Pearson's correlation r, from its formula¶

The slope — two formulas that must agree¶

The intercept, and the fitted line¶

Residuals — what the line misses¶

R² — two ways that must agree¶

The final check: against `np.polyfit`¶

Draw it — the fit, and the residual plot¶

4) The idiom — vectorized vs. a loop¶

5) Refactor into a reusable function¶

6) Now you code¶

What you learned¶

	area	price
0	1656	215000
1	896	105000
2	1329	172000
3	2110	244000
4	1629	189900

The code behind it: fit a regression line by hand¶

1) The goal¶

2) The data in code¶

Selecting columns, masking rows, and missing values¶

3) Build it step by step¶

Pearson's correlation r, from its formula¶

The slope — two formulas that must agree¶

The intercept, and the fitted line¶

Residuals — what the line misses¶

R² — two ways that must agree¶

The final check: against np.polyfit¶

Draw it — the fit, and the residual plot¶

4) The idiom — vectorized vs. a loop¶

5) Refactor into a reusable function¶

6) Now you code¶

What you learned¶

The final check: against `np.polyfit`¶