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 ptdf = 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.xarray([-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()| 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.410931.132899999999996