Likelihood, take 02

likelihood for linear regression

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")
y, X = pt.dmatrices("body_mass_g ~ flipper_length_mm", data = df)
X[:5, :]
array([[  1., 181.],
       [  1., 186.],
       [  1., 195.],
       [  1., 193.],
       [  1., 190.]])
rng = np.random.default_rng()
beta = rng.normal(size = 2)
np.sum(beta * X[0])
array(-315.02421429)
yhat = np.sum(beta * X, axis = 1)
np.sum((y[:, 0] - yhat) ** 2)
array(7.31618449e+09)
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}
o = minimize(ll_linear_regression, rng.normal(size = 2), 
             args = (data,), 
             method = "L-BFGS-B")
o.x
array([-5781.12418231,    49.68701617])

Check answers against StatsModels OLS.

fit = sm.OLS.from_formula("body_mass_g ~ flipper_length_mm", data = df).fit()
fit.summary()
OLS Regression Results
Dep. Variable: body_mass_g R-squared: 0.759
Model: OLS Adj. R-squared: 0.758
Method: Least Squares F-statistic: 1071.
Date: Tue, 21 Apr 2026 Prob (F-statistic): 4.37e-107
Time: 13:47:05 Log-Likelihood: -2528.4
No. Observations: 342 AIC: 5061.
Df Residuals: 340 BIC: 5069.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -5780.8314 305.815 -18.903 0.000 -6382.358 -5179.305
flipper_length_mm 49.6856 1.518 32.722 0.000 46.699 52.672
Omnibus: 5.634 Durbin-Watson: 2.190
Prob(Omnibus): 0.060 Jarque-Bera (JB): 5.585
Skew: 0.313 Prob(JB): 0.0613
Kurtosis: 3.019 Cond. No. 2.89e+03


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