1.2  Regression

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
N = 1200

# ── Simulate House Price Dataset ──────────────────────────────
area     = np.random.uniform(45, 200, N)
beds     = np.random.randint(1, 6, N).astype(float)
dist_km  = np.random.uniform(1, 25, N)
age_yr   = np.random.uniform(0, 50, N)
noise    = np.random.normal(0, 15, N)

price = (2.4*area + 18*beds - 4.2*dist_km - 0.8*age_yr + 80 + noise)  # $k
df = pd.DataFrame({'area':area,'beds':beds,'dist_km':dist_km,'age_yr':age_yr,'price':price})

X = df[['area','beds','dist_km','age_yr']].values
y = df['price'].values

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ── Train OLS ─────────────────────────────────────────────────
lr = LinearRegression()
lr.fit(X_tr, y_tr)
y_pred = lr.predict(X_te)

r2   = r2_score(y_te, y_pred)
rmse = np.sqrt(mean_squared_error(y_te, y_pred))
print("=== Linear Regression — House Price ===")
print(f"R²   = {r2:.4f}")
print(f"RMSE = {rmse:.2f} $k")
print(f"Intercept: {lr.intercept_:.3f}")
for name, coef in zip(['area','beds','dist_km','age_yr'], lr.coef_):
    print(f"  {name:10s}: β = {coef:+.4f}")

# ── 5-fold CV ─────────────────────────────────────────────────
cv_r2 = cross_val_score(lr, X, y, cv=5, scoring='r2')
print(f"\nCV R²: {cv_r2.mean():.4f} ± {cv_r2.std():.4f}")

# ── Residual diagnostics (print summary) ─────────────────────
residuals = y_te - y_pred
print(f"\nResidual stats — Mean:{residuals.mean():.3f}, Std:{residuals.std():.3f}")
print(f"  Skewness: {pd.Series(residuals).skew():.3f}")

# ── Agriculture: Crop Yield ───────────────────────────────────
np.random.seed(7)
n2 = 500
rain = np.random.uniform(400, 1000, n2)
fert = np.random.uniform(50, 250, n2)
temp = np.random.uniform(12, 28, n2)
yield_t = 0.35*(rain/100) + 0.012*fert - 0.08*temp + 2.5 + np.random.normal(0,.3,n2)
Xa = np.column_stack([rain, fert, temp])
ya = yield_t
lr_agr = LinearRegression().fit(Xa, ya)
print(f"\nCrop Yield R² = {r2_score(ya, lr_agr.predict(Xa)):.4f}")

import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.metrics import r2_score, mean_squared_error

np.random.seed(42)
N = 300

# ── Simulate Fertiliser Response (Quadratic) ──────────────────
nitrogen = np.random.uniform(0, 350, N)
yield_t  = 1.8 + 0.07*nitrogen - 0.0002*nitrogen**2 + np.random.normal(0, 0.25, N)

X = nitrogen.reshape(-1, 1)
y = yield_t

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ── Compare degrees via cross-validation ─────────────────────
print("Degree | CV R²     | CV RMSE")
print("-------|-----------|--------")
for deg in range(1, 6):
    pipe = Pipeline([
        ('poly', PolynomialFeatures(degree=deg, include_bias=False)),
        ('lr',   LinearRegression())
    ])
    cv_r2   = cross_val_score(pipe, X_tr, y_tr, cv=5, scoring='r2')
    cv_rmse = np.sqrt(-cross_val_score(pipe, X_tr, y_tr, cv=5,
                                        scoring='neg_mean_squared_error'))
    print(f"  d={deg}  | {cv_r2.mean():.4f}±{cv_r2.std():.3f} | {cv_rmse.mean():.4f}")

# ── Best model (degree 2) ─────────────────────────────────────
best_pipe = Pipeline([
    ('poly', PolynomialFeatures(degree=2, include_bias=True)),
    ('lr',   LinearRegression())
])
best_pipe.fit(X_tr, y_tr)
y_pred = best_pipe.predict(X_te)
print(f"\nTest R² = {r2_score(y_te, y_pred):.4f}")
print(f"RMSE    = {np.sqrt(mean_squared_error(y_te, y_pred)):.4f}")

lr_step = best_pipe.named_steps['lr']
print(f"\nCoefficients: intercept={lr_step.intercept_:.4f}, β₁={lr_step.coef_[1]:.6f}, β₂={lr_step.coef_[2]:.8f}")

# ── Predict optimal N ─────────────────────────────────────────
n_grid = np.linspace(0, 350, 1000).reshape(-1,1)
y_grid = best_pipe.predict(n_grid)
optimal_N = n_grid[np.argmax(y_grid)][0]
print(f"\nOptimal Nitrogen dose: {optimal_N:.1f} kg/ha → Yield: {y_grid.max():.2f} t/ha")

# ── Medicine: Pharmacokinetics ────────────────────────────────
np.random.seed(3)
time  = np.random.uniform(0.25, 24, 100)
conc  = 400 * time * np.exp(-0.35*time) + np.random.normal(0, 10, 100)
pipe_pk = Pipeline([('poly', PolynomialFeatures(3)), ('lr', LinearRegression())])
pipe_pk.fit(time.reshape(-1,1), conc)
print(f"\nPharmacokinetics Polynomial R² = {pipe_pk.score(time.reshape(-1,1), conc):.4f}")

import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge, RidgeCV, LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score, mean_squared_error

np.random.seed(42)
N, P = 300, 40

# ── Simulate Highly Correlated Financial Factors ──────────────
base      = np.random.randn(N, 5)
X_raw     = base @ np.random.randn(5, P) + 0.3*np.random.randn(N, P)
true_beta = np.random.randn(P) * np.random.choice([0,1], P, p=[0.6,0.4])
y         = X_raw @ true_beta + np.random.normal(0, 1.5, N)

X_tr, X_te, y_tr, y_te = train_test_split(X_raw, y, test_size=0.25, random_state=42)
scaler = StandardScaler()
X_tr_sc = scaler.fit_transform(X_tr); X_te_sc = scaler.transform(X_te)

# ── OLS vs Ridge ──────────────────────────────────────────────
ols = LinearRegression().fit(X_tr_sc, y_tr)
ols_r2 = r2_score(y_te, ols.predict(X_te_sc))

alphas = np.logspace(-3, 4, 100)
ridge_cv = RidgeCV(alphas=alphas, cv=5, scoring='r2')
ridge_cv.fit(X_tr_sc, y_tr)
best_alpha = ridge_cv.alpha_
ridge_r2 = r2_score(y_te, ridge_cv.predict(X_te_sc))

print("=== Ridge Regression — Portfolio Factors ===")
print(f"OLS  Test R²   = {ols_r2:.4f}")
print(f"Ridge Test R²  = {ridge_r2:.4f}  (best λ = {best_alpha:.4f})")

# ── Compare coefficient shrinkage ─────────────────────────────
print("\nCoefficient shrinkage (first 6 features):")
print(f"{'Feature':10} | {'OLS β':>10} | {'Ridge β':>10}")
for i in range(6):
    print(f"{'X'+str(i):10} | {ols.coef_[i]:10.4f} | {ridge_cv.coef_[i]:10.4f}")

# ── Ridge Path ────────────────────────────────────────────────
print("\nRidge path (λ → coef norm):")
for lam in [0.01, 1, 10, 100, 1000]:
    r = Ridge(alpha=lam).fit(X_tr_sc, y_tr)
    print(f"  λ={lam:6.2f}  |coef|₂ = {np.linalg.norm(r.coef_):.4f}")

import numpy as np
import pandas as pd
from sklearn.linear_model import Lasso, LassoCV
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score

np.random.seed(42)
N, P = 300, 80

# ── Simulate: 80 features, only 12 truly relevant ────────────
X = np.random.randn(N, P)
true_coef = np.zeros(P)
true_coef[:12] = np.random.randn(12) * 2  # only first 12 matter
y = X @ true_coef + np.random.normal(0, 1.5, N)

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
sc = StandardScaler()
X_tr_sc = sc.fit_transform(X_tr); X_te_sc = sc.transform(X_te)

# ── LassoCV finds optimal λ ───────────────────────────────────
lasso_cv = LassoCV(cv=5, alphas=np.logspace(-4, 2, 100), max_iter=10000, random_state=42)
lasso_cv.fit(X_tr_sc, y_tr)
best_alpha = lasso_cv.alpha_

print("=== Lasso Regression — Feature Selection ===")
print(f"Best λ (alpha) = {best_alpha:.6f}")
print(f"Test R²        = {r2_score(y_te, lasso_cv.predict(X_te_sc)):.4f}")

# ── Count selected features ───────────────────────────────────
n_nonzero = np.sum(lasso_cv.coef_ != 0)
print(f"Non-zero coefficients: {n_nonzero} / {P}")

# ── Lasso Regularisation Path ─────────────────────────────────
print("\nLasso Path:")
for alpha in [0.001, 0.01, 0.1, 1.0, 10.0]:
    m = Lasso(alpha=alpha, max_iter=10000).fit(X_tr_sc, y_tr)
    r2 = r2_score(y_te, m.predict(X_te_sc))
    nz = np.sum(m.coef_ != 0)
    print(f"  λ={alpha:6.3f}  features={nz:3d}  Test R²={r2:.4f}")

# ── Compare true vs recovered coefficients ────────────────────
print("\nRecovery of true coefficients (first 15):")
for i in range(15):
    print(f"  β{i:02d}: true={true_coef[i]:+.3f}  lasso={lasso_cv.coef_[i]:+.3f}")

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error

np.random.seed(42)
N = 1000

# ── Simulate Insurance Premium Dataset ───────────────────────
age    = np.random.randint(18, 70, N).astype(float)
bmi    = np.random.normal(27, 6, N)
smoker = np.random.binomial(1, 0.2, N).astype(float)
claims = np.random.poisson(0.5, N).astype(float)
region = np.random.randint(0, 4, N).astype(float)

premium = (250*age + 500*bmi + 12000*smoker + 2000*claims +
           300*region + 1000 + np.random.normal(0, 800, N))
premium = np.maximum(premium, 500)  # floor at 500

X = np.column_stack([age, bmi, smoker, claims, region])
feature_names = ['Age','BMI','Smoker','PriorClaims','Region']
y = premium

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ── Single Tree vs Random Forest ─────────────────────────────
dt = DecisionTreeRegressor(max_depth=6, random_state=42)
dt.fit(X_tr, y_tr)
dt_rmse = np.sqrt(mean_squared_error(y_te, dt.predict(X_te)))
dt_r2   = r2_score(y_te, dt.predict(X_te))

rf = RandomForestRegressor(n_estimators=200, max_features='sqrt',
                            oob_score=True, n_jobs=-1, random_state=42)
rf.fit(X_tr, y_tr)
rf_rmse = np.sqrt(mean_squared_error(y_te, rf.predict(X_te)))
rf_r2   = r2_score(y_te, rf.predict(X_te))

print("=== Random Forest vs Single Tree ===")
print(f"{'':18} | {'R²':>8} | {'RMSE':>10}")
print(f"{'Decision Tree':18} | {dt_r2:8.4f} | {dt_rmse:10.2f}")
print(f"{'Random Forest':18} | {rf_r2:8.4f} | {rf_rmse:10.2f}")
print(f"{'RF OOB R²':18} | {rf.oob_score_:8.4f} |")

# ── Feature Importances ───────────────────────────────────────
print("\nFeature Importances:")
for name, imp in sorted(zip(feature_names, rf.feature_importances_),
                         key=lambda x: -x[1]):
    bar = '█' * int(imp*40)
    print(f"  {name:15s}: {imp:.4f}  {bar}")

# ── Number of trees effect ────────────────────────────────────
print("\nEffect of number of trees on OOB R²:")
for n in [10, 50, 100, 200, 500]:
    rf_n = RandomForestRegressor(n_estimators=n, oob_score=True, n_jobs=-1, random_state=42)
    rf_n.fit(X_tr, y_tr)
    print(f"  n_trees={n:4d}  OOB R²={rf_n.oob_score_:.4f}")