-
Notifications
You must be signed in to change notification settings - Fork 8
ALM
The Augmented Linear Model (alm in R, ALM in Python) is the core estimator in greybox. It extends the standard linear model by supporting 26 distributions and 7 loss functions, enabling regression for a wide variety of data types — continuous, count, binary, bounded, and zero-inflated.
In R, alm() follows R's formula-based interface: alm(y ~ x1 + x2, data, distribution="dnorm").
In Python, ALM follows the scikit-learn estimator pattern: construct with ALM(distribution="dnorm"), then call model.fit(X, y). The formula module converts R-style formulas to design matrices.
library(greybox)
# Fit model
model <- alm(mpg ~ ., mtcars, distribution="dlnorm", subset=c(1:30))
summary(model)
# Predict
predict(model, newdata=tail(mtcars, 2), interval="prediction", level=c(0.9,0.95))import numpy as np
from greybox import ALM, formula, mtcars
# Parse formula and fit model
y_vec, X_mat = formula("mpg ~ .", mtcars.head(30))
model = ALM(distribution="dlnorm")
model.fit(X_mat, y_vec)
print(model.summary())
# Predict
_, X_new = formula("mpg ~ .", mtcars.tail(2))
model.predict(X_new, interval="prediction", level=[0.9,0.95])| Parameter | Type | Default | Description |
|---|---|---|---|
formula |
formula | — | R formula, e.g. y ~ x1 + x2
|
data |
data.frame | — | Data containing the variables |
subset |
vector | NULL |
Subset of observations to use |
na.action |
function | — | How to handle NAs |
distribution |
character | "dnorm" |
Distribution (see table below) |
loss |
character | "likelihood" |
Loss function (see table below) |
occurrence |
character | "none" |
Occurrence model: "none", "plogis", "pnorm"
|
scale |
formula | NULL |
Scale model formula |
orders |
vector | c(0,0,0) |
ARIMA orders (p,d,q) |
parameters |
vector | NULL |
Initial parameter values |
fast |
logical | FALSE |
Use fast estimation (OLS only) |
... |
Additional: alpha, size, nu, shape, lambda, lambdaBC, FI
|
| Parameter | Type | Default | Description |
|---|---|---|---|
distribution |
str | "dnorm" |
Distribution (see table below) |
loss |
str | "likelihood" |
Loss function (see table below) |
occurrence |
str | "none" |
Occurrence model: "none", "plogis", "pnorm"
|
scale_formula |
array-like | None |
Scale model formula |
orders |
tuple | (0, 0, 0) |
ARIMA orders (p, d, q). AR order p, differencing d, MA q (MA not implemented; must be 0). Example: (1,1,0) for ARI(1,1). |
alpha |
float | None |
Asymmetric Laplace quantile level (0–1) |
shape |
float | None |
Generalized Normal shape parameter |
lambda_bc |
float | None |
Box-Cox lambda (0–1) |
size |
float | None |
Negative Binomial / Binomial size |
nu |
float | None |
Student's t / Chi-squared degrees of freedom |
trim |
float | 0.0 |
Trim proportion for ROLE loss |
lambda_l1 |
float | None |
L1 regularization (LASSO) |
lambda_l2 |
float | None |
L2 regularization (RIDGE) |
nlopt_kargs |
dict | {} |
Optimizer settings (see Optimizer section) |
verbose |
int | 0 |
Verbosity level |
All 26 distributions supported by ALM:
| Code | Name | Link | Extra Param | Use Case |
|---|---|---|---|---|
dnorm |
Normal | identity | — | General continuous data |
dlaplace |
Laplace | identity | — | Heavy-tailed, robust to outliers |
ds |
S (half-Laplace) | identity | — | Light-tailed data |
dgnorm |
Generalized Normal | identity | shape |
Flexible tail weight |
dlogis |
Logistic | identity | — | Heavy tails, symmetric |
dt |
Student's t | identity | nu |
Heavy tails, small samples |
dalaplace |
Asymmetric Laplace | identity | alpha |
Quantile regression |
| Code | Name | Link | Extra Param | Use Case |
|---|---|---|---|---|
dlnorm |
Log-Normal | identity (log scale) | — | Positive, right-skewed data |
dllaplace |
Log-Laplace | identity (log scale) | — | Positive, heavy-tailed |
dls |
Log-S | identity (log scale) | — | Positive, light-tailed |
dlgnorm |
Log-Generalized Normal | identity (log scale) | shape |
Positive, flexible tails |
| Code | Name | Link | Extra Param | Use Case |
|---|---|---|---|---|
dfnorm |
Folded Normal | identity | — | Absolute/non-negative data |
drectnorm |
Rectified Normal | identity | — | Zero-inflated non-negative |
dbcnorm |
Box-Cox Normal | identity | lambda_bc |
Non-normal, power transform |
dlogitnorm |
Logit-Normal | identity (logit scale) | — | Data in (0, 1) |
dbeta |
Beta | identity | — | Data in (0, 1), two-part model |
| Code | Name | Link | Extra Param | Use Case |
|---|---|---|---|---|
dinvgauss |
Inverse Gaussian | log | — | Positive, right-skewed |
dgamma |
Gamma | log | — | Positive, right-skewed |
dexp |
Exponential | log | — | Time between events |
| Code | Name | Link | Extra Param | Use Case |
|---|---|---|---|---|
dpois |
Poisson | log | — | Count data (mean ≈ variance) |
dnbinom |
Negative Binomial | log | size |
Overdispersed count data |
dbinom |
Binomial | log | — | Binary count data |
dgeom |
Geometric | log | — | Trials until first success |
dchisq |
Chi-squared | squared | nu |
Sum of squared normals |
| Code | Name | Link | Extra Param | Use Case |
|---|---|---|---|---|
plogis |
Logistic CDF | logistic | — | Binary classification |
pnorm |
Normal CDF (probit) | probit | — | Binary classification |
| Loss | Description |
|---|---|
likelihood |
Maximum likelihood estimation (default) |
MSE |
Mean Squared Error |
MAE |
Mean Absolute Error |
HAM |
Half Absolute Moment |
LASSO |
Likelihood + L1 penalty (lambda_l1) |
RIDGE |
Likelihood + L2 penalty (lambda_l2) |
ROLE |
Robust Likelihood Estimation (trimmed, trim) |
| Attribute | Description |
|---|---|
coefficients |
Estimated coefficients |
fitted.values / fitted
|
Fitted values |
residuals |
Model residuals |
mu |
Location parameter values |
scale |
Scale parameter |
distribution |
Distribution used |
logLik |
Log-likelihood |
loss |
Loss function type |
lossValue |
Final loss function value |
df.residual |
Residual degrees of freedom |
df |
Model degrees of freedom |
other |
Extra parameter value (alpha, shape, etc.) |
timeElapsed |
Computation time (seconds) |
FI |
Fisher Information matrix |
| Property | Type | Description |
|---|---|---|
coef |
ndarray | Slope coefficients (excluding intercept) |
coefficients |
ndarray | Same as coef
|
intercept_ |
float | Intercept value |
fitted |
ndarray | Fitted values |
residuals |
ndarray | Model residuals |
actuals |
ndarray | Original response values |
scale |
float | Scale parameter |
sigma |
float | Residual standard error |
loglik / log_lik
|
float | Log-likelihood |
aic |
float | Akaike Information Criterion |
aicc |
float | Corrected AIC |
bic |
float | Bayesian Information Criterion |
bicc |
float | Corrected BIC |
loss_value |
float | Final loss function value |
time_elapsed |
float | Computation time (seconds) |
df_residual_ |
int | Residual degrees of freedom |
distribution_ |
str | Distribution name |
loss_ |
str | Loss function name |
nobs |
int | Number of observations |
nparam |
int | Number of parameters |
n_param |
dict | {"number": nparam, "df": df_residual} |
formula |
str | Formula string (if provided) |
other_ |
dict | Extra parameter values |
data |
ndarray | Alias for actuals
|
arima_polynomial_ |
dict or None | AR+I combined polynomial coefficients keyed by lag name, e.g. {"yLag1": -0.8}. None if no ARIMA terms. |
arima_string_ |
str or None | Human-readable ARIMA order string, e.g. "ARIMA(1,0,0)". None if no ARIMA terms. |
| R Generic | Python Method | Description |
|---|---|---|
summary(model) |
model.summary(level=0.95) |
Summary with coefficients, std errors, p-values, significance |
predict(model, newdata, interval, level, side) |
model.predict(X, interval, level, side) |
Predictions with optional intervals |
confint(model, parm, level) |
model.confint(parm, level) |
Confidence intervals for parameters |
vcov(model) |
model.vcov() |
Variance-covariance matrix |
coef(model) |
model.coef (property) |
Extract coefficients |
sigma(model) |
model.sigma (property) |
Residual standard error |
logLik(model) |
model.loglik (property) |
Log-likelihood |
nobs(model) |
model.nobs (property) |
Number of observations |
nparam(model) |
model.nparam (property) |
Number of parameters |
| — | model.score(X, y, metric) |
Score using metric: "likelihood", "MSE", "MAE", "R2"
|
| — | model.get_params() |
Get constructor parameters (sklearn compat) |
| — | model.set_params(**params) |
Set constructor parameters (sklearn compat) |
The predict() method returns a PredictionResult object.
| Slot | Type | Description |
|---|---|---|
mean |
ndarray | Point forecasts |
lower |
ndarray or None | Lower interval bounds |
upper |
ndarray or None | Upper interval bounds |
level |
float or list | Confidence level(s) used |
variances |
ndarray or None | Variance estimates per observation |
side |
str |
"both", "upper", or "lower"
|
interval |
str |
"none", "confidence", or "prediction"
|
| Method / Property | Description |
|---|---|
to_dataframe() |
Convert to pandas DataFrame |
columns |
Column names of DataFrame representation |
shape |
Shape of DataFrame representation |
index |
Index of DataFrame representation |
values |
Values of DataFrame representation |
len(pred) |
Number of observations |
pred["mean"] |
Column access by name |
In R, predict() returns a list with $mean, $lower, $upper, $level, $variances.
# R
model <- alm(mpg ~ wt + hp, mtcars, distribution="dnorm")
summary(model)# Python
import numpy as np
from greybox import ALM, formula
data = {"y": np.array([21,22,23,20,18,15,14,24,23,20]),
"x1": np.array([2.6,2.8,2.3,3.2,3.4,3.8,4.0,2.2,2.5,3.0]),
"x2": np.array([110,110,93,150,150,245,245,62,95,123])}
y, X = formula("y ~ x1 + x2", data)
model = ALM(distribution="dnorm")
model.fit(X, y)
print(model.summary())# R — robust to outliers
model <- alm(y ~ x, data, distribution="dlaplace")# Python
model = ALM(distribution="dlaplace")
model.fit(X, y)# R — prediction interval
pred <- predict(model, newdata=new_data, interval="prediction", level=0.95)
pred$mean
pred$lower
pred$upper
# One-sided upper interval
pred <- predict(model, newdata=new_data, interval="prediction", level=0.95, side="upper")# Python — prediction interval
pred = model.predict(X_new, interval="prediction", level=0.95)
print(pred.to_dataframe())
# One-sided upper interval
pred = model.predict(X_new, interval="prediction", level=0.95, side="upper")
# Multiple confidence levels
pred = model.predict(X_new, interval="prediction", level=[0.80, 0.95])Note: Prediction and confidence interval bounds are computed using distribution-specific quantile functions to match R's
predict.alm(). For example,dlnormintervals are asymmetric and always positive;dpois/dnbinomintervals are integer-valued;plogis/pnormbounds lie in (0, 1). Symmetric t-distribution bounds are used only fordnormanddt.
# Python
summary = model.summary()
print(summary)
# Output:
# Response variable: y
# Distribution used in the estimation: Normal
# Loss function type: likelihood
# Coefficients:
# Estimate Std. Error Lower 2.5% Upper 97.5%
# (Intercept) 5.0123 1.2345 2.5926 7.4320 *
# x1 2.0045 0.1234 1.7626 2.2464 *
#
# Error standard deviation: 3.0123
# ...ALM supports autoregressive and integrated (differencing) components via the
orders=(p, d, q) parameter. Currently p (AR) and d (I) are supported;
q (MA) must be 0. Lagged response values are appended automatically to the
design matrix during fitting and prediction.
# R
ts_data <- c(2.1, 3.4, 3.0, 4.2, 5.1, 4.8, 5.5, 6.0, 5.8, 6.5,
7.1, 6.8, 7.5, 8.0, 7.8)
model_ar <- alm(y ~ 1, data.frame(y=ts_data), distribution="dnorm", orders=c(1,0,0))
summary(model_ar)
predict(model_ar, h=3, interval="prediction", level=0.9)# Python
import numpy as np
from greybox import ALM
from greybox.formula import formula
ts = np.array([2.1, 3.4, 3.0, 4.2, 5.1, 4.8, 5.5, 6.0, 5.8, 6.5,
7.1, 6.8, 7.5, 8.0, 7.8])
y, X = formula("y ~ 1", {"y": ts}) # intercept-only design matrix
model_ar = ALM(distribution="dnorm", orders=(1, 0, 0))
model_ar.fit(X, y)
print(model_ar.summary())
# Response variable: y
# Distribution used in the estimation: Normal
# Loss function type: likelihood
# ARIMA(1,0,0)
# Coefficients:
# Estimate Std. Error Lower 2.5% Upper 97.5%
# (Intercept) 1.xxxx x.xxxx x.xxxx x.xxxx *
# yLag1 -x.xxxx x.xxxx -x.xxxx -x.xxxx *
# ...
print(model_ar.arima_string_) # 'ARIMA(1,0,0)'
print(model_ar.arima_polynomial_) # {'yLag1': -0.xx}
# Forecast 3 steps ahead with 90% prediction intervals
X_new = np.ones((3, 1)) # intercept column only (lags appended automatically)
pred = model_ar.predict(X_new, interval="prediction", level=0.9)
print(pred.to_dataframe())
# mean lower upper
# 0 x.xx x.xx x.xx
# 1 x.xx x.xx x.xx
# 2 x.xx x.xx x.xx# Python
import numpy as np
from greybox import ALM
from greybox.formula import formula
np.random.seed(42)
n = 40
x = np.linspace(0, 4, n)
y_raw = 1.5 * x + np.cumsum(np.random.randn(n) * 0.5) + 10.0
y2, X2 = formula("y ~ x", {"y": y_raw, "x": x})
model_ari = ALM(distribution="dnorm", orders=(1, 1, 0))
model_ari.fit(X2, y2)
print(model_ari.summary())
print(model_ari.arima_string_) # 'ARIMA(1,1,0)'
print(model_ari.arima_polynomial_) # {'yLag1': ..., 'yLag2': ...}
# Out-of-sample forecast (supply new x values; y placeholder not used)
x_new = np.array([4.1, 4.2, 4.3])
_, X_new2 = formula("y ~ x", {"y": np.zeros(3), "x": x_new})
pred2 = model_ari.predict(X_new2, interval="prediction", level=0.95)
print(pred2.to_dataframe())ALM uses nlopt for optimization (nloptr in case of R). Settings can be customized via the nlopt_kargs parameter (Python) or ... arguments (R).
| Key | Type | Default | Description |
|---|---|---|---|
algorithm |
str | "NLOPT_LN_NELDERMEAD" |
Optimization algorithm |
maxeval |
int | 40 * n_params |
Maximum function evaluations |
maxtime |
float | 600 |
Maximum time in seconds |
xtol_rel |
float | 1e-6 |
Relative tolerance on parameters |
xtol_abs |
float | 1e-8 |
Absolute tolerance on parameters |
ftol_rel |
float | 1e-4 |
Relative tolerance on function value |
ftol_abs |
float | 0 |
Absolute tolerance on function value |
print_level |
int | 0 |
Verbosity (0=none, 3=full) |
| Algorithm | Description |
|---|---|
NLOPT_LN_NELDERMEAD |
Nelder-Mead simplex (default, derivative-free) |
NLOPT_LN_SBPLX |
Subplex (improved Nelder-Mead) |
NLOPT_LN_COBYLA |
COBYLA (constrained) |
NLOPT_LN_BOBYQA |
BOBYQA (quadratic approx) |
NLOPT_LN_NEWUOA |
NEWUOA (unconstrained quadratic) |
NLOPT_LN_PRAXIS |
PRAXIS (gradient-free) |
NLOPT_LD_LBFGS |
L-BFGS (gradient-based) |
NLOPT_LD_SLSQP |
SLSQP (gradient-based, constrained) |
NLOPT_LD_MMA |
MMA (gradient-based, constrained) |
# Example: custom optimizer settings
model = ALM(
distribution="dnorm",
nlopt_kargs={
"algorithm": "NLOPT_LN_SBPLX",
"maxeval": 1000,
"xtol_rel": 1e-8,
"print_level": 1,
}
)| Aspect | R | Python |
|---|---|---|
| Interface | alm(formula, data) |
ALM().fit(X, y) |
| Formula | Built-in | Via formula() → X, y |
| Coefficients |
coef(model) or model$coefficients
|
model.coef property |
| Intercept | Included in coef()
|
Separate model.intercept_
|
| Predict | predict(model, newdata=df) |
model.predict(X) |
| Summary | summary(model) |
model.summary() |
| Optimizer | nloptr |
nlopt |
fast param |
Available (TRUE for OLS) |
Not implemented |
| sklearn compat | — |
get_params(), set_params(), score()
|
- Svetunkov, I. (2023). Statistics for Business Analytics. https://openforecast.org/sba/
- Svetunkov, I. (2022). Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM). https://openforecast.org/adam/