Integrated Numerical Methods and Advanced Machine Learning for Road Accident Analysis in Bangladesh

1. Introduction

Road-traffic crashes are difficult to model because they do not arise from a single cause. Exposure, road geometry, vehicle mix, driver behavior, enforcement, emergency response, reporting quality, and seasonal variation all influence the number and severity of crashes observed in a given place and time. This complexity has pushed road-safety research well beyond simple descriptive counts. Current work includes multidimensional severity prediction [1], IoT- and vision-based emergency-response systems [2], simulation-assisted accident reconstruction [3], motorcycle-specific risk analysis [4], multi-road-user safety simulation [5], rural–urban temporal scaling [6], large-language-model-assisted investigation of crash causation [7], and GIS-based hotspot detection in Bangladesh [8].

Machine learning has widened the range of tools available for such problems, particularly where nonlinear relationships and interactions are difficult to specify in advance. Reviews of artificial intelligence in crash-severity prediction show a rapid increase in the use of ensemble learning, deep learning, and hybrid approaches [9,11]. Seasonal and demographic analyses nevertheless continue to show that time structure matters [10], while rule-based generative-AI systems are beginning to appear in traffic-accident detection and interpretation [12]. More focused studies have used machine learning for accident prediction [32], compared statistical and machine-learning approaches for severity classification [33], developed AI-based pedestrian crash models [34], applied machine learning to crash and hotspot prediction [35], and proposed improved learning algorithms for severity estimation [36]. These studies confirm the value of flexible predictive models, but they also highlight a familiar problem: high predictive complexity does not by itself ensure that a model is numerically stable, interpretable, or reliable outside the sample on which it was trained.

For that reason, mathematical and statistical modelling remains central to road-safety research. Mathematical formulations have been used to study accident mechanisms and safety interventions [13], while data-driven models have been applied to broader crash-prediction tasks [14]. Time-series methods such as SARIMA and Prophet have been used to represent temporal dependence and changepoints [15], and ensemble methods have been explored for crash-severity prediction [16]. In Bangladesh, mathematical modelling has been used to examine factors associated with the growth of road accidents [17], and Poisson regression has been applied to crash-frequency modelling [18]. Other studies have considered operational traffic optimization under accident conditions [19], SEIR-type representations of accident severity [20], severity classification [21], alternative crash-prediction formulations [22], mathematical approaches to road-safety management [23], and data-based estimation of vehicle-crash parameters and dynamic crash response [24,25]. This body of work shows that road-safety analysis is not a single modelling problem; it may involve approximation, likelihood estimation, numerical linear algebra, dynamic state estimation, and supervised prediction.

The reliability of any of these methods also depends on the data on which they are built. Modern road-safety databases may combine administrative records, geospatial information, sensor data, connected-vehicle data, and other digital sources, each with its own limitations [26]. Descriptive studies have shown the importance of contextual and catastrophic crash factors [27], while conventional accident analyses remain useful for identifying broad burden patterns [28]. Risk-indicator studies further demonstrate that geographic priorities can change when different exposure denominators are used [29]. Data-mining research also shows that preprocessing, harmonization, missing-value treatment, and feature construction can materially affect the final results [30]. For Bangladesh, the authors' earlier work combined nonlinear deterministic, exposure-adjusted count, and seasonal time-series models [31]. That study provided a useful foundation, but it also pointed to the need for a framework that gives greater attention to numerical stability and out-of-time predictive performance.

Several gaps therefore remain. Numerical conditioning is seldom reported in road-safety studies, even though coefficient estimates may become unstable when predictors are poorly scaled or nearly dependent. Matrix rank, singular values, condition numbers, and perturbation sensitivity are often omitted from otherwise detailed modelling studies. Regularization is also commonly treated only as a prediction tool, rather than as a way of stabilizing a numerical system. In count modelling, Poisson and negative-binomial formulations are frequently chosen without fully examining whether the assumed mean–variance relationship matches the data under the selected exposure definition. Temporal studies, meanwhile, often rely on a single forecasting or smoothing technique rather than comparing interpolation, numerical differentiation, latent-state filtering, and spectral structure. A further concern is validation: many predictive studies still use random train–test splits, although neighbouring months and years are not independent and such splits can overstate generalization performance. These issues are especially relevant in Bangladesh, where mathematical, statistical, spatial, and machine-learning studies have generally developed along separate lines [8,17,20,31].

The present study addresses these gaps by combining numerical analysis and advanced machine learning within one reproducible framework. Numerical stability is examined using matrix rank, singular-value decomposition, condition numbers, normal equations, QR factorization, SVD pseudoinverse, Tikhonov regularization, generalized cross-validation, and controlled perturbation analysis. Polynomial trend estimation is accompanied by explicit examination of Vandermonde conditioning and time rescaling. Monthly fatality dynamics are studied using natural cubic splines, finite differences, Savitzky–Golay smoothing, Kalman filtering with Rauch–Tung–Striebel smoothing, and Fourier decomposition. Conditional fatality counts are modelled using Poisson and NB2 likelihoods with accident frequency incorporated as an exposure offset. Alongside these numerical and statistical methods, fifteen predictive models are evaluated, including linear and regularized regression, robust regression, count generalized linear models, kernel and neighbour methods, tree ensembles, boosting algorithms, and two Adam-trained artificial neural networks.

The aim of the study is to develop a numerically reliable and temporally credible framework for analysing road-accident burden in Bangladesh from 2023 to 2025. The analysis examines how accidents, fatalities, and injuries vary across years and divisions after population standardization; whether the numerical systems used for inference and approximation are well conditioned; what count, smoothing, state-space, and spectral methods reveal about fatality dynamics; and how numerical and machine-learning models perform when an entire year is withheld from training. The emphasis is therefore not only on which model fits best, but on which conclusions remain credible after numerical diagnostics, exposure adjustment, and out-of-time validation are taken into account.

2. Materials and Methods

2.1. Study Design, Data Sources, and Analytical Units

This retrospective analytical study covers January 2023 through December 2025. The primary analytical unit is the division-month, comprising eight administrative divisions observed over 36 consecutive months and yielding a balanced panel of 288 observations. Aggregation across divisions produces a 36-point national monthly series used for interpolation, numerical differentiation, state-space filtering, and spectral analysis. Four synchronized datasets are employed: 2022 division-level population statistics, the principal monthly road-accident panel, fatalities by vehicle type and division, and recorded vehicle counts by vehicle type and division. The overall analytical workflow, from data sources and preprocessing through numerical/statistical modelling, advanced machine learning, time-aware validation, model comparison, and road-safety interpretation, is summarized in Figure 1. The processed-data provenance is documented in the authors’ earlier Bangladesh road-safety study, which identifies Bangladesh Road Transport Authority road-accident statistics and Bangladesh Bureau of Statistics population statistics as the underlying public sources [31].

The complete analytical workflow was implemented in Python using NumPy, pandas, SciPy, Matplotlib, scikit-learn, and openpyxl. A fixed random seed of 42 was used for stochastic procedures, and bootstrap uncertainty was estimated using 5,000 resamples. The script also exports cleaned tables, a consolidated Excel workbook, solver diagnostics, cross-validated predictions, model metrics, a cleaning log, and reproducibility metadata to support transparent and reproducible analysis.

As summarized in Table 1, the main panel is complete for the variables used in modelling, whereas the vehicle files contain sparse missingness. Known division and vehicle-label inconsistencies are harmonized to canonical labels while retaining source values in the audit workflow. Numeric fields are validated, invalid values are converted to missing, analytical keys are checked for duplicates, and a first-of-month date index is constructed. Temporal ordering is preserved throughout.

2.2. Outcomes, Exposure Definitions, and Leakage Control

Let A_it, F_it, and I_it denote accidents, fatalities, and injuries for division i in month t, and let P_i denote the corresponding 2022 population. For a Y-year observation window, the annualized fatality rate per 100,000 persons is defined by

Rᶠᵢ = (100000 / Y) [ (Σₜ Fᵢₜ) / Pᵢ ]. (1)

The population denominator is intentionally fixed at 2022 values. Thus, the rate is a standardized three-year comparison rather than an estimate that incorporates annual population change. Division uncertainty is quantified from 5,000 nonparametric bootstrap resamples of monthly rates and percentile 95% confidence intervals.

Predictive leakage is controlled by excluding variables algebraically derived from each target. Fatalities-per-accident is therefore excluded when fatalities are predicted, and injuries-per-accident is excluded for injuries. Observed accident count is retained as an exposure-like covariate for conditional fatality and injury modelling; consequently, those predictions are interpreted as consequence models conditional on accident frequency, not as pre-crash forecasts.

Continuous predictors entering scale-sensitive models are standardized as

zⱼ = (xⱼ - x̄ⱼ) / sⱼ, (2)

and categorical predictors are one-hot encoded within each training fold so that preprocessing parameters are learned without using the held-out year.

2.3. Numerical Linear Algebra and Solver Audit

For an n x p design matrix X, singular-value decomposition (SVD) is written as

X = U Σ Vᵀ, (3)

where U and V have orthonormal columns and Σ contains non-negative singular values. Numerical sensitivity is summarized by the 2-norm condition number

κ₂(X) = σmax(X) / σmin(X). (4)

Four coefficient solvers are compared under the same linear model. The normal-equation estimator is

β̂_NE = (XᵀX)⁻¹ Xᵀy, (5)

whereas a thin QR factorization X = QR gives

β̂_QR = R⁻¹ Qᵀy. (6)

The SVD pseudoinverse solution is

β̂_SVD = V Σ⁺ Uᵀ y, (7)

where Σ+ reciprocates nonzero singular values subject to numerical tolerance. Tikhonov regularization uses

β̂_λ = (XᵀX + λD)⁻¹ Xᵀy, (8)

with D equal to the identity except for a zero intercept entry. Generalized cross-validation (GCV) selects λ using the smoothing matrix H_λ = X(XᵀX+λD)⁻¹Xᵀ:

GCV(λ) = n ||(I - H_λ)y||₂² / [n - tr(H_λ)]². (9)

The solver audit evaluates residual error, condition number, and sensitivity to controlled perturbations. This distinction is important because normal equations square the singular-value ratio, so κ₂(XᵀX) is approximately κ₂(X)² for a full-rank system.

2.4. Polynomial Approximation, Spline Interpolation, and Numerical Differentiation

To avoid the severe scaling of raw powers of month index t, time is mapped to [-1,1]:

xₜ = 2(t - tmin)/(tmax - tmin) - 1. (10)

For degree d, the Vandermonde matrix V_d = [1, x, x², ..., xᵈ] is fitted by least squares:

β̂(d) = arg min_β ||V_d β - y||₂². (11)

Degrees 1-8 are evaluated by leave-one-out cross-validation (LOOCV) on the 36-month national fatality series. Both raw- and scaled-basis condition numbers are reported to separate approximation error from numerical instability.

A natural cubic spline S(t) with zero second derivative at the two boundaries interpolates the national fatality series. Because exact interpolation can amplify local noise in derivatives, the analytical spline derivative is compared with second-order finite differences,

y'(t) ≈ [y(t+h) - y(t-h)] / (2h), h = 1 month, (12)

and with Savitzky-Golay local-polynomial smoothing. For a symmetric window j = -m,...,m and polynomial degree p, the local coefficients solve

minₐ Σⱼ₌₋ₘᵐ [ yₜ₊ⱼ - Σᵣ₌₀ᵖ aᵣ jʳ ]², ŷₜ = a₀. (13)

Agreement in derivative sign across methods is treated as stronger evidence of a turning direction than agreement in peak amplitude.

2.5. Exposure-Adjusted Poisson and Negative-Binomial Likelihoods

Fatality counts are modelled conditionally on accident frequency using log(A_it) as an offset. The Poisson specification is

Fᵢₜ ~ Poisson(μᵢₜ), log μᵢₜ = log Aᵢₜ + xᵢₜᵀβ. (14)

The covariate vector contains standardized year index, sine and cosine month terms, population, population density, dependency ratio, literacy rate, and seasonal indicators. Division fixed effects are omitted from this inferential specification because the division-level demographic predictors are time-invariant and would become exactly collinear with a complete division indicator set.

For the NB2 alternative, the conditional mean remains μ while variance is

Var(Fᵢₜ | xᵢₜ) = μᵢₜ + α μᵢₜ², α ≥ 0. (15)

Using r = 1/α, the NB2 probability mass function is

P(F=f) = Γ(f+r)/[Γ(r)Γ(f+1)] [r/(r+μ)]ʳ [μ/(r+μ)]ᶠ. (16)

Poisson and NB2 parameters are estimated by direct numerical maximum likelihood. Model support is compared using log-likelihood, AIC, BIC, and Poisson Pearson dispersion. For a standardized predictor coefficient β_j, the incidence-rate ratio is

IRRⱼ = exp(βⱼ). (17)

2.6. Local-Linear-Trend Kalman Filter and RTS Smoother

The national fatality series is also represented as a noisy observation of a latent local level and slope. With state vector s_t = [ℓ_t, b_t]^T, the transition equation is

sₜ = G sₜ₋₁ + wₜ, G = [[1,1],[0,1]], wₜ ~ N(0,Q), (18)

and the observation equation is

yₜ = H sₜ + vₜ, H = [1,0], vₜ ~ N(0,R). (19)

The Kalman prediction and correction equations are

sₜ⁻ = Gsₜ₋₁; Pₜ⁻ = GPₜ₋₁Gᵀ+Q; Kₜ = Pₜ⁻Hᵀ(HPₜ⁻Hᵀ+R)⁻¹; sₜ = sₜ⁻ + Kₜ(yₜ-Hsₜ⁻). (20)

Variance parameters are estimated by maximizing the Gaussian innovation likelihood. The Rauch-Tung-Striebel (RTS) backward recursion smooths filtered states using

Jₜ = PₜGᵀ(Pₜ₊₁⁻)⁻¹, sₜ|T = sₜ + Jₜ(sₜ₊₁|T - sₜ₊₁⁻). (21)

2.7. Fourier Decomposition

A linear trend is removed before fast Fourier transformation. If z_t denotes the detrended series, a K-harmonic reconstruction is

ŷₜ = a + bt + Σₖ₌₁ᴷ Aₖ cos(2π fₖ t + φₖ). (22)

Spectral components are ranked by power, periods are reported as 1/f_k months, and the four strongest non-zero harmonics are retained for reconstruction. With only 36 monthly observations, periodicity is interpreted as evidence of recurring structure rather than as a precise long-horizon frequency estimate [6,15].

2.8. Time-Aware Statistical and Machine-Learning Models

Fifteen predictive model families are compared under identical preprocessing: ordinary least squares, Ridge, Elastic Net, Huber robust regression, Poisson GLM, Tweedie compound-Poisson GLM, RBF support-vector regression, distance-weighted KNN, two- and three-hidden-layer neural networks trained with Adam, random forest, Poisson random forest, Extra Trees, Huber gradient boosting, and Poisson histogram gradient boosting. The expanded library spans linear, regularized, robust, count-aware, kernel, local-neighbor, neural, bagging, and boosting structures. Key computational settings from the reproducible analysis script are listed in Table 2.

For ordinary linear regression, the empirical squared-error objective is

β̂_OLS = arg min_β Σᵢ (yᵢ - xᵢᵀβ)². (23)

Ridge regression solves

β̂_R = arg min_β { Σᵢ (yᵢ - xᵢᵀβ)² + α||β||₂² }, (24)

while Elastic Net solves

β̂_EN = arg min_β { (1/2n)||y-Xβ||₂² + α[ρ||β||₁ + (1-ρ)||β||₂²/2] }. (25)

Huber regression minimizes a robust loss plus L2 regularization, with residual loss

Lδ(r) = 0.5r² if |r|≤δ; δ(|r|-0.5δ) otherwise. (26)

The Poisson GLM uses a logarithmic mean link

μᵢ = exp(xᵢᵀβ), (27)

whereas the Tweedie compound-Poisson model uses the variance function

Var(Yᵢ|xᵢ) = φ μᵢᵖ, p = 1.5. (28)

RBF support-vector regression uses the kernel K(x,x') = exp[-γ||x-x'||²] and the ε-insensitive primal objective

min 0.5||w||² + CΣᵢ(ξᵢ+ξᵢ*), subject to |yᵢ-f(xᵢ)| ≤ ε + slack. (29)

Distance-weighted KNN predicts from the k nearest training observations as

ŷ(x) = [Σⱼ∈Nₖ(x) d(x,xⱼ)⁻¹ yⱼ] / [Σⱼ∈Nₖ(x) d(x,xⱼ)⁻¹]. (30)

For the neural networks, hidden-layer propagation is

h⁽ˡ⁾ = ReLU(W⁽ˡ⁾h⁽ˡ⁻¹⁾ + b⁽ˡ⁾), ŷ = W⁽ᴸ⁾h⁽ᴸ⁻¹⁾ + b⁽ᴸ⁾. (31)

Adam updates first and second gradient moments and then the parameters:

mₜ=β₁mₜ₋₁+(1-β₁)gₜ; vₜ=β₂vₜ₋₁+(1-β₂)gₜ²; θₜ=θₜ₋₁-η m̂ₜ/(√v̂ₜ+ε). (32)

Bagged tree ensembles predict by averaging B regression trees,

ŷ_RF(x) = (1/B) Σᵦ₌₁ᴮ Tᵦ(x), (33)

with Poisson random forest changing the node-split criterion and Extra Trees increasing split randomization. Gradient boosting constructs an additive predictor

f_M(x) = f₀(x) + Σₘ₌₁ᴹ ν hₘ(x), (34)

where each h_m approximates the negative gradient of the selected loss; Huber and Poisson losses are used for the two boosting variants.

2.9. Temporal Validation and Performance Metrics

Validation uses GroupKFold with calendar year as the grouping variable, generating three leave-one-year-out folds: train on two years and predict the omitted year. No headline result is based on a random split. For observations y_i and out-of-fold predictions ŷ_i, root-mean-square error is

RMSE = √[(1/n) Σᵢ (yᵢ-ŷᵢ)²], (35)

mean absolute error is

MAE = (1/n) Σᵢ |yᵢ-ŷᵢ|, (36)

and the coefficient of determination is

R² = 1 - [Σᵢ(yᵢ-ŷᵢ)² / Σᵢ(yᵢ-ȳ)²]. (37)

Symmetric mean absolute percentage error is

SMAPE = (100/n) Σᵢ [2|yᵢ-ŷᵢ|/(|yᵢ|+|ŷᵢ|)], (38)

and normalized RMSE is

NRMSE(%) = 100 × RMSE / ȳ. (39)

For positive count predictions, mean Poisson deviance is

D_P = (2/n) Σᵢ [ yᵢ log(yᵢ/ŷᵢ) - (yᵢ-ŷᵢ) ], (40)

with the convention y log(y/ŷ)=0 at y=0. Spearman rank correlation is also reported. Model selection within each target is based on minimum out-of-fold RMSE, with NRMSE used for cross-outcome comparisons.

2.10. Vehicle-Exposure Screening Ratio

Vehicle fatality and recorded-vehicle files are matched by year, month, division, and harmonized vehicle type. Only non-missing records with strictly positive recorded-vehicle denominators are eligible. For vehicle type v,

VFRᵥ = 100 × [Σ Fᵥ / Σ Vᵥ]. (41)

The quantity is explicitly described as fatalities per 100 recorded vehicles. Because the denominator is a source-defined recorded count rather than verified registered-fleet exposure or vehicle-kilometres travelled, VFR is a descriptive screening indicator and not a causal or actuarial fatality probability.

3. Results

3.1. National Burden and Annual Change

Across 2023-2025, the synchronized panel contains 17,046 accidents, 15,994 fatalities, and 20,389 injuries. Table 3 shows that the three outcomes do not move in parallel. Accidents increased by 6.57% in 2024 and then decreased by 2.75% in 2025; fatalities increased by 9.08% in 2024 and remained almost unchanged in 2025 (+0.18%); injuries fell by 13.68% in 2024 and by another 0.71% in 2025.

Figure 2 visualizes the divergence. The 2024 increase in accident frequency and fatal burden occurred simultaneously with a marked injury reduction, while the 2025 reduction in accidents did not produce a reduction in fatalities. This supports modelling accidents, fatalities, and injuries as related but non-interchangeable outcomes [16,21,22].

3.2. Spatial Disparity and Bootstrap Uncertainty

Population standardization materially changes geographic prioritization. Table 4 ranks the eight divisions using fixed 2022 population exposure. Barishal has the highest annualized fatality rate (5.44 per 100,000), followed by Sylhet (4.44) and Rajshahi (3.83), whereas Rangpur (2.23) and Khulna (2.29) are lowest.

Figure 3 presents the ranked annualized rates and confirms that high raw-count divisions are not necessarily the highest-rate divisions. The figure therefore emphasizes the distinction between scale and standardized burden.

Figure 4 adds bootstrap uncertainty. Barishal's mean monthly rate is 0.453 per 100,000 (95% CI 0.383-0.528), compared with 0.186 (0.161-0.213) in Rangpur. The separation is consistent with risk-indicator research showing that prioritization depends on both spatial unit and denominator [29].

3.3. Correlation, Rank, and Numerical Conditioning

Figure 5 shows strong dependence among several safety and demographic variables. Accidents and fatalities are highly correlated (r=0.969), accidents and injuries are positively correlated (r=0.762), population correlates with accidents (r=0.735) and fatalities (r=0.743), and population density correlates with fatalities (r=0.595). These relationships motivate multivariable modelling while also making numerical diagnostics necessary.

Table 5 summarizes the core conditioning diagnostics. The inferential fatality design is full rank (11/11), with κ2(X)=16.73 and κ2(XᵀX)=279.97. The largest and smallest singular values are 26.517 and 1.585. In contrast, the raw cubic Vandermonde matrix is severely ill-conditioned until the time coordinate is scaled.

Figure 5 displays the singular-value spectrum. Its gradual decline, rather than collapse toward a nearly zero terminal singular value, supports the full-rank diagnosis and indicates that the inferential coefficient system is not dominated by one near-null direction.

Figure 6 evaluates ridge stabilization. The normal-equation system begins near a condition number of 279.97 and reaches a minimum near 4.79 around λ≈355.4 before rising at very large λ because the intercept is deliberately unpenalized. The non-monotonic path demonstrates why regularization should be inspected numerically rather than assumed to improve conditioning without limit.

Figure 7 isolates polynomial-basis conditioning. At the LOOCV-selected degree 3, scaling the month coordinate from its raw integer form to [-1,1] reduces κ(V) from 6.21×10^4 to 7.77. At degree 8, the raw condition number reaches approximately 6.21×10^12, whereas the scaled basis remains several orders of magnitude better conditioned.

3.4. Polynomial Trend, Smoothing, and Numerical Derivatives

LOOCV selects a cubic polynomial for the national fatality series. The degree-3 fit RMSE is 78.13 and LOOCV RMSE is 84.80. Increasing degree beyond three provides only modest in-sample improvement while worsening out-of-sample error and rapidly inflating the raw-basis condition number. Selected numerical time-series diagnostics are summarized in Table 6.

Figure 8 compares observed fatalities with the exact natural cubic spline, the LOOCV-selected cubic polynomial, and the Savitzky-Golay smoothed trend. The spline passes through every observation, the cubic polynomial captures broad low-frequency movement, and Savitzky-Golay smoothing preserves major turning behaviour while attenuating local noise.

Figure 9 compares first derivatives from the exact spline, second-order finite differences, and Savitzky-Golay smoothing. Around the pronounced mid-2024 decline, the approximate minima are -323, -206, and -88 fatalities per month, respectively. All methods agree on direction and timing, but the exact interpolant produces the largest amplitude; consequently, derivative sign is more robust than peak magnitude.

3.5. Time-Aware Predictive Performance

Predictive ranking is strongly outcome-specific under leave-one-year-out validation. Figure 10 compares NRMSE across the complete statistical and machine-learning library. Regularized and robust linear models remain competitive for accidents and fatalities, whereas a count-aware random forest is preferred for injuries. Neural networks converge successfully but do not dominate out-of-time performance.

Table 7 reports the best machine-learning/statistical model for each outcome. Ridge regression is best for accidents (RMSE=24.20; R²=0.467), Huber robust regression is best for fatalities (RMSE=7.97; R²=0.937), and Poisson random forest is best for injuries (RMSE=29.89; R²=0.633). The fatality result should be interpreted conditionally on observed accident frequency.

Figure 11, Figure 12 and Figure 13 provide observed-versus-predicted diagnostics for the three best target-specific models. Figure 11 shows that accident predictions follow the one-to-one trend but shrink toward the centre, consistent with the moderate R².

Figure 12 shows a substantially tighter relationship for fatalities under Huber robust regression, consistent with R²=0.937 and Spearman ρ=0.971.

Figure 13 shows the injury predictions from Poisson random forest. The model captures the principal ordering and much of the outcome spread, although extreme observations retain larger residuals.

3.6. Vehicle-Exposure Screening Patterns

Among matched observations with positive recorded-vehicle denominators, motorcycles have the highest overall exposure-style fatality ratio at 85.24 fatalities per 100 recorded vehicles, followed by auto rickshaws (81.68), other vehicles (70.01), easy bikes (66.06), and vans (61.47). Table 8 lists the eight highest categories.

Figure 14 visualizes the overall vehicle ranking. The concentration of high values among motorcycles and paratransit categories is consistent with literature emphasizing motorcycle and vulnerable-road-user risk [4,27], but the source-defined denominator requires cautious interpretation.

Figure 15 demonstrates substantial division-by-vehicle heterogeneity. Some cells exceed 100 fatalities per 100 recorded vehicles; this does not represent a probability greater than one because the denominator is an aggregated source-defined recorded count rather than a fixed cohort at risk. The heat map therefore strengthens the denominator-sensitivity warning.

3.7. Poisson Versus NB2 Conditional Fatality Inference

Both direct maximum-likelihood optimizations converge. Table 9 shows that Poisson attains log-likelihood -953.50, AIC 1928.99, and BIC 1969.28, whereas NB2 produces log-likelihood -953.55, AIC 1931.09, and BIC 1975.05. The NB2 overdispersion parameter is only 7.33×10^-5 and Poisson Pearson dispersion is 0.969, indicating near-equidispersion under the accident-offset formulation.

Table 10 reports the principal Poisson incidence-rate ratios. A one-standard-deviation increase in year index is associated with a 1.9% increase in conditional fatality incidence (IRR=1.019, 95% CI 1.003-1.035, p=0.019), while a one-standard-deviation increase in population density is associated with a 5.2% increase (IRR=1.052, 95% CI 1.022-1.083, p<0.001). Other listed covariates are not significant at the 5% level.

Figure 16 shows observed and fitted fatalities for Poisson and NB2. The fitted curves are visually almost indistinguishable, exactly as expected when the estimated NB2 dispersion collapses toward zero. This supports diagnosing overdispersion for the specific outcome and exposure definition rather than transferring it mechanically from other crash-count formulations [18,22,31].

3.8. State-Space and Spectral Fatality Dynamics

Figure 17 presents the local-linear-trend Kalman/RTS estimate. The latent level is substantially smoother than the observed monthly series, while the state uncertainty band expands in portions of the series with larger observational departures. The state-space view treats month-to-month variation as a combination of latent evolution and measurement/process noise rather than interpreting every observed fluctuation as a structural shift.

Figure 18 shows the detrended Fourier power spectrum. The strongest non-zero component has a 12-month period (power 23,947; amplitude 51.58), followed by 9- and 6-month components. The annual component is therefore the clearest recurring frequency in the 36-month sample.

Figure 19 reconstructs the national series from the linear trend and four strongest harmonics. The reconstructed series has RMSE 52.46 and MAE 41.49, and the estimated linear component increases by approximately 1.14 fatalities per month over the study window. The reconstruction captures broad periodic structure without claiming precise long-horizon seasonality from only three annual cycles.

3.9. GCV-Selected Tikhonov Regularization and ANN Optimization

Figure 20 displays the GCV criterion together with the associated conditioning path for Tikhonov regularization. The selected regularization represents a compromise between fit and coefficient-system stability, and it becomes the strongest numerical baseline in the direct numerical-versus-ML comparison reported below.

Figure 21 shows the training-loss histories for the three-hidden-layer Adam networks across targets. The curves demonstrate successful optimization and early-stopping behaviour, but optimization convergence should not be confused with out-of-time predictive superiority.

Figure 22, Figure 23 and Figure 24 show ANN-3H-Adam observed-versus-predicted diagnostics for accidents, fatalities, and injuries. The neural networks learn nontrivial relationships, especially for fatalities and injuries, yet the accident predictions exhibit substantial shrinkage and dispersion.

3.10. Direct Comparison of Numerical Solvers and Machine Learning

The direct family comparison in Figure 25 uses NRMSE so that the three outcomes can be compared on the same relative scale. Tikhonov GCV is the best numerical solver for all three outcomes, whereas Ridge regression, Huber robust regression, and Poisson random forest are the best machine-learning/statistical models for accidents, fatalities, and injuries. Table 11 quantifies the margins.

Figure 26 compares annual out-of-time predictions from the best numerical and best machine-learning/statistical model. The plot makes clear that model-family differences are not constant across years; some years are captured similarly by both families, while others account for most of the relative-performance gap.

Figure 27 places all numerical and machine-learning models in a single NRMSE heat map. Table 12 reproduces the displayed values for auditability. The result is not a simple numerical-versus-ML dichotomy: Tikhonov GCV is nearly tied with Ridge for accidents and with Huber/Elastic Net for fatalities, while the larger advantage of Poisson random forest appears for injuries. ANN-3H-Adam produces NRMSE values of 79.0%, 23.9%, and 43.3% for accidents, fatalities, and injuries, respectively; successful neural optimization therefore does not guarantee superior temporal generalization.

4. Discussion

4.1. Principal Findings and Numerical Interpretation

This study produces four principal findings.

First, road-safety burden is outcome-specific and spatially heterogeneous: accidents, fatalities, and injuries follow different annual trajectories, and population standardization places Barishal and Sylhet above the larger divisions in annualized fatality burden.

Second, the main inferential design is numerically usable but not immune to sensitivity: κ2(X)=16.73 is moderate, whereas κ2(XᵀX)=279.97 illustrates the amplification created by normal equations. The polynomial subsystem provides the sharper warning because scaling reduces the cubic Vandermonde condition number by almost four orders of magnitude.

Third, the conditional fatality counts are close to equidispersed, so Poisson is more parsimonious than NB2 under the accident-offset specification.

Fourth, out-of-time performance is model- and outcome-specific: no single family dominates all targets.

The polynomial result is particularly relevant to applied mathematics. Increasing model degree can reduce in-sample residuals while simultaneously degrading both conditioning and cross-validated reconstruction. This is a concrete example of why approximation quality must be considered jointly with numerical sensitivity. The SVD spectrum, Ridge path, Tikhonov GCV curve, and raw-versus-scaled Vandermonde comparison together show that numerical diagnostics are not ancillary plots; they determine whether estimated coefficients and trends are reproducible under finite precision and small data perturbations [19,23,24,25].

4.2. Count Inference, Exposure, and Denominator Sensitivity

The Poisson-versus-NB2 comparison also illustrates why distributional assumptions should be evaluated against the exact modelling question. Count data are often assumed to require negative-binomial overdispersion, but here fatalities are modelled conditional on observed accident exposure. With Pearson dispersion 0.969 and NB2 alpha close to zero, the extra NB2 parameter is not supported by AIC or BIC. The significant positive association of population density with conditional fatality incidence is consistent with the idea that urban exposure and consequence mechanisms deserve attention, but the model remains associational rather than causal [18,22].

The vehicle-exposure analysis requires even greater care. Motorcycle and paratransit categories rank highly, consistent in direction with motorcycle-risk and catastrophic-factor studies [4,27], yet the denominator is a source-defined recorded-vehicle count. A ratio above 100 in some division-vehicle cells is mathematically possible under this aggregate denominator and should not be interpreted as a probability. This is exactly the type of denominator problem highlighted by general risk-indicator and data-source research [26,29].

4.3. Temporal Dynamics and Model Generalization

The agreement between Fourier analysis and the smoothed time-series representations provides complementary evidence of recurring annual structure. The 12-month component is the strongest non-zero spectral frequency, while the Kalman/RTS model separates a slowly evolving latent level from short-term observation noise. These approaches answer different questions: Fourier analysis characterizes recurrent periodic components, whereas state-space estimation asks how the latent level evolves when observations are noisy. The result aligns with prior seasonal forecasting and scaling studies [6,10,15].

The predictive results reinforce the importance of temporal validation. Ridge regression for accidents, Huber regression for fatalities, and Poisson random forest for injuries outperform more complex alternatives on their respective targets. ANN-2H and ANN-3H models achieve stable training convergence, yet their NRMSE values do not establish superiority. This finding is consistent with broader reviews warning that high-capacity AI can underperform when samples are small, temporally structured, or weakly informative [9,11]. The practical implication is that model sophistication should be justified by out-of-time evidence rather than by architecture alone.

The direct numerical-versus-ML comparison adds an applied-mathematical perspective that is often missing from road-safety benchmarking. Tikhonov GCV remains within 0.2-0.3 NRMSE percentage points of the best data-driven model for fatalities and accidents, respectively, while the advantage of Poisson random forest becomes more visible for injuries. This pattern suggests that regularized numerical solvers can serve as strong transparent baselines. Machine learning should be interpreted as an incremental modelling choice whose value depends on outcome structure, not as a universal replacement for numerical methods.

4.4. Policy Relevance

For road-safety management in Bangladesh, three implications follow. First, geographic prioritization should report both absolute burden and population-standardized rates because these lead to different rankings. Second, conditional fatality models can identify covariate associations and consequence patterns after accounting for accident exposure, but they should not be described as crash-occurrence models. Third, operational model selection should use future-period or leave-period-out validation whenever deployment will occur in a later year. GIS hotspot analysis in Bangladesh [8], data-mining frameworks [30], and the authors' prior national modelling [31] provide complementary tools that could be integrated with the present numerical diagnostics in future decision-support systems.

4.5. Limitations and Future Work

The study has several limitations. The panel is aggregated at division-month level and spans only three years, limiting the identification of rare events, local roadway mechanisms, and long-period spectral components. Population is fixed at 2022 values, so rates are standardized comparisons rather than dynamic demographic rates. Vehicle-count denominators are source-defined and are not equivalent to registered fleet, traffic volume, or vehicle-kilometres travelled. Accident count is used as an exposure-like predictor for fatality and injury models, which improves conditional consequence prediction but prevents interpretation as a pre-crash forecast. Finally, no model can correct under-reporting, source-definition changes, or omitted roadway/weather/driver variables that are absent from the administrative data.

Future research should extend the panel temporally, incorporate crash-location geometry, roadway design, weather, enforcement, traffic volume, and emergency-response variables, and test geographic transfer across finer spatial units. Hierarchical Bayesian or spatiotemporal count models could explicitly represent spatial dependence; causal designs could address intervention effects; and external validation on a later, completely unseen year would provide a stronger deployment test. For machine learning, nested temporal tuning, calibrated uncertainty, explainability, and model-drift monitoring should accompany any move from retrospective benchmarking to operational prediction [7,9,11,12].

5. Conclusions

An integrated numerical and machine-learning analysis of Bangladesh road accidents during 2023-2025 shows that numerical reliability and temporal generalization are as important as predictive flexibility. The fatality design matrix is full rank with moderate conditioning, but raw polynomial time bases become severely ill-conditioned; scaling the cubic Vandermonde system reduces κ(V) from 6.21×10^4 to 7.77. Conditional fatality counts are nearly equidispersed, and Poisson is preferred to NB2 by AIC/BIC. The national fatality series contains a dominant 12-month component, while spline, finite-difference, Savitzky-Golay, and Kalman/RTS methods provide complementary views of local change and latent trend.

Under leave-one-year-out validation, Ridge regression is best for accidents, Huber robust regression for fatalities, and Poisson random forest for injuries. Tikhonov GCV remains a competitive transparent numerical baseline, especially for accidents and fatalities, whereas the ANN models demonstrate that successful optimization does not guarantee best out-of-time performance. Population-standardized burden is highest in Barishal at 5.44 fatalities per 100,000. Overall, the study supports a road-safety workflow in which exposure definitions, numerical conditioning, likelihood adequacy, and temporal validation are checked before policy conclusions are drawn.