Delivery Time Cost Prediction using Bayesian Regression in ML

Logistics companies need to forecast the total cost of each delivery—before dispatch—using early indicators such as delivery distance, package weight, number of stops, and time of day. Delivery costs exhibit nonlinear scale effects (e.g., per‑mile fuel discounts, driver‑hour premiums) and are subject to uncertainty from traffic conditions and variable labour rates. A single point‑estimate model ignores this uncertainty, risking budget overruns or under‑quoting. By applying Bayesian Regression, we obtain both:

1. A point estimate of expected delivery cost.

2. A credible interval that quantifies our uncertainty—enabling data‑driven pricing, resource allocation, and risk management.

Libraries Required

import pandas as pd                              # data loading & manipulation  
import numpy as np                               # numerical operations  

import matplotlib.pyplot as plt                  # plotting  
import seaborn as sns                            # visualization  

import pymc3 as pm                               # Bayesian modeling  
import arviz as az                               # posterior analysis  

from sklearn.model_selection import train_test_split  
from sklearn.preprocessing import StandardScaler  
from sklearn.metrics import mean_absolute_error  

Dataset

Cost Prediction for Logistic Company (ML)

Step-by-Step Code Implementation

Import Libraries & Load Data

import pandas as pd

# Load training data
df = pd.read_csv("data/train.csv")

# Preview key columns
df[['distance_km','weight_kg','num_stops','is_peak_hour','cost']].head()

Preprocessing & Train/Test Split

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# Select predictors and target
X = df[['distance_km','weight_kg','num_stops','is_peak_hour']].values
y = df['cost'].values  # cost in USD

# Chronological split if timestamp exists; else random split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Standardize numeric features for stable MCMC
scaler = StandardScaler().fit(X_train[:, :3])
X_train_s = X_train.copy()
X_train_s[:, :3] = scaler.transform(X_train[:, :3])
X_test_s  = X_test.copy()
X_test_s[:, :3]  = scaler.transform(X_test[:, :3])

Define & Fit the Bayesian Regression Model

Priors:

• α ∼ Normal(0, 100): broad intercept prior for baseline cost.
• β ∼ Normal(0, 50): moderate uncertainty on each predictor—distance, weight, stops, and peak‑hour flag.
• σ ∼ HalfNormal(50): positive residual noise scale.

Model Structure:

• Linear predictor μ = α + β·X_standardized captures how each feature drives cost.
• Observations cost ∼ Normal(μ, σ).

Inference:

• We run MCMC with 2,000 posterior draws (after 1,000 tuning) at target_accept=0.9 to ensure stable convergence.
• Posterior predictive sampling yields full predictive distributions for new inputs.

import pymc3 as pm

with pm.Model() as model:
    # Priors
    α = pm.Normal("α", mu=0, sigma=100)
    β = pm.Normal("β", mu=0, sigma=50, shape=X_train_s.shape[1])
    σ = pm.HalfNormal("σ", sigma=50)

    # Linear predictor
    μ = α + pm.math.dot(X_train_s, β)

    # Likelihood
    Y_obs = pm.Normal("Y_obs", mu=μ, sigma=σ, observed=y_train)

    # Sample posterior
    trace = pm.sample(
        draws=2000,
        tune=1000,
        target_accept=0.9,
        return_inferencedata=True
    )

Posterior Analysis & Point Predictions

Posterior means of α and β provide point forecasts; MAE on held‑out data quantifies average prediction error.

import arviz as az

# Posterior summary
az.summary(trace, round_to=2)

# Posterior predictive sampling
with model:
    ppc = pm.sample_posterior_predictive(trace, var_names=["Y_obs"])

# Extract posterior means
α_post = trace.posterior["α"].mean().item()
β_post = trace.posterior["β"].mean(dim=["chain","draw"]).values

# Compute point predictions on the test set
y_pred = α_post + X_test_s.dot(β_post)

# Evaluate MAE
from sklearn.metrics import mean_absolute_error
mae = mean_absolute_error(y_test, y_pred)
print(f"Test MAE: ${mae:.2f}")

Visualise Predictions & Credible Intervals

By sweeping one feature (distance) and holding others constant, we plot the posterior mean cost curve and its 94% Highest Posterior Density interval—illustrating both expected cost scaling and uncertainty.

import numpy as np
import matplotlib.pyplot as plt

# Vary distance_km; fix other features at their median
dist_grid = np.linspace(X_train_s[:,0].min(), X_train_s[:,0].max(), 100)
grid = np.tile(np.median(X_train_s, axis=0), (100,1))
grid[:,0] = dist_grid

with model:
    pm.set_data({"X": grid})
    ppc_grid = pm.sample_posterior_predictive(trace, var_names=["Y_obs"])

preds     = ppc_grid["Y_obs"]
mean_pred = preds.mean(axis=0)
hpd       = az.hdi(preds, hdi_prob=0.94)

# Convert distance back to original scale
dist_orig = scaler.inverse_transform(
    np.column_stack([grid[:,0], grid[:,1], grid[:,2]])[:, :3]
)[:,0]

plt.figure(figsize=(8,5))
plt.plot(dist_orig, mean_pred, label="Posterior mean")
plt.fill_between(dist_orig, hpd[:,0], hpd[:,1], alpha=0.3,
                 label="94% credible interval")
plt.scatter(
    scaler.inverse_transform(X_test_s[:, :3])[:,0], y_test,
    color="k", alpha=0.5, label="Test data"
)
plt.xlabel("Distance (km)")
plt.ylabel("Predicted Delivery Cost (USD)")
plt.title("Bayesian Regression: Cost vs. Distance")
plt.legend()
plt.tight_layout()
plt.show()

Summary

This Bayesian Regression workflow for delivery‐cost forecasting provides:

1. Point estimates of expected delivery cost from early indicators (distance, weight, stops, peak‑hour).

2. Credible intervals quantifying uncertainty—enabling risk‑aware pricing and resource planning.

3. Actionable insights: logistics teams can quote with both an expected cost and uncertainty bounds, optimize route planning, and allocate workforce under uncertainty.