import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
from scipy.optimize import minimize
import pandas as pd
import statsmodels.api as sm
import patsy as pt
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/penguins.csv")
f = "body_mass_g ~ island + flipper_length_mm"
y, X = pt.dmatrices(f, data = df)
X[:, :]
array([[  1.,   0.,   1., 181.],
       [  1.,   0.,   1., 186.],
       [  1.,   0.,   1., 195.],
       ...,
       [  1.,   1.,   0., 193.],
       [  1.,   1.,   0., 210.],
       [  1.,   1.,   0., 198.]])
def ll_linear_regression(beta, data):
    y = data["y"]
    X = data["X"]
    yhat = np.sum(beta * X, axis = 1)
    return np.sum((y[:, 0] - yhat) ** 2)
data = {"y": y, "X": X}
rng = np.random.default_rng()
o = minimize(ll_linear_regression, rng.normal(size = 4), 
             args = (data,), 
             method = "L-BFGS-B")
o.x
array([-4625.50244192,  -260.00380406,  -187.43335478,    44.54360838])
fit = sm.OLS.from_formula("body_mass_g ~ island * flipper_length_mm", data = df).fit()
fit.summary()
OLS Regression Results
Dep. Variable: body_mass_g R-squared: 0.786
Model: OLS Adj. R-squared: 0.783
Method: Least Squares F-statistic: 246.5
Date: Thu, 23 Apr 2026 Prob (F-statistic): 4.55e-110
Time: 13:43:25 Log-Likelihood: -2508.3
No. Observations: 342 AIC: 5029.
Df Residuals: 336 BIC: 5052.
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -5463.9397 431.358 -12.667 0.000 -6312.442 -4615.438
island[T.Dream] 3550.7000 969.105 3.664 0.000 1644.422 5456.978
island[T.Torgersen] 3217.8277 1679.621 1.916 0.056 -86.071 6521.726
flipper_length_mm 48.5438 2.052 23.653 0.000 44.507 52.581
island[T.Dream]:flipper_length_mm -19.4038 4.938 -3.929 0.000 -29.117 -9.690
island[T.Torgersen]:flipper_length_mm -17.4109 8.730 -1.994 0.047 -34.584 -0.238
Omnibus: 6.268 Durbin-Watson: 2.557
Prob(Omnibus): 0.044 Jarque-Bera (JB): 6.128
Skew: 0.324 Prob(JB): 0.0467
Kurtosis: 3.102 Cond. No. 1.81e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.81e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

Learning about reductions in numpy.

A = rng.normal(size = (3, 4, 6, 16))
A.shape
(3, 4, 6, 16)
np.sum(A, axis = 0).shape
(4, 6, 16)

Towards Non-Linear Regression

Re slopes

ndf = pd.DataFrame({"flipper_length_mm": [200, 201], 
                    "island": ["Torgersen", "Torgersen"]})
np.diff(fit.predict(ndf))
array([31.13288034])
48.5438-17.4109
31.132899999999996